Exposures from Oil and Gas Development and Childhood Leukemia Risk in Colorado: A Population-Based Case–Control Study Open Access

In communities near upstream oil and natural gas (O&G) development, childhood leukemia is of particular interest. Childhood leukemia may be an early health outcome from environmental exposure to leukemogens emitted from upstream O&G sites (1) due to its relatively short latency and the particular vulnerability of children (2, 3). Among children aged 0 to 14 years, acute lymphoblastic leukemia (ALL), one of the most commonly diagnosed cancers, has a peak incidence occurring at 2 to 5 years of age (4).

Herein, “upstream O&G development” refers to the activities related to the extraction of O&G resources from onshore below-ground resources. The rapid expansion of upstream O&G development in the 21st century (5) has resulted in an extensive dispersion of O&G wells across populated areas in the United States. Approximately 1.37 to 2.7 million children less than 5 years of age live within 1.6 km (1 mile) of at least one O&G well (6, 7). In areas of intense upstream O&G development, there can be hundreds of wells within 1.6 km of a home, and many wells (e.g., 3–40) can be co-located on a single site or well pad (8). Some of the most intensive upstream O&G activity is located in Colorado, where 35,290 O&G wells were drilled between 2003 and 2016 (9), and more than 378,000 people live within 1.6 km of at least one O&G well (8). This population is growing at a faster rate and may be at an economic disadvantage compared with Colorado’s general population (8). For example, in Colorado’s Denver-Julesburg Basin, a larger proportion of low-value homes are located within 152 m of an O&G well site (8).

An early study using an ecological county-level design and limited by an inadequate latency period did not observe more cases of childhood leukemia than would be expected in counties where natural gas wells were present (10, 11). Our registry-based case–control study in rural Colorado found that children aged 5 to 24 years diagnosed with ALL were 3.2 to 4.6 times as likely to live in areas with active O&G wells than children diagnosed with nonhematologic cancers, and the association between ALL and residential density of O&G wells increased monotonically from the lowest to highest inverse distance weighted (IDW) well count categories (12). A population-based case–control study in Pennsylvania found that children aged 2 to 7 years diagnosed with ALL were 1.74 times more likely to have at least one O&G well within 2 km of their birth residence (13). Although the two case–control studies improve on the early ecological study by evaluating exposure and outcomes at the individual level, full latency periods, and specific leukemia types, limitations remain. These include a low number of cases near upstream O&G development and exposure misclassification introduced by a lack of temporal and spatial specificity in simple IDW well counts.

As described elsewhere (14), upstream O&G development is a complex, multiphase process that continuously emits leukemogenic air pollutants potentially affecting all individuals residing nearby (15, 16). Benzene is the most documented leukemogen emitted from upstream O&G development activities (17), with ambient concentrations exceeding the US Environmental Protection Agency 0.102 parts per billion by volume risk-based screening level for benzene in residential air (18). The lifetime excess cancer risk is estimated to be 41 times higher for people living within 152 m and three times higher for people living within 610 m of Colorado’s upstream O&G facilities than for people living further than 1.6 km, with benzene contributing more than 95% of the risk estimates (19). These cancer risk estimates considered increased susceptibility to early life benzene exposure. Additionally, a biomarker study observed a fivefold increase in urinary trans-muconic acid (a metabolite of benzene) levels in pregnant women living near upstream O&G development (20).

Additional leukemogens, such as 1,3-butadiene, formaldehyde, and polycyclic aromatic hydrocarbons (PAH), are also present in diesel exhaust emitted from the thousands of trucks, generators, compressor engines, and other equipment supporting the upstream development of O&G resources (1). Concentrations of two indicators of diesel exhaust particulate matter less than 2.5 microns and black carbon almost doubled during the development of a multiwell O&G site in Colorado (21).

A growing body of evidence demonstrates associations between childhood ALL and maternal exposure to benzene (22, 23) as well as prenatal exposure to ambient levels of several leukemogens, including benzene (24, 25), PAHs (26), and traffic-related pollution (27, 28). These lines of evidence, coupled with the peak incidence of childhood ALL at 2 to 5 years of age (4), suggest that perinatal and early life exposure to leukemogens emitted during upstream O&G development could contribute to childhood ALL pathogenesis.

