Association of Glyphosate Exposure with Blood DNA Methylation in a Cross-Sectional Study of Postmenopausal Women

The study population consisted of postmenopausal women residing in southern California between the ages of 45 and 66 y ($N=392$). All participants provided informed consent to participate, and the study was approved by the University of California, Irvine institutional review board, HS #2016-3,127. Recruitment, questionnaires, and specimen collection were described in detail previously.⁶⁵ Briefly, women were recruited from 2017 to 2019 via two recruitment arms: targeted mailings to members of the Athena Breast Health Network (a cohort of women receiving screening mammograms at University of California, Irvine health facilities), and untargeted recruitment from the broader community via flyers, social media, events, and word of mouth. The study was designed to identify epigenetic markers for factors potentially related to breast cancer risk; thus, women with a personal history of breast cancer or mastectomy were excluded from the study.

Specimens were collected and processed as previously described.⁶⁵ Briefly, participants self-collected first-morning urine samples on 2 d within a 10-d period not known to them in advance to capture their typical behaviors and exposures. Urine samples were stored in a freezer ($-20°\mathrm{C}$) until they could be transported to the laboratory. At an appointment attempted to be scheduled within 10 d of the second urine sample collection (median 1 d, range 0–24 d), a certified phlebotomist collected peripheral blood in three glass Vacutainer blood collection tubes: one containing ethylenediaminetetraacetic acid (EDTA), one containing acid citrate dextrose (ACD), and one with no anticoagulants. DNA was extracted from the buffy coat (from EDTA and ACD tubes, pooled) within 6 h of collection using the QIAamp DNA Blood Maxi Kit (Cat. No. 51,194, QIAGEN) and stored at $-80°\mathrm{C}$ for later analysis.

Urinary glyphosate and AMPA were measured using liquid chromatography–tandem mass spectrometry (LC-MS/MS) at the Collaborative Center for Translational Mass Spectrometry (CCTMS) (Tgen) using a Vanquish UHPLC coupled to TSQ Altis triple quadrupole mass spectrometer (Thermo Scientific) as previously described.¹² Briefly, glyphosate and AMPA assay validation was performed in a commercially available urine pool from LeeBio Solutions prior to analysis. The calibration curves for glyphosate and AMPA were prepared over a linear range of $0–5\text{\hspace{0.17em}ng}/\mathrm{mL}$ (coefficient of determination $R^2>0.99$) by spiking variable concentrations of glyphosate and AMPA and their respective isotopically labeled internal standards, ${}^{13}{\text{C}}_2$${}^{15}\mathrm{N}\text{-Glyphosate}$ and ${\mathrm{D}}_2{}^{13}\mathrm{C}{}^{15}\mathrm{N}\text{-AMPA}$ (Sigma-Aldrich) at a fixed concentration of $6.25\text{\hspace{0.17em}ng}/\mathrm{mL}$. Data acquisition and processing were performed using Xcalibur 4.1.50 and QuanBrowser, 4.1.50 (Thermo Scientific). Both assays were linear (coefficient of determination $R^2>0.99$) over a range $0–5\text{\hspace{0.17em}ng}/\mathrm{mL}$ (Figure S1). The limits of detection (LOD) for glyphosate and AMPA were 0.014 and $0.013\text{\hspace{0.17em}ng}/\mathrm{mL}$, and the limits of quantitation (LOQ) were 0.041 and $0.040\text{\hspace{0.17em}ng}/\mathrm{mL}$, respectively (Table S1). Creatinine was measured using the DetectX urinary creatinine detection kit (Arbor Assays, K002-H5) according to the manufacturer’s instructions. The data were tested for batch effects using the Kruskal-Wallis test, and batch correction was performed for glyphosate values using the removeBatchEffects function in limma.⁶⁶ Four samples with implausibly low creatinine ($<10\mathrm{mg}/\text{dL}$) were excluded from subsequent analysis. Measurements were averaged for the two urine samples for each analyte after replacing values $<\text{LOD}$ with $\text{LOD}/\surd 2$.⁶⁷

