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.

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 ozoneinduced 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 . 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 (FEV 1 ) 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).

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, NO 2 , SO 2 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). sure, 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 Y i obs (w i = 0)-ignoring the crossover time period, and the seventeen ones observed after exposure to ozone, denoted by Y i obs (w i = 1), are shown in Fig. 1. We denote the mean participant methylation at site k after clean air and ozone exposures by m 0,k = ( 1 17 ) Σ i=1:17 Y i,k obs (w i = 0), and m 1,k = ( 1 17 ) Σ i=1:17 Y i,k obs (w i = 1), respectively. The full set of K = 484,531 CpG sites includes locations near promoter regions. Although the empirical distribution of m 0,k was bimodal among the overall set of 484,531 CpGs, the empirical distribution of m 0,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{Y i,k (w i = 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.  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(m 0,k , m 1,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(m 1,k ) k=1:484,531 versus f(m 0,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, Y i obs (w i =2) vs. Y i obs (w i =0), with the largest deviation for Participant 4 (Fig. 2).

Methods
The overall statistical strategy we follow comprises five stages: 1. "Global analyses" assess the overall effect of ozone exposure on the DNA methylome, 2. "Regional analyses" assess whether ozone exposure changes DNA methylation of genomic regions, 3. "Local analyses" estimate the ozone effect on DNA methylation measured at each CpG site, 4. "Responsiveness analyses" assess whether the seventeen participants exhibit differential epigenomic responses to ozone exposure, 5. "Enrichment analyses" identify for each participant whether CpG sites are over-represented among the set of responsive CpG sites.

Global tests for the O 3 effect on the epigenome. Distributional tests. The mean participant
methylations at site k after clean air and ozone exposures can be re-expressed as m 0, where N CA-O3 and N O3-CA are the number of participants who were first exposed to clean air (CA) or ozone (O 3 ), 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 2 17 = 131,072 possible allocations. Using the Fisherian procedure to test the sharp null hypothesis (H 00 ) stating for each participant i, for each CpG site k, Y i,j=1,k (w i,j=1 = 0) = Y i,j=1,k (w i,j=1 = 1) and Y i,j=2,k (w i,j=1 = 1, w i,j=2 = 0) = Y i,j=2,k (w i,j=1 = 0, w i,j=2 = 1), as described by Bind and Rubin 21 Table 3. Science table of the crossover randomized experiment.

Participant i
Visit j

Epigenomic potential outcomes
Unit-level causal effect Y i,j,k (w i,j = 0) Y i,j,k (w i,j = 1)  Causal estimands and estimators. We define three causal estimands, none of which, as with all causal estimands, is directly observable: 1. the i-j unit causal effect (UCE) for site k, τ i,j,k , as the unobservable contrast τ i,j,k = Y i,j,k (w i,j = 1) − Y i,j,k (w i,j = 0) (see Table 3), 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 = ( 1 2 ) (τ i,1,k + τ i,2,k ), and 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 = ( 1 17 ) Σ 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, d i,k , is: if participant i is exposed to clean air first (i: 1, …, N CA-O3 ), and Y i,j=1,k (w i,j=2 = 1) − Y i,j=2,k (w i,j=2 = 0) if participant i is exposed to ozone first (i: 1, …, N O3-CA ) (see Table 4). We estimated the K = 484,531 average participant causal effects (APCEs) using d k , the average of the seventeen d i,k 's.
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. f(τ low ) is the density of the k low (= 197,754)-dimension vector of APCEs such that m 0,k ≤ 0.2%5mC, 2. f(τ med ) is the density of the k med (= 112,190)-dimension vector of APCEs such that 0.2%5mC < m 0,k < 0.8%5mC, 3. f(τ high ) is the density of the k high (= 174,587)-dimension vector of APCEs such that 0.8%5mC ≤ m 0,k . Table 5. Global Fisher-exact tests. 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. www.nature.com/scientificreports/ 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 17 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 (m 0,k ) versus the participant methylation mean after ozone exposure (m 1,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., 1. quantile-quantile of observed − log 10 (p-value) versus "expected" − log 10 (p-value), 2. − log 10 (p-value) versus chromosome number, and, 3. − log 10 (p-value) versus APCEs, respectively.
We compared (2) and (3)    www.nature.com/scientificreports/ "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.   www.nature.com/scientificreports/

Results
Global tests of the O 3 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. The histogram of the estimated average participant causal effects (APCEs) for all loci, f(d k ) 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(d k,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(d k,med ) k=1:112,190 , appears to follow a bellshaped distribution (Fig. 4, middle panel), and (3) for loci with high methylation (> 0.8%5mC), f(d k,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 www.nature.com/scientificreports/ 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.
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).
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).  www.nature.com/scientificreports/ p-value: cg00605859, cg21964391, cg25265081, cg20129242, and cg11797346 (see Table 8 for a description of these CpG sites). 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 − log 10 (p value ) compared to the "expected" − log 10 (p value ) (Supplemental Fig. 7). The deviation from the "expected" − log 10 (p value ) could be explained by the CpG dependence, which is not taken into account for the construction of the "expected" − log 10 (p value ). 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/2 17 (or on the logarithmic scale, − log 10 (2/2 17 ) = 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).

Local analysis of the average participant causal effects and associated p-values. Supplemental
"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". 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 www.nature.com/scientificreports/ 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 NO 2 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.