Epigenome-Wide DNA Methylation and Pesticide Use in the Agricultural Lung Health Study.

BACKGROUND: Pesticide exposure is associated with many long-term health outcomes; the potential underlying mechanisms are not well established for most associations. Epigenetic modifications, such as DNA methylation, may contribute. Individual pesticides may be associated with specific DNA methylation patterns but no epigenome-wide association study (EWAS) has evaluated methylation in relation to individual pesticides.
OBJECTIVES: We conducted an EWAS of DNA methylation in relation to several pesticide active ingredients.
METHODS: The Agricultural Lung Health Study is a case-control study of asthma, nested within the Agricultural Health Study. We analyzed blood DNA methylation measured using Illumina's EPIC array in 1,170 male farmers of European ancestry. For pesticides still on the market at blood collection (2009-2013), we evaluated nine active ingredients for which at least 30 participants reported past and current (within the last 12 months) use, as well as seven banned organochlorines with at least 30 participants reporting past use. We used robust linear regression to compare methylation at individual C-phosphate-G sites (CpGs) among users of a specific pesticide to never users.
RESULTS: Using family-wise error rate (p<9×10-8) or false-discovery rate (FDR<0.05), we identified 162 differentially methylated CpGs across 8 of 9 currently marketed active ingredients (acetochlor, atrazine, dicamba, glyphosate, malathion, metolachlor, mesotrione, and picloram) and one banned organochlorine (heptachlor). Differentially methylated CpGs were unique to each active ingredient, and a dose-response relationship with lifetime days of use was observed for most. Significant CpGs were enriched for transcription motifs and 28% of CpGs were associated with whole blood cis-gene expression, supporting functional effects of findings. We corroborated a previously reported association between dichlorodiphenyltrichloroethane (banned in the United States in 1972) and epigenetic age acceleration. DISCUSSION: We identified differential methylation for several active ingredients in male farmers of European ancestry. These may serve as biomarkers of chronic exposure and could inform mechanisms of long-term health outcomes from pesticide exposure. https://doi.org/10.1289/EHP8928.

Introduction

Pesticides can have both acute and long-term health effects. The mechanisms for acute health effects (e.g., cough, skin or eye irritation, headache, respiratory depression) are generally well understood. Pesticide exposure has also been associated with numerous long-term adverse effects in humans, including neurologic diseases and cancers (Mostafalou and Abdollahi 2017). The potential mechanisms underlying some of these long-term health effects are better understood than others. For example, organophosphates kill insects by inhibiting acetylcholinesterase, an enzyme critical for neurotransmission, and thus, leads to neurotoxicity (Costa et al. 2008). Organophosphates have been associated with neurologic diseases in humans, including Parkinson’s and Alzheimer’s (Mostafalou and Abdollahi 2017). However, mechanisms for other pesticides and long-term health effects are less well understood. For example, the International Agency for Research on Cancer (IARC) has concluded that there is sufficient evidence that lindane causes non-Hodgkin’s lymphoma in humans (IARC 2018). The exact mechanism is not clear, but the IARC has determined there is strong evidence that lindane is immunosuppressive and induces oxidative stress (IARC 2018).

Epigenetics may contribute to the molecular mechanisms underpinning some long-term health effects of pesticides. Epigenetics refer to DNA modifications that do not change the DNA sequence. The best-studied epigenetic mechanism in humans is DNA methylation. The advent of stable platforms has enabled researchers to measure DNA methylation at individual C-phosphate-G sites (CpGs) with reasonable genome-wide coverage. Agnostic epigenome-wide methylation studies of smoking and alcohol intake have identified specific differentially methylated CpGs that are reproducible and can serve as signatures of chronic exposure (Joehanes et al. 2016; Liu et al. 2018; Maas et al. 2019). Differential methylation at these reproducible CpGs could also contribute to the biological mechanisms underlying chronic exposure and health effects, supporting the use of epigenome-wide studies to identify signatures for other chronic exposures such as pesticides.

Several studies have reported associations between pesticide exposure and global DNA methylation (e.g., Alu, LINE-1) (Benitez-Trinidad et al. 2018; Itoh et al. 2014; Kim et al. 2010; Lind et al. 2013). Only three epigenome-wide methylation studies of pesticides have been published and each used the Illumina 450K array, which measures methylation at . In a population-based cohort in the Netherlands, van der Plaat et al. (2018) used a single exposure metric for all occupational pesticide exposures combined and no differential methylation was observed. In a population-based study of individuals living in an agricultural region in California, ambient exposure to organophosphates and pyrethroids was assessed from residential and workplace addresses. Differential methylation was reported in relation to the two pesticide classes, organophosphates (Paul et al. 2018) and pyrethroids (Furlong et al. 2020). Because a given chemical class contains many different active ingredients, each with a unique chemical structure, individual active ingredients may have unique methylation signals not shared among the class as a whole. To our knowledge, no epigenome-wide association studies (EWASs) have examined methylation in relation to specific active ingredients.

The advent of stable platforms to measure DNA methylation has also led to the development of epigenetic clocks composed of CpGs correlated with age and have subsequently been shown to predict mortality and age-related morbidities, capturing the variability in biologic, as opposed to chronologic, age (Bell et al. 2019). A previous study reported a positive association between levels of a metabolite of dichlorodiphenyltrichloroethane (DDT) and epigenetic aging (Lind et al. 2018), suggesting that pesticide exposure may accelerate one’s aging process. No studies have confirmed these findings nor have they examined other pesticides.

We conducted an EWAS to examine differential methylation related to different pesticide active ingredients, using the newer EPIC array, which has greater genomic coverage (). Additionally, we examined the association between these active ingredients and epigenetic aging. Analyses were conducted among men of European ancestry from a U.S. farming population.

Methods

Study Population

The Agricultural Lung Health Study (ALHS) is a case–control study of current asthma nested within the Agricultural Health Study (AHS) (data version P3REL201209.00). Details of the AHS and ALHS have been previously described (Alavanja et al. 1996; House et al. 2017). Briefly, the AHS enrolled licensed pesticide applicators (mostly farmers) from Iowa or North Carolina and their spouses between 1993 and 1997 (Phase 1). Participants were followed up with a computer-assisted telephone interview between 1999–2003 (Phase 2) and 2005–2010 (Phase 3). Information about asthma from the Phase 3 interview was used to enroll farmers or their spouses into the ALHS between 2009 and 2013. From the 2,363 asthma cases identified in the AHS, 1,223 cases were enrolled in the ALHS, and 2,078 noncases were randomly selected among the AHS participants who did not meet our definition of asthma. All participants provided informed consent for the ALHS. The institutional review board at the National Institutes of Health and its contractors approved this study.