Genomic DNA was bisulfite converted using the Zymo EZ DNA methylation kit (Zymo Research), and then DNA methylation at over 850,000 CpG sites was measured using the Illumina Infinium MethylationEPIC BeadChip (Illumina) at the University of Southern California Molecular Genomics Core. Laboratory staff were blinded to glyphosate and AMPA measurements. Methylation data preprocessing was performed according to recommended steps,⁶⁸ including probe filtering, normalization, and batch correction. All data analysis was performed in R, (version 3.6.2; R Development Core Team). Probes with a detection $p>0.05$ were considered missing, and other low-quality probes were removed as follows: missing in at least 20% of samples ($n=648$); had SNPs with global minor allele frequency $>1\%$ within 5 base pairs of the target sequence or mapping problems with the probe sequence ($n=\mathrm{99,109}$);⁶⁹ hybridized to multiple locations ($n=15$);⁷⁰ or located on the X or Y chromosome ($n=\mathrm{16,927}$). All samples passed quality control, including efficiency of bisulfite conversion, verification of reported sex and identity, and having $<1\%$ failed probes (maximum value: 0.13%), and were included in the final data analysis. Values were normalized with noob normalization as implemented in the minfi package (version 1.32.0)⁷¹ for background correction and dye bias adjustment followed by beta mixture quantile normalization (BMIQ) to correct for type II probe bias.⁷² Except where noted, filtering and normalization was completed using the ChAMP package (version 2.12.4).⁷³ Finally, DNA methylation measurements were corrected for batch and position on chip using ComBat⁷⁴ as implemented in sva version 3.30.1.⁷⁵ We used a reference-based method⁷⁶ to estimate the proportions of six WBC types in our samples, which did not significantly differ with glyphosate or AMPA tertile (Table S2). We used the methylation M value (the base 2 logarithm of the ratio of methylated to unmethylated intensities) for all analyses; in some cases, we also reported $\mathrm{\beta }$ values (percent methylation) for ease of interpretation.⁷⁷

Relevant covariates, including age, self-reported race/ethnicity, height, weight, smoking status (current, former, or never smoker), alcohol consumption, organic eating habits, physical activity, and recent herbicide use, were collected via questionnaires. Race and ethnicity were self-reported in response to two questions: “What is your racial background?” (response choices: “Black or African American,” “White,” “Asian,” “American Indian or Alaska Native,” “Native Hawaiian or other Pacific Islander,” “Some other race,” “Do not know,” or “Prefer not to answer”) and “Are you of Hispanic, Latino, or Spanish origin or ancestry?” (yes/no; “yes” responses were prompted to further self-identify as “Mexican, Mexican American, or Chicano,” “Puerto Rican,” “Cuban,” and/or “Other Hispanic, Latino, or Spanish origin”). Responses to these questions were pooled into a combined race/ethnicity measure with categories non-Hispanic White, non-Hispanic Asian, non-Hispanic Black, Hispanic, or Other. All those who answered “yes” to the Hispanic/Latino/Spanish origin question were placed into the Hispanic category, regardless of response to the race question. Multiple selections were allowed; individuals who selected more than one race ($n=7$) were placed into the “Other” category. Those who answered “American Indian or Alaska Native,” “Native Hawaiian or other Pacific Islander,” or “Some other race” were placed into the “Other” category due to small numbers in the study cohort (total $n$ including $>1$ $\text{race}=10$). Race/ethnicity was included in the analysis due to known associations with DNA methylation.^78,79,80

Self-reported height and weight were used to calculate body mass index (BMI) by dividing the weight in kilograms by the height in meters squared. Organic eating habits were self-reported as “seldom or never,” “sometimes,” or “often or always.”⁸¹ Physical activity was reported as typical frequency and duration of moderate and vigorous physical activity. These responses were converted to minutes per week of moderate exercise, with each minute of vigorous exercise equivalent to 2 min of moderate exercise.⁸² Diet quality was estimated by using up to three unannounced ASA24 24-h dietary recalls,⁸³ completed during the 10 d prior to blood collection, to calculate the Healthy Eating Index-2015 (HEI) averaged across all recalls. The HEI is scored on a 0–100 scale, with higher values indicating greater adherence to the Dietary Guidelines for Americans.⁸⁴ Some participants were asked whether they had used herbicides at home or work within the past 7 d. This question was added to the questionnaire partway through study recruitment, so data on recent herbicide use is missing for 37% ($n=144$) of participants.

Variability and correlation of urinary glyphosate and AMPA levels between and within samples were characterized with the intraclass correlation (ICC) and Pearson’s correlation, respectively. The association between all covariates and the natural logarithm of glyphosate and AMPA concentration (adjusted for urinary creatinine) was assessed using linear regression. In all linear regression analyses, two-sided $p<0.05$ were considered statistically significant. In the DNA methylation analyses, missing data for covariates (from $n=23$ individuals) were substituted with the median (numeric variables) or mode (categorical variables).

