Methylome-wide association study provides evidence of particulate matter air pollution-associated DNA methylation.

BACKGROUND: DNA methylation (DNAm) may contribute to processes that underlie associations between air pollution and poor health. Therefore, our objective was to evaluate associations between DNAm and ambient concentrations of particulate matter (PM) ≤2.5, ≤10, and 2.5-10 μm in diameter (PM2.5; PM10; PM2.5-10).
METHODS: We conducted a methylome-wide association study among twelve cohort- and race/ethnicity-stratified subpopulations from the Women's Health Initiative and the Atherosclerosis Risk in Communities study (n = 8397; mean age: 61.5 years; 83% female; 45% African American; 9% Hispanic/Latino American). We averaged geocoded address-specific estimates of daily and monthly mean PM concentrations over 2, 7, 28, and 365 days and 1 and 12 months before exams at which we measured leukocyte DNAm in whole blood. We estimated subpopulation-specific, DNAm-PM associations at approximately 485,000 Cytosine-phosphate-Guanine (CpG) sites in multi-level, linear, mixed-effects models. We combined subpopulation- and site-specific estimates in fixed-effects, inverse variance-weighted meta-analyses, then for associations that exceeded methylome-wide significance and were not heterogeneous across subpopulations (P < 1.0 × 10-7; PCochran's Q > 0.10), we characterized associations using publicly accessible genomic databases and attempted replication in the Cooperative Health Research in the Region of Augsburg (KORA) study.
RESULTS: Analyses identified significant DNAm-PM associations at three CpG sites. Twenty-eight-day mean PM10 was positively associated with DNAm at cg19004594 (chromosome 20; MATN4; P = 3.33 × 10-8). One-month mean PM10 and PM2.5-10 were positively associated with DNAm at cg24102420 (chromosome 10; ARPP21; P = 5.84 × 10-8) and inversely associated with DNAm at cg12124767 (chromosome 7; CFTR; P = 9.86 × 10-8). The PM-sensitive CpG sites mapped to neurological, pulmonary, endocrine, and cardiovascular disease-related genes, but DNAm at those sites was not associated with gene expression in blood cells and did not replicate in KORA.
CONCLUSIONS: Ambient PM concentrations were associated with DNAm at genomic regions potentially related to poor health among racially, ethnically and environmentally diverse populations of U.S. women and men. Further investigation is warranted to uncover mechanisms through which PM-induced epigenomic changes may cause disease.

Introduction

Ambient particulate matter (PM) air pollution is a modifiable exposure that has been consistently associated with morbidity and mortality (Cohen et al., 2017; Di et al., 2017; Miller et al., 2007) attributed to cardiovascular disease (Brook et al., 2004, 2010), respiratory disease (Dominici et al., 2006; Gan et al., 2013; Laumbach and Kipen, 2012), and lung cancer (Pope et al., 2002; Raaschou-Nielsen et al., 2013). Despite the ubiquity of air pollution exposure and the continued population burden of PM (Cohen et al., 2017), the causal mechanisms underlying PM associations with poor health have not been adequately investigated.

One such mechanism could involve methylation of deoxyribonucleic acids (DNAm), conventionally measured at Cytosine-phosphate-Guanine (CpG) sites. DNAm is a heritable, but dynamic epigenetic modification that can influence gene expression without altering the DNA sequence (Clouaire and Stancheva, 2008; Neidhart, 2016) and may be central to mediation of PM-associated disease risk (Baccarelli et al., 2010; Bollati and Baccarelli, 2010; Zhong et al., 2016). Indeed, PM exposure has been implicated in whole blood DNAm near candidate genes involved in inflammation, oxidative stress, coagulation and vasoconstriction (Bellavia et al., 2013; Chen et al., 2015, 2016; Tarantini et al., 2009, 2013), abnormalities of which have established associations with cardiovascular and respiratory disease. A few studies have agnostically evaluated DNAm associations with PM on a methylome-wide scale (de F.C. Lichtenfels et al., 2018; Panni et al., 2016; Plusquin et al., 2017), but none have done so in large, sociodemographically and environmentally diverse, well-characterized populations of adult women and men.

The present study therefore examined methylome-wide associations between DNAm and ambient concentrations of PM ≤ 2.5, ≤10, and 2.5–10 μm in diameter (PM2.5, PM10, and PM2.5–10) within the Women’s Health Initiative (WHI) and the Atherosclerosis Risk in Communities study (ARIC) cohorts, and their replication in subpopulations of the Cooperative Health Research in the Region Augsburg (KORA) study.

Methods

Study design and populations

The study included 8397 consenting participants from subpopulations within the WHI and ARIC cohorts who had available peripheral blood leukocyte DNA.

