Heterogeneous ozone effects on the DNA methylome of bronchial cells observed in a crossover study.

We used a randomized crossover experiment to estimate the effects of ozone (vs. clean air) exposure on genome-wide DNA methylation of target bronchial epithelial cells, using 17 volunteers, each randomly exposed on two separated occasions to clean air or 0.3-ppm ozone for two hours. Twenty-four hours after exposure, participants underwent bronchoscopy to collect epithelial cells whose DNA methylation was measured using the Illumina 450 K platform. We performed global and regional tests examining the ozone versus clean air effect on the DNA methylome and calculated Fisher-exact p-values for a series of univariate tests. We found little evidence of an overall effect of ozone on the DNA methylome but some suggestive changes in PLSCR1, HCAR1, and LINC00336 DNA methylation after ozone exposure relative to clean air. We observed some participant-to-participant heterogeneity in ozone responses.

RCT Entities:   Population Interventions Outcomes

Introduction

Ozone is a ubiquitous air pollutant that has been studied more extensively than perhaps any other environmental toxicant. Both short-term[1,2] and long-term[3,4] exposures have been linked with increased mortality and morbidity[5,6]. Numerous randomized controlled human exposure studies have shown that ozone exposure causes decrements in lung function and increased pulmonary inflammation (reviewed in US EPA air quality criteria for ozone and related photochemical oxidants) with considerable heterogeneity across people, and a randomized study has described cardiovascular changes in humans exposed to ozone[7]. Multiple studies have reported ozone-induced changes in pro-inflammatory cytokines in the lung of humans exposed to ozone[7,8], as well as changes in cytokine mRNA expression in bronchoalveolar lavage cells obtained following exposure of humans to ozone[9]. In addition, ozone has been shown to alter upstream mitogen-activated protein (MAP) kinase pathways that control cytokine gene expression[10,11].

It has been proposed that adverse health effects caused by exposure to air pollutants may be mediated by epigenetic changes[12]. Exposure to air pollution has been associated in several recent non-randomized studies with epigenetic changes, especially DNA methylation. In one study[13], exposure to traffic-related pollutants was associated with reduced lung function in elderly men and epigenetic changes in genes related to inflammation and immunity were proposed to modify the air pollution-lung function association. Increased air pollution was significantly associated with changes in DNA methylation of several genes, including those involved in inflammation[14]. A relationship between particulate air pollution and hypomethylation of the inducible nitric oxide synthesis gene was reported in a cohort of children[15].

Because all these studies are non-randomized cohort studies in which changes in ambient pollutant levels are associated with epigenetic changes, it is difficult to claim a causal relationship between air pollutant exposure and epigenetic changes. Furthermore, these studies measured epigenetic changes in blood cells rather than the primary targets of inhaled pollutants (i.e., respiratory tract cells). To show a direct causal link between air pollutant exposure and biologic changes in humans, randomized controlled human exposure studies are considered to be the “gold standard”. A recent randomized controlled human exposure study measured DNA methylation changes in T helper cells found in blood and reported changes in methylation in genes involved in mitochondrial oxidative energy metabolism in subjects exposed to particulate air pollution[16]. Here, we report the results of a random exposure study in which human volunteers were randomly exposed to either 0.3-ppm ozone or clean air on two occasions separated by several weeks. Bronchoscopy was performed 24 h after each exposure to remove cells directly targeted by ozone; i.e., airway epithelial cells. DNA methylation following each person’s clean air exposure was then compared with DNA methylation following that person’s ozone exposure. Because of the small sample size of the experiment, we eschew standard asymptotic inference assuming Student’s t-distributions under the null hypothesis, and instead capitalize on exact randomization-based inference.

In our study, we go beyond the current state of knowledge on the effect of ozone on DNA methylation. In observational studies, we do not know whether the estimated associations between ozone and the health outcomes are confounded by other environmental exposures correlated with ozone (e.g., NO2, temperature). In this randomized setting, we can directly estimate, with uncertainties expressed by the results of the randomization tests, the effect of ozone on the DNA methylome of bronchial cells (another “gold standard” for the field of air pollution epigenetics).

The crossover experiment

Study participants

The study population and exposure design have been described in detail previously by Devlin et al.[7], who estimated the causal effect of ozone versus clean air on the cardiovascular system. A total of seventeen healthy individuals (see their characteristics in Table 1) were recruited to participate in this study under a contract with Westat Corporation. Participants were excluded if they were smokers, pregnant, had any previous cardiopulmonary disease or allergies (as determined by their medical history and physical examination), or had a forced vital capacity (FVC) or forced expiratory volume in the first second of expiration (FEV1) of less than 80% predicted from their height and age. Prior to enrollment, all participants were informed of the study procedures and potential risks, and all provided a written informed consent. The consent forms and protocol were approved by the University of North Carolina School of Medicine and the US Environmental Protection Agency. The study was registered on ClinicalTrials.gov (NCT01492517).