We calculated epigenetic age for each sample according to three epigenetic clocks: Hannum,⁵² Horvath,⁵³ and Levine,⁵⁴ and one epigenetic measure of pace of aging, DunedinPoAm.⁵⁵ The Hannum and Horvath clocks are designed to correlate with chronological age, whereas the Levine clock also incorporates phenotypic measures, including mortality and physical functioning. DunedinPoAm is intended to capture the current rate of change of biological age. Epigenetic age was calculated from the model coefficients provided by each author using normalized methylation beta values. Values for missing probes ($n=9\text{\hspace{0.17em}out of\hspace{0.17em}}71$ for Hannum, $n=25\text{\hspace{0.17em}out of\hspace{0.17em}}353$ for Horvath, $n=10\text{\hspace{0.17em}out of\hspace{0.17em}}513$ for Levine, $n=12\text{\hspace{0.17em}out of\hspace{0.17em}}46$ for DunedinPoAm) were imputed using k-nearest neighbors imputation. DunedinPoAm was calculated using code provided by its developers.⁵⁵ Epigenetic age acceleration (residuals from a model regressing chronological age on epigenetic age from each clock) and epigenetic pace of aging from DunedinPoAm was examined for association with glyphosate and AMPA, first in univariate linear regression models adjusted for urinary creatinine concentration and then adjusted for age, race/ethnicity, BMI, smoking status, alcohol consumption, organic eating, and HEI, in addition to batch and position on chip. A third model also included estimated WBC type proportions predicted according to Houseman’s method.⁷⁶ These covariates were selected based on known relationships with either DNA methylation, glyphosate/AMPA, or both; creatinine was included to account for differences in urine concentration.

For probe- and region-specific analyses, the study samples were divided into training ($n=332$) and validation ($n=60$) sets stratified by glyphosate tertile. Demographic and dietary variables were compared for the training and validation sets (Table S3). All variables except smoking were not significantly associated with the randomly assigned set; there was a larger proportion of former smokers in the validation set ($p=0.03$). Because all the former smokers in the validation set reported cessation $>20\mathrm{y}$ prior to blood sample collection, we proceeded with the planned analysis. Using the training set, we identified differentially methylated probes (DMPs) and differentially methylated regions (DMRs) according to a resampling-based method.⁸⁵

We selected a random subsample consisting of 90% of the training set ($n=299$). Using this subsample, we identified candidate DMPs using limma (version 3.44.3),⁶⁶ to fit linear models and candidate DMRs using DMRcate (version 2.0.7),⁸⁶ with a kernel bandwidth of 1,000 base pairs and scaling factor of 2. Each model was adjusted for the same variables as in the final epigenetic aging model. Models without adjustment for dietary variables (organic eating and HEI) did not show substantially different results (Pearson’s $R>0.99$ for test statistics from the two versions of the model for both analytes), nor did models that used creatinine-standardized glyphosate and AMPA concentrations instead of including creatinine as a covariate in the model (Pearson’s $R=0.987$ for glyphosate and 0.985 for AMPA). A false discovery rate (FDR) q-value of $<0.1$ was considered statistically significant. This process was repeated in a series of 1,000 subsamples of 90% of the training set and sampling distributions for all summary statistics constructed. Probes and regions that were selected as candidate DMPs or DMRs in $>90\%$ of subsamples were considered differentially methylated and carried forward for further analysis. This approach results in a more stable list of DMPs and DMRs, because DNA methylation microarray results are known to be sensitive to small differences in the study cohort.^87,88,89,90 The median FDR q-value was used to rank the relative significance of DMPs. For DMRs, overlapping regions were combined, and then significance was ranked after combining median q-values for the probes within the region using Stouffer’s method. A traditional epigenome-wide association analysis with the entire training set⁷⁴ was also done to determine whether findings were robust to changes in the analysis approach.