Data Collection

In the ALHS, field technicians visited individuals’ homes to collect blood samples. A computer-assisted telephone interview collected information on smoking history that was used to update participants’ smoking status (never, former, current) and pack years reported on earlier questionnaires from the parent AHS cohort.

DNA Methylation Measurement, Quality Control, and Estimation of Cell Type Proportions

To minimize the effects of population stratification, we measured DNA methylation in 2,391 ALHS participants of European ancestry. Ancestry was determined from principal component analysis based on genome-wide genetic data, and participants with principal components suggestive of non-European ancestry did not have methylation measured (Sikdar et al. 2021). DNA from whole blood was bisulfite converted using the EZ-96 DNA Methylation kit (Zymo Research Corporation). Illumina’s MethylationEPIC BeadChip (EPIC array) was used to assess DNA methylation (Illumina, Inc.).

We excluded samples with of CpGs with detection , intensity values less than three standard deviations below the mean bisulfite control intensities, or samples with sex mismatch. We removed 31,533 CpGs with of samples with detection . Using the ENmix R package, background and dye bias were corrected using Relic (Xu et al. 2016, 2017), data were normalized using interarray (quantile) normalization (Xu et al. 2016), and probe-type bias adjustment was conducted using Rcp (Niu et al. 2016). Batch effects (i.e., plate effects) were corrected using ComBat from the sva R package (Johnson et al. 2007). Extreme methylation outliers were replaced with winsorized values () using ewaff.handle.outliers from the ewaff R package (Suderman et al. 2019). We filtered out 66,353 probes where methylation may not be accurately targeted or measured: probes with a single nucleotide polymorphism (SNP) in the extension base that causes a color channel switch, probes masked for extension base inconsistent with specified color channel, probes with low quality mapping, and probes that start with “ch” or are SNPs (Zhou et al. 2017). After excluding 19,681 CpGs from the sex chromosomes, we analyzed 755,489 CpGs on the autosomes. After methylation quality control, DNA methylation was available in 2, 288 participants (Figure S1).

We estimated cell type proportions for seven estimated cell types (monocytes, neutrophils, eosinophils, natural killer cells, B cells, cells, and cells) using the Houseman method (Houseman et al. 2012) with the Reinius reference panel (Reinius et al. 2012).

Pesticide Exposure Assessment

Because we focused on direct exposure to pesticide, we restricted analyses to the 1,170 males because the licensed pesticide applicator was nearly always the male spouse. Additionally, the frequency of pesticide use reported among female spouses was lower. In the AHS, participants provided the names of chemicals they had ever used in their lifetime on the Phase 1 questionnaire (1993–1997), in the most recent farming season in Phase 2 (1999–2003), and since the last study contact in Phase 3 (2005–2010). In the ALHS (2009–2013), participants reported names of the pesticides they had used within the past 12 months.

For each pesticide, responses in the AHS questionnaires were combined to classify participants’ past use of an active ingredient (yes/no). Participants who indicated use of an active ingredient from any of the AHS questionnaires had past use. Participants who reported no use of the active ingredient in both Phase 1 and Phase 3 had no past use. Those who had missing information in Phase 1 or Phase 3 were set to missing. Responses in the ALHS were used to classify participants’ current use (within the last 12 months) of an active ingredient (yes/no). Analyses of smoking show that methylation reverts to unexposed levels at many loci with time since quitting smoking (Joehanes et al. 2016). Therefore, we compared participants who reported use of a given active ingredient in any prior AHS questionnaires and use within the 12 months before the ALHS blood draw (i.e., past and current use) with those who had never used that active ingredient. Thus, individuals who reported using the active ingredient only on a previous AHS questionnaire (i.e., past use only) or only within the past 12 months (i.e., current use only) were not included in analyses of that active ingredient. To reduce false-positive results, epigenome-wide analyses were restricted to active ingredients with at least 30 pesticide users. We analyzed nine active ingredients for which at least 30 men reported both past and current use: dicamba, picloram, mesotrione, acetochlor, metolachlor, glyphosate, 2,4-dichlorophenoxyacetic acid (2,4-D), atrazine, and malathion. Although agricultural use of most organochlorine pesticides was banned by the 1980s, these pesticides have long half-lives and bioaccumulate. Therefore, we also analyzed seven banned active ingredients related to organochlorine with at least 30 men reporting past use: aldrin, chlordane, DDT, dieldrin, heptachlor, lindane, and toxaphene. Thus, a total of 16 active ingredients were analyzed.

Lifetime years and average days per year of use were calculated using information collected in the three AHS questionnaires (Phase 1, Phase 2, and Phase 3) for 15 of the 16 active ingredients. We estimated lifetime days of use for all active ingredients except for picloram because years of use were not available for this chemical. Farmers who used the active ingredient were dichotomized at the median lifetime days.

Statistical Analyses

Epigenome-wide analyses were conducted for each active ingredient. We used robust linear regression to compare DNA methylation among farmers who reported use of an active ingredient (both past and current use for currently marketed pesticides, past use for banned pesticides) with those who never used the active ingredient. We calculated robust standard errors using White’s estimator in the sandwich R package. Untransformed methylation beta values, ranging from 0 (unmethylated) to 1 (methylated), were the outcome. All epigenome-wide analyses were adjusted for the study selection factor (asthma status), age at blood collection, smoking status (never, former, current), pack years (continuous), state of residence (Iowa or North Carolina), and seven estimated cell type proportions (monocytes, neutrophils, eosinophils, natural killer cells, B cells, CD4T cells, and CD8T cells—all continuous). For each separate active ingredient, significance was assessed using both a family-wise error rate (FWER) () (Mansell et al. 2019) and the Benjamini and Hochberg false-discovery rate (FDR) () (Benjamini and Hochberg 1995). For all FDR-significant CpGs, we visually inspected the distribution for potentially polymorphic CpGs and removed CpGs with a bimodal distribution (suggestive of influence of a SNP) from the functional enrichment and gene expression analyses.