Table 1

Characteristics of the seventeen participants.

	Min	25th	Mean	75th	Max
Age	21	23	25.4	27.5	33
Body mass index	17	23	25	27	29
Systolic blood pressure (mmHg)	106	113	121	130	136
Diastolic blood pressure (mmHg)	62	71	74	78	90
Heart rate (beats/min)	52	63	70	75	91
Gender
Male			88%
Female			12%
Race
White			88%
Non-white			12%

Study design

We conducted a randomized, single-blind, crossover study in which each participant was exposed twice, for 2 h to clean air (i.e., PM, CO, NO2, SO2 concentrations were below detection limit) or 0.3-ppm ozone (2015 US 8-h ozone standard: 0.07 ppm). The exposures were only for 2 h, in contrast to the standard, which is for 8 h. The actual number of ozone molecules to which a person was exposed in this study was about what they would have breathed for 8 h at the current standard. The exposure sessions were separated by a minimum of 13 days in an attempt to avoid carry-over effects associated with the first exposure. The exposure chamber and pollutant generation system are described elsewhere[7]. During the 2-h exposure, participants alternated between exercising for 15 min on a cycle ergometer and resting for 15 min. Minute ventilation was measured during each exercise session, and exercise levels were adjusted to obtain a minute ventilation of approximately 25 L min−1 m−2 body surface area. All exposures were conducted at the same time of the day and on the same day of the week for all participants between July 2010 and June 2011. The meteorological conditions in Chapel Hill during the time of the study[17] are presented in Table 2. As with any empirical studies, one can always object and find limitations in the study design: the volunteers were not housed and exposed to clean air or ozone for days, and some factors were not controlled (e.g., air pollution exposure preceding bronchoscopy).

Table 2

Meteorological conditions in Chapel Hill during the study time.

	Min	25th percentile	Mean	Median	75th percentile	Max
Ozone (ppm)	0.002	0.019	0.026	0.025	0.033	0.063
PM2.5 (µg/m3)	1.0	7.8	10.9	9.9	13.2	28.2
Temperature (°C)	− 7.6	7.3	15.1	15.6	23.3	32.8
Relative humidity (%)	26.5	57.9	67.9	70.0	78.4	96.3
Barometric pressure (hPa)	1,001.4	1,014.0	1,018.2	1,018.0	1,021.9	1,036.4

Measurement and pre-processing of DNA methylation

Twenty-four hours following each exposure, the subjects underwent a research bronchoscopy with brush biopsy to obtain bronchial epithelial cells as described previously[18]. The cytology brushes containing epithelial cells tips were placed in a 1.5 mL tube with 200 µL Trizol Buffer (ThermoFisher Scientific), DNA extracted using the Gentra Puregene Buccal Cell Kit (Qiagen, Inc.). DNA samples were stored frozen at − 80 °C until analysis. DNA extracted from the bronchial epithelial cells was sent to a commercial laboratory (Expression Analysis, Durham, NC) for DNA methylation assessment using the Illumina HumanMethylation 450 K BeadChip platform. The extracted DNA samples were placed on four chips (see positions in Supplemental Fig. 1). Participants labelled 1, 3, 4, 6, 7, and 13 had both extracted DNA, i.e., after exposures to clean air and ozone, positioned on the same Chip.

We performed background correction using noob[19], dye bias corrections, and corrected for probe design bias arising from Type I and Type II probes with the Beta-mixture quantile normalization method (BMIQ) method[20]. Although the modes of probe-type DNA methylation distributions were different, they were well aligned after the BMIQ procedure (see Supplemental Fig. 2). Methylation values are reported as the methylation beta value, i.e., the ratio of methylated probe signal intensities to methylated and unmethylated probe signal intensities. A small offset of 100 is added to the denominator to prevent inflated values when the sum of methylated and unmethylated probe signal intensities is low. After this pre-processing, the epigenomic analyses considered a set of 484,531 CpG sites, which we index by k.

Epigenomic data

The seventeen DNA methylome distributions observed after exposure to clean air, denoted by Yiobs(wi = 0)—ignoring the crossover time period, and the seventeen ones observed after exposure to ozone, denoted by Yiobs(wi = 1), are shown in Fig. 1. We denote the mean participant methylation at site k after clean air and ozone exposures by m0,k = () Σi=1:17 Yi,kobs(wi = 0), and m1,k = () Σi=1:17 Yi,kobs(wi = 1), respectively. The full set of K = 484,531 CpG sites includes locations near promoter regions. Although the empirical distribution of m0,k was bimodal among the overall set of 484,531 CpGs, the empirical distribution of m0,k was unimodal among 91,795 CpGs sites near promoter regions (see Supplemental Fig. 3), which supports the hypothesis that methylation sites near promoter regions tend to be hypomethylated compared to other CpG sites. Similar modes occur for the empirical participant median densities after exposure to clean air, median{Yi,k(wi = 0)}, among both sets of methylation sites (results not shown). The boundaries displayed in Supplemental Fig. 3 (i.e., at levels 0.2%5mC and 0.8%5mC) define three levels of methylation (i.e., low, medium, and high) throughout the paper. These values were chosen based on the bimodal empirical distributions of the participant methylation mean and median among the overall set of CpGs after exposure to clean air.