DMPs and DMRs were annotated using the Illumina manifest (version B4) to identify associated genes and genomic context. Probes were also mapped to ChromHMM data from ENCODE⁹¹ to determine predicted chromatin state; data from the GM12878 lymphoblastic cell line was used because it is the tissue most similar to leukocytes. Chromatin states were pooled into six categories as follows: promoter (active, weak, or poised promoters), enhancer (strong or weak enhancers), transcribed (transcriptional transition, transcriptional elongation, or weak transcribed), polycomb-repressed, inactive (heterochromatin and repetitive regions), and insulators. DMPs were examined for enrichment of certain locations relative to CpG islands or chromatin states, in comparison with the background of all probes included on the array, using Fisher’s exact test. We used the “gometh” function from the missMethyl package (version 1.16.0)⁹² to examine enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and gene ontology terms in genes associated with DMPs and assessed enrichment of pathways using Ingenuity Pathway Analysis (IPA) software (version 62089861; release Date: 17 February 2021). Terms or pathways with $p<0.01$ and at least three genes associated with DMPs were considered enriched.

We developed methylation indices to predict urinary glyphosate and AMPA levels using blood DNA methylation at the glyphosate- and AMPA-associated DMPs. Methylation ($\mathrm{\beta }$) values for each DMP were passed to an elastic net ($\text{alpha}=0.5$) regression model with 10-fold cross-validation repeated 100 times to estimate lambda (the penalty coefficient) and regression coefficients for each site. The final lambda was selected as the largest lambda within one standard error of the minimum prediction error. This approach was implemented using the glmnet package (version 2.0-18).⁹³ The performance of the methylation indices were assessed by calculating the Pearson correlation between predicted and actual glyphosate and AMPA concentration in the 60 validation samples. The indices’ relationship with glyphosate and AMPA tertile was also tested using analysis of variance (ANOVA), and the area under the receiver operating characteristic curve (AUC) was computed to characterize the discriminatory accuracy for the highest vs. lowest tertile.

A total of 96% of study participants had detectable ($>\text{LOD}$) glyphosate, and 82% had detectable AMPA in at least one urine sample. Table 1 describes glyphosate and AMPA concentrations stratified by various demographic, dietary, and lifestyle factors. The median glyphosate concentration (averaged between two urine samples for each participant) was $0.12\text{\hspace{0.17em}ng}/\mathrm{mL}$ (range $<\text{LOD}-1.65\text{\hspace{0.17em}ng}/\mathrm{mL}$), and the median AMPA concentration was $0.06\text{\hspace{0.17em}ng}/\mathrm{mL}$ (range $<\text{LOD}-1.36\text{\hspace{0.17em}ng}/\mathrm{mL}$). Both glyphosate and AMPA concentrations were strongly right-skewed and thus transformed with the natural logarithm for all further analyses (Figure S2). Between-sample agreement was moderate for both glyphosate ($\text{ICC}=0.53$) and AMPA ($\text{ICC}=0.34$), as was the within-sample correlation between glyphosate and AMPA measurements (Pearson’s $R=0.48$, $p<0.001$) (Figure S3).

Women who reported “often or always” eating organic food had lower glyphosate [median (interquartile range, IQR) 0.10 (0.06, 0.20) ng/mL compared with 0.14 (0.07, 0.23) ng/mL for “seldom or never,” $p=0.05$] and AMPA [median (IQR) 0.04 (0.01, 0.11) ng/mL compared with 0.07 (0.03, 0.13) ng/mL for “seldom or never,” $p=0.01$]. Lower diet quality was associated with elevated AMPA [median (IQR) 0.07 (0.02, 0.12) ng/mL for highest quartile of diet quality vs. 0.04 ($<\text{LOD}$, 0.12) ng/mL for lowest quartile, $p=0.03$], but not glyphosate ($p=0.99$). Those who met the Physical Activity Guidelines for Americans ($\ge 150\text{\hspace{0.17em}minutes}$ of moderate-intensity exercise per week) had lower glyphosate [median (IQR) 0.10 (0.05, 0.24) ng/mL for those meeting guidelines vs. 0.13 (0.07, 0.22) for those not meeting, $p=0.03$]. Former smokers had marginally lower AMPA concentrations in comparison with never smokers [median (IQR) 0.04 (0.01, 0.09) ng/mL vs. 0.06 (0.02, 0.13) ng/mL, $p=0.04$]. Women in the “Other” race/ethnicity category had lower concentrations of both glyphosate [median (IQR) 0.07 (0.04, 0.13) vs. 0.12 (0.06, 0.24) ng/mL, $p=0.004$] and AMPA [median (IQR) 0.02 (0.01, 0.05) vs. 0.05 (0.02, 0.12) ng/mL, $p=0.005$] in comparison with White women, although there were only 20 women in the “Other” category. None of the other variables examined (age, BMI, alcohol consumption, or herbicide use) were associated with urinary glyphosate or AMPA.