The WHI is a multicenter prospective study of risk factors for cardiovascular disease (CVD), cancer, osteoporotic fractures, and other causes of morbidity and mortality among postmenopausal women (Anderson et al., 2003; NIH, 1998). Between 1993 and 1998, women aged 50–79 years from forty WHI clinical centers throughout the United States (US) were enrolled in the Clinical Trials (CT) (n = 68,132) or Observational Study (OS) (n = 93,676). All WHI participants completed a screening visit (SV). CT participants also completed an annual visit (AV) at one, three, six, and nine years after randomization (AV1, AV3, AV6, AV9), and OS participants three years after enrollment (AV3). An additional visit of CT and OS participant subsets occurred between 2011 and 2012 (ranging from 14 to 19 years after enrollment) as part of the WHI Long Life Study (LLS) (Anderson and LaCroix, n.d.).

For the current study, WHI participants were drawn from three ancillary studies: Epigenetic Mechanisms of PM-Mediated CVD Risk (WHI-EMPC) (Whitsel, n.d.), Broad Agency Announcement 23 (WHI-BAA23) (Assimes et al., n.d.) and Ancillary Study 311 (WHI-AS311) (Jordahl et al., 2018). WHI-EMPC is a study of epigenetic mechanisms underlying associations between ambient PM air pollution and CVD within the WHI CT. From this population, DNAm was measured in 2200 randomly selected participants (stage 1: SV, AV3, or AV6), remeasured in 200 participants at a second visit (stage 2: AV3 or AV6), and remeasured again in 43 participants at a third visit among those who participated in the WHI Long Life Study (stage 3: LLS), yielding 2443 total observations. WHI-BAA23, also known as Integrative Genomics and Risk of CHD and Related Phenotypes in the Women’s Health Initiative, is a case-control study of coronary heart disease within the WHI CT (n = 1546) and OS (n = 442). By design, WHI-BAA23 oversampled African Americans and Hispanic/Latino Americans and required all participants to have undergone genome-wide genotyping and profiling of seven cardiovascular disease biomarkers. DNAm was measured in blood collected at the SV, before the incidence of coronary heart disease. WHI-AS311 is a matched case-control study of bladder cancer among women within the WHI CT (n = 405) and OS (n = 455). Bladder cancer cases were matched to controls based on enrollment year, age at enrollment, follow-up time, and DNAm extraction method. DNAm was measured in blood collected at the SV, before the incidence of bladder cancer.

ARIC is a community-based prospective study of atherosclerosis and its clinical outcomes in four US communities: Washington County, Maryland; Forsyth County, North Carolina; selected suburbs of Minneapolis, Minnesota; and Jackson, Mississippi (ARIC Investigators, 1989). Enrollment in 1987–1989 (Visit 1) was followed by five sub-sequent visits (Visits 2–6) between 1990 and 2017. The present study included all 2796 African Americans from Forsyth County or Jackson (ARIC-AA) with DNA and 1139 European Americans from Forsyth County or Minneapolis (ARIC-EA) with cerebral magnetic resonance imaging data (Mosley et al., 2005), all at Visits 2 (1990–1992) or 3 (1993–1995).

Replication involved up to 2176 participants from two studies of the population-based KORA cohort: F3 (n = 464) and F4 (n = 1712). KORA F3 (2004–2005) and F4 (2006–2008) are follow-up studies of the KORA S3 and S4 cohort participants, including German nationals aged 25–74 years from Augsburg, Germany (Holle et al., 2005; Wichmann et al., 2005).

Particulate matter exposure estimation

The study focuses on three ambient particulate matter (PM) air pollutants, including two (PM2.5 and PM10) that are regulated under the Clean Air Act by the US Environmental Protection Agency (EPA) according to its National Ambient Air Quality Standards (NAAQS) (EPA, 2017).

PM exposures were estimated at all geocoded WHI and ARIC participant addresses (Whitsel et al., 2004, 2006) in the contiguous US since the baseline examinations using two exposure modeling approaches, both based on US EPA Air Quality System (AQS) monitoring data for PM10 (since 1987) and PM2.5 (since 1999). In the WHI, the median distance from geocoded participant addresses to PM10 and PM2.5 EPA monitors was 7.8 and 7.6 km. In ARIC, it was 4.8 and 7.2 km. Geocoded address-specific daily mean PM10 concentrations (μg/m3) were spatially estimated using national-scale, log-normal ordinary kriging. Exposure measurement error using kriging methods may yield misclassification and increase variance or bias associations (Alexeeff et al., 2014; Lee et al., 2012), therefore validity of the estimation was assessed, using standard cross-validation statistics: average prediction error (PE), standardized prediction error (SPE), root mean square standardized (RMSS), and standard error (SE). Observed values of PE and SPE near zero, RMSS near one, and RMS near SE have provided evidence of model validity (Liao et al., 2006, 2007).