Figure 1

Participant distributions of DNA methylation after exposures to clean air (red) and ozone (blue).

Here, we perform statistical analyses on four sets of CpG sites: the epigenomic set of K = 484,531 CpG sites, and the sets of CpG sites with low, medium, and high methylation.

Throughout, we use f(D) to denote the empirical distribution of the generic K-component vector D. The joint bivariate epigenomic distribution of: (1) the participant mean methylation distribution under clean air, and (2) the mean methylation distribution under ozone exposure, f(m0,k, m1,k)k=1:484,531, also suggests two modes, around low and high methylation levels (data not shown). The Q–Q plot between the empirical epigenomic participant mean methylation distribution after ozone versus after clean air, f(m1,k)k=1:484,531 versus f(m0,k)k=1:484,531 suggests that these two distributions are similar (results not shown). However, the participant-specific Q-Q plots show some deviations between the empirical epigenomic participant methylation distributions after ozone versus after clean air, Yiobs(wi=2) vs. Yiobs(wi=0), with the largest deviation for Participant 4 (Fig. 2).

Figure 2

Participant-specific Q-Q plots for the empirical epigenomic participant methylation distributions after ozone versus after clean air.

Methods

The overall statistical strategy we follow comprises five stages:

“Global analyses” assess the overall effect of ozone exposure on the DNA methylome,

“Regional analyses” assess whether ozone exposure changes DNA methylation of genomic regions,

“Local analyses” estimate the ozone effect on DNA methylation measured at each CpG site,

“Responsiveness analyses” assess whether the seventeen participants exhibit differential epigenomic responses to ozone exposure,

“Enrichment analyses” identify for each participant whether CpG sites are over-represented among the set of responsive CpG sites.

Global tests for the O3 effect on the epigenome

Distributional tests

The mean participant methylations at site k after clean air and ozone exposures can be re-expressed as m0,k = (CA-O3) Σi,j:{wi,j=1}=0 Yi,j=1,k(wi,j=1 = 0) + (O3-CA) Σi,j:{wi,j=2}=0 Yi,j=2,k(wi,j=2 = 0), and m1,k = (CA-O3) Σi,j:{wi,j=2}=1 Yi,j=2,k(wi,j=2 = 1) + (O3-CA) Σi,j:{wi,j=1}=1 Yi,j=1,k(wi,j=1 = 1), where NCA-O3 and NO3-CA are the number of participants who were first exposed to clean air (CA) or ozone (O3), respectively. The stochasticity of the exposure assignment mechanism resides in the order each participant is exposed to clean air or ozone, which, if we assume a Bernoulli (with probability 1/2) assignment mechanism, implies the existence of 217 = 131,072 possible allocations. Using the Fisherian procedure to test the sharp null hypothesis (H00) stating for each participant i, for each CpG site k, Yi,j=1,k(wi,j=1 = 0) = Yi,j=1,k(wi,j=1 = 1) and Yi,j=2,k(wi,j=1 = 1, wi,j=2 = 0) = Yi,j=2,k(wi,j=1 = 0, wi,j=2 = 1), as described by Bind and Rubin[21], we performed four randomization tests that compared the empirical participant mean methylation distributions under ozone versus clean air exposures for the four sets of CpG sites: (1) f(m0,k)k=1:484,531 versus f(m1,k)k=1:484,531, (2) f(m0,k)k=1:197,754 versus f(m1,k)k=1:197,754, (3) f(m0,k)k=1:112,190 versus f(m1,k)k=1:112,190, and (4) f(m0,k)k=1:174,587 versus f(m1,k)k=1:174,587. We chose the Kolmogorov–Smirnov (KS) distance as the test statistic to compare the two mean distributions of clean air versus ozone. Then, we calculated the KS distance for all possible allocations of the N × 2 matrix of exposure assignment w = (wj=1,wj=2) assuming the sharp null hypothesis (H00). Here, note that because the test statistic is symmetric, the number of unique values of KS distances under H00 is 217/2 = 65,536 (assuming no ties).

Causal estimands and estimators

We define three causal estimands, none of which, as with all causal estimands, is directly observable:

