Short- and intermediate-term exposure to ambient fine particulate elements and leukocyte epigenome-wide DNA methylation in older men: the Normative Aging Study.

BACKGROUND: Several epigenome-wide association studies (EWAS) of ambient particulate matter with aerodynamic diameter ≤ 2.5 µm (PM2.5) have been reported. However, EWAS of PM2.5 elements (PEs), reflecting different emission sources, are very limited.
OBJECTIVES: We performed EWAS of short- and intermediate-term exposure to PM2.5 and 13 PEs. We hypothesized that significant changes in DNAm may vary by PM2.5 mass and its elements.
METHODS: We repeatedly collected blood samples in the Normative Aging Study and measured leukocyte DNA methylation (DNAm) with the Illumina HumanMethylation450K BeadChip. We collected daily PM2.5 and 13 PEs at a fixed central site. To estimate the associations between each PE and DNAm at individual cytosine-phosphate-guanine (CpG) sites, we incorporated a distributed-lag (0-27 d) term in the setting of median regression with subject-specific intercept and examined cumulative lag associations. We also accounted for selection bias due to loss to follow-up and mortality prior to enrollment. Significantly differentially methylated probes (DMPs) were identified using Bonferroni correction for multiple testing. We further conducted regional and pathway analyses to identify significantly differentially methylated regions (DMRs) and pathways.
RESULTS: We included 695 men with 1,266 visits between 1999 and 2013. The subjects had a mean age of 75 years. The significant DMPs, DMRs, and pathways varied by to PM2.5 total mass and PEs. For example, PM2.5 total mass was associated with 2,717 DMPs and 10,470 DMRs whereas Pb was associated with 3,173 DMPs and 637 DMRs. The identified pathways by PM2.5 mass were mostly involved in mood disorders, neuroplasticity, immunity, and inflammation, whereas the pathways associated with motor vehicles (BC, Cu, Pb, and Zn) were related with cardiovascular disease and cancer (e.g., "PPARs signaling").
CONCLUSIONS: PM2.5 and PE were associated with methylation changes at multiple probes and along multiple pathways, in ways that varied by particle components.

Introduction

Exposure to ambient fine particulate matter (an aerodynamic diameter ≤2.5 µm; PM2.5), have been linked with death (Kloog et al. 2013) and multi-systemic diseases (Brook et al. 2010; Gu et al. 2020; Hunt et al. 2003; Nassan et al. 2021; Turner et al. 2011), in short-, intermediate-, and long-term time windows. Moreover, PM2.5 is a complex mixture of many particulate elements (PEs) that differ in their physicochemical, toxicological properties (Council 2004), reflecting specific emission sources. Emerging epidemiological studies have found that PM2.5-related health effects vary by different PEs (Dai et al. 2016; Franklin et al. 2008; Wang et al. 2020b; White et al. 2019). However, the underlying molecular alterations caused by PM2.5 and its PEs have not been adequately investigated.

DNA methylation (DNAm), a chemical modification of DNA with a methyl group addition predominantly at a cytosine-phosphate-guanine (CpG) site (Robertson and Jones 2000), has been reported to be associated with human health, such as inflammation (Gonzalez-Jaramillo et al. 2019), oxidative stress (Menezo et al. 2016), aging (Gonzalo 2010), cardiovascular disease (Pepin et al. 2019), and cancer (Wang and Lei 2018). Meanwhile, DNAm has been linked with PM2.5 across different time windows (Eze et al. 2020; Gondalia et al. 2019; Gruzieva et al. 2019; Jiang et al. 2014; Li et al. 2018; Panni et al. 2016; Plusquin et al. 2017). For instance, Gondalia et al. assessed epigenome-wide association studies (EWAS) of monthly mean PM2.5 before health examinations across twelve elderly cohorts and found three suggestive significant CpGs by random-effects meta-analysis (Gondalia et al. 2019). Panni et al. studied EWAS of PM2.5 in both short- and intermediate-term time windows across three elderly cohorts. The random-effects meta-analysis observed 12 significant CpGs (Panni et al. 2016).

Since PM2.5 is a combination of different PEs with different characteristics, these PEs may be associated with DNAm at different sites. To date, only five studies have assessed the associations between PEs and DNAm, among which three (Baccarelli et al. 2009; Chen et al. 2015; Hou et al. 2014) focusing on short-term and one (Madrigano et al. 2011) focusing on intermediate-term effects of PEs on specific loci or repetitive elements. And only one (Dai et al. 2017) examined the long-term effects of PEs on DNAm in an epigenome-wide scope, which was conducted by our group previously. We investigated EWAS of one-year moving averages of PEs and observed 20 significant CpGs for iron, 8 for nickel, and 1 for vanadium (Dai et al. 2017). Pathway analysis also suggested specific effects on DNAm for each PE. For example, iron and nickel were associated with type II diabetes mellitus and insulin signaling pathways. However, we did not examine the associations of PM2.5 total mass and DNAm. And linear mixed-effect models that we used are not optimal for DNAm data, which is usually non-normally distributed. Additionally, studies examining whether the associations between PM2.5/PEs and DNAm vary by different time windows remain quite limited. To date, no EWAS of short- and intermediate-term exposure to PEs has been performed.