Also, geocoded address-specific monthly mean concentrations (μg/m3) were spatiotemporally estimated using generalized additive mixed models and geographic information system-based predictors. Because EPA AQS monitoring data for PM2.5 were not widely available until 1999, spatiotemporal estimation also involved the log-transformed ratio of PM2.5 to predicted PM10 between 1987 and 1999. A five- or ten-fold, out-of-sample cross-validation of the estimates in which the squared Pearson correlation between excluded monthly observations and model predictions (R2 = 0.68–0.77) indicated that estimation models performed well (Yanosky et al., 2014).

Daily mean concentrations of PM10 were averaged over the 2-, 7-, 28-, and 365-day periods ending on (including) the examination day. Monthly mean concentrations of PM2.5 and PM10 were averaged over the 12-month period ending on (including) the calendar month of ex-amination. Finally, coarse PM (PM2.5–10) concentrations for each averaging duration were calculated as differences between PM10 and PM2.5 concentrations.

DNA methylation

Peripheral blood leukocytes were isolated from visit-specific, fasting blood drawn from study participants. DNA was extracted from the peripheral blood leukocytes and then DNAm was measured on a methylome-wide scale at 485,577 CpG sites using the Illumina 450K Infinium Methylation BeadChip (Illumina Inc.; San Diego, CA, USA). Methylation was quantitatively represented by beta, the proportion of methylated cytosines over the sum of methylated and unmethylated cytosines across the same loci. The data from all studies were quality controlled (Table S1), Beta Mixture Quantile (BMIQ)-normalized to adjust for probe bias (Teschendorff et al., 2013), and in WHI-EMPC, ComBat-adjusted for stage and plate using empirical Bayes methods (Johnson et al., 2007). Otherwise, technical covariates (assay plate, chip, and row) were available to control for batch effects; and leukocyte proportions (CD8+ T cell, CD4+ T cell, B cell, natural killer cell, monocyte, and granulocyte) to account for leukocyte composition (Houseman et al., 2012). Among ARIC-AA participants, missing lymphocyte, monocyte, neutrophil, eosinophil, and basophil proportions were imputed based on measured proportions. Analyses excluded CpG sites at which DNAm distributions were multi-modal (Andrews et al., 2016) in at least one study.

Multiple imputation

To avoid potential for selection bias in complete-data analysis when data are missing at random (Hernan et al., 2004), multivariate imputation by chained equations (MICE) (Azur et al., 2011; Stuart et al., 2009) as implemented in SAS 9.3 (Cary, NC) was used to impute in-frequently missing PM2.5, PM10, and PM2.5–10 concentrations (missing range: 3.3%, 3.5%) and other covariates (missing range: 0%, 10.4%), excluding methylome-wide DNAm. Binary and categorical data were imputed using the logistic and discriminant functions whereas interval-scale data were imputed using predictive means matching with a k-nearest neighbor (k = 5) approach.

Statistical analysis

All analyses were stratified by cohort and race/ethnicity (African-, European-, and Hispanic/Latino-American) and adjusted for age (years) at blood draw, education (high school education or lower, more than high school), smoking status (current, former, never), alcohol use (current, former, never), physical activity (metabolic equivalent of task [MET-hours/week]), body mass index (BMI, kg/m2), neighborhood socioeconomic status (Roux et al., 2001), mean temperature (°C), mean dew point (°C), mean barometric pressure (kPa), season, and methylation-related variables, which included ten principal components (PCs) for genetic ancestry (when available), leukocyte proportions, and technical covariates. Analyses additionally controlled for cohort-specific covariates, including binary sex (male, female) in ARIC; randomly assigned treatment group (CT subpopulations of WHI-AS311, WHI-BAA23, WHI-EMPC); case-control status (WHI-AS311, WHI-BAA23); and control matching criteria (WHI-AS311).

In each subpopulation, covariate-adjusted, multi-level, linear, mixed-effects models (LMMs) were used to estimate DNAm-PM associations. In WHI-EMPC, three-level, longitudinal models had a random intercept for examination at the participant level, a random intercept and slope for PM at the WHI center level, and a random intercept for chip, as given by

In WHI-BAA23 CT & OS, and WHI-AS311 CT & OS, two-level cross-sectional models had a random intercept and slope for PM at the WHI center level and a random intercept for plate and chip, as given by

In ARIC-AA and ARIC-EA, one-level cross-sectional models had a random intercept for plate and chip, as given by

Above, i, j and k denote the i examination of the j participant in the k center; DNAm is the CpG site-specific beta value; β0 is the intercept; PM is the 2-, 7-, 28-, 365-day, or 1- or 12-month mean of PM2.5, PM10, or PM2.5–10; and Z is a vector of covariates. The terms are a random intercept and a random slope for PM at the center level, is a random intercept for examination at the participant level, are random intercepts for technical covariates, and is the random error at the examination level. Measures of association (β1) and their 95% confidence intervals (β1 ± 1.96 × standard error) were reported as an absolute percentage change in DNAm per 10 μg/m3 increase in PM.