We used lifetime days of use to dichotomize pesticide users at the median and conducted a trend test to evaluate for dose response. Lifetime days was analyzed as an ordinal variable where never users were coded as 0, lower than the median was 1, and above the median was 2. For these analyses, significance was assessed nominally (). All analyses were conducted in R (version 3.4.0; R Development Core Team).

Gene Annotation

We used Hypergeometric Optimization of Motif EnRichment software (HOMER; version 4.10.3; Heinz et al. 2010) and the human genome database (hg19; version 5.10) provided by HOMER (Heinz et al. 2010) to annotate the CpGs to genes using location information on the genome. Each CpG was mapped to the closest transcription start site based on the RefSeq annotation, using the annotatePeaks.pl script with the default setting. CpGs from the transcription start site were annotated to the nearest gene name from the Gene Symbol column.

Functional Enrichment

We used eFORGE (version 2.0; Breeze et al. 2019) to examine whether the FDR-significant CpGs were enriched for DNase I hypersensitive sites, 15 chromatin states, and 5 histone marks in specific tissues and cells. eFORGE-TF (Breeze et al. 2019) was used to identify enrichment of specific transcription factor motifs.

Association between Significant CpGs and Gene Expression

We assessed whether FDR-significant CpGs correlated with the expression of nearby genes (cis-expression). Because we found no gene expression studies with methylation data from the EPIC array, we evaluated only our significant CpGs present on the array. Four cohorts from the Biobank-based Integrative Omics Study (BIOS) Consortium (Leiden Longevity Study, LifeLines Study, Rotterdam Study, and Netherland Twin Study) had methylation data from the Illumina 450K array and gene expression data from RNA-seq in blood in 3,075 adults (Bonder et al. 2017). We examined the association between the CpG and each gene expression transcript within of the CpG. Within each cohort, gene expression was regressed on methylation (based on M-values) using limma, a common statistical R package used for analyzing gene expression data (Ritchie et al. 2015). Analyses were adjusted for age, sex, lymphocyte proportion, monocyte proportion, and RNA flow cell number. To reduce potential false positives, inflated test statistics were corrected using bacon (van Iterson et al. 2017). The results from each cohort were meta-analyzed using a fixed-effects model with inverse variance weights in METAL (Willer et al. 2010). Significance was assessed at .

Look Up of Implicated Genes in Experimental Studies