The present study, therefore, sought to investigate the associations between short- and intermediate-term exposure to PM2.5/PEs and DNAm by conducting EWAS and quantile regression analyses, using whole-blood samples in the same population. We hypothesized that the changes in DNAm varied by PM2.5 and the PEs.

Methods

Study population

The NAS is a closed and ongoing cohort established in 1963 by the U.S. Veterans Administration in the Greater Boston Area (Bell et al. 1966). The participants were aged 21–82 years and were free of any known chronic diseases at enrollment. They have undergone health examinations in a clinical center, including blood collection, every 3–5 years.

In the present study, we included subjects who lived in the Greater Boston Area during the study period and had visits with DNA samples collected starting in 1999 (Wang et al. 2020b). We dropped non-white participants (~3%) to diminish study heterogeneity that may be introduced by diverse genetic ancestry (Wang et al. 2020b). All study participants provided written informed consent before enrollment and sample collection. This study was approved by the Harvard T.H. Chan School of Public Health and the Institutional Review Boards of the Department of Veterans Affairs.

DNAm measures

DNA samples were extracted using the IQAamp DNA Blood Kit (Qiagen, CA, U.S.) from the buffy coat of the whole blood collected between 1999 and 2013. We measured DNAm by the Illumina Infinium Human Methylation450K BeadChip (450 K; Illumina Inc., San Diego, CA, U.S.), which provides information on ~ 485,000 CpG sites. To minimize batch effects, we randomized the samples across 450 K BeadChip and 96-well plates based on a two-stage age-stratified algorithm so that age distributed similarly across plates (Dai et al. 2017).

We preprocessed DNAm data via the ewastools package in Github (Heiss and Just 2018; 2019). At the sample level, we finally included 1,559 high-quality samples (Heiss and Just 2018). We then corrected dye-bias using a regression on the logarithm of internal control probes described in Xu et al. (Xu et al. 2017). At the probe level, we first dropped the probes with a detection p-value > 0.01 based on an estimation of the background distribution using non-specific fluorescence (Heiss and Just 2019). We then excluded probes with outliers or low quality (Xu et al. 2016). We also dropped sex-chromosome probes, non-CpG probes, probes with a single nucleotide polymorphism (SNPs) at the CpG site (minor allele frequency (MAF) ≥ 5%), probes with SNPs at single-base extension (MAF ≥ 5%), and probes containing an SNP (MAF ≥ 5%) (Cardenas et al. 2019). In addition, we excluded probes with any annotated SNPs within ten base pairs of the CpG site as provided in the Illumina annotation regardless of MAF (The 1000 Genomes Projects Consortium 2012), cross-reactive probes (Chen et al. 2013), and probes that were non-unimodal in the dip test (Hartigan and Hartigan 1985). As a result, 360,272 high-quality probes were included in this study. The steps for probe cleaning are shown in more detail in Table S1.

We normalized DNAm data by controlling for the normalization factors in the outcome regression instead of using other commonly used approaches, such as beta-mixture quantile normalization (Teschendorff et al. 2013). This normalization approach ensures a better adjustment for batch effects as their impact often varies across probes, and we have applied it previously (Wang et al. 2020a; Wang et al. 2020b). Specifically, we used an elastic-net regularized generalized linear model. We chose the mixing parameter to be 0.5 and did the grid search over lambda because it was generally a good choice and led to good results. We then extracted the top five important experimental covariates (i.e., Non-polymorphic Red, Specificity I Red, Bisulfite Conversion I Red, Bisulfite Conversion II, Extension Red) which has the fewest number of coefficients equal to 0 (BeadArray controls report 2015).

DNAm level was expressed as the ratio of methylated cytosines over the sum of the methylated and unmethylated cytosines at the position and then multiplied by 100 (mean %5-methylcytosine, i.e., %5-mC). Thus, the DNAm level ranged from 0- to 100 %5-mC.

PM2.5 and its elements (PEs) measures

