Introduction

Ozone is a ubiquitous air pollutant that has been studied more extensively than perhaps any other environmental toxicant. Both short-term1,2 and long-term3,4 exposures have been linked with increased mortality and morbidity5,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 ozone7. Multiple studies have reported ozone-induced changes in pro-inflammatory cytokines in the lung of humans exposed to ozone7,8, as well as changes in cytokine mRNA expression in bronchoalveolar lavage cells obtained following exposure of humans to ozone9. In addition, ozone has been shown to alter upstream mitogen-activated protein (MAP) kinase pathways that control cytokine gene expression10,11.

It has been proposed that adverse health effects caused by exposure to air pollutants may be mediated by epigenetic changes12. Exposure to air pollution has been associated in several recent non-randomized studies with epigenetic changes, especially DNA methylation. In one study13, 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 inflammation14. A relationship between particulate air pollution and hypomethylation of the inducible nitric oxide synthesis gene was reported in a cohort of children15.

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 pollution16. 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.
Full size table

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 elsewhere7. 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 study17 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.
Full size table

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 previously18. 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 noob19, 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) method20. 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 = (\(\frac {1} {17}\)) Σi=1:17 Yi,kobs(wi = 0), and m1,k = (\(\frac {1} {17}\)) Σ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
figure 1

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

Full size image

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
figure 2

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

Full size image

Methods

The overall statistical strategy we follow comprises five stages:

  1. 1.

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

  2. 2.

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

  3. 3.

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

  4. 4.

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

  5. 5.

    “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 = (\(\frac {1} {N}\)CA-O3) Σi,j:{wi,j=1}=0 Yi,j=1,k(wi,j=1 = 0) + (\(\frac {1} {N}\)O3-CA) Σi,j:{wi,j=2}=0 Yi,j=2,k(wi,j=2 = 0), and m1,k = (\(\frac {1} {N}\)CA-O3) Σi,j:{wi,j=2}=1 Yi,j=2,k(wi,j=2 = 1) + (\(\frac {1} {N}\)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 Rubin21, 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:

  1. 1.

    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),

  2. 2.

    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 = (\(\frac {1} {2}\)) (τi,1,k + τi,2,k), and

  3. 3.

    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 = (\(\frac {1} {17}\)) Σi=1:17 τi,k; this is estimable under fewer assumptions than either UCE or PCE is.

Table 3 Science table of the crossover randomized experiment.
Full size table

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.
Full size table

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:

  1. 1.

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

  2. 2.

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

  3. 3.

    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 autocorrelation22. We calculate the univariate Fisher-exact p-values based on the standard paired statistic and on the 2\(\frac {17} {2}\) possible randomizations:

$$\begin{aligned} &{\text{T}}_{{{\text{k}},{\text{ paired}}}} = \left| {\left( {\frac{1}{17}} \right) \, \Sigma_{{{\text{i}} = {1}:{17}}} {\text{d}}_{{{\text{i}},{\text{k}}}} } \right|/\left[ {{\text{ s}}_{{\text{D}}} /{ 17}^{{{1}/{2}}} } \right], \hfill \\ &{\text{where s}}_{{\text{D}}}^{{2}} = \, \left[ {{1}/{16}} \right] \, \left[ {\Sigma_{{{\text{i}} = {1}:{17}}} \left( {{\text{d}}_{{{\text{i}},{\text{k}}}} - \, \left( {{1}/{17}} \right) \, \Sigma_{{{\text{i}} = {1}:{17}}} {\text{d}}_{{{\text{i}},{\text{k}}}} } \right)^{{2}} } \right]. \hfill \\ \end{aligned}$$

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 DMRcate23. 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 Rubin21. We then constructed Q–Q plots of observed versus “expected” p-values, Manhattan24 and Volcano25 plots, i.e.,

  1. 1.

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

  2. 2.

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

  3. 3.

    − 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 distance26 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.
Full size table

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
figure 3

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

Full size image
Figure 4
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]).

Full size image

Regional tests of the ozone effect on the epigenome

Comb-p

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.
Full size table

DMRcate

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.
Full size table

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.
Full size table

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
figure 5

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

Full size image
Figure 6
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.

Full size image
Figure 7
figure 7

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

Full size image

“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
figure 8

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

Full size image

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 health1,2 or epigenetic13,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 studies21. 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 accumulation27. Other strategies that focus on computing accurate and robust p-value adjustments should also be explored28.

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 NO229. 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 exposure30. Particulate matter-associated amino acids were also strongly associated with markers of inflammation and lung function31.

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 angiogenesis32 and is related to wound healing33. 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 growth34. 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 lung35.

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.