the i–j unit causal effect (UCE) for site k, τi,j,k, as the unobservable contrast τi,j,k = Yi,j,k(wi,j = 1) − Yi,j,k(wi,j = 0) (see Table 3),

Table 3

Science table of the crossover randomized experiment.

Participant i	Visit j	Epigenomic potential outcomes		Unit-level causal effect
		Yi,j,k(wi,j = 0)	Yi,j,k(wi,j = 1)
1	1	Y1,1,k(w1,1 = 0)	Y1,1,k(w1,1 = 1)	τ1,1,k = Y1,1,k(w1,1 = 1) − Y1,1,k(w1,1 = 0)
…	…	…	…	…
17	1	Y17,1,k(w17,1 = 0)	Y17,1,k(w17,1 = 1)	τ17,1,k = Y17,1,k(w17,1 = 1) − Y17,1,k(w17,1 = 0)
1	2	Y1,2,k(w1,2 = 0)	Y1,2,k(w1,2 = 1)	τ1,2,k = Y1,2,k(w1,2 = 1) − Y1,2,k(w1,2 = 0)
…	…	…	…	…
17	2	Y17,2,k(w17,2 = 0)	Y17,2,k(w17,2 = 1)	τ17,2,k = Y17,2,k(w17,2 = 1) − Y17,2,k(w17,2 = 0)

the participant i (also unobservable) causal effect (PCE) for site k, τi,k, as the average of the two unit-level causal effects, τi,1,k, and τi,2,k, i.e., τi,k = () (τi,1,k + τi,2,k), and

the finite population, or average participant causal effect (APCE) for site k, τk, as the average of the N individual-level causal effects τi,k’s, τk = () Σi=1:17 τi,k; this is estimable under fewer assumptions than either UCE or PCE is.

The observed participant difference due to ozone versus clean air for site k, di,k, is: Yi,j=2,k(wi,j=2 = 1) − Yi,j=1,k(wi,j=1 = 0) if participant i is exposed to clean air first (i: 1, …, NCA-O3), and Yi,j=1,k(wi,j=2 = 1) − Yi,j=2,k(wi,j=2 = 0) if participant i is exposed to ozone first (i: 1, …, NO3-CA) (see Table 4). We estimated the K = 484,531 average participant causal effects (APCEs) using dk, the average of the seventeen di,k’s.

Table 4

Observed table of the crossover randomized experiment.

Participant i	Visit j	wi,j	Epigenomic potential outcomes		Observed participant difference di,k
			Yi,j,k(wi,j = 0)	Yi,j,k(wi,j = 1)
1	1	w1,1 = 0	Y1,1,k(w1,1 = 0)	?	d1,k = Y1,2,k(w1,2 = 1) − Y1,1,k(w1,1 = 0)
1	2	w1,2 = 1	?	Y1,2,k(w1,2 = 1)
…	…	…	…	…	…
17	1	w17,1 = 1	?	Y17,1,k(w17,1 = 1)	d17,k = Y17,1,k(w17,1 = 1) − Y17,2,k(w17,2 = 0)
17	2	w17,2 = 0	Y17,2,k(w17,2 = 0)	?

Average participant causal effects (APCEs) along the epigenome

We consider the marginal average participant causal effect (APCE) along the entire epigenome, and three conditional APCEs defined by boundaries of mean participant methylation at site k after clean air exposure, m0,k: low methylation; m0,k ≤ 0.2%5mC; medium methylation, 0.2%5mC < m0,k < 0.8%5mC; and high methylation, 0.8%5mC ≤ m0,k ≤ 0.2%5mC).

In particular, f(τ) denotes the density of the K(= 484,531)-dimension vector of APCEs. Slightly abusing notation, we write τ = (τlow, τmed, τhigh) as a permutation of τ such that:

f(τlow) is the density of the klow(= 197,754)-dimension vector of APCEs such that m0,k ≤ 0.2%5mC,

f(τmed) is the density of the kmed(= 112,190)-dimension vector of APCEs such that 0.2%5mC < m0,k < 0.8%5mC,

f(τhigh) is the density of the khigh(= 174,587)-dimension vector of APCEs such that 0.8%5mC ≤ m0,k.

We examined the characteristics of the empirical distributions of the APCEs for all sites, that is, f(τ), and for the three subsets of CpG sites (low, medium, and high), that is, (1) f(τlow), (2) f(τmed), and (3) f(τhigh).

Regional epigenomic analyses

Comb-p

To identify potential differentially methylated regions (DMRs), we used Comb-p (Python), which groups neighboring Fisher-exact p-values that are small and weighted according to their observed autocorrelation[22]. We calculate the univariate Fisher-exact p-values based on the standard paired statistic and on the 2 possible randomizations:

We use these Fisher-exact p-values to identify DMRs associated with ozone exposure with a minimum distance of 500 base pairs and a required p-value of 10−3 to start a region. For each CpG site within regions suggested by the Comb-p analysis, we report the estimated APCE and the associated Fisher-exact p-value.

