Mixture Effects of Commonly Applied Herbicides on County Level Obesity Rates in the United States: An Exploratory Ecologic Study (2013–2018)

TL;DR

This study explores how mixtures of herbicides may be linked to higher obesity rates in U.S. counties, especially in rural areas.

Contribution

The study is novel in examining the combined effects of multiple herbicides on population-level obesity rates.

Findings

Glyphosate and other herbicides showed significant positive associations with county-level obesity rates.

Herbicide mixtures were strongly linked to higher obesity rates, with rural counties showing stronger associations.

The study highlights the importance of considering cumulative herbicide exposure rather than individual chemicals.

Abstract

Metabolic disorders such as obesity have increased globally in recent decades and are a major public health concern. Previous research suggests that herbicide exposures may contribute to metabolic dysfunction, but few studies have examined mixture effects of multiple herbicides on obesity at a population level. Using county-level data from 2013 to 2018, we examined the associations between obesity rates and the application of 13 commonly applied herbicides in the U.S. We first conducted adjusted single-pollutant mixed-effects models and then used quantile-based g-computation mixture modeling to assess combined herbicide mixture effects on county-level obesity rates. Models were adjusted for demographic and socioeconomic covariates and accounted for geographic clustering. Significant positive associations were identified between county-level obesity rates and applications of glyphosate, 2,4-D, atrazine, acetochlor, metolachlor, and several other herbicides in adjusted single-pollutant models. Glyphosate showed one of the strongest individual associations (β = 0.29 per standard deviation increase, 95% CI: 0.21–0.36). Increases in herbicide mixture were significantly associated with higher obesity rates (Psi = 0.71 per quantile exposure mixture, 95% CI: 0.65–0.76) from mixture modeling. Inclusion of significant interaction terms did not appreciably increase the mixture effect. Glyphosate, 2,4-D, metolachlor, dimethenamid-P, and glufosinate contributed most strongly to the weighted mixture effect. Mixture effects varied by rurality, with stronger associations observed in rural counties, particularly in micropolitan regions. Our findings highlight the importance of considering cumulative herbicide mixture exposures rather than individual chemicals in isolation. The observed rural–urban disparities emphasize the need for targeted public health interventions and policy actions in rural communities, which may be particularly vulnerable to the adverse metabolic impacts of herbicide mixtures.

Figures

Full text

1. Introduction

Recent decades have seen a substantial increase in human metabolic disorders—including obesity, diabetes, and metabolic syndrome—which now represent a major public health challenge in the United States and globally [1,2]. Agricultural practices have evolved in tandem with the widespread application of pesticides. Of these, herbicides are by far the most heavily applied class of pesticides, accounting for ~88% of agricultural herbicide mass during 2013–2017 in the United States, whereas insecticides and fungicides each contributed ~6%. Moreover, by 2016 glyphosate alone comprised ~44% of all herbicide mass applied to crops [3,4]. Our prior study demonstrated a significant association between glyphosate exposure, measured in urine, and metabolic syndrome using nationally representative NHANES data, with marked differences observed across racial groups [5]. These findings suggest that herbicide exposure may be an underappreciated environmental co-factor of metabolic dysfunction. However, given that agricultural herbicide exposures rarely occur in isolation, a comprehensive analysis that includes multiple compounds is warranted to elucidate interactive and mixture effects on metabolic health.

In the present study, we extend our investigation with an ecologic study of herbicide application and its associations with county-level obesity rates by incorporating a panel of thirteen herbicides—glyphosate, herbicide 2,4-D (2,4-D), atrazine, dicamba, trifluralin, acetochlor, dimethenamid-P, glufosinate, metolachlor, metolachlor + metolachlor-S, pendimethalin, paraquat—and one fungicide—chlorothalonil—that are most extensively applied in U.S. agriculture and have been the subject of extensive toxicologic and epidemiologic evaluation due to their widespread application and potential health and environmental impacts [6]. In addition glyphosate applications have increased dramatically since 1996 (15-fold by 2016) following the development of “Roundup^TM^-ready” crops [7]. Enthusiastic adoption of “Roundup^TM^-ready” crops has resulted in increased applications of dicamba and herbicide 2-4D to manage herbicide-resistant weeds. The selection of the fourteen compounds studied here was informed by multiple criteria: (1) their high prevalence in agricultural use in the United States, and documented potential for environmental persistence; (2) emerging toxicological evidence implicating them in endocrine disruption, oxidative stress, and perturbations in glucose and lipid metabolism [8,9,10]; and (3) epidemiological observations that link these exposures to adverse metabolic outcomes [11]. Importantly, many of these herbicides are applied together in modern cropping systems, potentially leading to synergistic and mixture effects that are not adequately captured by single-pollutant studies.