Our study objective is to further advance understanding of the association between upstream O&G development and childhood ALL by accounting for temporal and spatial variation in activity at O&G well sites in a large, population-based case–control study of childhood ALL in Colorado. This allows us to provide new information on how distance from, density of, and activity on well sites affect ALL risk in local and regional populations. Colorado is particularly well suited for this study because many people live within 1.6 km of an O&G well, and the state maintains a publicly available O&G well dataset. Our hypothesis is that upstream O&G activity increases the risk of ALL among children living near O&G well sites.

We conducted a population-based case–control study of 3,157 children born in Colorado between 1992 and 2019. We used birth and cancer registry data provided by the Colorado Department of Public Health and Environment’s (CDPHE) Center for Health Environmental Data and the Colorado Central Cancer Registry (CCCR). The data were deidentified to maintain the confidentiality of registry records. The Colorado birth registry includes information from birth certificates for all live births in Colorado. The CCCR uses active and passive surveillance to identify more than 99% of incident cancer cases in Colorado.

Our study population includes all children with a birth date between 1992 and 2019 in Colorado’s birth registry. From this cohort, CDPHE staff selected cases and controls as described in the next section and shown in Fig. 1.

CDPHE staff selected all children aged 2 to 9 years meeting the following criteria: (i) a confirmed ALL diagnosis in the CCCR between 2002 and 2021 and (ii) no previous diagnosis of cancer of any kind reported to the CCCR. We selected these criteria to (i) include the two most prominent age ranges of ALL incidence (29) and exclude children with infant ALL that differs in etiology from childhood ALL (30), (ii) exclude ALL associated with previous cancer treatment (31), and (iii) focus on the growth of unconventional O&G development, which began in 2002 in Colorado. Unconventional O&G development is characterized by hydraulic fracturing and/or directional drilling technologies (5). CDPHE staff next matched these ALL cases to Colorado’s birth registry and excluded ALL cases without a Colorado birth certificate. Based on these criteria, 451 children with ALL were selected as cases.

For each case, CDPHE staff randomly selected six controls per case without a diagnosis of cancer of any type in the CCCR from Colorado’s birth registry. Controls were matched to cases based on birth month, year, and ethnicity (Hispanic and not Hispanic) and assigned a reference date corresponding to the ALL diagnosis date of their matched case. Using these criteria, 2,706 children were selected as controls.

We then applied intensity-adjusted IDW (IA-IDW; ref. 13) proximity measures for O&G wells and other air pollution sources to estimate exposures as described in the following paragraphs.