DMRcate

We also examined averaged regional DNA methylation changes after ozone exposure using the R Bioconductor package DMRcate[23]. Briefly, this algorithm for regional DNA methylation analyses fits linear regression models via limma for each CpG and subsequently applies a Gaussian kernel smoothing function to test-statistics grouping significant probes based on genomic proximity. We modeled the mean DNA methylation difference between the clean air and ozone exposure at each CpG with an intercept only. For the discovery of candidate regions, given the multiple comparison issue, we chose the 10−6 threshold for statistical significance for the asymptotic, meaningless, p-values provided by limma.

Local analysis of the average participant causal effects (APCEs) and associated meaningful p-values

For each CpG k, we compared the participant methylation mean after clean air exposure (m0,k) versus the participant methylation mean after ozone exposure (m1,k) using univariate Fisher-exact two-sided tests, similarly as in Bind and Rubin[21]. We then constructed Q–Q plots of observed versus “expected” p-values, Manhattan[24] and Volcano[25] plots, i.e.,

quantile–quantile of observed − log10(p-value) versus “expected” − log10(p-value),

− log10(p-value) versus chromosome number, and,

− log10(p-value) versus APCEs, respectively.

We compared (2) and (3) plots for p-values calculated with the Fisher-exact and approximating asymptotic paired Student’s tests.

“Epigenomic responsiveness” to ozone

We estimated whether each participant i “responded” to ozone (with respect to the epigenome) by calculating the Cook’s distance[26] of the estimated individual causal effect for the 484,531 CpG sites. Participant i was defined as an ozone “epigenomic responder” at CpG k if the Cook’s distance of its estimated individual causal effect was greater than four times the mean of the Cook’s distances of the estimated individual causal effects across the 17 participants. For each participant, we compare “responsiveness” counts, which are bounded by 484,531, the total number of CpG sites. As sensitivity analyses, we also increased the “epigenomic responsiveness” threshold from four to five and six times the mean of the Cook’s distances.

Enrichment analyses for “responsive” CpGs

Based on the Cook’s distances described in the previous section, we identified (1) the individual(s) with a large number of responsive CpGs and describe them as “most-extreme responder(s)”, and (2) the individual(s) with a low number of responsive CpGs, the “least responder(s)”. We compared them in enrichment analyses. For each responder, we took the list of CpGs that were responsive in that individual, then removed from that list all the CpGs that were also responsive in the non-responders and performed enrichment analysis on the final gene list. Thus, enrichment analyses were performed for each responder based on an inter-individual comparison of responsive CpGs for responder(s) versus non-responder(s). We examined enrichment for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using the missMethyl R package while adjusting for varying gene size to avoid bias from large genes being overrepresented among the CpGs. We used the corrected p-values using the false discovery rate (FDR) to identify possible enriched pathways at a significance level of 0.05.

All methods were performed in accordance with the relevant guidelines and regulations.

Results

Global tests of the O3 effect on the epigenome

According to the four randomization tests comparing the participant mean methylation distributions under clean air versus ozone, we found some evidence against one or two of the four sharp null hypotheses (0.07 ≤ p-values ≤ 0.21, Table 5). The randomization distributions assuming the sharp null hypothesis are presented in Supplemental Fig. 4, suggesting that the asymptotic KS p-values, all less than 10−10, are deceptive.

Table 5

Global Fisher-exact tests.

Null hypothesis: no effect of ozone versus clean air	Kolmogorov–Smirnov distance	Fisher-exact p-value
f(m0,k)1:484,531 versus f(m1,k)1:484,531	D = 0.0083	0.21
f(m0,k)k=1:197,754 versus f(m1,k)k=1:197,754	D = 0.0204	0.17
f(m0,k)k=1:112,190 versus f(m1,k)k=1:112,190	D = 0.0138	0.12
f(m0,k)k=1:174,587 versus f(m1,k)k=1:174,587	D = 0.0124	0.07

The asymptotic Kolmogorov-tests led to highly significant p-values (< 10−10) compared to the ones obtained using the actual randomization of the crossover experiment.