The biology linking these compounds to metabolic disturbances is supported by several mechanistic pathways. For instance, glyphosate has been implicated in disrupted insulin signaling and altering the gut microbiome—both of which are critical regulators of energy homeostasis and metabolic function [12,13,14,15]. Similarly, atrazine and 2,4-D have been associated with endocrine disruption, which may interfere with hormonal regulation of adipogenesis and glucose metabolism [10,16]. The neurotoxic potential of paraquat [17,18] further raises concerns about systemic oxidative stress, a contributor to insulin resistance and beta-cell dysfunction [19,20]. In integrating diverse lines of evidence, our study seeks to clarify whether cumulative exposures to herbicides contribute to metabolic disorders at a population level.

After assessing associations between county-level obesity rates and single pollutants using mixed-effects models with a random term for county FIPS, we use a mixture model approach to assess mixture effects, using a county-level random effect to account for location. This approach enables us to accommodate geographic variations in herbicide usage and to control for spatially specific local demographic and socioeconomic factors that may modify exposure–outcome associations. Finally, by stratifying analyses on National Center for Health Statistics (NCHS) urban–rural designation, we aim to reveal potential environmental health disparities that could inform targeted public health interventions.

Mixture effects. Mixture effects occur when multiple environmental pollutants interact in ways that amplify their individual health effects, and may be more common when compounds are in related families [21,22]. Pollutant mixtures may have joint effects [23], even at low doses [24,25], and if there are interactions between components the total mixture effect can be amplified [26]. As pollutant mixture effects are a long-recognized source of health risk, methods to assess the effects of pollutant mixtures, primarily for air pollutant mixtures, were in use as early as 2000 [27]. Despite some focus on acute herbicide exposure (poisoning), the effects of chronic ambient environmental exposure herbicide mixtures are little studied in view of increasing global application of herbicide mixtures for everything from weed control in cropping to crop desiccation, aquatic or wildland weed control, and urban area management [28,29].

In summary, our study addresses critical gaps in the current understanding of environmental risk factors for metabolic disorders. By examining a spectrum of herbicides with documented endocrine and metabolic perturbations, this study provides a nuanced assessment of the complex interplay between agricultural chemical exposures and metabolic health outcomes as a function of exposure setting and magnitude. Such insights are vital for informing regulatory policies and for designing future studies that can elucidate the causal pathways underlying these associations.

2. Materials and Methods

2.1. Study Design

Our study uses aggregated county-level exposure and outcome data to explore relationships between herbicide exposure and health outcomes. Because both exposure and outcome comprise county-level measurements, there is an inherent study assumption that agricultural herbicide applications are correlated with local population exposure. This is biologically and empirically plausible because urinary glyphosate concentrations rise during the agricultural spray season among residents living near fields [30,31], indicating that local application intensity tracks human exposure. Diet may contribute, but provides a smaller contribution to exposure for subjects living in proximity [32]. Moreover, our prior NHANES study demonstrated an association between urinary glyphosate measurements and metabolic syndrome score, increasing confidence that population-level exposure can serve as a meaningful proxy for human exposure in this context [5]. The U.S. Geological Survey (USGS) [33] provides estimates of agricultural pesticide and herbicide application in kilograms, aggregated by FIPS code and year, which we employ alongside county-level obesity rates from the Behavioral Risk Factor Surveillance System (BRFSS) and demographic information from the American Community Survey (ACS) for study years 2013 to 2018. Sources are described further below.

2.2. Data Sources, Curation, and Integration

We curated a county–year analytic file by integrating multiple public, nationally representative sources and joining records by five-digit FIPS and calendar year. Agricultural herbicide use estimates were obtained from the U.S. Geological Survey (USGS) Pesticide National Synthesis Project (e-Pest), for which peer-reviewed evaluations of data generation and occurrence are available [3,4]. County-level adult obesity prevalence was drawn from CDC small-area estimation based on BRFSS (PLACES), which uses multilevel regression and post-stratification and has been validated for local health indicators [25,26]. Socio-demographic and behavioral covariates (e.g., smoking, uninsured, unemployment, education, age structure, race/ethnicity) were curated from the University of Wisconsin’s County Health Rankings, which aggregate the American Community Survey and other federal series under a documented framework [34,35]. This approach mirrors recent ecologic studies that explicitly curate and integrate PLACES and other public datasets for population-level analyses [36,37,38].