We collected daily ambient PM2.5 and 13 PEs at a stationary monitoring site located at the roof of the Harvard University Countway Library that is less than 1 km from the health examination clinical center. We measured PM2.5 total mass using a Harvard Impactor Sampler (Koutrakis et al. 1993). PEs included black carbon (BC), iron (Fe), lead (Pb), zinc (Zn), copper (Cu), aluminum (Al), calcium (Ca), silicon (Si), potassium (K), nickel (Ni), vanadium (V), sodium (Na), and sulfur (S). We focused on these PEs because their concentrations were mostly above detection limits and representative of different sources of emission (Dai et al. 2017; Wang et al. 2020b) – S is a regional pollutant primarily from power plants and regional pollution; K is from wood burning; BC, Cu, Pb, and Zn are components from motor vehicles (Pb is no longer in gasoline, but can be emitted from brake wear and ground down tire weights, and is present in re-suspended dust from motor vehicles); Ca and Fe are tracers of road dust; Al and Si are mostly from soil; Ni and V originate from oil combustion; Na is a major component of sea salt. We measured BC using an Aethalometer (Magee Scientific Inc., Berkeley CA, U.S.) (Lepeule et al. 2014) and the other 12 PEs using the Energy Dispersive X-ray Fluorescence Spectrometers (Epsilon 5, PAN-alytical, Almelo, Netherlands) (Wang et al. 2020b).

To estimate the short- and intermediate-term effects of PM2.5 and its 13 PEs on DNAm, we considered a single time window from day 0 (the current day of the event) to day 27 (27 days prior to the event) a priori based on previous studies (Gondalia et al. 2019; Mehta et al. 2015; Panni et al. 2016).

Covariates assessment

The covariates were selected a priori based on the relevant literature (Dai et al. 2016; Wang et al. 2020b): 1) sociodemographic factors including age (years), years of education; 2) lifestyle factors including smoking status (ever/never), cigarette pack-years, alcohol consumption (<2 or ≥ 2 drinks/day); body mass index (BMI, kg/m2); the metabolic equivalent of task (hours/week); 3) biological factors: the estimated cell-type compositions (CD4 + T lymphocytes, CD8 + T lymphocytes, natural killer cells, B cells, and monocytes) by the Houseman method (Houseman et al. 2012); 4) technical factors such as batch effects and five normalization factors; and 5) weather variables including season (March-May; June-August; September-November; December-February), and ambient temperature and relative humidity measured at the Boston Logan Airport weather station, which is 8 km away from the clinical center (Dai et al. 2016). For the covariates mentioned above, only years of education had four missing. The bootstrapping-based algorithm was used to impute the missing values (Honaker et al. 2011).

Statistical analyses

We used three different approaches of EWAS to analyze the accumulative effects of short- and intermediate-term exposure to PM2.5 and its 13 PEs on DNAm: site-by-site, regional, and pathway analyses.

Site-by-site analyses

To identify statistically differentially methylated probes (DMPs) due to PM2.5 and PEs, we performed site-by-site analyses. Since most of the DNAm levels were not normally distributed, statistical inferences based on classical least squares regression may be inappropriate. Therefore, we used median regression, which is influenced neither by outliers nor skewness in the distribution of the dependent variable (Koenker and Hallock 2001), to analyze the associations between particles and DNAm. In the median regression, we included in the models constrained polynomial distributed lag terms, which used a 3-degree polynomial structure to constrain how the size of the particle effects vary with the lags between exposure and DNAm measurement from the same day up to lags of 27 days (Schwartz 2000). Only the overall effect across the 28-d time window was reported. The distributed lag terms account for the latency of the effects of particles and the overall effects minimize the unstable estimate due to the high degree of collinearity between lags at neighboring days (Schwartz 2000).

We fitted median regression for longitudinal data by the Koenker et al. method (Koenker 2004) because approximately 80% of the participants had repeated DNAm measures. Briefly, this method enables one to fit fixed-effects and correlated random-effects within the same subject and uses Bootstrap inference. In addition, we adjusted for two selection biases: 1) healthier men were more likely to return for subsequent exams after 1999 and 2) mortality that occurred prior to the first DNAm visit (1999). We applied inverse probability weighting (Hogan and Lancaster 2004; Robins et al. 1995) via logistic regression to calculate the two probabilities given chronological age, education, BMI, blood pressure, smoking status, cigarette pack-years, alcohol consumption, C-reactive protein, asthma, chronic bronchitis, and emphysema at previous visit (Wang et al. 2020a, 2020b; McCracken et al. 2010). We then multiplied the two inverse probability weights. Therefore, the visits in this study remained representative of the original population.

We fitted the following models to assess the cumulative effects of ambient PM2.5 total mass only without considering its elements (as shown in formula (1)) and the cumulative effects of each of the 13 PEs one at a time, adjusting for PM2.5 total mass in the model (as shown in formula (2)). where Y is the DNAm level of subject at visit . β in formula (1), and φ in formula (2) are the lag-specific coefficients of PM2.5 total mass, with l ∈ [0, 27] lag days starting from the day of the visit to the previous 27 days. Σβl is the cumulative effect of PM2.5 total mass on DNAm over 28 days (lags 0–27). Similarly, φ in formula (2) are the lag-specific co-efficients of each of the PEs, with l ∈ [0, 27] lag days. Σφ is the cumulative effect of each of the PE’s exposure on DNAm over 28 days (lags 0–27). X are the covariates we mentioned above (Covariates Assessment) for subject at visit . δ in formula (1) and ω in formula (2) are the random intercepts for participant . ∊ in formula (1) and ε in formula (2) are the residuals. Median regressions with distributed lags, correlated random-effects, and inverse probability weighting were performed using the rqpd package in R (Bind et al. 2016).