Given the focus on fixed effects, LMMs were fit with maximum likelihood using the MixedModels package (Bates, 2017) in Julia v0.6 (Bezanson et al., 2017). Stratum-specific results were combined using fixed-effects, inverse-variance weighted meta-analysis. Homogeneity of associations was assessed using Cochran’s Q test statistic (Cochran, 1954). A P < 0.10 and Bonferroni-corrected threshold of P < 1 × 10−7 (i.e. assuming 500,000 independent CpG tests) were used to identify significant CpG associations. The threshold of suggestive significance was P < 1 × 10−5.

Examination of stratified and meta-analyzed results included reviewing quantile-quantile (QQ) plots of the observed −log10-transformed P values for each CpG site against the expected values from a theoretical χ2 distribution and estimating the associated genomic inflation factor (λ), where λ is defined as the ratio of the observed to expected median −log10P values (Devlin et al., 2001).

Technical validation

In a random subset of 200 WHI-EMPC participants, bisulfite pyrosequencing was used to validate the Illumina 450K measures of DNAm at ten PM10- or PM2.5-sensitive CpG sites (P < 1 × 10−5). CpG sites with poor next generation sequencing data or situated in CpG-rich, repetitive element, or low sequence complexity regions of the genome were not candidates for pyrosequencing. Site-specific comparisons of DNAm measures were based on mean Illumina 450K minus bisulfite pyrosequencing differences (Δ), Pearson correlation coefficients (r), and Deming regression estimates of their intercepts (α) and slopes (β) (Cornbleet and Gochman, 1979). When the two measures are nearly identical, Δ, r, α, and β approach values of 0, 1, 0, and 1, respectively.

Functional annotation

Published genotype-phenotype associations for variants annotated to or within 100 kilobases of genes containing statistically significant PM-sensitive CpG sites were identified in the National Human Genome Research Institute (NHGRI) Genome-Wide Association Study (GWAS) Catalog (Welter et al., 2014). Tissue-specific gene expression was assessed using the Genotype-Tissue Expression (GTEx) database (Lonsdale et al., 2013) and associations between DNAm and gene expression in human blood cells were obtained from a study of approximately 400,000 CpG sites and > 13,000 transcripts in the Multi-Ethnic Study of Atherosclerosis (MESA) and Grady Trauma Project (GTP) cohorts (Kennedy et al., 2018). PM-sensitive CpG sites (P < 1 × 10−5) were functionally characterized using experimentally derived Functional element Overlap analysis of ReGions from EWAS (eFORGE) v2.0 (Breeze et al., 2016) with data from the Encyclopedia of DNA elements (ENCODE) (Consortium, 2012), Roadmap Epigenomics Project (Bernstein et al., 2010), and BLUEPRINT (Stunnenberg et al., 2016). Overlap of CpG site-specific PM sensitivity, histone modification, and DNase I hypersensitivity were evaluated in eFORGE with a false discovery rate (FDR) threshold of 0.05.

Replication

Significant CpG sites that were not heterogeneous across sub-populations (P < 1.0 × 10−7; P > 0.10) underwent replication and meta-analyses in KORA F3 and F4. Pollutant-and averaging duration-specific replication thresholds were Bonferroni-corrected by dividing the conventional alpha level (0.05) by the number of CpG sites carried into replication.

Results

The study consisted of twelve ARIC and WHI subpopulations, collectively representing 8397 participants, of whom 45.8% were African American, 8.4% were Hispanic/Latino American, and 83.0% were female (Table 1). Participants were on average 61.3 years of age and contributed methylation data at ≥461,014 CpG sites. One-month mean concentrations of PM10, PM2.5, and PM2.5–10 were 20.9, 13.2, and 7.7 μg/m3; varied by subpopulation and race/ethnicity (Tables 1 and S2); and did not exceed NAAQS in place at the time of data collection. Between-pollutant Pearson correlation coefficients depended on size fraction and averaging duration (Table 2). Overall, the median (range) was 0.35 (−0.14, 0.79) and among 2-, 7-, 28, and 365-day mean PM10 concentrations, it was 0.64 (0.43, 0.79). Correlations between PM10 and PM2.5 concentrations were 0.73 and 0.64 when they were averaged over 1 and 12 months.

Table 1

Characteristics of the study participants, by subpopulation.

Subpopulation			Race/ethnicity	n	% female	Age, yrs	Maximum CpGs	PM (μg/m3), 1 mo x¯ (SD)
						x¯ (SD)		PM10	PM2.5	PM2.5–10