The histogram of the estimated average participant causal effects (APCEs) for all loci, f(dk)k=1:484,531, appears to follow roughly a double-exponential distribution (Fig. 3). Similarly, (1) for loci with low methylation (≤ 0.2%5mC), f(dk,low)k=1:197,754, appears to follow roughly a double-exponential distribution (Fig. 4, top panel), (2) for loci with medium methylation (> 0.2%5mC and ≤ 0.8%5mC), f(dk,med)k=1:112,190, appears to follow a bell-shaped distribution (Fig. 4, middle panel), and (3) for loci with high methylation (> 0.8%5mC), f(dk,high)k=1:174,587, appears to follow a roughly double-exponential distribution (Fig. 4, bottom panel). The characteristics of these distributions for different boundaries (e.g., at 0.1%5mC and 0.9%5mC or at 0.3%5mC and 0.7%5mC) do not alter these conclusions (results not shown). We chose the exponential distribution to model the empirical distributions across CpG sites of the: (1) positive and (2) negative estimated APCEs and estimated the rate using the “fitdist” R function. We constructed a Q–Q plot between random draws from an exponential distribution (with estimated rate equal to 190) and the positive APCEs across 484,531 CpG sites (Supplemental Fig. 5, left graph). We proceed similarly for the negative APCEs and constructed a Q–Q plot between random draws from an exponential distribution (with estimated rate equal to 186) and the absolute value of the negative APCEs across 484,531 CpG sites (Supplemental Fig. 5, right graph). We observed that the negative APCEs can be fit with exponential distributions, although it appears more complex to find the approximating distributions for the positive APCEs.

Figure 3

Empirical distribution of the estimated average participant causal effects (K = 484,531 APCEs).

Figure 4

Empirical distributions of the estimated average participant causal effects (APCEs) for low, medium, and high DNA methylation sites (i.e., low: [0.0%5mC; 0.2%5mC], medium: ]0.2%5mC–0.8%5mC[, and high: [0.8%5mC; 1.0%5mC]).

Regional tests of the ozone effect on the epigenome

Although we found limited evidence of a global effect of ozone on the epigenome, using regional analyses, we observed differences in DNA methylation of two CpGs in the Phospholipid Scramblase 1 (PLSCR1) gene, nine CpGs in the Hydroxycarboxylic Acid Receptor 1 (HCAR1) gene, and eight CpGs in the Long Intergenic Non-Protein Coding RNA 336 (LINC00336) gene after exposure to ozone compared to clean air (Table 6).

Table 6

Differentially methylated regions (DMRs) and associated CpGs found using Comb-p analysis.

Genomic location	Number of CpGs	Nearby Gene	CpG site	Estimated average participant causal effect (APCE)	Fisher exact p-valuea
Chromosome 3 (146262175–146,262,422)	2	PLSCR1 (software p-value = 0.0026)	cg10679755	0.0252 − 0.0212 = 0.0040	0.0187
			cg23486067	0.0258 − 0.0194 = 0.0064	0.0154
Chromosome 12 (123215011–123215554)	9	GPR81 alias HCAR1 (software p-value = 0.0019)	cg22534509	0.3834 − 0.3560 = 0.0274	0.0135
			cg05290737	0.3522 − 0.3328 = 0.0194	0.0153
			cg22972858	0.4241 − 0.4025 = 0.0216	0.0068
			cg20566840	0.4294 − 0.3998 = 0.0296	0.0428
			cg01306688	0.2559 − 0.2356 = 0.0203	0.1680
			cg19975916	0.2122 − 0.1987 = 0.0135	0.2715
			cg19328828	0.1503 − 0.1251 = 0.0252	0.1268
			cg23505823	0.3863 − 0.3537 = 0.0327	0.0499
			cg13702536	0.4898 − 0.4543 = 0.0355	0.0138
Chromosome 6 (33560954–33561450)	8	C6orf227 Alias LINC00336 (software p-value = 0.0084)	cg04329454	0.2349 − 0.2018 = 0.0331	0.0039
			cg01392313	0.1600 − 0.1315 = 0.0284	0.0015
			cg06289138	0.1711 − 0.1535 = 0.0176	0.0949
			cg08301503	0.1479 − 0.1276 = 0.0203	0.0015
			cg07873320	0.1926 − 0.1852 = 0.0073	0.1752
			cg05602975	0.1069 − 0.0978 = 0.0091	0.0877
			cg19869469	0.1632 − 0.1454 = 0.0178	0.0007
			cg00536532	0.6582 − 0.6683 =  − 0.0101	0.2623

aBased on 217 replications.

The genes detected by the Comb-p approach were also detected by the limma-based approach DMRcate: six CpGs in the PLSCR1 gene, four CpGs in the HCAR1 gene, and four CpGs in the LINC00336 gene (Table 7).

Table 7

Differentially methylated regions (DMRs) and associated CpGs found using DMRcate with a threshold p = 10−6.

Genomic location	Number of CpGs	Nearby Gene	CpG site	Estimated average participant causal effect (APCE)	Fisher exact paired p-valuea
chromosome3 (146261941–146262421)	6	PLSCR1	cg05452836	0.0284 − 0.0250 = 0.0034	0.0083
			cg10679755	0.0252 − 0.0212 = 0.0040	0.0187
			cg12662193	0.0346 − 0.0257 = 0.0089	0.0004
			cg14795253	0.0287 − 0.0226 = 0.0061	0.0012
			cg15437043	0.2160 − 0.2143 = 0.0017	0.8154
			cg23486067	0.0258 − 0.0194 = 0.0064	0.0154