To further validate our findings, we searched three databases with data on how pesticide exposure affects gene expression in experimental studies: Illumina’s BaseSpace Correlation Engine (Kupershmidt et al. 2010) (formerly known as NextBio) (Illumina 2020, the Comparative Toxicogenomics Database (Davis et al. 2019; CTD 2020), and the Chemical Dashboard of the U.S. Environmental Protection Agency (EPA) (Williams et al. 2017; U.S. EPA 2020). The BaseSpace Correlation Engine combines published and unpublished data sets from public, private, and proprietary sources; the Comparative Toxicogenomics Database relies on published literature; and the Chemical Dashboard is based on in vitro studies on a limited number of genes. In each database, we extracted the genes that have been shown to have altered expression after exposure to the active ingredient. The resulting genes identified using these databases were compared with the genes annotated to the FDR-significant CpGs.

Look Up of Our Significant CpGs in an Independent Study

We identified one study with information on active ingredients and blood DNA methylation measured in adults. The Parkinson’s Environment and Genes (PEG) Study is a case–control study of Parkinson’s disease (Paul et al. 2018). Briefly, participants were recruited from the agricultural regions of the California Central Valley from 2000 to 2007 (Wave 1) and 2010–2015 (Wave 2). Using a geographic information system, the researchers combined pesticide use records from California and land-use maps to estimate ambient exposures to active ingredients at residential and workplace addresses. For each pesticide, the pounds applied per year and per acre within a buffer of each residential and occupational address were summed. For this study, participants were considered exposed to an active ingredient if they were exposed to the pesticide within the last 5 y prior to the year of blood draw. Methylation was measured in peripheral whole blood using the Illumina Infinium 450K array in 789 adults (469 males). The methylation values were preprocessed using background normalization from Genome Studio. Only significant CpGs () in the ALHS were analyzed. Data from both PEG Study waves were analyzed together. Analyses were conducted for active ingredients with at least 15 exposed participants and were adjusted for age, sex, current smoking status (yes/no), six estimated cell types (monocytes, CD4T cells, CD8T cells, B cells, natural killer cells, granulocytes), study wave, and the study selection factor (Parkinson’s disease). Because the PEG Study included participants of European and non-European descent, genetic ancestry based on ancestry informative markers was also included in the model to account for population stratification.

Look Up of Previously Reported CpGs Related to Pesticide in the ALHS

The PEG Study previously reported 70 CpGs associated with organophosphates (Paul et al. 2018) at () and four differentially methylated CpGs associated with pyrethroids at (Furlong et al. 2020). We used the active ingredients that the PEG Study used to define organophosphates and pyrethroids to classify exposure in ALHS participants. ALHS participants were considered organophosphate users if they reported use of any organophosphate active ingredient in the AHS and within the past 12 months. Similarly, ALHS participants were considered pyrethroid users if they reported use of any pyrethroid active ingredient in the prior AHS questionnaires and within the past 12 months from blood draw (Excel Table S1). In the ALHS, comparing users to never users, we ran epigenome-wide analyses using robust linear regression and adjusted for the same covariates included in our active ingredient analyses. These analyses enabled the ALHS to look up the 74 CpGs reported in the PEG Study.

Epigenetic Age Acceleration

We explored whether using specific active ingredients contributes to epigenetic aging by calculating participants’ epigenetic ages based on four epigenetic clocks: the Horvath (Horvath 2013), Hannum (Hannum et al. 2013), Skin and Blood (Horvath et al. 2018), and PhenoAge clocks (Levine et al. 2018). The Horvath, Hannum, and Skin and Blood clocks use specific sets of CpGs to estimate chronological age as a surrogate of biological age. The PhenoAge clock estimates phenotypic aging based on a different set of CpGs related to a composite of clinical measures. Measures of epigenetic aging were calculated with the Normalize Data and Advanced Analysis options using the New Methylation Age Calculator (DNA Methylation Age Calculator 2020). We subtracted chronological age from epigenetic age to obtain , a measurement of epigenetic age acceleration. We used linear regression to evaluate the association between epigenetic age acceleration and each of the 16 active ingredients, adjusting for the same covariates used in the epigenome-wide analyses. Significance was assessed using a nominal .

Results

The median age of the 1,170 male farmers at ALHS enrollment was 63 y (Table 1). The majority of the farmers enrolled in the ALHS were from Iowa (70%) and a small percentage were current smokers (6%). Some participants had missing past use information for certain specific active ingredients and thus could not be categorized into one of the four categories of users (median missing: 34; range: 0 to 74). The number and percentage of past-only or current-only users for each currently marketed active ingredient is provided in Table 2 (median: 54.2%).

Table 1

Demographic characteristics of 1,170 male farmers of European ancestry in the Agricultural Lung Health Study, 2009–2013.

Characteristic	n (%) or median (IQR)
Smoking status
Never	679 (58.0)
Former	427 (36.5)
Current	64 (5.5)
State
Iowa	820 (70.1)
North Carolina	350 (29.9)
Asthma
Yes	438 (37.4)
No	732 (62.6)
Age at blood draw (y)	63 (55–72)
Pack years	14 (3–32)a
Estimated cell type proportion
Monocyte	0.08 (0.07–0.10)
Neutrophil	0.59 (0.52–0.66)
Eosinophil	0 (0–0.001)b
Natural killer cell	0.06 (0.04–0.10)
B cell	0.05 (0.03–0.07)
CD4T cell	0.16 (0.11–0.20)
CD8T cell	0.05 (0.02–0.08)

Note: IQR, interquartile mean.

Among ever smokers.

Median (3rd–97th percentile).

Table 2

Summary of epigenome-wide analyses of DNA methylation and pesticide use, enriched transcription factor motifs, and associations with nearby gene expression.

Pesticide	Users (n)a	Never users (n)	Past-only or current-only users [(n) (%)]b	CpGs (n) p<9×10^−8c	CpGs (n) FDR<0.05c	λ d	Enriched TF motifs (n)e	CpGs associated with cis-gene expression (n) FDR<0.05f
Herbicide
Acetic acid
Dicamba	43	424	629 (57.4)	3	6	0.94	3	3
Picloram	32	979	159 (13.6)	8	48	1.05	35	40
Anilide
Mesotrione	39	1,020	111 (9.5)	12	72	1.11	27	20
Chloroacetanilide
Acetochlor	41	968	161 (13.8)	6	21	0.98	24	11
Metolachlor	53	454	600 (54.2)	5	6	1.05	11	1
Organophosphorus
Glyphosate	523	69	573 (49.2)	1	1	0.96	1	—
Phenoxy
2,4-D	325	141	691 (59.7)	0	0	1.20	—	—
Triazine
Atrazine	172	191	799 (68.8)	1	1	0.96	0	1
Insecticide
Organochlorine
Aldrin	286	845	—	0	0	1.04	—	—
Chlordane	382	750	—	0	0	1.03	—	—
DDT	375	765	—	0	0	1.07	—	—
Dieldrin	108	1,019	—	0	0	0.94	—	—
Heptachlor	242	889	—	3	6	1.60	11	3
Lindane	328	804	—	0	0	1.08	—	—
Toxaphene	215	915	—	0	0	1.22	—	—
Organophosphate
Malathion	36	196	912 (79.7)	1	1	1.04	0	—

Note: —, not applicable; BIOS, Biobank-based Integrative Omics Study; CpGs, C-phosphate-G sites; DDT, dichlorodiphenyltrichloroethane; FDR, false discovery rate; TF, transcription factor; 2,4-D, 2,4-dichlorophenoxyacetic acid.

Users reported both past and current use of the pesticide, except for organochlorines. Only past use information was available for organochlorines.

The sum of never users, past-only users, current-only users, and past and current users may not add up to 1,170 for all pesticides where use of that specific activate ingredient in the Agricultural Health Study questionnaires is missing.

Model adjusted for selection factor (asthma status), age at blood collection, pack years, smoking, state of residence, and seven estimated cell types.

Genomic inflation statistic.

Significant CpGs were pasted into eFORGE-TF (https://eforge-tf.altiusinstitute.org/) for enrichment.

Association between significant CpGs on the Illumina 450K and expression of nearby genes () in the BIOS consortium.

A total of 162 CpGs were differentially methylated ( or ) in relation to eight of the nine currently marketed active ingredients (dicamba, picloram, mesotrione, acetochlor, metolachlor, glyphosate, atrazine, and malathion) and one banned organochlorine (heptachlor) (Table 2). None of 162 differentially methylated CpGs overlapped across the active ingredients. Methylation was higher in farmers who had used the specific active ingredient compared with those who never used it at approximately half of the significant CpGs (83/162) (Tables S2–S10). The active ingredients with the highest number of differentially methylated CpGs () were picloram (48 CpGs; 32 users) and mesotrione (72 CpGs; 39 users). The CpGs differentially methylated at the FWER () in relation to a specific active ingredient are presented in Table 3; all CpGs differentially methylated at an are shown in Excel Tables S2–S10. The full epigenome-wide DNA methylation results for each of the 16 pesticides are available at https://zenodo.org/record/5156585#.YS-gd4hKg2w. There was minimal evidence of systematic biases or genomic inflation across the analyses of active ingredients (range of : 0.94–1.60). Manhattan, volcano, and QQ plots for each of the analyses of the 16 active ingredients are presented in Figures S2–S17. Because DDT was banned in the U.S. market at the end of 1972, younger farmers in our study would not have had access to DDT. Repeating the DDT epigenome-wide analyses and restricting it to the 818 farmers who were at least 18 years of age in 1972 did not yield any significant CpGs (data not shown).

Table 3

Differentially methylated CpGs at a family-wise error rate () in relation to use of active ingredients, Agricultural Lung Health Study, 2009–2013.

Active ingredient	Probe	Chr	Positiona	Users (n)b	Never users (n)	Beta	SE	p-Valuec	Gene name [symbol (full name)]d
Herbicide
Dicamba	cg08004045	3	127797258	43	424	0.008	0.001	1.05 ×10^–10	RUVBL1-AS1 (RuvB Like AAA ATPase 1 antisense RNA 1)
	cg14830371	8	25991602	43	423	−0.006	0.001	1.62 ×10^–9	EBF2 (early B cell factor transcription factor 2)
	cg23238147	1	15439013	43	424	−0.006	0.001	2.17 ×10^–8	TMEM51-AS1 (transmembrane protein 51 antisense RNA 1)
Picloram	cg06216309	7	101603359	32	976	0.029	0.004	9.27 ×10^–12	CUX1 (cut like homeobox 1)
	cg04852443	1	2990064	32	974	−0.004	0.001	5.16 ×10^–10	PRDM16 (PR domain containing 16)
	cg09671989	16	88330828	32	976	0.006	0.001	7.15 ×10^–9	LINC02182 (long intergenic non-protein coding RNA 2182)
	cg21760990	6	121656542	32	979	0.029	0.005	1.19 ×10^–8	TBC1D32 (TBC1 domain family member 32)
	cg04815626	8	9912509	32	979	−0.001	0.0003	1.44 ×10^–8	MSRA (methionine sulfoxide reductase A)
	cg16188532	13	47471090	32	979	−0.001	0.0002	5.76 ×10^–8	HTR2A (5-hydroxytryptamine receptor 2A)
	cg16273766	11	28131811	32	979	−0.001	0.0001	6.03 ×10^–8	METTL15 (methyltransferase-like 15)
	cg02034311	13	112547999	32	978	−0.056	0.010	7.10 ×10^–8	LINC00354 (long intergenic non-protein coding RNA 354)
Mesotrione	cg18025430	11	60048004	39	1,018	−0.011	0.002	5.04 ×10^–9	MS4A4A (membrane spanning 4-domains A4A)
	cg02156591	3	130480267	39	1,019	−0.033	0.006	8.86 ×10^–9	PIK3R4 (phosphoinositide-3-kinase regulatory subunit 4)
	cg01919701	13	112996346	39	1,020	0.051	0.009	1.12 ×10^–8	LINC01044 (long intergenic non-protein coding RNA 1044)
	cg22053500	6	111580379	39	1,017	−0.002	0.0003	1.25 ×10^–8	MFSD4B (major facilitator superfamily domain containing 4B)
	cg21270074	20	49412129	39	1,020	−0.001	0.0001	2.50 ×10^–8	BCAS4 (breast carcinoma amplified sequence 4)
	cg03584506	11	39689828	39	1,019	0.057	0.010	3.46 ×10^–8	LRRC4C (leucine rich repeat containing 4C)
	cg16258223	5	134377058	39	1,020	0.054	0.010	4.34 ×10^–8	C5orf66-AS1 (esophagus epithelial intergenic associated transcript)
	cg11223552	3	158984799	39	1,020	0.044	0.008	4.37 ×10^–8	SCHIP1 (schwannomin interacting protein 1)
	cg14870603	8	31986050	39	1,018	0.059	0.011	4.39 ×10^–8	NRG1-IT1 (neuregulin 1 intronic transcript 1)
	cg24255843	2	192140815	39	1,020	0.004	0.0007	5.22 ×10^–8	MYO1B (myosin IB)
	cg15263667	2	29118655	39	1,019	−0.001	0.0002	6.51 ×10^–8	WDR43 (WD repeat domain 43)
	cg09849314	3	195283427	38	1,018	−0.067	0.012	7.87 ×10^–8	PPP1R2 (protein phosphatase 1 regulatory inhibitor subunit 2)
Acetochlor	cg16722983	14	91244421	41	968	0.031	0.005	1.57 ×10^–9	TTC7B (tetratricopeptide repeat domain 7B)
	cg08165462	19	652493	41	968	0.013	0.002	1.47 ×10^–8	RNF126 (ring finger protein 126)
	cg01869283	7	37209906	40	964	0.009	0.002	1.77 ×10^–8	ELMO1-AS1 (engulfment and cell motility 1 antisense RNA 1)
	cg04108706	1	180126520	41	967	0.017	0.003	3.72 ×10^–8	QSOX1 (quiescin sulfhydryl oxidase 1)
	cg10252897	17	38275458	41	968	0.008	0.002	5.56 ×10^–8	MSL1 (male-specific lethal 1)
	cg00919585	3	51975770	41	968	−0.002	0.0004	7.90 ×10^–8	RRP9 (ribosomal RNA processing 9, U3 small nucleolar RNA binding protein)
Metolachlor	cg17001333	1	213843678	53	454	−0.079	0.014	5.31 ×10^–9	LINC00538 (long intergenic non-protein coding RNA 538)
	cg20390515	2	945697	53	454	−0.002	0.0003	9.07 ×10^–9	SNTG2 (syntrophin gamma 2)
	cg18085807	6	30313188	51	439	−0.002	0.0003	2.51 ×10^–8	RPP21 (ribonuclease P/MRP subunit p21)
	cg14219242	5	179051847	53	454	0.002	0.0004	4.06 ×10^–8	HNRNPH1 (heterogeneous nuclear ribonucleoprotein H1)
	cg16702144	9	131218309	52	448	0.0004	0.00007	8.18 ×10^–8	ODF2 (outer dense fiber of sperm tails 2)
Glyphosate	cg06950346	7	123978181	523	69	−0.027	0.005	1.21 ×10^–8	LOC101928211 (uncharacterized LOC101928211)
Atrazine	cg01867395	11	31839628	172	187	0.033	0.006	1.99 ×10^–8	PAX6 (paired box 6)
Insecticide
Heptachlor	cg18284070	14	57736207	242	889	−0.0004	0.00007	1.19 ×10^–8	EXOC5 (exocyst complex component 5)
	cg04265842	3	156807650	242	889	−0.007	0.001	5.40 ×10^–8	LINC00881 (long intergenic non-protein coding RNA 881)
	cg16310513	18	55471075	242	889	0.026	0.005	7.08 ×10^–8	ATP8B1 (ATPase phospholipid transporting 8B1)
Malathion	cg26210302	3	42913809	36	196	0.045	0.008	5.13 ×10^–9	CYP8B1 (cytochrome P450 family 8 subfamily B member 1)

Note: All significant findings at are shown in Excel Tables S2–S10. Chr, chromosome; CpGs, C-phosphate-G sites; FDR, false-discovery rate; SE, standard error.

Build hg19/GRCh37.

Except for heptachlor, users reported both past and current use of the pesticide. Heptachlor users only had past use.

Analyses conducted using robust linear regression, adjusting for the selection factor (asthma status), age at blood collection, pack years, smoking (never, former, current), state of residence, and seven estimated cell types.

Annotated using HOMER (version 4.10.3; Heinz et al. 2010).

Of the nine active ingredients with differentially methylated CpGs, lifetime days of use data were available for all active ingredients except picloram. Of the 114 FDR-significant CpGs across these eight active ingredients, the ordinal test for trend gave a for all but one CpG (Table S11)—6/6 CpGs with dicamba: cg08004045 [RuvB Like AAA ATPase 1 antisense RNA 1 (RUVBL1-AS1)], cg14830371 [early B cell factor transcription factor 2 (EBF2)], cg23238147 [transmembrane protein 51 antisense RNA 1 (TMEM51-AS1)], cg03033508 [ring finger protein 6 (RNF6)], cg01422293 long intergenic non-protein coding RNA 920 (LINC00920), cg09230290 (LINC01399); 72/72 CpGs with mesotrione (top 6 listed): cg20421702 [integrin subunit alpha 6 (ITGA6)], cg13780053 [semaphorin 4D (SEMA4D)], cg21270074 [breast carcinoma amplified sequence 4 (BCAS4)], cg14063331 [collagen beta(1-O)galactosyltransferase 1 (COLGALT1)], cg01536905 [microRNA 4802 (MIR4802)], cg04461219 [nuclear receptor subfamily 2, group F, member 2 (NR2F2)]; 21/21 CpGs with acetochlor (top 6 listed): cg08165462 [ring finger protein 126 (RNF126)], cg15300856 [cystathionine beta-synthase (CBS)], cg24578004 [DNA fragmentation factor subunit beta (DFFB)], cg18677834 [activating signal cointegrator 1 complex subunit 3 (ASCC3)], cg11702745 [long intergenic non-protein coding RNA 1422 (LINC01422)], cg17840355 [tRNA methyltransferase 44 homolog (TRMT44)]; 6/6 CpGs with metolachlor: cg17001333 [long intergenic non-protein coding RNA 538 (LINC00538)], cg20390515 [syntrophin gamma 2 (SNTG2)], cg18085807 [ribonuclease P/MRP subunit p21 (RPP21)], cg14219242 [heterogeneous nuclear ribonucleoprotein H1 (HNRNPH1)], cg16702144 [outer dense fiber of sperm tails 2 (ODF2)], cg19350141 [leucine rich repeat containing 49 (LRRC49)]; 1/1 CpG with glyphosate: cg06950346 [uncharacterized LOC101928211 (LOC101928211)]; 1/1 CpG with atrazine: cg01867395 [paired box 6 (PAX6)]; 5/6 CpGs with heptachlor: cg05088513 [tubulin beta 6 class V (TUBB6)], cg11691038 [chromosome 12 open reading frame 40 (C12orf40)], cg18284070 [exocyst complex component 5 (EXOC5)], cg07627706 [chromosome 3 open reading frame 80 (C3orf80)], cg16310513 [ATPase phospholipid transporting 8B1 (ATP8B1)]; and 1/1 CpG with malathion: cg26210302 [cytochrome P450 family 8 subfamily B member 1 (CYP8B1)].

One CpG (cg10701801, associated with acetochlor) had a bimodal distribution and was removed from the functional enrichment analyses. Using eFORGE, we observed only tissue- and cell-specific enrichment with picloram. The 48 FDR-significant CpGs were enriched for DNase I hypersensitive sites in blood and fetal brain tissues (Figure S18). Across all active ingredients with significant CpGs, none were enriched for any of the chromatin states or histone markers. Using eFORGE-TF, we observed enriched transcription factor motifs for seven of the nine active ingredients with differentially methylated CpGs (number of transcription factor motifs: 1–35) (Table 2; Excel Table S12).

Association between Significant CpGs and Nearby Gene Expression

Of the 162 differentially methylated CpGs identified across the active ingredients, 89 were available in BIOS (Table 2) and 28 (31%) were associated with expression in cis. The 28 CpGs were associated with expression of nearby gene transcripts for 79 CpG-transcript pairs, of which 36 (46%) had an inverse relationship (i.e., higher methylation was associated with lower gene expression) (Excel Table S13).

We looked up the genes implicated by our differential methylation findings in databases with experimental studies. In the BaseSpace Correlation Engine (Kupershmidt et al. 2010), altered gene expression has been reported for the following pesticides: heptachlor (2/6 implicated genes) and malathion (1/1 implicated gene) (Excel Table S14). In that database, there is evidence that exposing mouse liver cells to heptachlor increases expression of TUBB6 and decreases expression of ATP8B1. Similarly, in the same database, there is evidence that malathion increases CYP8B1 expression in mouse lung tissue and decreases CYP8B1 expression in mouse liver tissue. None of the genes implicated in our study had been reported to be induced by the active ingredient in the Comparative Toxicogenomics Database or the U.S. EPA’s Chemical Dashboard.

We attempted to look up our FDR-significant results in the PEG Study (Paul et al. 2018), which has methylation data on 789 individuals assessed using the Illumina 450K array. Descriptive characteristics of the PEG Study are available in Excel Table S15. However, glyphosate was the only active ingredient with at least 15 exposed participants in the PEG Study; thus, we evaluated replication for only glyphosate. The CpG associated with glyphosate in the ALHS, cg06950346 (; ), was not nominally significant in the PEG Study, but the inverse direction of the association was consistent (128 exposed in the PEG Study, ; ).

The PEG Study previously reported differentially methylated CpGs in relation to organophosphate, 70 CpGs at () (Paul et al. 2018), and pyrethroid, 4 CpGs at (Furlong et al. 2020). We briefly describe how organophosphate and pyrethroid exposures were defined in these previous publications. For each specific organophosphate (from 1974 to the year of the blood draw) or pyrethroid pesticide (5 y prior to the blood draw), PEG Study participants were dichotomized on the basis of the median level of exposure estimated within a buffer around both residential and workplace addresses. Participants could have been exposed at neither, one, or both locations. For organophosphates, the number of specific organophosphates above the median was summed, and participants were exposed to an average of eight organophosphates (range: 0 to 46). For pyrethroids, the summed variable was further dichotomized to indicate any pyrethroid exposure above the median level (45 exposed).

We looked up the PEG Study CpGs in the ALHS. Comparing 92 past and current organophosphate users to 79 nonusers, we were able to look up 58 of the 70 CpGs reported as associated with organophosphate in the PEG Study. In the ALHS, of these 58 CpGs, only 3 (cg03655023, cg15083522, and cg19952704) had a and the same direction of association (Excel Table S16). For pyrethroid, we compared 78 past and current users to 666 nonusers and looked up 3 of 4 CpGs identified in the PEG Study; none had a . The epigenome-wide analyses of organophosphates and pyrethroids as a class in the ALHS did not identify any FWER- or FDR-significant CpGs.

Evaluating associations of the 16 active ingredients with the four epigenetic clocks, we observed only two significant associations (). Past use of DDT was associated with positive epigenetic age acceleration calculated using the Hannum clock { [95% confidence interval (CI): 0.34, 1.96]} (Hannum et al. 2013) and the Skin and Blood clock [ (95% CI: 0.07, 1.55)] (Horvath et al. 2018) (Table 4). DDT was banned in the United States in 1972. When we repeated the analyses and restricted it to 818 farmers who were in 1972, the magnitude of the association was slightly stronger across all four clocks but was only significant in the Hannum and Skin and Blood clocks (Table 4).

Table 4

Association between past DDT use and epigenetic age acceleration among male farmers of European ancestry in the Agricultural Lung Health Study, 2009–2013.

Epigenetic clock	Alla			At least 18 years of age in 1972b
	β	95% CI	p-Value	β	95% CI	p-Value
Horvath	0.44	−0.39, 1.26	0.30	0.67	−0.23, 1.58	0.14
Hannum	1.15	0.34, 1.96	0.006	1.26	0.36, 2.15	0.005
Skin and blood	0.81	0.07, 1.55	0.03	1.14	0.33, 1.95	0.006
PhenoAge	0.77	−0.26, 1.80	0.14	1.08	−0.05, 2.22	0.06

Note: Analyses conducted using linear regression, adjusting for the selection factor (asthma status), age at blood collection, pack years, smoking (never, former, current), state of residence, and seven estimated cell types. CI, confidence interval; DDT, dichlorodiphenyltrichloroethane.

Past use/never use () was 375/765.

Past use/never use () was 365/453.

Discussion

We observed differentially methylated CpGs in relation to nine pesticide active ingredients. All the differentially methylated CpGs were unique to an active ingredient, suggesting that there are pesticide-specific methylation signals. For example, dicamba and picloram are acetic acid pesticides, but none of the differentially methylated CpGs () identified in relation to dicamba or picloram overlapped. Similarly, we identified six differentially methylated CpGs with heptachlor but none with the other organochlorines we analyzed. Past use of DDT was associated with positive age acceleration that was based on the Hannum and Skin and Blood clocks.

As a validation of our findings, we looked up the implicated genes in pesticide–gene interaction databases. The most comprehensive database used was the BaseSpace Correlation Engine, which combines published and unpublished data sets from public, private, and proprietary sources (Kupershmidt et al. 2010). Across the nine active ingredients for which we observed differential methylation, there were up to 17 data sets available in this database. Given the paucity of specific pesticide–gene interaction studies, we were unable to confirm most implicated genes because they had not been previously examined. The Comparative Toxigenomics Database relies on published literature (Davis et al. 2019), and the U.S. EPA’s Chemical Dashboard contains only in vitro studies on a limited number of genes (Williams et al. 2017), which most likely explains why none of the implicated genes were confirmed in those two databases compared with the BaseSpace Correlation Engine.

In BIOS, about 30% of the differentially methylated CpGs were associated with the expression of nearby genes in blood. The proportion of CpGs associated with nearby gene expression is smaller than that previously reported in a study of differential methylation in relation to air pollution, a ubiquitous ambient exposure (Lee et al. 2019). The four cohorts in BIOS were population based and did not select for farmers. For genes that are inducible by pesticide exposure, levels of expression might be low or vary little in a nonfarming population that has no or minimal pesticide exposure. Thus, correlations between gene expression and methylation may be difficult to detect.

We also observed enriched transcription factor motifs across most pesticides; these may inform a biological linkage between pesticide exposure and health outcomes. For example, the EPA and the IARC considers heptachlor to be a possible human carcinogen. Heptachlor use has been associated with an increased risk of bladder and pancreatic cancer (Andreotti et al. 2009; Koutros et al. 2016). The differentially methylated CpGs identified with heptachlor were enriched for two transcription factor motifs related to cancer: HIC1 (a tumor repressor gene on chromosome 17) and MYB (an oncogene on chromosome 6). The United States banned the commercial sale of heptachlor in 1988, but heptachlor can linger in the environment. These findings need to be interpreted cautiously given that we had no suitable data to replicate these findings.

Because the ALHS is a case–control study of asthma, we included asthma in our model as a selection factor. Using a weighting strategy instead of adjusting for the selection factor can be preferred if the selection factor is a confounder or intermediate variable (Richardson et al. 2007). Using logistic regression to examine the association between each of the 16 active ingredients as defined in this study and asthma status, lindane was the only active ingredient where the -value was nominally significant (, data not shown), suggesting that asthma is not on the causal path between pesticide use and DNA methylation. Notably, we did not observe any significant CpGs with lindane. Furthermore, we ran an EWAS of asthma vs. no asthma, and looked up the 162 FDR-significant CpGs that were identified across 9 active ingredients (data not shown). Only one of the 162 CpGs differentially methylated in relation to an active ingredient was FDR significant in the asthma EWAS, which is far less than expected by chance (). These results strongly suggest that asthma is neither a confounder nor an intermediate variable.

We identified only one study with specific pesticide information and blood DNA methylation assessed in adults (Paul et al. 2018). That study, the PEG Study, was conducted in the Central Valley of California. Replicating our significant CpGs was challenging because of differences in the agricultural crops—and thus the pesticides—used, as well as pesticide exposure assessment. For example, common crops in Iowa and North Carolina include wheat, corn, sugarcane, and soybean, whereas the California Central Valley grows a variety of produce—including tomatoes, grapes, tree fruits, and nuts. The PEG Study had too few participants exposed to active ingredients commonly used by ALHS farmers to look up most of the CpGs we found. Additionally, ALHS farmers likely had more direct and higher exposure to pesticides compared with PEG Study participants given that ALHS farmers stored, mixed, and applied the pesticides. If exposure is not persistent, methylation at some CpGs might revert to unexposed levels, which has been well described for smoking (Joehanes et al. 2016). The CpG associated with glyphosate (cg06950346) was not significant in the PEG Study, but the direction of the association was consistent (). This CpG annotates to a long noncoding RNA, a class of RNAs that are thought to be highly inducible (Karlsson and Baccarelli 2016).

Farmers use multiple pesticides, and it may be of interest to understand the impact of pesticide mixtures. Analysis of chemical mixtures is a rapidly evolving research area. Current methods to analyze chemical mixtures such as Bayesian kernel machine regression and Bayesian g-computation are computationally intensive, even for a single outcome, and would not be practical to run epigenome-wide. Future development of mixtures methods and convenient code suitable for epigenome-wide analysis would be of interest. Thus, although we did not conduct a mixtures analysis, we attempted to replicate ( in the ALHS and the same direction of association) the differentially methylated CpGs that the PEG Study previously reported to be associated with organophosphate and pyrethroid, two common pesticide classes (Furlong et al. 2020; Paul et al. 2018). We replicated just 3 of the 70 organophosphate-related CpGs and none of their pyrethroid-related CpGs. It may be difficult to replicate the CpGs reported in the PEG Study because of differences in exposure type (occupational vs. community exposures), the exposure assessments methods used in these studies, and the mixtures of active ingredients being different for each study population. The low replication rate suggests that signals related to pesticide class may be difficult to identify and replicate because mixtures of active ingredients may likely be study specific.

Our findings suggest potential research questions to be explored. For example, cg01867395 was implicated with atrazine and annotated to PAX6. PAX6 is involved in the development of the central nervous system, eyes, and pancreas (Blake and Ziman 2014). In a systematic review, Goodman et al. (2014) identified four epidemiological studies that examined the association between maternal atrazine exposure and birth defects related to the nervous system. All four studies had null findings. Although there are no epidemiological studies of atrazine and eye defects, animal studies have reported abnormal eye development with exposure to atrazine (Scahill 2008; Wang et al. 2015). PAX6 is important for glucose metabolism, and disruption of PAX6 expression increases the risk of diabetes (Panneerselvam et al. 2019). Atrazine has been associated with gestational diabetes in the AHS (Saldana et al. 2007). There is also evidence of a positive, nonsignificant association with over 100 lifetime days of atrazine use and incident diabetes (Montgomery et al. 2008). We are unaware of other studies that have examined the association between atrazine and diabetes. Atrazine may be associated with eye defects and diabetes in humans through PAX6, but further research is warranted.

Another potential research question suggested by our results is the potential association between picloram and neurological outcomes. Using eFORGE, the differentially methylated CpGs associated with picloram were enriched in fetal brain tissue. Some of the genes implicated in relation to picloram play roles in neuronal biology or neurologic disease. For example, CUX1 (cut-like homeobox 1), glial cell line-derived neurotrophic factor family receptor alpha 1 (GFRA1), and protocadherin gamma subfamily B, 3 () are involved in neuron differentiation, survival, or dendritic arborization (Cubelos et al. 2010; Garrett et al. 2012; Irala et al. 2016), EN2 (engrailed homeobox 2) has been associated with autism (Chelini et al. 2019), HTR2A (5-hydroxytryptamine receptor 2A) encodes for a serotonin receptor and may be associated with obsessive-compulsive disorder (Noh et al. 2017). Additionally, the 48 CpGs identified with picloram were also enriched for transcription factor motifs involved in brain development [e.g., BRF1 (RNA polymerase III transcription initiation factor 90 kDa subunit) (Borck et al. 2015), ZFP423 (zinc finger protein 423) (Massimino et al. 2018), RFX4 (regulatory factor X4) (Loeliger et al. 2020), EOMES (eomesodermin) (de la Torre-Ubieta et al. 2018), and OTX1 (orthodenticle homeobox 1) (Huang et al. 2018)]. Notably, the top enriched transcription factor motif was BRF1, a constituent of RNA polymerase III, which suggests that picloram exposure may disrupt transcription. There is some evidence in mice that exposure to picloram damages primary neurons (Reddy et al. 2011), but we are unaware of epidemiologic studies on picloram exposure in relation to any cognitive outcome.

Our finding that past use of DDT was associated with positive epigenetic age acceleration confirms results from a previous study (Lind et al. 2018). That study measured plasma levels of ′-dichlorodiphenyldichloroethylene (DDE; the main metabolite of DDT) in a nonfarming Swedish population of 70-y-old adults and reported a positive association between ′-DDE and epigenetic age acceleration using the Hannum clock (Lind et al. 2018). Similar to our results, that study did not observe an association using the Horvath clock (Lind et al. 2018). They did not examine the association using the Skin and Blood (Horvath et al. 2018) or PhenoAge clocks (Levine et al. 2018) because their paper was published the same year as those two clocks. Together, these results suggest that DDT, banned in the United States in 1972 and in Sweden in 1970, might have global and lasting impacts across the epigenome that contribute to the acceleration of an individual’s biological age and, potentially, to adverse health outcomes associated with accelerated aging.

Our study had some limitations. There are few studies with DNA methylation and pesticide information to replicate our results. Because we restricted our analyses to male farmers of European ancestry from North Carolina and Iowa, our results may not be generalizable to women if sex-specific effects exist, men of other ancestries, or the general population. However, we analyzed only the men because they were nearly always the licensed pesticide applicator. Pesticide use was self-reported via questionnaire. AHS applicators have been shown to provide valid responses to decade of first use of a specific pesticide (Hoppin et al. 2002). In addition, a previous study showed that AHS applicators can reliably recall ever/never use of specific pesticides 1 y later (Blair et al. 2002). Together, these two publications suggest pesticide applicators can accurately report their exposure to specific pesticides. Self-reported data can capture long-term exposure better than current biological measurements, which are suitable only for measuring short-term exposure for most nonpersistent chemicals. However, misclassification of exposure and recall bias cannot be ruled out. We would expect this to be nondifferential with respect to methylation and thus generally be a source of bias toward the null. Our analyses did not consider the amount of chemical use, use of protective equipment, number of years of pesticide use, nor application methods. Given the limited number of differentially methylated CpGs identified, we did not perform pathway analyses. Although we ran several epigenome-wide analyses, significance was assessed for each active ingredient and not across all models, an approach that is commonly practiced in genome-wide association studies of multiple related phenotypes.

In this EWAS of blood DNA methylation and pesticide use, we identified specific differentially methylated CpGs in relation to nine pesticide active ingredients. If replicated, these findings may contribute to the development of biomarkers of chronic exposure to these different pesticides and improve our understanding of the specific biological impacts that these pesticides might have.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.