ARIC			AA	2664	63%	56.6 (5.9)	463,431	20.5 (4.6)	13.2 (3.1)	7.3 (2.1)
			EA	1100	58%	59.9 (5.4)	462,543	23.2 (5.3)	15.4 (4.3)	7.8 (3.5)
WHI	AS311	CT	EA	351	100%	64.7 (7.1)	461,136	19.8 (6.6)	11.9 (3.82)	7.9 (4.6)
		OS	EA	395	100%	66.2 (6.9)	461,136	19.9 (5.7)	12.0 (3.9)	7.9 (4.1)
	BAA23	CT	AA	371	100%	61.8 (6.3)	461,014	22.6 (6.2)	14.3 (4.2)	8.3 (3.8)
			EA	926	100%	67.8 (6.2)	461,014	19.7 (5.7)	11.7 (3.7)	8.0 (4.4)
			HLA	220	100%	60.7 (6.4)	461,014	21.4 (8.1)	10.3 (4.1)	11.1 (5.7)
		OS	AA	259	100%	62.8 (6.8)	461,014	22.3 (5.9)	14.0 (4.0)	8.3 (4.2)
			HLA	174	100%	62.8 (7.3)	461,014	23.0 (8.1)	11.0 (4.2)	11.9 (6.4)
	EMPC[a]		AA	553	100%	62.7 (6.9)	463,916	22.2 (6.2)	15.2 (5.1)	7.0 (4.7)
			EA	1072	100%	64.6 (7.1)	463,916	19.4 (6.0)	13.0 (5.0)	6.4 (5.2)
			HLA	312	100%	61.5 (6.1)	463,916	21.9 (7.1)	12.8 (6.3)	9.1 (5.3)
All			AA (45.8%)
			HLA (8.4%)	8397	83%	61.3 (7.4)	463,916	20.9 (5.8)	13.2 (4.3)	7.7 (4.0)
			EA (45.8%)

Abbreviations: AA, African American; ARIC, Atherosclerosis Risk in Communities; AS311, Ancillary Study 311; BAA23, Broad Agency Award 23; CpG, Cytosine-phosphate-Guanine; CT, Clinical Trial; EA, European American; EMPC, Epigenetic Mechanisms of PM-Mediated CVD Risk; HLA, Hispanic/Latino American; mo, month; OS, Observational Study; PM10, PM < 10 μm in diameter; PM2.5, PM < 2.5 μm in diameter; PM2.5–10, PM > 2.5 and < 10 μm in diameter; SD, standard deviation; WHI, Women’s Health Initiative; , mean.

At the 1st visit. Methylation data also were available among 185 & 43 WHI-EMPC participants @ the 2nd & 3rd visits

Table 2

Particulate matter concentration (μg/m3) means and Pearson correlations in the total population (n = 8397).

		PM10	PM10	PM10	PM10	PM10	PM10	PM2.5	PM2.5	PM2.5–10	PM2.5–10
		2 d	7 d	28 d	365 d	1 mo	12 mo	1 mo	12 mo	1 mo	12 mo
	x¯ (SD)	31.9 (12.1)	31.1 (9.2)	30.9 (7.1)	31.2 (5.1)	20.9 (5.8)	20.9 (4.0)	13.2 (4.3)	13.2 (3.0)	7.7 (4.0)	7.8 (3.1)
PM10	2 d	1.00
PM10	7 d	0.74	1.00
PM10	28 d	0.58	0.79	1.00
PM10	365 d	0.43	0.56	0.70	1.00
PM10	1 mo	0.39	0.48	0.54	0.27	1.00
PM10	12 mo	0.15	0.18	0.24	0.35	0.62	1.00
PM2.5	1 mo	0.29	0.36	0.41	0.17	0.73	0.39	1.00
PM2.5	12 mo	0.11	0.12	0.15	0.23	0.40	0.64	0.66	1.00
PM2.5–10	1 mo	0.25	0.31	0.35	0.21	0.67	0.48	−0.02	−0.13	1.00
PM2.5–10	12 mo	0.08	0.12	0.17	0.23	0.41	0.67	−0.14	−0.14	0.74	1.00

Abbreviations: d, day; mo, month; PM, particulate matter; PM10, PM < 10 μm in diameter; PM2.5, PM < 2.5 μm in diameter; PM2.5–10, PM > 2.5 and < 10 μm in diameter; SD, standard deviation; , mean.

QQ plots (Fig. 1) based on the trans-ethnic, fixed-effects, inverse variance-weighted meta-analyses provided little evidence of inflation across pollutants and averaging durations: median (range) λ = 1.01, (0.89–1.07). Manhattan plots (Fig. 2) show three significant (P < 1 × 10−7) and 55 suggestively significant (1 × 10−5 < P < 1 × 10−7) PM-sensitive CpG sites (Tables 3 and S3). The three significant CpG sites (cg19004594; cg24102420; cg12124767) were neither within ten base pairs of single nucleotide polymorphisms (minor allele frequency > 1%) nor previously identified as cross-reactive probes (Chen et al., 2013).