Chromosome 12 (123215010–123215308)	7	HCAR1	cg00357958	0.4270 − 0.4037 = 0.0234	0.0034
			cg01306688	0.2559 − 0.2356 = 0.0203	0.1680
			cg05290737	0.3522 − 0.3328 = 0.0194	0.0153
			cg19975916	0.2122 − 0.1987 = 0.0135	0.2715
			cg20566840	0.3998 − 0.4294 = 0.0296	0.0428
			cg22534509	0.3834 − 0.3560 = 0.0274	0.0135
			cg22972858	0.4241 − 0.4025 = 0.0216	0.0068
Chromosome 6 (33560953–33561449)	9	LINC00336	cg00536532	0.6582 − 0.6683 =  − 0.0100	0.2623
			cg01392313	0.1600 − 0.1315 = 0.0284	0.0015
			cg04329454	0.2349 − 0.2018 = 0.0331	0.0039
			cg05602975	0.1069 − 0.0978 = 0.0091	0.0877
			cg06289138	0.1711 − 0.1535 = 0.0176	0.0949
			cg07873320	0.1926 − 0.1852 = 0.0073	0.1752
			cg08301503	0.1479 − 0.1276 = 0.0203	0.0015
			cg19735538	0.0763 − 0.0644 = 0.0118	0.1153
			cg19869469	0.1632 − 0.1454 = 0.0178	0.0007

aBased on 217 replications.

Local analysis of the average participant causal effects and associated p-values

Supplemental Fig. 6 displays the distribution of the univariate Fisher-exact p-values. Five CpG sites achieved the minimum p-value: cg00605859, cg21964391, cg25265081, cg20129242, and cg11797346 (see Table 8 for a description of these CpG sites).

Table 8

Characteristics of CpG sites that achieve the minimum Fisher-exact p-value.

CpG site	Chromosome	Location	UCSC gene	UCSC region	Relation to CpG island	CpG island
cg00605859	15	73089234			Island	73089132–73089545
cg21964391	1	218253368
cg20129242	6	75915754	COL12A1	TSS200	Island	75914705–75916387
cg11797346	8	55367611			Island	55366180–55367628
cg25265081	15	29346097	APBA2	Body

Even though the ratios of the median of the empirically observed distribution of the test statistic to the “expected” median are close to λ = 1, we observed some deviation of the observed − log10(pvalue) compared to the “expected” − log10(pvalue) (Supplemental Fig. 7). The deviation from the “expected” − log10(pvalue) could be explained by the CpG dependence, which is not taken into account for the construction of the “expected” − log10(pvalue). The Manhattan and Volcano plots are presented in Figs. 5 and 6, respectively. These plots were constructed using the Fisher-exact and approximating asymptotic p-values. We observed that the minimum possible p-value 2/217 (or on the logarithmic scale, − log10(2/217) = 4.8) was subceeded by some univariate asymptotic p-values. The Manhattan plot also indicates that the Bonferroni adjustment (with a threshold around 6 on the logarithmic scale) offers no power in this high-dimensional setting with small sample size. The asymptotic and Fisher-exact (two-sided) p-values do not agree in this crossover experiment (Fig. 7).

Figure 5

Manhattan plots of the univariate asymptotic p-values (top) and Fisher exact p-values (bottom).

Figure 6

Volcano plots of the univariate p-values versus the average causal effect. Top plot: asymptotic p-values, bottom top: Fisher-exact p-values, and red line: maximum of − log10(p-value) that can be reached in this crossover experiment.

Figure 7

Asymptotic versus Fisher-exact (two-sided) p-values.

“Responsiveness” to ozone with respect to the epigenome

We observed a difference in “epigenomic responsiveness” to ozone exposure: Participants 4, 5, 10, and 14 had larger “epigenomic responses” to ozone than other participants, whereas Participant 6 had lower “epigenomic responses” (Fig. 8). That is, Participants 4, 5, 10, and 14 were defined as “most-extreme responders” and Participant 6 as “least responder”.

Figure 8

Responsiveness across subjects based on detection with Cook’s distance (threshold = 4, 5, and 6, from left to right).

Enrichment analysis

When comparing the CpGs present among the most-extreme responders, there was no CpG that “responded” in all four of the responsive individuals. In contrast, at the pathway level, there was overlap between the pathways that were enriched among the responsive CpGs for Participants 4, 5, 10, and 14. Using Cook’s distance greater than four to define responsive CpGs, thirty pathways were enriched in the responsive CpGs for Individuals 4, 5, 10, and 14, but were not enriched in Individual 6, the least responder (see Supplemental Table). Among these thirty pathways, twenty-one were enriched in all four responders and only two were enriched in just one responder. Pathways enriched in all four responders were linked to fatty acid and lipid metabolism, amino acid metabolism, and hormone/signaling molecule metabolism/biosynthesis.