Each of the three epigenetic clocks (Hannum, Horvath, and Levine) was well-correlated with chronologic age (Pearson’s $r=0.71$, 0.64, and 0.54, respectively), although the Hannum and Levine clocks were poorly calibrated in this sample (Figure S4). Glyphosate was not significantly associated with differences in epigenetic age acceleration according to any of the epigenetic clocks (Table 2). AMPA was associated with increased epigenetic age acceleration according to the Hannum clock [adjusted coefficient 0.36 y; 95% confidence interval (CI): 0.03, 0.70, $p=0.04$]. AMPA was also associated with increased epigenetic age acceleration according to the Horvath clock in the adjusted model (adjusted coefficient 0.41; 95% CI: 0.02, 0.80, $p=0.04$) but not the univariate model ($p=0.24$); the opposite was true for the Levine clock (univariate coefficient 0.61; 95% CI: 0.04, 1.18, $p=0.04$; adjusted $p=0.32$). Epigenetic pace of aging according to DunedinPoAm was not significantly associated with glyphosate or AMPA concentration (Table 2).

Twenty-four probes associated with urinary glyphosate concentration and two with urinary AMPA were identified in $>90\%$ of subsamples and considered differentially methylated (Table 3). Seventeen of the 24 probes (71%) were hypomethylated with higher glyphosate (Figure 1), whereas both AMPA-associated probes were hypermethylated (Table 3). Glyphosate-associated probes with the smallest median $p$-values were located within the PEX26, SF3B2, and CHMP1A genes. The largest methylation difference between glyphosate tertiles were observed in probes within the VMO1, KCP, and QARS1 genes. Our findings were robust to changes in the analysis approach: All glyphosate- and AMPA-associated DMPs were statistically significant (FDR $q<0.05$) using a traditional epigenome-wide association analysis with the entire training set, except for one AMPA-associated DMP (cg23261070, $q=0.06$).

Biological and functional analyses were performed on only glyphosate-associated probes because there were only two AMPA-associated probes (Figure 2). In comparison with all other probes on the array, a greater proportion of glyphosate-associated DMPs were located within CpG islands (33.3% vs. 19.2%), although this difference was not statistically significant (Fisher $p=0.21$). Glyphosate-associated DMPs were significantly enriched for enhancer regions (29.2% vs. 12.6%, Fisher $p=0.02$) and depleted for inactive/heterochromatin regions (4.2% vs. 36.8%, Fisher $p=0.0004$). Three gene ontology terms were enriched in the 20 annotated genes containing glyphosate-associated probes: DNA metabolic process (biological process), intracellular organelle part (cellular component), and organelle part (cellular component), all with $p<0.01$. No KEGG or IPA pathways were enriched in glyphosate-associated genes.

Four regions for glyphosate and none for AMPA were significant in $>90\%$ of subsamples and considered differentially methylated (Table 4). Three regions were significantly associated with AMPA at a relaxed threshold of $>60\%$ of subsamples and are also included in Table 4. The glyphosate-associated regions were all located within gene promoters. Three were hypomethylated with greater glyphosate (MSH4, KCNA6, and ABAT), and the other was hypermethylated (NDUFAF2/ERCC8). The top AMPA-associated regions were all hypomethylated. Two were located within gene bodies (RNF39, TRIM31) and one within a gene promoter (ESR1).

All 24 DMPs were selected by the elastic net for inclusion in the final glyphosate methylation index (Table S4). The index was significantly correlated with urinary glyphosate in the training (Pearson’s $R=0.62$; 95% CI: 0.55, 0.68, $p<0.0001$; Figure 3A) and validation ($R=0.35$; 95% CI: 0.10, 0.55, $p=0.006$; Figure 3B) sets. In the validation set, the methylation index was significantly associated with glyphosate tertile [median (IQR) $-2.4$ ($-2.7$ to $-2.2$) for lowest tertile vs. $-2.0$ ($-2.3$ to $-1.6$) for highest tertile, ANOVA $p=0.009$; Figure 3C] and showed excellent discrimination between the top and bottom tertiles of urinary glyphosate ($\text{AUC}=0.74$, 95% CI: 0.57, 0.90; Figure 3D). A modified index trained using only probes included on the HumanMethylation450 BeadChip (the previous version of the Illumina methylation array), which selected 14 probes in the final model, achieved similar performance (Figure S5). The AMPA index was not significantly correlated with urinary AMPA in either the training or validation set (Figure S6).