Fig. 1.

Quantile-quantile (QQ) plot of observed vs. expected −log10 P-value of each CpG site from trans-ethnic, fixed-effects meta-analyses of 2-, 7-, 28-, and 365-day PM10 and 1- and 12-month PM10 and PM2.5. The red diagonal line references the methylome-wide significance threshold (P < 1.0 × 10−7). Lambda (λ) is the inflation factor. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2.

Manhattan plot of −log10 P-value vs. chromosomal position of each CpG site from trans-ethnic, fixed-effects meta-analyses of 2-, 7-, 28-, and 365-day PM10 and 1- and 12-month PM10 and PM2.5. The red line references the methylome-wide significance threshold (P < 1.0 × 10−7). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 3

Findings from trans-ethnic, fixed-effects meta-analyses (P < 1 × 10−7, P > 0.10).

Chr	Position[a]	CpG	Exposure	%Δ (95% CI)[b]	P	nobs	Gene
20	43,926,884	cg19004594	PM10, 28 d	0.3 (0.2, 0.4)	3.33 × 10−8	8622	MATN4
3	35,785,890	cg24102420	PM10, 1 mo	−0.5 (−0.7, −0.3)	5.84 × 10−8	8575	ARPP21/miR128–2
7	117,299,297	cg12124767	PM2.5–10, 1 mo	−0.5 (−0.7, −0.3)	9.96 × 10−8	8577	CFTR

Abbreviations: Δ, change; Chr, chromosome; CI, confidence interval; CpG, Cytosine-phosphate-Guanine; d, days; mo, month; PM10, PM < 10 μm in diameter; PM2.5, PM < 2.5 μm in diameter; PM2.5–10, PM > 2.5 and < 10 μm in diameter.

Build 37.

Absolute percentage point per 10 μg/m3 increase in PM10.

On chromosome 20 within an exonic CpG island of MATN4, a 10 μg/m3 increase in 28-day mean PM10 was associated with a 0.3% (95% confidence interval [CI]: 0.2, 0.4) higher DNAm at cg19004594 (P = 3.33 × 10−8; Fig. 3A). On chromosome 3 intronic to ARPP21, a 10 μg/m3 increase in 1-month mean PM10 was associated with a 0.5% (95% CI: 0.3, 0.7) lower DNAm at cg24102420 (P = 5.84 × 10−8; Fig. 3B). Cg24102420 is approximately 200 base pairs upstream from the transcriptional start site for microRNA 128–2 (miR128–2). On chromosome 7 intronic to CFTR, a 10 μg/m3 increase in 1-month mean PM2.5–10 was associated with a 0.5% (95% CI: 0.3, 0.7) lower DNAm at cg12124767 (P = 9.86 × 10−8; Fig. 3C). Furthermore, PM associations with cg19004594, cg24102420, and cg12124767 were similar across race/ethnic strata (Fig. S1). Complete annotations for all PM-sensitive CpG sites (P < 1 × 10−7) are available in Excel Table S1.

Fig. 3.

Forest plots of PM-CpG associations (95% confidence intervals) for A) cg19004594, B) cg2410240, and C) cg12124767 with a 10 μg/m3 increase in PM by subpopulation and overall after fixed-effects meta-analysis.

Overall, bisulfite pyrosequencing and Illumina 450K-based DNAm measures were similar (Table S4). The medians (interdecile ranges) of Δ, r, α and β were: 0.01 (−0.06, 0.07), 0.73 (0.20, 0.83), 0.04 (−0.27, 0.24), and 0.98 (0.09, 1.62). Corresponding estimates (95% CIs) for cg24102420 were −0.04 (−0.04, −0.03), 0.79 (0.73, 0.83), −0.16 (−0.38, 0.07) and 1.13 (0.88, 1.39). Cg19004594 and cg12124767 were not pyrosequenced.

MATN4 is highly expressed in the pancreas, reproductive tract, and skin (Fig. S2), but variants of this gene have not been significantly associated (P < 5 × 10−8) with any phenotypes in prior GWAS. ARPP21 is primarily expressed in the brain (Fig. S3), is significantly associated with neuroticism and severe H1N1 influenza, and suggestively associated (5 × 10−8 < P < 5 × 10−6) with entorhinal cortical thickness and childhood-onset asthma in prior GWAS. CFTR is expressed in various tissues, including the pancreas, colon, minor salivary gland, digestive tract, and lung (Fig. S4). CFTR polymorphisms are associated with cystic fibrosis (CF), Barrett’s esophagus/esophageal carcinoma, and coronary artery disease.