Discussion

We found some modest evidence for regional changes in the DNA methylome of bronchial cells after ozone exposure compared to clean air exposure. In contrast to observational studies examining the health[1,2] or epigenetic[13,14] effects of short-term exposure to ozone, these findings (1) are designed to be unconfounded, (2) are based on DNA methylation of respiratory tract cells, primary targets of inhaled ozone, and (3) do not assume approximating asymptotic distribution that can be inaccurate in small studies[21]. Airway epithelial cells present in where we obtained the biopsies consist of at least three cell types: ciliated cells, Clara cells (mucin secreting), and basal cells. In this study, we cannot determine the contribution of each cell type to the methylation results.

The thorough description of the characteristics of: (1) the DNA methylation distributions after exposures to clean air and ozone, (2) the empirical distributions of the average participant causal effects for the sites with low, medium, and high methylations, should motivate future statistical work. Traditional epigenome-wide analyses rely mostly on approximate asymptotic p-values, which is an issue in human experimental settings with small samples. Here, we perform statistical tests of the ozone effect on the DNA methylome and on site-specific tests using randomization-based inference to compute Fisher-exact p-values. Global tests did not reveal substantial changes in DNA methylome measured after ozone versus clean air, but there were some regional differences, as well as inter-individual differences. By comparing inter-individual differences, we were able to discover pathways that are potentially related to the “epigenomic responsiveness” to short-term ozone exposure.

To our knowledge, this is the first time that personalized epigenomic responsiveness is suggested. It is interesting that healthy participants are not responding homogeneously to short-term ozone exposure. New statistical methods should be developed to investigate responsiveness, in particular in high-dimensional settings. A study should further examine the estimated individual causal effects, especially whether the observed heterogeneity can be explained by background covariates such sex, age, race, DNA methylation under clean air. Future statistical research should also focus on the understanding of the high-dimensional characteristics of estimated average causal effects and should explore whether these distributions approximately follow double-exponential distributions, which has been used to model noise accumulation[27]. Other strategies that focus on computing accurate and robust p-value adjustments should also be explored[28].

Using the number of CpGs per individual with substantial responsiveness to ozone based on their Cook’s distance, we were able to identify individuals who had an apparent excess, or deficit, of responsive CpGs. Four individuals were highly responsive to ozone from the perspective of DNA methylation loci changes, whereas one individual was “under-responsive”. Although there was no overlap in the individual CpGs among the responders, there was substantial overlap in the pathways enriched for altered CpGs among the most-extreme responders. Among “responsive” CpGs identified in the responders, 70% were present in all “ozone responders” and linked to pathways related to lipid and fatty acid metabolism and amino acid metabolism, among others. Long-chain fatty acids have been previously shown to be responsive to short-term exposure to NO2[29]. In the same participants of our study, fatty acid metabolism and amino acid metabolism were found to be enriched after ozone exposure compared to before the exposure[30]. Particulate matter-associated amino acids were also strongly associated with markers of inflammation and lung function[31].

At the pathway level, there were three genes (i.e., HCAR1, PLSCR1, and LINC00336) that were uncovered in both the DMRcate and Comb-p analyses. HCAR1 is located in an ozone exposure DMR on chromosome 12. HCAR1 is an important lactic acid receptor molecule, which itself has a variety of known functions. Some functions are mediated by HCAR1, including inhibition of pro-inflammatory and cytotoxic responses. In the brain, HCAR1 mediates angiogenesis[32] and is related to wound healing[33]. These functions have not been examined in the lung but may be a key function linking HCAR1 and ozone given the tissue damage induced by ozone exposure. PLSCR1 is related to epidermal growth factor (EGF) and its receptor (EGFR) signaling pathway. This gene is closely linked with apoptosis and cellular growth. In a mouse model, cells with elevated PLSCR1 expression showed an eightfold decrease in growth[34]. Cellular growth is an important component of epithelial-layer revival after exposure to ozone. LINC00336 is a long non-coding RNA that is involved in ferroptosis (i.e., type of programmed cell death dependent on iron and characterized by the accumulation of lipid peroxides) in the lung[35].

Conclusion

Short-term exposure to ozone may alter the epigenome in young, healthy individuals who had epigenomes with variable responsiveness to an identical exposure to ozone. The most responsive individuals had several shared enriched pathways that mirrored effects of short-term exposure to ozone on fatty acid and amino acid metabolism as seen in metabolomics data from the same individuals and independent cohorts. Future studies are needed to further expand on the effects observed here and to integrate data to obtain a multiomic understanding of short-term ozone exposure effects.

Supplementary Information.