We reported our results as the median difference in DNAm (%5-mC) per one interquartile range (IQR) increase in PM2.5 and its elements after 28 days (0–27) exposure. The results from individual CpGs were adjusted for multiple testing by Bonferroni threshold p-value < 4.96 × 10−9 (0.05/(360,272 × 28) for PM2.5 total mass and < 3.81 × 10−10 (0.05/(360,272 × 28 × 13) for each of the 13 PEs.

Regional and pathway analyses

In addition to examining the associations of PM2.5/PEs and DNAm at the individual probe level, we also investigated statistically differentially methylated regions (DMRs) and pathways in relation to PM2.5/PEs. To identify the DMRs, we applied the combp function from the ENmix package in R Bioconductor (Xu et al. 2020b) because the comb-p tool has the best sensitivity and highest control of false-positive rate compared to other DMR tools (Mallik et al. 2019). The input contained individual CpG sites with their p-values from site-by-site analyses as well as information on chromosomal locations (Xu et al. 2020b). We defined a significant DMR as one with three or more probes and its Sidak p-value < 0.05 (Li et al. 2020).

Ingenuity Pathway Analysis (IPA) database (QIAGEN Inc.) was used to identify canonical pathways significantly enriched with genes located around the corresponding top100 CpGs associated with PM2.5 total mass and its elements. We calculated permutation p-values using the results of 10,000 random shuffles of association p-values for the CpGs on the 450 K array (Xu et al. 2020a). We defined significant pathways if p-value < 0.05 and gene set contains ≥ 3genes.

Sensitivity analyses

We conducted two sensitivity analyses in the site-by-site analyses to check the robustness of our results. First, instead of adjusting the top five important experimental covariates in the median regression models, we additionally adjusted the following two experimental covariates (i.e., Bisulfite Conversion I Green and Hybridization High Medium). Second, we only considered the selection bias due to healthier men being more likely to return for subsequent exams after 1999 when DNA samples began being collected. We compared the effect sizes and p-values of the top one probe for PM2.5 total mass and each PE from the main analyses with the ones from the two sensitivity analyses.

Results

Population description

We included 695 men with 1,266 visits. The characteristics of the study population are presented in Table 1. Eighty-two percent of the participants had more than one visit. The mean age [standard deviation (SD)] at the first and the last visits were 73 (7) and 76 (5) years, respectively. The men were generally highly educated, with a mean of 15 years (SD = 3) of education. Thirty-one percent of subjects never smoked, and 80% had fewer than two drinks per day across all visits.

Table 1

Characteristics of elderly white men from the Normative Aging Study, 1999–2013.

Variables	First visit (N = 695)	Second visit (N = 432)	Third visit (N = 136)	Fourth visit (N = 3)	All visits (N = 1,266)
Age, Mean ± SD	73 ± 7	76 ± 7	80 ± 6	76 ± 5	75 ± 7
BMI (kg/m2), Mean ± SD	28.14 ± 4.15	27.90 ± 4.14	27.81 ± 4.33	30.02 ± 7.45	28.03 ± 4.17
Years of education, Mean ± SD	15.03 ± 3.00	15.27 ± 3.09	15.64 ± 3.07	16.67 ± 5.51	15.18 ± 3.05
Smoking status, n (%)
Never	220 (31.65%)	133 (30.79%)	48 (35.29%)	1 (33.33%)	402 (31.75%)
Ever	475 (68.35%)	299 (69.21%)	88 (64.71%)	2 (66.67%)	864 (68.25%)
Pack-years, Mean ± SD	20.71 ± 24.62	20.64 ± 24.17	17.65 ± 22.11	5.00 ± 8.66	20.32 ± 24.19
Alcohol consumption (drinks/day), n (%)
<2	564 (81.15%)	354 (81.94%)	110 (80.88%)	2 (66.67%)	1,030 (81.30%)
≥ 2	131 (18.85%)	78 (18.06%)	26 (19.12%)	1 (33.33%)	236 (18.70%)
Season, n (%)
Spring (March-May)	169 (24.31%)	108 (25.00%)	34 (25.00%)	0 (0%)	311 (23.85%)
Summer (June-August)	174 (25.04%)	115 (26.62%)	31 (22.79%)	0 (0%)	320 (25.28%)
Fall (September-November)	238 (34.15%)	137 (31.71%)	41 (30.15%)	1 (33.33%)	417 (32.88%)
Winter (December-February)	114 (16.50%)	72 (16.67%)	30 (22.06%)	2 (66.67%)	218 (17.27%)
Estimated cell type (%), Mean ± SD
Granulocytes	57.20 ± 9.46	58.58 ± 9.63	62.78 ± 9.63	59.27 ± 4.39	58.27 ± 9.46
Monocytes	10.79 ± 2.77	10.19 ± 3.00	10.13 ± 2.65	9.88 ± 2.87	10.52 ± 2.85
B cells	1.56 ± 2.61	1.17 ± 2.02	1.31 ± 3.50	0.83 ± 0.73	1.40 ± 2.54
CD4 + T lymphocytes	8.53 ± 4.37	8.44 ± 4.63	7.10 ± 3.99	9.17 ± 5.56	8.35 ± 4.44
CD8 + T lymphocytes	8.69 ± 3.08	8.62 ± 3.19	7.16 ± 2.79	7.51 ± 0.77	8.50 ± 3.12
Natural killer cells	7.21 ± 3.80	7.54 ± 4.40	7.06 ± 3.92	7.94 ± 3.27	7.31 ± 4.02

Abbreviations: BMI, body mass index; SD, standard deviation.

Concentrations of PM2.5 and 13 PEs concentrations

Table 2 describes of daily concentrations of PM2.5 and its elements during the study period (1999–2013). The mean (SD) concentration of daily PM2.5 was 10.3 (6.6) µg/m3, with an IQR of 7.3 µg/m3. Among the examined PEs, S (power plants and regional pollution) accounted for the largest proportion of PM2.5 total mass (10.4%). The major components of motor vehicles (BC, Pb, Zn, Cu) accounted for 8.0%, although other motor vehicle components not measured would increase the total from this source. The PEs from sea salt (Na) and soil (Al, Si) were responsible for 1.9% and 1.3%, respectively. Tracers of road dust (Ca, Fe), wood-burning (K), and oil combustion (Ni and V) took up 1.1%, 0.35%, and 0.1% of PM2.5 total mass, respectively.

Table 2

Distributions of the daily concentrations of PM2.5, its 13 PEs, and meteorological data during the study period from the Normative Aging Study, 1999–2013.

Particles	Mean (SD)	Min	Median (25%, 75%)	Max	IQR
PM2.5 (µg/m3)	10.32 (6.56)	1.07	8.45 (5.60, 12.91)	52.72	7.31
S (ng/m3)	1070.95 (883.31)	0.79	776.24 (485.44, 1398.32)	6663.98	912.88
BC (ng/m3)	801.68 (413.05)	187.38	707.91 (481.64, 1040.87)	2426.88	559.23
Pb (ng/m3)	5.81 (6.04)	0.01	5.15 (3.22, 7.53)	94.59	4.31
Zn (ng/m3)	12.69 (10.04)	0.01	10.07 (6.49, 15.60)	87.6	9.11
Cu (ng/m3)	3.55 (3.04)	0.01	3.33 (1.82, 4.86)	53.08	3.04
Na (ng/m3)	192.52 (139.77)	0.01	161.26 (97.20, 255.54)	1270.69	158.34
Al (ng/m3)	51.86 (39.91)	0.01	44.26 (30.73, 63.31)	601.71	32.58
Si (ng/m3)	78.39 (67.65)	1.10	64.59 (41.87, 94.48)	914.02	52.6
Ca (ng/m3)	34.34 (22.03)	0.37	30.76 (21.91, 39.69)	269.33	17.78
Fe (ng/m3)	74.97 (56.51)	5.36	65.90 (48.10, 88.74)	699.72	40.64
K (ng/m3)	36.02 (21.90)	0.42	30.95 (22.05, 44.88)	215.31	22.82
Ni (ng/m3)	3.03 (3.16)	0.01	2.06 (1.05, 4.07)	26.32	3.02
V (ng/m3)	3.50 (3.08)	0.01	2.61 (1.26, 4.73)	23.14	3.47
Temperature (℃)	12.33 (8.85)	−18.89	12.78 (5.60, 19.40)	31.10	13.80
Relative humidity (%)	68.13 (15.94)	19.92	67.83 (56.38, 81.17)	99.33	24.79

Abbreviations: PM2.5, particulate matter with an aerodynamic diameter ≤ 2.5 µm; S, sulfur; BC, black carbon; Pb, lead; Zn, zinc; Cu, copper; Na, sodium; Al, aluminum; Si, silicon; Ca, calcium; Fe, iron; K, potassium; Ni, nickel; V, vanadium; IQR, interquartile range; SD, standard deviation.

EWAS of PM2.5 and 13 PEs

Significantly differentially methylated probes

In the site-by-site analyses, we observed multiple significant DMPs for PM2.5 total mass and its 13 PEs across 28 days of exposure (see Table 3). For example, 3,173 significant DMPs were observed due to Pb exposure; 164 significant DMPs were identified by Fe (see Table 3). We presented the significant DMPs ranked by p-value with their annotated genes in Tables S2–15. We also compared the significant DMPs by PM2.5 total mass and the PEs. We found that among the 2,717 significant DMPs by PM2.5 total mass, there were 53 that were also associated with at least one of the 13 PEs (see Fig. 1). Among the 53 DMPs, 45 were associated with PEs from motor vehicles, and three (i.e., cg26201011, cg15415507, and cg07776285) were associated with PEs from more than one emission source.

Table 3

The number of differentially methylated probes (DMPs) of cumulative effects (lag 0–27d) from exposure to PM2.5 and its elements in the site-by-site analyses.

Particles		Number of DMPs
Sources	PM2.5 mass	2,717
Motor vehicles	BC	394
	Cu	615
	Pb	3,173
	Zn	281
Oil combustion	Ni	24
	V	20
Road dust	Ca	2
	Fe	164
Soil	Al	31
	Si	19
Wood burning	K	343
Power plants	S	270
Sea salt	Na	1,707

Abbreviations: PM2.5, particulate matter with an aerodynamic diameter ≤ 2.5 µm; BC, black carbon; Cu, copper; Pb, lead; Zn, zinc; Ni, nickel; V, vanadium; Ca, calcium; Fe, iron; Al, aluminum; Si, silicon; K, potassium; S, sulfur; Na, sodium; DMP, differentially methylated probes.

Fig. 1.

The heatmap of 53 differentially methylated probes that were associated with both PM2.5 mass and at least one of its 13 PEs. Abbreviations: PM2.5, particulate matter with an aerodynamic diameter ≤ 2.5 µm; BC, black carbon; Cu, copper; Pb, lead; Zn, zinc; Ni, nickel; V, vanadium; Ca, calcium; Fe, iron; Al, aluminum; Si, silicon; K, potassium; S, sulfur; Na, sodium.

We then extracted the top three DMPs associated with PM2.5 total mass and PEs (N = 42; i.e., 3 × 14). We found that among the 42 DMPs, 6 were located in the genes’ promoter, with all related with PEs, not PM2.5 total mass (see Table 4). For example, 28-day cumulative exposure to BC was associated with an increase in DNAm at cg24899205 (0.47 %5-mC, p = 0.00 × 10−12), which was annotated to the TXNRD2 gene’s promoter; whereas 28-day cumulative to S was associated with a decrease in DNAm at cg22056044 (−2.32 %5-mC, p = 0.00 × 10−12), which was mapped in the TM9SF3 gene’s promoter. For these 6 probes, we compared the effect sizes for their related PEs in Table 4 with the effect size for age and found that the absolute values of the effects sizes per IQR increase for all the PEs were greater than those for a one-year increase of age (Table S16).

Table 4

The top three DMPs associated with one IQR increase in PM2.5 and its PEs if they are located in the genes’ promoter in the site-by-site analyses.

Sources	Particles	CpG	UCSC gene [a]	Chr [a]	Position	p-value	Difference in DNAm [b]
Motor vehicles	BC	cg24899205	TXNRD2	22	19,929,184	0.00 × 10−16	0.47
	Cu	cg12298562	RBP4	10	95,361,132	0.00 × 10−16	−1.68
	Cu	cg24912904	ACTN3	11	66,314,217	0.00 × 10−16	−0.97
	Zn	cg01910869	EDNRB	13	78,493,182	0.00 × 10−16	0.75
Oil combustion	Ni	cg19821536	OR6B3	2	24,098,376	1.33 × 10−12	−1.71
Power plants	S	cg22056044	TM9SF3	10	98,347,055	0.00 × 10−16	−2.32

Note: BC, black carbon; Cu, copper; Zn, zinc; Ni, nickel; S, sulfur.

The University of California Santa Cruz UCSC database (GRCh37/hg19)

Difference in DNAm level (% 5-mC) per one IQR increase in particles

Significant regions and pathways

We identified multiple DMRs for PM2.5 mass and its PEs across 28 days of exposure (see Table S17–30). For example, 10,470 DMRs were observed due to PM2.5 total mass; 11,293 DMRs were found by Na. We compared the top 10 DMRs and their annotated genes from PM2.5 total mass and PEs and found no common DMRs or annotated genes.

In the pathway analyses, we also identified multiple significant pathways in relation to PM2.5 total mass and PEs: 17 for PM2.5, 45 for BC, 1 for Cu, 1 for Pb, 3 for Zn, 0 for Ni, 1 for V, 6 for Ca, 5 for Fe, 1 for Al, 1 for S, 4 for Si, 1 for K, and 0 for Na. We presented the significant pathways ranked by p-value in Tables S31–42. We also compared the enriched pathways by PM2.5 total mass and its PEs and found there were 12 common pathways identified by PM2.5 and some of its elements, for example, “Gg alpha signaling” and “RhoGDI signaling” (see Table 5).

Table 5

Common enriched pathways by PM2.5 and its elements.

Common pathways	Particles
Calcium Signaling	PM2.5	Zn
Dopamine-DARPP32 Feedback in cAMP Signaling	PM2.5	Ca
Glucocorticoid Receptor Signaling	PM2.5	Al
Gq alpha Signaling	PM2.5	BC
Aryl Hydrocarbon Receptor Signaling	BC	Al
Thrombin Signaling	BC	Al
Xenobiotic Metabolism Signaling	BC	Ca
Production of Nitric Oxide and Reactive Oxygen Species in Macrophages	BC	Ca
Human Embryonic Stem Cell Pluripotency	BC	Si
P2Y Purigenic Receptor Signaling Pathway	BC	Si
Wnt/β-catenin Signaling	Zn	Fe
RhoGDI Signaling	Al	Si

Note: PM2.5, particulate matter with an aerodynamic diameter ≤ 2.5 µm; Al, aluminum; Ca, calcium; Si, silicon; BC, black carbon; Fe, iron; Zn, zinc.

Sensitivity analyses

We extracted the top one probes from PM2.5 and its 13 PEs in the main analyses and compared their effect sizes and p-values with that from two sensitivity analyses. It showed no substantive changes (see Table S43). The effect sizes in the main analyses were close to the ones in sensitivity analyses #1. And the effect sizes were almost the same in the main analyses as in the sensitivity analyses #2.

Discussion

To our knowledge, this is the first EWAS of short- and intermediate-term exposure to PM2.5 and PEs. We identified multiple DMPs, DMRs, and pathways cumulatively associated with these particles across lag 0–27 days. Moreover, the identified DMPs, DMRs, and pathways by PM2.5 were different from those by its elements. Specifically, the significant pathways suggest that PM2.5 total mass was related to DNAm involved in mood disorders, neuroplasticity, immunity, and inflammation, whereas the pathways associated with motor vehicle emissions (BC, Cu, Pb, and Zn) were mostly annotated to cardiovascular disease and cancer pathways.

We start our discussion from the significant pathways associated with PM2.5 total mass. Pathways of “dopamine-DARPP32 feedback in cAMP signaling”, “synaptic long-term depression”, and “synaptic long-term potentiation” were involved in mood disorders (Amare et al. 2017; Bliss and Cooke 2011) and neuroplasticity (Morimoto et al. 2014). Other pathways such as “CD28 signaling in T helper cells”, “role of NFAT in regulation of the immune response”, and “Phospholipase C signaling” were associated with the immunity and inflammatory system (Dong and Flavell 2000; Zhu et al. 2018). Pathways associated with PEs were different from those associated with PM2.5 total mass. More specifically, elements from motor vehicle emissions (BC, Cu, Pb, Zn) were associated with pathways related to cardiovascular disease, such as “Cardiac β-adrenergic signaling”, “relaxin signaling”, “PPARs signaling”, “endothelin-1 signaling”, “leptin signaling in obesity”, “blood coagulation”, “eNOS signaling” (Johnson and Faraci 2011; Orekhov et al. 2016; Sandoval et al. 2014; St-Louis and Massicotte 1985; Yang and Barouch 2007), and cancer, for instance, “small cell lung cancer signaling”, “ceramide signaling”, “p53 signaling”, “aryl hydrocarbon receptor signaling”, “glioma signaling”, and “P2Y purinergic receptor signaling pathway”, and “Wnt/β-catenin signaling” (Feng et al. 2013; MacDonald et al. 2009; Schulien et al. 2020; Sheridan and Ogretmen 2021). The pathways associated with road dust (Ca and Fe) and soil (Al and Si) were involved in cancer, for example, “HIPPO signaling” (Zygulska et al. 2017) and “RhoGDI signaling” (Harding and Theodorescu 2010). These findings help us understand that specific PEs may affect health via different pathways compared with the PM2.5 total mass and suggest the necessity to examine the associations between DNAm and different PM2.5 elements. We note that our PEs do not include organic carbon, nitrate, or ammonium particles, which may explain why we see PM2.5 total mass associated with pathways (and DMR and DMP) not associated with the particle elements we had.

Many toxicological studies have linked particles and the pathways that we identified in this study. For example, Jiang et al. have examined exposure to PM2.5 in Nrf2−/− mice and found a reduction in PPARα (Jiang et al. 2020) in the liver. Zheng et al. have suggested that exposure to PM2.5 reduces the expression of PPARγ in both liver tissues and hepatic stellate cells (Zheng et al. 2013; Zheng et al. 2015). PPARs are ligand-inducible transcription factors and consist of three nuclear receptor isoforms, i.e., PPARγ, PPARα, and PPARδ. All three isoforms have a wide spectrum of beneficial effects on the prevention and treatment of cardiovascular disease (Orekhov et al. 2016). Pan et al. have assessed Pb-induced toxicity mechanisms via the skin permeability in mice and shown the overexpression of Rho GDP-dissociation inhibitor 2 (RhoGDI2) (Pan et al. 2010). Although we did not find the association between Pb and RhoGDI2 in the present study, we identified “RhoGDI signaling” by Al and Si. This pathway has been reported to mediate several processes during tumorigenesis and cancer progression, such as pancreatic and bladder cancer (Harding and Theodorescu 2010). Therefore, the intervention of this pathway has been proposed as a novel cancer therapeutic strategy (Harding and Theodorescu 2010).

It is worth mentioning that inflammation is involved in most of the pathways by PM2.5. We also find the pathways of inflammation-related diseases such as arthritis (due to BC exposure) in this study. These associations have been reported in both mechanistic and epidemiological studies. For example, in an experimental study, Liu et al. have found that PM2.5 might increase the risk of osteoarthritis via the production of IL-6 synthesis (Liu et al. 2021). Adami et al. have shown a positive association between air pollution and rheumatoid arthritis severity in a case-crossover study (Park et al. 2021). We note that these pathways are not necessarily associated with poor health endpoints. But our findings do suggest that particles may influence health by triggering different biological changes, and elements from motor vehicles are related with DNAm in more pathways compared with elements from other emission sources.

To date, only a limited number of research studies investigated associations between PM2.5 elements and DNAm and most of them focused on DNAm at specific loci or repeated elements. Our previous study was the only one examining PEs and DNAm in an epigenome-wide but focusing on a one-year moving average time window (Dai et al. 2017). The number of the identified DMPs and pathways was much less in our previous work compared with this present study (We did not perform the regional analyses in our previous study). Our previous work identified a total of 29 DMPs (20 for Fe, 8 for Ni, and 1 for V) and 9 significant pathways for Fe, Ni, and V. We did not find any common DMPs, but a few overlapping pathways in our two studies, such as pathways in cancer by Fe. The different findings (especially the significant DMPs) between the two studies indicate that the exposure windows are crucial when we examine the particles ~ DNAm relationships.

This study has some limitations: 1) The concentrations of ambient PM2.5 total mass and PEs were estimated from monitors at a stationary site because there were not enough monitors in the study area and AOD from MISR could only help with S and organic carbon prediction. The misclassification from a single ambient site occurs because of geographic differences in ambient exposure. In addition, the subjects with more outdoor activity would have higher exposure to ambient particles. All these may lead to the inevitable exposure misclassification. Our group is working on models predicting short-term metals’ exposures, which will be used in future. 2) We do not have data on gene expression, and thus we are not able to determine the regulation directions between DNAm and the coded protein. 3) This study only included elderly white men. The generalizability of our findings to other age groups, sex, and races are uncertain although studies that assess the modification effects of age, sex, and races in the association between particles and DNAm are limited (Ding et al. 2017; Ladd-Acosta et al. 2019; Wu et al. 2021). 4) IPA database is built based on multiple tissues; the role of the identified pathways in leukocytes may not be relevant for all diseases. 5) As in other cohort studies, this study is subject to selection bias such as index event bias and depletion of susceptible participants (Dudbridge et al. 2019; Fireman et al. 2020; Stovitz et al. 2018). However, we addressed this issue by using the inverse probability weighting to account for two main selection biases, i.e., we controlled for key predictors of mortality and drop-out. 6) It is quite challenging to determine potential confounders since this is an agnostic study (without pre-specifying any outcomes within the DNAm sites covered by the analytical platform we used). In our follow-up studies with specific health outcomes or pathways, we will be able to deal with these confounders in our models.

The notable strengths of this study include: 1) It was a relatively large EWAS to assess short- and intermediate-term effects of PM2.5 total mass and its elements. 2) Statistical inference using median regression does not require the normally distributed assumption for the residuals. 3) Repeated measurements of methylation give us a variation of the exposure and outcome within-subject. Hence, it increases statistical power. Also, since the visits are approximately 4 years apart, there is a large amount of variation in exposure between-subjects. This reduces the dependence of the analysis on between-subject variation in exposure.

Conclusions

In summary, this EWAS of short- and intermediate-term exposure to PM2.5 and its 13 PEs indicate that PM2.5 mass alone is not a sufficient metric when understanding the health effects of particles. PM2.5 total mass was associated with pathways in mood disorders, neuroplasticity, immunity, and inflammation whereas elements from motor vehicles (BC, Cu, Pb, and Zn) were mostly involved in cardiovascular disease and cancer. More EWAS of PEs with diverse study populations from different areas, especially from low-income, high-pollution regions are needed, to enrich present findings.