Differential methylation at cg19004594, cg24102420, or cg12124767 was not associated with gene expression in blood cells at any of the > 13,000 transcripts evaluated (P > 10−5) in the MESA/GTP cohorts. Although genomic regions around PM-sensitive CpG sites were associated with tri-methylation of histone 3 at lysine 9 (H3K9me3) in natural killer cells, derived mesenchymal stem cells, the fetal adrenal gland, fetal lung fibroblasts, and foreskin fibroblasts (FDR < 0.05; Fig. 4), they were not associated with mono- or tri-methylation of histone 3 at lysine 4, 27, or 36 (H3K4me1, H3K4me3, H3K27me3, or H3K36me3) or DNase I hypersensitivity in any tissues catalogued by eFORGE.

Fig. 4.

Enrichment of PM-sensitive CpG sites in regions overlapping H3K9me3 using Roadmap data.

The three statistically significant, non-heterogeneous PM-sensitive CpG sites (cg19004594; cg24102420; cg12124767) did not replicate in KORA F3/F4 (Table S5).

Discussion

This methylome-wide association study (MWAS) discovered three CpG sites at which higher levels of monthly mean ambient particulate matter air pollution concentrations were associated with DNAm. The DNAm-PM associations at all three CpG sites were homogeneous across the twelve subpopulations and each site was annotated to a neurological, pulmonary, endocrine, or cardiovascular disease-related gene (MATN4, ARPP21 or CFTR). Although a recent MWAS also implicated cigarette smoking in DNA methylation at ARPP21 and CFTR (Joehanes et al., 2016)—two genes that may underlie epigenetically mediated responses to inhalable environmental exposures—the CpG sites discovered herein are in different regions of ARPP21 and CFTR, suggesting varied responses to particulate exposures, and none of them were associated with gene expression of blood cells in MESA/GTP.

Methylation of cg19004594 (exon of MATN4) was positively associated with 28-day mean PM10 concentrations. MATN4 encodes Matrilin 4, a von Willebrand factor A domain-containing protein, which contributes to cardiac remodeling (Barallobre-Barreiro et al., 2012) and inhibits the proliferation of hematopoietic stem cells at rest. Additionally, environmental stressors trigger expression of the CXCL12- encoded chemokine (SDF1) (Liberda et al., 2010) and activation of its G protein-coupled receptor (CXCR4), leading to inhibition of Matrilin 4 and subsequent expansion of hematopoietic stem cell pools (Uckelmann et al., 2016). SDF1-activated CXCR4 also inhibits beta-adrenergically activated calcium influx through myocardial L-type calcium ion channels (Pyo et al., 2006), a process that may affect PM10-associated ventricular action potential and electrocardiographic QT interval duration (Gondalia et al., 2017). Methylation of MATN4 may therefore underlie commonly observed hematological and electrocardiographic of effects of PM10.

Methylation at cg24102420 (intron of ARPP21) was positively associated with 1-month mean PM10 concentrations. ARPP21 encodes a neuronal cAMP-regulated phosphoprotein, a regulator of calmodulin signaling (RCS) that is highly enriched in medium spiny neurons within the basal ganglia, cerebral cortex, and other regions of the brain (Rakhilin et al., 2004), with dual evidence of expression in cardiac tissues (Kahr et al., 2011; Kirchhof et al., 2011; Mathar et al., 2013). Variants of ARPP21 have been associated with entorhinal cortical thickness (Furney et al., 2010). Calmodulin signaling (O’Day et al., 2015), entorhinal cortical thickness (Velayudhan et al., 2013), and PM air pollution (Cacciottolo et al., 2017) are all associated with Alzheimer’s disease progression, suggesting a potential epigenetic mechanism of PM10-related neuropathology.

Indeed, ARPP21 and miR128–2, a microRNA within ARPP21, are both regulators of dendritic growth (Rehfeld et al., 2018). In a study of rats, exposure to ammonium sulfate, a major component of PM2.5, was associated with diminished dendritic complexity in hippocampal neurons (Cheng et al., 2017). Additionally, miR128 expression in peripheral blood of steel plant workers increased with increases in PM exposure, as was confirmed by an in vitro study of PM-treated pulmonary tissue (Bollati et al., 2015). Additional roles of miR128 include the inhibition of ABCA1 and ABCG1, adenosine triphosphate-binding cassette (ABC) transporter genes also involved in homeostasis of cholesterol (Adlakha et al., 2013), an established risk factor for stroke, myocardial infarction, and other common forms of cardiovascular disease.