2.3. Designation of Rurality

Rural–urban setting designations are from the 2013 National Center for Health Statistics (NCHS) Urban-Rural Classification Scheme for counties and include six categories: “large central metro”, “large fringe metro”, “medium metro”, “small metro”, “micropolis” and “non-core.” Details for each category are shown in Table 1. The 2013 NCHS scheme is based on the 2010 census and the February 2013 Office of Management and Budget (OMB) delineation of metropolitan statistical areas (MSAs) and micropolitan statistical areas, where a micropolitan statistical area is defined as a population of 10,000 to 49,999. There were 374 metropolitan statistical areas and 581 micropolitan statistical areas defined as of 2013. These designations were recently updated to a 2023 standard using 2022 census data, but our data comprise years 2013 to 2018; thus, we employ the 2013 system. Note that there were minimal changes to county assignments between 2013 and 2023.

2.4. Exposure

Primary exposure variables include the estimated annual agricultural use of atrazine, dicamba, trifluralin, acetochlor, chlorothalonil, dimethenamid-P, glufosinate, glyphosate, herbicide 2,4-D, metolachlor, metolachlor-S, metolachlor + metolachlor-S, pendimethalin, and paraquat. Herbicides were selected based on application levels across U.S. counties as follows: Exposure estimates were downloaded from the USGS e-pest (Pesticide National Synthesis Project [33]) site for years 2013 through 2018. E-pest estimates are reported as “low” and “high” values (in kilograms) for 448 unique compounds over the time period for most states. California supplies only one comprehensive number per county, taking data directly from its comprehensive Pesticide Use Reporting (PUR) database. For the remaining states, we computed averages of low and high values (in kilograms) for each exposure, then omitted exposures with a total mean value of less than 1400 kg/hectare, irrespective of missing units, resulting in 43 candidate exposures. Columns with more than 30% missing units were dropped, resulting in the final 14 exposures listed above. Remaining missingness varied from <1% to at most 20%. State and county FIPS numbers were joined to create five-digit FIPS values for each county, and exposures were merged with county-level outcome and demographic data based on study year and county FIPS number.

2.5. Outcome

Adult (ages 18 and older) obesity data are compiled by the CDC from the Behavioral Risk Factor Surveillance System (BRFSS). The BRFSS is an annual state-level random digit dial survey used to assess health and risk-related behaviors. From 2016 onwards, the CDC employed a multilevel modeling approach to estimate obesity, along with other health conditions, based on telephone survey responses and respondent age, sex, and race/ethnicity, combined with county-level poverty and other relevant county- and state-level features [39,40]. For counties where there is insufficient data, the approach borrows data from the entire BRFSS sample as well as old census estimates, using a parametric bootstrap to produce standard errors. A companion study will examine associations between county-level obesity, hypertension and hypercholesterolemia rates (%), and herbicide exposures within larger Metropolitan Statistical Areas (MSAs), using data from CDC Places for 2013 to 2018. The CDC uses a similar multilevel modeling approach to estimate CDC Places health outcomes.

2.6. Covariates