Using information available in the publicly accessible Colorado Oil and Gas Information System (https://ecmc.state.co.us/data.html), we built and geocoded a dataset that contains the American Petroleum Institute well identification number, latitude, longitude, and monthly status (development, producing, shut-in, and abandoned) of all unconventional and conventional O&G wells (32) in Colorado between 1992 and 2019. From the years 2000 to 2019, we applied our IA-IDW model to estimate the monthly relative intensity of emissions at each O&G well site in Colorado. These relative intensities are based on measured phase-specific emissions of air pollutants, incorporating well pad information including the number of wells, the phase of development, and the amounts of O&G produced. A more in-depth description of this model has been published previously (21). Because the Colorado Oil and Gas Information System only had binary information on well activity (active and not active) prior to the year 2000, we estimated IA-IDW for each case and control prior to the year 2000 by developing a model based on O&G wells coming online between January 2000 and December 2001, as described in Supplementary Fig. S1.

To this dataset, we added the latitude and longitude coordinates for other air pollution sources in Colorado. We obtained information for the other air pollution sources from (i) the U.S. Environmental Protection Agency’s Toxics Release Inventory (TRI) program (https://www.epa.gov/toxics-release-inventory-tri-program) and Enforcement and Compliance History Online database (https://echo.epa.gov/facilities/facility-search), (ii) the U.S. Geological Survey National Mines Information Center, (iii) the CDPHE’s Concentrated Animal Feeding Operations, and (iv) CDPHE’s Composting, Solid Waste, and Wastewater Treatment Facility data. In addition to the latitude and longitude coordinates, we added the number of animal units for CAFOs, the site emissions for composting facilities, the flow (millions of gallons per day) for wastewater treatment facilities, and the annual reported air emissions for TRI facilities. Finally, we used ArcGIS Pro to calculate the average annual daily traffic within 152 m of each address on major roads and highways in 2015, according to the Colorado Department of Transportation.

We provided the final dataset, containing monthly estimates of O&G well activity and annual estimates of other air pollution activity, to CDPHE staff.

CDPHE staff determined the conception date based on the gestational age recorded on the birth certificate and then identified all O&G wells and other air pollution sources within a 10-mile (16 km) buffer of the maternal residence provided on the birth certificate for each month from conception through the diagnosis (case) or reference (control) date for each case and control. Based on the associations observed with adverse health outcomes in previous studies, a 16-km buffer represents a conservative geographic area of interest that could plausibly affect exposure (12, 33–36). It also allows for the evaluation of smaller buffers. CDPHE staff also computed distances between the maternal residence on the birth certificate and each O&G well and other air pollution source between conception and the diagnosis or reference date and returned to us a deidentified dataset.

We used our IA-IDW model (21) to estimate the monthly relative intensity of O&G well site activity in the 16-km buffer around the maternal residence of each case and control for each month from conception to the beginning of the latent period. We defined the induction and latent periods as described in Hicks and colleagues (37). Based on the induction period between exposure and benzene, known to be emitted from O&G well sites (17), and leukemia initiation occurring between 1 and 10 years (38–40), we defined the beginning of the latent period as 1 year prior to diagnosis. To evaluate if there was a specific developmental period when children were more susceptible to exposure, we then computed the total sum of monthly IA-IDW counts for six exposure windows: conception to the beginning of latency, birth to the beginning of latency, and conception to birth, as well as the first, second, and third trimesters of pregnancy.

Because the O&G wells included in our IA-IDW metric are weighted by the distance between the well and the residence, a well that is closer to the individual will contribute more to that individual’s metric than a well with the same intensity that is further away. Our IA-IDW metric differs from other methods that define an individual as exposed if they have a well within a given buffer without adjustment for the phase of well development or the intensity of operations that occur at the well site (13).

We assigned concentrated animal feeding operations an intensity of 1 to 5 based on the percentile of animal units, mining sources an intensity of 1 to 4 based on the type of mine, and TRI and solid waste and wastewater treatment facility sources an intensity of 1 to 5 based on the percentile of annual air emissions and waste flow, respectively. We assigned Composting and Enforcement and Compliance History sources an intensity of 1 because information on emissions and activity was not available.

We used unconditional logistic regression (41) to model the relationships between variables with case or control status. The sum of IA-IDW well count ≥1 was log-transformed and divided into quintiles for subsequent logistic regression (Supplementary Fig. S2). Each quintile was compared with the referent group of an IA-IDW well count of <1.

We then estimated the crude OR for the association between ALL and each quintile of summed IA-IDW well counts from conception to 1 year prior to the diagnosis (cases) or reference (controls) date. Because we observed points of departure in the OR for <20 and >80 quintiles and similar ORs between the 20 and 80 quintiles (Supplementary Table S1), we grouped IA-IDW well counts into three groups (low: <20 percentile, medium: 20–80 percentile, and high: >80 percentile) for subsequent statistical analysis to improve the robustness of the statistical analysis and avoid small cell counts.

We further investigated associations between ALL and IA-IDW groups from conception to the beginning of the latent period by adjusting for maternal and child covariates obtained from the birth registry, based on a priori knowledge of their association with both exposure and outcome. In our first adjusted analysis (model 1), we specifically considered maternal age greater than 35 years (yes, no; ref. 42), child’s biological sex (4), birth weight >4,000 g (yes, no; ref. 43), and birth address in a rural area (yes, no; ref. 12). We next added IA-IDW counts of other air pollution sources, UV exposure (44–46), and distance to the nearest highway <152.4 m (yes, no) to model 1 for our second analysis (model 2). Likelihood ratio tests indicated that IA-IDW groups had a statistically significant association with ALL, adjusting for confounders, at a significance level of 0.05. Evaluation of residual deviance between our crude and adjusted models indicated no overdispersion in the adjusted models. We also performed a stratified analysis by maternal and child covariates (sex, rurality, Hispanic ethnicity, and child’s age) to evaluate the potential for effect modification and investigated associations between ALL and IA-IDW groups between (i) conception and birth and (ii) birth and the beginning of the latent period, as well as IA-IDW groups for each trimester of pregnancy, as defined by the American College of Obstetricians and Gynecologists (47).

Because regulators typically consider the distance between an O&G well and a home (known as the setback) when establishing health-protective O&G well siting policies (48, 49), we compared the distance between the nearest O&G well site and the percentage of children in each IA-IDW group (low, medium, and high) for conception to latency, as well as the referent group (Fig. 2). For 100% of children in the high and medium IA-IDW groups, the nearest O&G well site to the maternal residence was within 5 and 13 km, respectively. Based on this evaluation and previous studies that observed increased ALL risk within 5 km of O&G wells and increased concentrations of air pollutants within 3 km of O&G wells (13, 50), we calculated summed IA-IDW between conception and the beginning of latency; determined the <20, 20 to 80, and >80 percentile exposure groups; and analyzed the association between ALL and IA-IDW for three additional buffers (3-, 5-, and 13-km) as described for the 16-km buffer. We also conducted tests to evaluate linear trends in binomial proportions with increasing IA-IDW by treating the categorical IA-IDW variable as ordinal and used the Wald χ² parameter to test for statistical significance (51).

Because data for several potential confounders were incomplete, we analyzed subsets of subjects with complete data for these covariates in the birth registry: no evidence of maternal smoking during pregnancy on the birth certificate (42), full-term birth (gestational age >36 weeks; refs. 52, 53), and socialeconomic status (SES) index (Supplementary Table S2) between 1 and 4 (4). We also analyzed subsets of children with an IA-IDW for other air pollution sources less than the 90th percentile, UV exposure less than the 90th percentile, and no highway within 500 feet to evaluate the effect of extreme exposure to other environmental stressors on our results. To evaluate potential exposure misclassification that could be introduced if the child changed residence over the induction period, we analyzed a subset of children whose mothers could be found in a commercially available database and who had the same address at conception and diagnosis/reference date (54). To evaluate whether our model for estimating well activity prior to the year 2000 affected our results, we analyzed the subset of children born after 1999.

We reported results for each outcome associated with the IA-IDW group (low, medium, and high) compared with the referent (IA-IDW <1 within each buffer) with 95% confidence intervals (CI). We considered the direction and magnitude of individual ORs, based on American Statistical Association guidance (55). All statistical analyses were conducted using SAS software version 9.5 (SAS Institute Inc.) or R Statistical Software (version 4.3.1; R Core Team 2023). The Colorado Multiple Institutional Review Board determined our study protocol to be exempt from human subject research per category 4 (COMIRB Protocol 21-3249).

The data that support the findings of this study are available from the Colorado Department of Public Health and Environment. Restrictions apply to the availability of these data, which were used under license for this study. Data are available from the author(s) with the permission of the Colorado Department of Public Health and Environment.

Mothers of children with ALL (cases) were more likely to be >35 years of age, were living in a zip code with a low SES score, were less likely to be Hispanic, have smoked during pregnancy, and were living in a rural area than mothers of controls. Cases were more likely to be male, were aged 3 to 4 years, have a birth weight <4,000 g, and were less likely than controls to have a highway within 152 m of their birth residence. Cases also had slightly lower UV exposure and slightly higher IA-IDW counts of other air pollution sources (Table 1).

For all windows of exposure and buffer sizes, a smaller proportion of cases than controls are in the referent group (IA-IDW <1). As the exposure window or buffer size decreases, the proportion of both cases and controls in the referent group increases, and the difference between the proportions of cases and controls in the referent group increases (Supplementary Tables S3 and S4).

Figure 3 presents associations between childhood ALL and IA-IDW from conception through the beginning of latency within the 16-km buffer. The unadjusted analysis suggests a 14% (OR = 1.14; 95% CI, 0.908–1.43) and 30% (OR = 1.30; 95% CI, 0.940–1.79) increase in ALL risk for the medium and high IA-IDW groups, respectively, compared with the referent group, as well as a 26% (OR = 0.742; 95% CI, 0.510–1.08) decrease in ALL risk for the low IA-IDW from conception to the beginning of latency group, compared with the referent group. After adjusting for maternal and child characteristics in model 1, we observed a 22% (OR = 1.22; 95% CI, 0.967–1.53) and 90% (OR = 1.90; 95% CI, 1.35–2.68) increase in ALL risk for the medium and high IA-IDW groups, respectively, compared with the referent group and a similar potential decrease (OR = 0.719; 95% CI, 0.493–1.05) in ALL risk for the low IA-IDW group compared with the referent group. Following further adjustment for environmental exposures in model 2, we observed a monotonic increase in risk (P for trend test <0.05) with a 32% (OR = 1.32; 95% CI, 1.04–1.69) and 118% (OR = 2.18; 95% CI, 1.51–3.15) increase in ALL risk for the medium and high IA-IDW groups, respectively, compared with the referent group and little change (OR = 0.756; 95% CI, 0.517–1.10) in ALL risk for the low exposure group.

Figure 4 presents associations between childhood ALL and IA-IDW within the 16-km buffer for the other five exposure windows, after adjustment for mother and child characteristics and coexposures (model 2). For each window, we observed a monotonic increase in risk (P for trend test <0.05). From birth to the beginning of latency, we observed a 31% (OR = 1.31; 95% CI, 1.02–1.67) and 127% (OR = 2.27; 95% CI, 1.56–3.27) increase in ALL risk for the medium and high IA-IDW groups, respectively, and possibly a decrease (OR = 0.821; 95% CI, 0.554–1.19) in ALL risk for the low exposure group, compared with the referent group. From conception to birth, we observed a 51% (OR = 1.31; 95% CI, 1.17–1.95) and a 127% (OR = 2.27; 95% CI, 1.52–3.34) increase in ALL risk for the medium and high IA-IDW O&G groups, respectively, and no difference (OR = 1.02; 95% CI, 0.690–1.49) in ALL risk for the low exposure group compared with the referent group. We observed similar increases in risk in the medium and high IA-IDW groups for each trimester of pregnancy. We observed a 20% to 23% potential increase in ALL risk for the low IA-IDW group in each trimester exposure window.

Figure 5 presents associations between childhood ALL and IA-IDW from conception to the beginning of latency within the 3-, 5-, and 13-km buffers, after adjustment for mother and child characteristics and coexposures (model 2). For each buffer, we observed a monotonic increase in risk (P for trend test <0.05). In the 3-km buffer, we observed a 64% (OR = 1.64; 95% CI, 0.866–2.90), 101% (OR = 2.01; 95% CI, 1.41–2.81), and 107% (OR = 2.02; 95% CI, 1.08–3.74) increase in ALL risk for the low, medium, and high IA-IDW groups, respectively, compared with the referent group. In the 5-km buffer, we observed a 62% (OR = 1.62; 95% CI, 0.964–2.62), 84% (OR = 1.84; 95% CI, 1.35–2.48), and 100% (OR = 2.00; 95% CI, 1.14–3.37) increase in ALL risk for the low, medium, and high IA-IDW groups, respectively, compared with the referent group. In the 13-km buffer, we observed a 59% (OR = 1.59; 95% CI, 1.03–2.37), 40% (OR = 1.40; 95% CI, 1.09–1.80), and 164% (OR = 2.64; 95% CI, 1.80–3.86) increase in ALL risk for the low, medium, and high IA-IDW groups, respectively, compared with the referent group.

We observed effect modification between rural and urban residences in the low IA-IDW group (Supplementary Fig. S3; Supplementary Table S4). In the low IA-IDW group, we observed a potential increased ALL risk (OR = 1.41; 95% CI, 0.503–3.95) for rural residences compared with a potential decreased ALL risk (OR = 0.679; 95% CI, 0.450–1.02) for urban residences. Although we did not observe effect modification by ethnicity, we did observe higher ALL risks in children of non-Hispanic mothers compared with children of Hispanic mothers (Supplementary Fig. S3). We also observed differences by the child’s age, with children aged 4 to 5 years experiencing greater ALL risk than younger and older children (Supplementary Fig. S4; Supplementary Table S4). We did not observe effect modification by the child’s biological sex (Supplementary Fig. S3).

In secondary analyses for children born at full term (≥37 weeks of gestation), born after the year 1999, with mothers who did not smoke during pregnancy, mothers who we could confirm did not move, an SES index between 2 and 4, an IA-IDW for other air pollution sources less than the 90th percentile, no major highway within 152 m of their birth residence, and UV exposure less than the 90th percentile, we observed results similar to the whole population (Supplementary Fig. S5; Supplementary Table S4).

In this study, we found that children aged 2 to 9 years diagnosed with ALL between 2002 and 2019 were 1.40 to 2.64 times more likely to live within 13 km of an O&G well site than children without a cancer diagnosis after adjusting for the child’s sex and birth weight, the mother’s age and rurality of address, IA-IDW of other sources of air pollution, UV exposure, and proximity to a major highway. Our results also indicate potentially elevated measures of association for exposures occurring during the perinatal period although these results are not conclusive. We found the burden of ALL risk to differ with the distance of the birth residence from the nearest O&G well site, as well as the intensity of O&G well site activity around the birth residence in regions with active upstream O&G development.

Previous case–control studies of children in Pennsylvania and rural Colorado also indicated positive associations between ALL and proximity to O&G well sites (12, 13). In this study, we observed lower associations than the 3 to 4 OR reported for children aged 5 to 24 years in our previous study in rural Colorado (12). The difference may be due to a more rigorous study design (registry-based vs. population-based), the inclusion of urban and suburban areas, and improved specificity of the exposure metric (IA-IDW vs. IDW), as well as adjustment for other environmental exposures in this study. Clark and colleagues (13) reported that children diagnosed with ALL between the ages of 2 and 7 years were 1.74 times more likely to have at least one unconventional O&G well within 2 km of their home within a window of three months prior to conception through to diagnosis, with results for buffer sizes of 5 and 10 km elevated but attenuating toward the null. Although this is within the 1.64 to 2.01 range we observed within our 3-km buffer, we did not observe an attenuation of association in our 5- and 13-km buffers. Because of the low prevalence of children living in their 2-km (2%), 5-km (7%), and 10-km (15%) buffers, Clark and colleagues did not apply an IDW measure of exposure or consider the intensity of activity at well sites. Instead, they considered a dichotomized exposure (e.g., presence or nonpresence of an O&G well within each buffer), which may have increased imprecision in exposure measurement as buffer size increased and biased results toward the null. Clark and colleagues did not consider buffers greater than 10 km.

Our results indicate that in regions with upstream O&G development, it is necessary to consider the density of wells and the level of activity on well sites in addition to the distance from the nearest O&G well in policies and epidemiologic studies. We observed that children living within 5 km of an O&G well and with O&G activity above the 80th percentile bear the greatest ALL risk. However, children living within 13 km of an O&G well and with O&G activity between the 20th and 80th percentiles also bear an increased risk of ALL. Although children living within 16 km of one O&G well and with O&G activity below the 20th percentile do not seem to experience an increased ALL risk, caution is warranted in interpretation because there are fewer than 50 cases in the low exposure group (Supplementary Table S4). One hundred percent of children in the high IA-IDW group for the 16-km buffer had a birth residence within 5 km of at least one O&G well site (Figs. 2 and 3) and an IA-IDW well count above the 80th percentile (2,308 and 132,254 intensity of well activity/km²). Ninety-nine percent of children in the medium IA-IDW group for the 16-km buffer (Supplementary Table S5) had a birth residence within 13 km of at least one O&G well site (Figs. 2 and 3) and an IA-IDW between the 20th and 80th percentiles (3.96 and 2,308 intensity of well activity/km²). One hundred percent of the children in the low IA-IDW group for the 16-km buffer had a birth residence within 16 km of at least one O&G well site (Fig. 2) and an IA-IDW ≤20th percentile. The 60% increase in effect size for the high exposure group observed in the 13 km may indicate that limiting analysis to smaller buffer sizes (e.g., 3 and 5 km) could lead to the inclusion of exposed children in the referent group and bias results toward the null.

The Colorado 2,000-foot setback (48) as well as the California 1-km setback, which is the largest setback between O&G well sites and residences in the United States (49), may not be sufficient to protect children from an increased burden of ALL risk. The results of Clark and colleagues (13) indicate increased ALL risk for children living up to 2 km from the nearest O&G well site, regardless of the intensity of O&G well site activities. Our results indicate increased ALL risk up to 13 km from the nearest O&G well site and provide clear evidence that the intensity of activity and the density of well sites around the birth residence, in addition to the distance to the nearest O&G well, affect the burden of ALL risk. For example, the nearest O&G well was further than 1 km from the birth residence for 33% and 95% of children in the high and medium IA-IDW groups, respectively, for the 16-km buffer (Fig. 2).

Although the etiologic pathway underlying an association between O&G activity and childhood ALL is unknown, the leukemogenic potential of air pollutants, such as benzene and PAHs, released in the areas of O&G activity is a plausible mechanism. For example, benzene metabolites, including t,t-muconaldehyde (aka E,E-muconaldehyde), benzene oxide, and benzoquinone, readily react with peptides, proteins, and DNA, thereby interfering with cellular function (56, 57). Benzoquinone formation in the bone marrow is a key step in leukemogenesis (56, 57). At ambient levels, benzene is preferentially metabolized through a pathway producing the cytotoxic metabolite, t,t-muconaldehyde (58–60). Human cells treated with t,t-muconaldehyde exhibit increases in reactive oxygen species (59), which can damage cellular nucleic acids (61, 62) and result in hematotoxicity and carcinogenesis (63, 64). In rodent cells, t,t-muconaldehyde is a potent inhibitor of gap junction channels, which has been linked to carcinogenesis and abnormal hematopoietic development (65, 66). Leukemogenic PAHs contained in diesel exhaust can lead to genetic toxicity by covalently binding to DNA (67). If not repaired, the resulting DNA adducts may be mutagenic and result in the disruption of the cell’s microenvironment, enzyme inhibition, cell death, and carcinogenesis (68). Although most ALL cases arise from B-cell somatic mutations (69), the somatic mutations to expect in ALL caused by benzene are unclear. Future research with phenotyping of ALL would advance the understanding of the etiologic pathway.

Children are particularly vulnerable to risks from elevated levels of leukemogens in their homes because high cell proliferation and rapid DNA turnover in early life provide a greater opportunity for DNA damage (65). Furthermore, young children spend a greater portion of their time at home (65).

Our rigorous population-based case–control study design and large study population are sufficient to detect differences in ALL incidence between cases and controls and establish temporality. Because all our data were obtained from the CCCR and Colorado’s birth registry, neither recall bias nor measurement bias affected our results. Although the possibility of random errors does exist in the registries, they are unlikely to affect our findings (70). We focused on ages of diagnosis (2–9 years) that are inclusive of ages of peak childhood ALL incidence and exclusive of ages for infant ALL (71). We reduced the potential for exposure misclassification by employing our validated spatiotemporal industrial activity model to improve the assessment of exposure that incorporates region-specific, data-driven weights based on empirical measures of emissions from O&G sites to quantify the relative intensity of air pollution emissions across four distinct O&G activity phases (i.e., construction, drilling, completions, and production; refs. 28, 30). Our model’s predicted O&G activity is strongly correlated with measured hydrocarbon concentrations over all phases of well development (r_Spearman = 0.74, P < 0.001) and is able to distinguish high-intensity phases, such as well completions, from low-intensity activities, such as production. Additionally, our model yielded a 19 times greater dynamic range of estimates compared with the simple IDW model, demonstrating a greater ability to identify a contrast in exposure to upstream O&G development among individuals.

Although we considered exposures to other sources of air pollution and potential leukemogens, as well as several characteristics that could be associated with childhood ALL in our analysis, we were not able to consider other potential confounders. However, our secondary analyses for children born full term (≥37 weeks of gestation) with mothers who did not smoke during pregnancy or an SES index between 1 and 4 did not indicate substantial confounding in our ALL results. Our inability to adjust for other potential confounders, including early common infections, nutrition, family history of neoplasms, trisomy and other genetic factors, water source, and daycare attendance, as well as individual income, may have resulted in residual confounding.

Using IA-IDW well counts in the absence of information on meteorology and topography likely reduced the temporal and spatial specificity in the IA-IDW metric. However, for the monthly IA-IDW well counts we used, a monthly dominant wind direction would not be sufficient, as the Rocky Mountains create upslope and downslope winds that are dependent on local topology, the time of year, the time of day, and other local conditions, with winds regularly shifting directions throughout the day (72, 73). Temperature inversions, rather than simply wind direction, are often related to high pollution days (74). The overall effect of the resulting exposure misclassification for not accounting for these many factors is unknown.

To reduce exposure misclassification from residential mobility, we considered using reconstructed residential histories from a commercially available database. We found that we were more likely to find a record in the commercial database, an address at conception, and an address at the reference date for mothers of cases than for mothers of controls (54). We also observed that the residential address at conception and the reference date differed from the address on the birth certificate for 29% and 45% of mothers found in the database, respectively. Because these differences between cases and controls could introduce systematic bias into our analysis and our secondary analysis indicated the potential for bias (Supplementary Fig. S4), we did not incorporate the reconstructed residential histories in our analyses and assumed that the address in the birth registry represented a child’s residence over the entire exposure period (conception to latency). Another assumption inherent in using a residential history reconstruction for the child’s birth mother is that the child resides at the same address as the birth mother throughout the latency period, which is not always true (65). However, we did conduct a secondary analysis of 658 children of mothers found in the database for whom we could confirm that the mother’s residential address did not change over the exposure period. In this secondary analysis, we observed a larger measure of association (OR = 3.71; 95% CI, 1.79–7.69) in the high IA-IDW group than we observed in our main analysis (Supplementary Table S6).

Secondary analysis for children born after 1999 indicates that less precise exposure estimates prior to the year 2000 may have introduced a slight bias. Our stratified analysis by child age also indicates the possibility of increased exposure misclassification of IA-IDW in the low exposure group for longer latency periods, as evidenced by the increasing OR in the low exposure group as the child’s age increased (Supplementary Fig. S4).

Our evaluation of the pregnancy and third trimester exposure window is limited by possible differences in the time at risk for exposure between cases and controls. Differences in time at risk for exposure between cases and controls may have biased our results for the pregnancy and third trimester exposure windows.

This study advances understanding of the relationship between distance to, density of, and level of activity on O&G well sites and childhood ALL. It provides evidence that children living up to 13 km from O&G well site activities are at an increased risk for ALL, with children living within 5 km bearing the greatest risk. These results indicate that current regulatory setbacks between O&G well sites and residences may not be sufficient to protect the health of children living in areas of upstream O&G activities. We emphasize the pressing need to continue comprehensive and rigorous research on the health consequences of early life exposures to upstream O&G activities. Future research should address the limitations in this and previous studies and focus on populations with greater ALL risk, such as Hispanic children and specific age groups.

L.M. McKenzie reports grants from the American Cancer Society and nonfinancial support from the Colorado Department of Public Health and Environment during the conduct of the study. W.B. Allshouse reports grants from the American Cancer Society during the conduct of the study. M. Cockburn reports grants from the University of Colorado during the conduct of the study. No disclosures were reported by the other authors.

L.M. McKenzie: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. W.B. Allshouse: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. D.R. Johnson: Resources, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. C.C. DeVoe: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. M. Cockburn: Conceptualization, resources, funding acquisition, validation, visualization, methodology, writing–review and editing. D. Ghosh: Conceptualization, resources, funding acquisition, validation, visualization, methodology, writing–review and editing.

This work was supported by an award from the American Cancer Society, grant number RSF-21-011-01-CPPB, to L.M. McKenzie. Any opinions, findings, and conclusions or recommendations expressed are those of the authors and do not necessarily reflect the views of the American Cancer Society. These data were collected and provided by the CCCR participating in the National Program of Cancer Registries of the Centers for Disease Control and Prevention and the Colorado Vital Statistics Program. The Centers for Disease Control and Prevention, CCCR, Colorado Vital Statistics Program, and the CDPHE specifically disclaim responsibility for any analyses, interpretations, or conclusions they have not provided. We thank Kirk Bol and John Arend for providing us with a deidentified dataset and helpful comments on the manuscript.