Methylation at cg12124767 (intron of CFTR) was inversely associated with 1-month mean PM2.5–10 concentrations. CFTR encodes a transmembrane conductance regulator; specifically, an ABC transporter of chloride and thiocyanate ions. The CFTR-encoded ABC transporter controls fluid secretion and absorption in epithelial tissues (Saint-Criq and Gray, 2017). Its most common mutation impairs folding and trafficking of the encoded protein in pulmonary and pancreatic epithelia, causing CF and CF-related diabetes (Brennan et al., 2004). However, cigarette smoke and chronic inflammation also reduce CFTR chloride channel function (Rasmussen et al., 2014), a hypothesized molecular pathway underlying the development of chronic obstructive pulmonary disease (Rab et al., 2013). Furthermore, CFTR chloride channel currents in the myocardium shorten action potential and QT interval duration (Duan, 2013). Their activation by cAMP protein kinase A (PKA), protein kinase C (PKC), or extracellular adenosine triphosphate (ATP) through purinergic receptors (al-Awqati, 1995; Duan, 2013) can be arrhythmogenic (Cacciapuoti et al., 1991; Engler and Yellon, 1996; Leonard et al., 2017; Najeed et al., 2002; Yamazaki and Hume, 1997). Hypomethylation of CFTR at this site therefore highlights another epigenetic mechanism that may underlie PM10-related pulmonary and electrocardiographic manifestations of disease.

While the putative mechanisms described above are biologically plausible, analyses on which they are based are limited by their reliance on DNAm derived from leukocytes. Although other (e.g. heart, lung, nervous) tissues may be more appropriate for studying the role of DNAm on human disease, their collection is highly invasive (McCullough et al., 2017; Zhong et al., 2016); as such, leukocytes extracted from peripheral blood are widely used surrogate tissues (Zhong et al., 2016) with demonstrated consistency of DNAm patterns across relevant tissues types (Byun et al., 2009; Fan and Zhang, 2009; Ma et al., 2014). Still, DNAm at cg19004594, cg24102420, cg12124767 was not associated with gene expression of blood cells in GTP/MESA (Kennedy et al., 2018). Unlike DNAm patterns though, gene expression is highly variable by tissue type (Aguet et al., 2017), and MATN4, ARPP21 and CFTR are primarily expressed in other tissues.

The inability to replicate associations in KORA F3 and F4 participants is noteworthy. Although independent from the discovery populations, KORA represents a population of white, European men and women living in Augsburg, Germany, one distinct from that of the environmentally diverse, multi-racial/ethnic U.S. populations in the discovery. In addition, PM composition in ARIC and WHI (1990–2012) may differ from that in Augsburg during KORA F3 and F4 (2004–2006). Furthermore, PM concentrations in KORA were measured at community monitors, while those in WHI and ARIC were spatially or spatiotemporally estimated at participant geocoded addresses from monitoring networks in the 48 contiguous US states.

DNAm associations with PM2.5 – potentially the driver for PM-associated disease (Brook et al., 2010) – were not detected in this study. Inability to do so may be due to lower power to detect PM2.5 versus PM10 associations with DNAm given lower-variance PM2.5 exposure estimates, lack of short-duration PM2.5 data before 1999 when EPA AQS started monitoring it, and/or induction of PM2.5 health effects that are not epigenetically mediated.

The analyses also were limited by predominantly cross-sectional data, high multiple testing burden, small effect sizes, and residual need for functional characterization. However, repeated measures of PM and DNAm over time were leveraged in WHI-EMPC to increase statistical power. Among-pollutant correlations also were moderate in this context, so the multiple comparisons made were not strictly independent. Similarly, the Bonferroni-corrected threshold used herein (P < 1 × 10−7) is conservative because of methylome-wide correlations among CpG sites (Saffari et al., 2018; Tsai et al., 2012), decreasing the likelihood of false positives. Moreover, observed effect sizes were consistent with those seen in other epigenetic studies of particulate matter exposure (de F.C. Lichtenfels et al., 2018; Panni et al., 2016; Plusquin et al., 2017) and smoking (Joehanes et al., 2016). Further investigation is nonetheless needed to determine the clinical impact of CpG-specific changes in methylation although functional validation of epigenetic associations was outside the scope of presently funded work. Still, this is a well-powered study of geographically diverse, multi-racial/ethnic populations of women and men with methylome-wide DNAm and geocoded address-specific PM data, that leveraged multi-variate imputation to minimize selection-related biases otherwise known to affect epidemiologic associations in complete data analyses.

Conclusions

Findings from this large, racially/ethnically and environmentally diverse methylome-wide association study of women and men in EPA regions 1–10 suggest that ambient particulate matter air pollution affects DNAm at regions of the genome potentially related to neurological, pulmonary, endocrine, and cardiovascular disease. Although the discovered associations are biologically plausible, functional characterization in relevant tissues or animal models remains necessary to validate associations and elucidate putative epigenetic mechanisms of PM-associated disease.