Additional demographic and socioeconomic variables—including adult smoking prevalence, percent uninsured, percent unemployed, age structure (≥65 and <18 years), educational attainment (percent high school graduate; percent with some college), racial/ethnic composition, and county rurality (NCHS 2013)—were included a priori to address confounding. Percent uninsured comprises age-adjusted prevalence of subjects aged 18 to 64 who lack health insurance; percent smoking reflects the percentage of the county population who currently smoke and have smoked more than 100 cigarettes in their lifetime. Smoking was treated as a key confounder because rural smoking prevalence exceeds urban prevalence, including in New York State and even after accounting for poverty, and because smoking is associated with adiposity/metabolic risk behaviors. To partially address dietary confounding—for which no county-level intake data were available—we included food insecurity and educational attainment as proxies for diet quality/access and health behaviors; recent national evidence shows poorer diet quality in non-metropolitan/rural areas independent of income/education and food-desert status [41]. These county-level metrics are compiled by study year and made available for download by the University of Wisconsin Population Health Institute (https://www.countyhealthrankings.org, accessed on 24 September 2024) and are also available at the Centers for Disease Control.

2.7. Statistical Analysis

Due to the distributed nature of missing units across county-level herbicide estimates, missingness decreases our sample size from 18,382 samples to 7811 samples if we take a complete case analysis approach. To preserve sample size, data were thus first multiply imputed using a fully conditional chained equations approach implemented in R via the mice package, creating 10 imputed datasets, each given 20 iterations to allow convergence, and using a custom prediction matrix. Please see Appendix A for details of the imputation process and model.

2.8. Single Pollutant Models

After imputation and before multivariable analysis, we assessed potential collinearity between herbicide application data by computing and visualizing a Pearson correlogram (Appendix A, Figure A1), as multicollinearity will inform further modeling. We then examined single-pollutant associations using mixed-effects models adjusted for racial percentages, percentage age 65 and above, and percentage below age 18; county-level smoking rates; percentage uninsured, unemployed, and food-insecure; and NCHS rurality designation, with a random effect to account for county-level variability. Exposures were mean-standardized before analysis and the results are reported as per standard deviation exposure, and results were pooled over 10 imputed datasets using Rubin’s Rules [42]. p-values were adjusted for multiple comparisons via the False Discovery Rate (FDR) method [43]. Results are given in tabular format and visualized as forest plots. A second set of adjusted single-pollutant regressions with the exposure categorized by quintile was run to assess potential for nonlinear dose–response relationship.

2.9. Mixture Models

Our data showed correlation between the herbicides; thus, a standard multivariable model was inappropriate due to potential multicollinearity and associated variance inflation, unstable coefficient estimates, and difficulty assessing exposure importance. Correlation analysis warns us of potential multicollinearity in our models yet cannot detect which correlated exposures are driving the associations, or whether there are significant interactions between the pollutants; we rely upon other methods to assess this. As the size of our dataset makes application of flexible nonparametric kernel-based methods such as Bayesian Kernel Machine Regression (BKMR) [44,45] computationally impractical, we used quantile-based g-computation models to estimate the joint effects of herbicide applications at the county level on county-level obesity rates [46]. Based on our correlation analysis (see Appendix A, Figure A1), we dropped Metolachlor-S from the mixture components. We first ran an unadjusted model to examine crude mixture effects, and then a model adjusted for racial percentages, percentage age 65 and above, and percentage below age 18; county-level smoking rates; and percentages uninsured, unemployed, and food-insecure, using a grouping term to account for clustering at the county level and a bootstrap approach to correctly estimate variances. Optimal quantile settings were selected by examining Z-scores and were set at 15 quantiles for both crude and adjusted models. Quantile-based g-computation models were run using the gcomp package in the R programming language [46].

2.10. Interactions

Potential interactions between herbicide exposures were assessed for inclusion in g-computation models by running interaction forests, a variant of the random forest (RF) method that explicitly captures interaction effects in the bivariable splits performed by the decision trees in RF, using the diversity Forest package in R [47].

2.11. Stratified Mixture Models/Subgroups Analysis

Since application of herbicides in farming is associated with rurality, we perform mixture modeling stratified on NCHS Urban-Rural designation.

3. Results

3.1. Study Population

The study population resides in 3066 counties spread across the United States (Table 2). Only 13.5% of the counties are in large metropolitan or large fringe metropolitan areas. In total, 63% of counties are in either micropolitan or non-core areas. Our study population is largely (median 86.3) non-Hispanic White and has graduated high school (86%). Over half have attended college (56%). County-level smoking rates are about 19%, most people are employed, and only 16% are without health insurance, with only 14% food insecurity. County-level obesity rates across the study period have mean and median values of 31% and rates are normally distributed.

3.2. Herbicide Exposure and Health Outcomes by Rurality

Herbicides of interest in this study had median annual application levels ranging from 128 kg per hectare (chlorothanonil) to 15,360 kg per hectare (glyphosate). Maximum annual applications ranged from 70,608 kg (metolachlor) to 594,336 (glyphosate) kg per hectare (Appendix A, Table A1). Annual applications at the county level generally increased with increasing rurality, flattening out as rurality increased above 45%. Counties within micropolitan urban–rural designations (52.3% median rurality) had the highest median annual applications of metolachlor species, herbicide 2-4D, atrazine, acetochlor, dimenthenamid-P, and glufosinate, and the second highest applications of glyphosate (Appendix A, Table A2). Large metropolitan areas and metropolitan fringe areas had the lowest application levels. Mean/median obesity rates also increase with rurality, again flattening out at about 31% at rurality 45% and above (NCHS small metro, micropolitan and non-core regions).

3.3. Single Pollutant Models

Across all counties, the adjusted analysis revealed positive single-pollutant associations between glyphosate, metolachlor and metolachlor-S, acetochlor, atrazine, pendimethalin, herbicide 2-4D, Glufosinate, and Dimenthenamid-P applications (in decreasing magnitude), respectively, and county-level obesity rates (Figure 1, Table 3). Associations for glyphosate were strongest (+0.30 per SD exposure, 95% CI: 0.2, 0.4), followed by metolachlor + metolachlor-S (+0.24, 95%CI: 0.2, 0.3) and eight other compounds. A total of 10 of 14 county-level herbicide applications showed positive and significant associations with obesity rates, though those for dimenthenamid-P were marginal after correction for multiple comparisons. Trifluralin, paraquat, dicamba, and chlorothanonil applications show no significant adjusted linear associations with county-level obesity rates. Quantile regressions with exposure categorized by quintile show little evidence of a nonlinear dose–response relationship in most herbicides. Results categorized by quintile show that 13 of 14 herbicides are significantly positive and monotonically increasing, and all are significant for at least two of the four reported quintiles (Appendix A, Table A2). Again, the largest positive associations are shown for glyphosate, followed by metolachlor + metolachlor-S and 11 other herbicides.

3.4. Mixture Models

In unadjusted models, a one-quantile increase in herbicide application mixture was associated with a significant increase in obesity rates ( $[eqn]$ = 1.02, 95% CI: 0.94, 1.1; p < 0.0001, Table 4). Adjusted models showed slightly reduced but still significant mixture effects with a one-quantile increase in obesity rates ( $[eqn]$ = 0.71, 95%CI: 0.65, 0.76, p < 0.0001). Model weightings for the mixture effects showed large positive contributions from herbicide 2-4D, metolachlor, glyphosate, dimenthenamid-P, and glufosinate, and large negative contributions from paraquat, trifluralin, and pendimethalin (Figure 2). Here, positive weights indicate positive additive effects on obesity rate, whereas negative weights are negative contributions for a given component in the per-quantile mixture effect. For a forest plot of coefficients (beta values) from the underlying fitted mixture model, please see Appendix A, Figure A2. Surprisingly, the interactions assessed via interaction forests are few. The largest, a quantitative interaction between dicamba and glyphosate, had an Effect Importance Measure (EIM) strength of 0.06 (Appendix A, Figure A3). Adjusted mixture models incorporating these interactions, however, do not produce significantly different total mixture effects (Table 4).

3.5. Mixture Analysis Stratified by Rurality

As in unstratified models, mixture-only models (estimates and confidence intervals are shown in pale blue in Figure 3) generally had higher mixture effects per quantile exposure mixture than adjusted models. Large central metros are an exception, but the sample size is also much smaller (comprising fewer counties) than other levels. Mixture effects from adjusted models rise with increasing rurality from large fringe metropolitan (Psi ( $[eqn]$ ) = 0.6, 95%CI: 0.4, 0.8) to micropolitan areas (Psi ( $[eqn]$ ) = 0.9, 95%CI: 0.7, 1.0, see Appendix A, Table A3). Surprisingly, non-core regions have among the lowest mixture effects, though the intercepts for this area show the highest rates of baseline obesity (Table A4). Note however that median exposures for most herbicides peak in small metropolitan and micropolitan districts during the period of interest (see Figure 4 and Appendix A, Table A4).

4. Discussion

Using freely available public data on county-level obesity and application of the 14 most commonly used herbicides in a five-year period (2013–2018), we have found significant adjusted associations between county-level obesity rates, individual herbicides, and their mixture. We have shown that these associations roughly mirror pesticide applications, and that these are in part, but not fully, associated with rurality and setting.

Our analysis demonstrates that U.S. counties with greater herbicide use, including glyphosate and several other common agents, tend to exhibit higher obesity rates, consistent with prior evidence that herbicide exposure contributes to adverse metabolic outcomes [48,49]. We observed significant positive associations from single-pollutant models between county-level obesity prevalence and the application of multiple herbicides (10 of the 14 analyzed in standardized continuous data, 12 of 14 in quantile categorized), most notably glyphosate, acetochlor, metolachlor -S, metolachlor (including its S-isomer), atrazine, metolachlor, pendimethalin, 2,4-D, glufosinate, and dimenthenamid-P (Figure 1), as well as some evidence of a nonlinear dose–response relationship. Glyphosate, as the most heavily used herbicide (Figure 4), showed the strongest individual associations with obesity, followed by 2,4-D and atrazine, consistent with human studies linking these herbicides to metabolic syndrome, of which obesity is a key component [5,50,51,52,53,54]. Experimental studies provide plausible biological mechanisms by which these herbicides may influence metabolic outcomes like obesity. Glyphosate has been shown to alter gut microbiota composition, increase oxidative stress, and up-regulate expression of the NF-κB (Nuclear Factor kappa-light-chain-enhancer of activated B cells), all of which can impair insulin signaling and promote adipogenesis [14,55,56]. Atrazine and 2,4-D have been linked to endocrine disruption through interference with estrogenic and thyroid pathways, activation of PPARγ—a transcription factor central to fat storage and adipocyte differentiation—and mitochondrial dysfunction, leading to insulin resistance [50,57]. These mechanistic insights support the plausibility that chronic exposure to herbicide mixtures could dysregulate energy balance, lipid metabolism, and glucose homeostasis, thereby contributing to obesity risk.

Not all compounds had significant linear independent effects in continuous data (e.g., dicamba, paraquat, trifluralin, and chlorothalonil showed no positive single-pollutant association with obesity in continuous data). However, most showed partial positive effects (usually the highest levels, Q4 and Q5 relative to Q1) in quintile categorized data, highlighting that the health impacts of herbicides are not uniform across compounds, and that slight nonlinearities may mask significant associations when examined in continuous data. This heterogeneity likely reflects different toxicological mechanisms and exposure patterns which may influence metabolic outcomes in subtle and indirect ways.

Importantly, our mixture modeling analysis suggests that combined exposures to herbicides may better explain population obesity patterns than individual chemicals alone. Using quantile-based g-computation, we found that a one-quantile increase in the overall herbicide mixture was associated with a significant increase in county obesity rates, even after adjusting for socioeconomic covariates. In the adjusted mixture model, herbicides such as 2,4-D, glyphosate, metolachlor, glyphosate, dimethenamid-P, glufosinate, metolachlor (including its S-isomer), dicamba, acetochlor, and atrazine showed the largest positive weight contributions to this mixture effect, whereas others (paraquat, trifluralin, pendimethalin, and chlorothalonil) contributed negative weights. Notably, some herbicides that showed null or modest single-pollutant associations (e.g., paraquat) exhibited negative contributions in the mixture context, possibly reflecting complex interactions (see Appendix A, Figure A3) or spatial usage patterns that are not fully accounted for in the mixture model (for instance, areas with high paraquat use may have lower use of other obesogenic herbicides, diluting the overall mixture effect). On the other hand, glyphosate—while strongly associated on its own—did not singularly dominate the mixture effect, as several co-occurring herbicides also exerted substantial influence. These results underscore that glyphosate is not uniquely responsible for the observed obesity link; rather, it is one significant part of a total herbicide mixture effect. In practical terms, evaluating herbicides in isolation may bias risk estimates, whereas accounting for concurrent exposure mixtures provides a more realistic assessment of environmental influences on obesity. Our study is among the first ecological analyses to examine herbicide mixtures in relation to metabolic health, and it provides novel evidence that the cumulative effects of multiple agrichemicals are relevant to population obesity rates.

4.1. Contextualizing Rural–Urban Disparities

A notable finding from our study is the heightened impact of herbicide exposure in more rural areas. herbicide usage tends to be highest in counties of intermediate-to-high rurality (e.g., small metropolitan and micropolitan areas) and lowest in major urban centers (Figure 4). Consistently, we observed that the overall herbicide mixture effect on obesity was more pronounced in more rural counties: the estimated obesity increase per mixture quartile was smallest in large metropolitan areas and grew larger moving toward small metro and micropolitan counties. Interestingly, mixture effects in the most remote rural counties (“non-core” areas) did not follow a strictly linear trend—they showed the highest baseline obesity rates but a somewhat lower mixture effect estimate—suggesting that beyond a certain point, additional exposures may yield diminishing returns, or that these communities may face saturated obesity risk from many non-measured contributors. Overall, however, the pattern indicates that communities with intensive agriculture (often rural) may experience a double burden: greater chemical exposure alongside underlying vulnerabilities such as higher poverty, limited health care access, and other lifestyle risk factors [58]. These contextual factors can amplify the health impact of environmental exposures. Indeed, rural populations often have fewer resources to mitigate or treat chronic conditions, which could exacerbate the observed associations. Our findings therefore support calls for region-specific public health interventions and regulatory approaches. In practice, this could mean prioritizing cumulative risk assessments for heavily agricultural rural regions and tailoring obesity prevention programs to address both lifestyle and environmental factors in these communities. By recognizing that rurality can be a proxy for both elevated exposure and increased susceptibility, public health officials can better target efforts to reduce herbicide exposure and bolster health care support in high-risk counties. Recent work has underscored how socioeconomic and infrastructural disparities can modify the health impacts of environmental exposures [58].

4.2. Confounding by Smoking and Diet

We adjusted for county-level adult smoking prevalence because rural smoking remains higher than urban smoking in the U.S. after adjustment for socio-demographics, showing substantially higher smoking in rural/upstate counties even where poverty is high [59,60]. We adjusted for diet proxies, specifically food insecurity and educational attainment, because both are well-established county-level predictors of diet quality, and because counties with higher food insecurity and lower education consistently show poorer diet quality, including lower fiber intake, higher added sugar intake, and worse overall Healthy Eating Index scores [41,61,62]. Furthermore, although we could not directly account for food consumption patterns, national studies confirm that rural populations tend to consume ultra-processed diets, with less fiber and more added sugar, which could contribute to higher obesity prevalence independent of herbicide exposure [63,64]. Therefore, our adjustment of selected dietary proxies help reduce confounding by diet, which may be systematically poorer in rural areas, although cannot fully address the lack of direct county-level measures of diet and physical activity. Other exposures common in agricultural settings, including the use of other types of pesticides (e.g., insecticides, rodenticides, etc.), air pollution, and heavy metals, may also confound or interact with herbicide–obesity associations, and their omission represents an additional limitation.

4.3. New Dimensions in Exposure Assessment

Our study contributes methodologically by leveraging publicly available county-level herbicide application data to approximate population exposures. This ecological approach enabled us to probe potential exposure–response relationships on a national scale. We observed, for example, that counties with the highest herbicide application levels tended to have some of the highest obesity prevalence, with obesity rates plateauing around ~31% in the most agriculturally intensive areas. Although based on observational correlations, this pattern raises the hypothesis of threshold effects whereby metabolic health risks may accelerate once herbicide use (and by inference, community exposure) exceeds a certain level. Our analysis, which captured geographic variability and multiple chemicals simultaneously, provides a more comprehensive real-world exposure scenario than single-chemical studies. By accounting for dozens of U.S. states and a spectrum of herbicide compounds, we capture the complex environmental conditions under which human populations actually live. Furthermore, our ecological approach adds a new dimension to exposure science by marrying large-scale data with mixture modeling to yield insights that are not apparent when studying one chemical or one location at a time. This integrated exposure assessment is a significant step forward in bridging the gap between epidemiological observations and mechanistic toxicology.

4.4. Strengths, Limitations, and Future Directions

Strengths of our study include that we analyzed a large, nationally representative dataset covering six years (2013–2018) and over 3000 counties, which provides ample statistical power and broad generalizability. We integrated high-quality data from federal sources—including USGS agricultural herbicide use estimates and CDC-modeled county obesity prevalence—and adjusted for a range of demographic and socioeconomic covariates to reduce confounding. To address missing data, we employed multiple imputation methods customized to the dataset, maximizing the use of available information and limiting bias from incomplete records. Notably, we applied a rigorous mixture modeling approach (quantile g-computation) to estimate the joint effect of 13 herbicides on obesity. This method, developed specifically for epidemiologic mixtures, allowed us to calculate an overall effect estimate (Psi) for the herbicide blend while yielding weights for each component. The quantile-based approach is robust to distributional extremes and multicollinearity, enhancing our confidence that the observed mixture effect is not an artifact of one highly prevalent chemical. Taken together, the study design and analytical techniques provide a robust triangulation of evidence at the population level, complementing prior individual-level study.

Our study has several limitations. First, the study is ecologic and analysis was performed on data that is aggregated at county level, thus results cannot be interpreted as causal effects at the individual level, and indeed, spatiotemporal aggregation may obscure finer-scale relationships. However, a gold-standard randomized controlled trial (RCT) for these exposures and outcomes would establish causality at an individual level, but would be unethical. And, while a longitudinal study would strengthen causality, one of this size would be prohibitively expensive and, if observational, would not establish causality on its own. Thus, while it has stated limitations, our study leverages publicly available datasets to explore associations between a ubiquitous set of exposures and a chronic health outcome on national scales. Our study is based on an assumption that county-level herbicide application is correlated with local human exposure, and results may be influenced by unmeasured county-level factors such as diet, physical activity, or other correlated exposures. We address this by including county-level random effects, adjusting for rurality, and controlling for multiple socio-demographic variables. Aggregation at the county level may introduce a modifiable area unit problem. Our herbicide use metric is a surrogate for human exposure and does not account for individual behaviors or chemical drift/dynamics, which may lead to non-differential exposure misclassification and bias estimates toward the null. County obesity prevalence was obtained from a modeling method (BRFSS small-area estimates) that has its own uncertainty. However, by including year as a fixed effect and county random effects we attempted to account for differences associated with time and spatially associated variability. Fourth, while quantile g-computation is a powerful tool for mixtures, it computes additive effects of increasing all exposures by one quantile and de facto assumes linearity unless nonlinear terms are explicitly incorporated into a model. We found minimal evidence of pairwise interactions among herbicides in supplementary analyses, suggesting the additivity assumption was reasonable. Still, very high correlation between certain herbicides (e.g., metolachlor and metolachlor-S) required us to drop one variable to avoid collinearity, highlighting a general challenge in multi-pollutant studies. Our current dataset (county-level CDC obesity data with additional county-level covariates from the ACS) comprises population rates by county and thus does not allow stratification by sex, and does not include occupational proportions. Therefore, we were unable to stratify by sex or by rural occupational status, other than including rate of unemployment by county. Future studies with individual-level data could address these important sources of heterogeneity.

Finally, our focus was on obesity as outcome; we did not examine other metabolic outcomes (such as diabetes, hypercholesterolemia, or hypertension rates) in this analysis. It remains possible that herbicide mixtures impact various metabolic health indicators differently, an aspect a future study will explore. In view of these limitations, our study should be considered exploratory.

We recommend several directions for further research. Controlled longitudinal studies—for example, following birth cohorts or agricultural communities over time—are needed to test the temporality of herbicide exposure and obesity onset, which our ecologic design cannot establish. Incorporating individual-level exposure data (e.g., urinary or blood biomarkers of herbicides) would greatly strengthen causal inference by reducing exposure misclassification and allowing dose–response relationship assessment. Indeed, one recent longitudinal study of young adults (the CHAMACOS cohort) reported that cumulative glyphosate exposure was associated with elevated metabolic syndrome risk, illustrating the value of detailed exposure tracking [51]. Future studies might also consider experimental and mechanistic investigations, such as animal or cellular models of herbicide mixtures, to unravel how these chemicals might jointly disrupt metabolic regulation. For instance, do low-dose combinations of glyphosate, 2,4-D, and atrazine induce greater adipogenesis or insulin resistance in vivo than each compound alone? Questions like this remain unanswered, and toxicological research could elucidate potential synergistic or antagonistic interactions at the molecular level. Additionally, examining spatial patterns using finer geographic resolutions (e.g., census tract data, as in recent environmental determinant studies) could help identify localized “hotspots” of metabolic disease tied to herbicide use, thereby refining intervention targets.

5. Conclusions

In summary, this ecological analysis adds to our understanding of environmental influences on metabolic health. While it cannot prove causation, the alignment of our population-level findings with evidence from individual-based studies [5,15,51,52] strengthens inference that chronic herbicide exposures may be contributing to chronic obesity. Our study offers a broad spatiotemporal perspective on this issue, suggesting that areas of heavy herbicide application are associated with obesity hotspots in the community. This big-picture view complements mechanistic and epidemiologic research at the micro level, and underscores the importance of considering real-world chemical mixtures in public health analyses. Ultimately, tackling complex problems like obesity requires integrating information across scales—from molecules and individuals up to communities and ecosystems. By highlighting an environmental dimension to obesity that operates at the county scale, our findings aim to stimulate more nuanced investigations and preventive strategies. Caution is warranted in interpretation due to our study design, but the evidence presented here contributes to a growing consensus that environmental chemical exposures are an important piece of the metabolic health puzzle.