Effect of concentration and duration of particulate matter exposure on the transcriptome and DNA methylome of bronchial epithelial cells

Abstract Exposure to particulate matter (PM) from ambient air pollution is a well-known risk factor for many lung diseases, but the mechanism(s) for this is not completely understood. Bronchial epithelial cells, which line the airway of the respiratory tract, undergo genome-wide level changes in gene expression and DNA methylation particularly when exposed to fine (<2.5 µm) PM (PM2.5). Although some of these changes have been reported in other studies, a comparison of how different concentrations and duration of exposure affect both the gene transcriptome and DNA methylome has not been done. Here, we exposed BEAS-2B, a bronchial epithelial cell line, to different concentrations of PM2.5, and compared how single or repeated doses of PM2.5 affect both the transcriptome and methylome of cells. Widespread changes in gene expression occurred after cells were exposed to a single treatment of high-concentration (30 µg/cm2) PM2.5 for 24 h. These genes were enriched in pathways regulating cytokine–cytokine interactions, Mitogen-Activated Protein Kinase (MAPK) signaling, PI3K-Akt signaling, IL6, and P53. DNA methylomic analysis showed that nearly half of the differentially expressed genes were found to also have DNA methylation changes, with just a slightly greater trend toward overall hypomethylation across the genome. Cells exposed to a lower concentration (1 µg/cm2) of PM2.5 demonstrated a comparable, but more attenuated change in gene expression compared to cells exposed to higher concentrations. There were also many genes affected by lower concentrations of PM2.5, but not higher concentrations. Additionally, repeated exposure to PM2.5 (1 µg/cm2) for seven days resulted in transcriptomic and DNA methylomic changes that were distinct from cells treated with PM2.5 for only one day. Compared to single exposure, repeated exposure to PM2.5 caused a more notable degree of hypomethylation across the genome, though certain genes and regions demonstrated increased DNA methylation. The overall increase in hypomethylation, especially with repeated exposure to PM2.5, was associated with an increase in expression of ten–eleven translocation enzymes. These data demonstrate how variations in concentration and duration of PM2.5 exposure induce distinct differences in the transcriptomic and DNA methylomic profile of bronchial epithelial cells, which may have important implications in the development of both acute and chronic lung disease.


Introduction
Particulate matter (PM) remains one of the most harmful forms of air pollution, contributing to more than 3 million deaths per year worldwide (1). Traffic, diesel, and manufacturing are major sources of PM in industrialized cities and countries. Small particles, including those less than 2.5 mm in diameter (PM 2.5 ), are especially harmful because of their ability to travel deep in the respiratory tract and cause toxicity to cells lining the airway. Over the decades, public health efforts have curtailed levels of PM 2.5 in many countries throughout the world. However, recent studies have shown that low levels of exposure, even below that set by governmental regulatory agencies in many advanced countries, continue to cause harm in the general population (2,3).
Bronchial epithelial cells line the airways and are the first cells that are exposed to PM 2.5 . Many studies have shown that PM 2.5 causes oxidative stress and toxicity in bronchial epithelial cells in a concentration-dependent fashion (4)(5)(6)(7)(8). Several studies have also shown that PM 2.5 causes widespread changes in the transcriptome of these cells (9)(10)(11)(12)(13)(14)(15)(16). These studies, however, mostly employ PM 2.5 at high concentrations and examine the effects of PM 2.5 only after a single exposure and short time interval. How the transcriptome changes with lower concentrations and chronic exposure to PM 2.5 has not been extensively studied. Recognizing that many humans are exposed to low levels of chronic air pollution, many investigators have attempted to model this phenomenon in vitro by treating epithelial cells with low levels of PM 2.5 repeatedly over several days. These studies have shown that epithelial cells under these conditions persistently release inflammatory cytokines, develop epigenetic changes, and undergo epithelial-mesenchymal transition (17)(18)(19). Our laboratory has shown that repeated exposure to low-concentrations of PM 2.5 over seven days resulted in the upregulation of many genes that were not observed after a single exposure to PM 2.5 (20). However, a full transcriptomic analysis of this treatment protocol-along with a comparison with higher concentrations and short-term exposure-has not been reported. We hypothesize that even though a single exposure to high concentrations of PM 2.5 may induce oxidative stress and DNA damage, along with reciprocal anti-oxidant and anti-toxin response, exposure to lower concentrations, especially over time, may cause more subtle changes in gene expression that might differ from exposure to higher concentrations, but play equally important roles in disease development.
DNA methylation is a well-characterized epigenetic mechanism recognized for its ability to influence gene expression, often in a persistent, heritable way. DNA methylation is critical to normal development, and alterations in DNA methylation have been shown to contribute to cancer (21), autoimmune disease (22), and asthma (23). DNA methylation, especially within gene promoters and in CpG islands, is traditionally associated with suppression of gene expression (24). More recent studies, however, have shown that increased DNA methylation, especially within gene bodies and exons, are associated with increased transcription (25), emphasizing the importance of distinguishing where in the gene DNA methylation changes are occurring and how it might affect gene transcription. PM 2.5 has been shown in numerous in vitro and in vivo studies to cause changes in DNA methylation, either hyper-or hypo-, as measured in whole blood or in isolated cell types (26)(27)(28), though the significance of methylation changes at individual CpG loci is not always clear. Here, we examine the transcriptomic and DNA methylomic changes that occur in bronchial epithelial cells in response to PM 2.5 at different concentrations and duration of exposure. We observed that PM 2.5 causes widespread changes in the expression and DNA methylation of genes that are important in epithelial cell biology. That we observed distinct patterns of DNA methylomic and transcriptomic changes with different treatment protocols further highlight the importance of considering variables such as concentration and duration of exposure when trying to model and understand the pleiotropic effects of PM 2.5 .

Methods
Collection of PM 2.5 and Extraction from Filters PM 2.5 was isolated from ambient air pollution collected in Beijing, China from 19-21 January 2015 as previously described (20). Low volume manual samplers placed on the rooftop of Peking University School of Public Health collected PM 2.5 over a 24-h period on 90 mm Emfab filters, which are made of borosilicate fibers reinforced with woven cloth and bonded with polytetrafluoroethylene (TX40HI20WW, part #7234, Pall Company, Beijing Office, Beijing, China). Each filter was folded, wrapped in aluminum foil, and stored in À20 C until extraction. Twentyfour hours before extraction, filters were placed in sterile amber jars and equilibrated at fixed humidity and room temperature in a sterile biosafety containment hood, and weighed before wetting and extraction. To extract PM 2.5 , 20 ml of double distilled water was added to each amber jar and filters were sonicated (VWR, model no. 97043-968, VWR International, Radnor, Pennsylvania, USA) on ice at 15 min intervals for a total of 3 h. After sonication, filters were air-dried for 3 days in amber jars located in the same biosafety containment hood at constant humidity and room temperature before being weighed on a microbalance. The difference in weight (averaged from 3-5 measurements) before and after extraction was used to calculate the concentration (mg/ml). The extracted PM 2.5 was aliquoted and stored for future use at À80 C.

Cells and Treatments
BEAS-2B, a bronchial epithelial cell line, was cultured on collagen-coated tissue-culture plates in serum-free Bronchial Epithelial Growth Medium (BEGM; CC-3170, Lonza, Walkersville, MD), which consists of basal medium supplemented with standardized growth factors provided by the manufacturer (BEGM BulletKit CC-3171 and CC-4175; Lonza). Cells were maintained in a 37 C incubator with 5% CO 2 . Collagen-coated tissue-culture plates were prepared by coating plates with pre-made bovine collagen solution (PureCol-Type I Bovine Collagen Solution, Advanced BioMatrix, San Diego, CA) diluted to a concentration of 3 mg/ml with 0.1 N HCl overnight at 4 C. The liquid was then aspirated and plates were UV-irradiated for 30 min before they were washed three times with sterile water.
For experiments, cells were plated at a density of 5 Â 10 5 in six-well plates in BEGM and allowed to adhere overnight. Medium was replaced the next day with either fresh BEGM alone or fresh BEGM containing a low (1 mg/cm 2 ) or high (30 mg/ cm 2 ) concentration of PM 2.5 for 24 h. For experiments involving chronic, repeated exposure, cells were initially plated at a density of 2.5 Â 10 5 in six-well plates. Medium was replaced the next day with fresh BEGM either alone or with 1 mg/cm 2 of PM 2.5 . This was repeated each day, with all cells being washed with phosphate-buffered saline (PBS) between treatments, for a total of six days before cell lysates were collected on day 7 for RNA and DNA analysis.

RNA-Seq
RNA was initially isolated from cells using Trizol (Catalog Number 15596018, Invitrogen, Carlsbad, CA, USA); extracted RNA underwent an additional step for cleanup using the RNeasy kit (Qiagen, Germantown, MD). RNA was assessed for quality on the Agilent TapeStation and samples with RNA integrity number (RIN) > 7 were used for subsequent library preparation. Poly-A enrichment was used to select mRNA for library preparation, and samples were sequenced as a single-end 50 bp fragment on the Illumina HiSeq 4000, which was performed by the University of Michigan Advanced Genomics Core. Each experiment, including the single exposure protocol and repeated exposure protocol, was performed three independent times; a total of 15 samples were multiplexed and sequenced over two lanes.

Enhanced Reduced Representation Bisulfite Sequencing (eRRBS)
Assessment of genome-wide DNA methylation was performed by eRRBS, as previously described (29). DNA was isolated from cells using DNeasy kit (Qiagen) and submitted to the University of Michigan Epigenomics Core for library construction and bisulfite sequencing. Briefly, DNA was fragmented using the restriction enzyme Msp I, end repaired, A-tailed, and ligated (all using reagents from New England Biolabs, Ipswich, MA, catalog #s: R0106M, M0203L, M0210L, M0212L) to pair-end methylated adapters from Illumina. DNA was then resolved on agarose gel and fragments (150-250 bp and 250-450 bp) were excised and eluted. Samples then underwent bisulfite conversion using the EZDNA Methylation kit (Zymo Research, Inc., Irvine, CA) and DNA amplified using the Roche High Fidelity FastStart system. Final libraries were cleaned using Agencourt AMPure beads before sequencing on the Illumina HiSeq 4000. Each experiment, including the single exposure protocol and repeated exposure protocol, was performed three independent times; a total of 15 samples were multiplexed and sequenced over five lanes. Raw data for all RNA-Seq and eRRBS samples were uploaded into the National Center for Biotechnology Information Gene Expression Omnibus (GEO) database under accession number GSE155617.

Quantitative RT-PCR
RNA was reverse-transcribed to cDNA using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Carlsbad, CA, USA) according to the manufacturer's recommendations; quantitative real-time PCR was performed on cDNA using SYBR green PCR Master Mix (Applied Biosystems) on a StepOne Realtime PCR System (Applied Biosystems). The fold-change in expression of target genes (purchased as pre-designed TaqMan assays with primer-probes from Applied Biosystems) was calculated by the DDCt method relative to b-actin as the endogenous control. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an alternative endogenous control to verify the findings.

Data Analysis
For RNA-Seq analysis, FASTQ raw reads were aligned using STAR (Spliced Transcripts Alignment to a Reference) version 2.5.3 to the reference genome GRCh38. Reads were counted using featureCounts and normalized using edgeR via the Bioconductor package. Each experiment was performed three independent times and the data from all three experiments was averaged together; differential expression relative to vehicle control was performed using edgeR and analyzed by glmQLFit with correction for multiple testing. A linear fold-change > 1.5 (log 0.6) with adjusted P-value (adjusted for false discovery rate) < 0.05 was deemed significant. For eRRBS analysis, sequencing reads were assessed for quality by FastQC, and reads were trimmed using TrimGalore to trim low-quality bases (quality score lower than 20), adapter sequences, and end-repair bases from the 3 0 end of reads. Bismark was used for methylation calling and alignment (30) to the reference genome using Bowtie2 (31). The methylSig R package was used to identify sites of differential methylation (32). An absolute difference of > 10% methylation with adjusted P-value (adjusted for false discovery rate) < 0.05 was deemed significant.
Pathway analysis was performed using Advaita iPathway (Advaita Corportation, Ann Arbor, MI) analysis. Statistically significant pathways were defined by adjusted P-value < 0.05.

Results
Transcriptomic Changes in BEAS-2B Cells after Single Exposure of PM 2.5 Other studies have reported the results of compositional analysis of PM 2.5 from Beijing, China taken on the same or similar dates (January through February, 2015) as the PM 2.5 that we used in our study (33,34). In those analyses, PM 2.5 was composed of anions SO 2À 4 and NO À 3 , elemental metals such as Ag, As, Cd, Cu, Hg, Pb, Se, and Zn, and polyaromatic hydrocarbons including Benz(a)anthracene, Chrysene, and 1,8-Naphthalic anhydride (33,34). The largest source of this PM 2.5 was identified as coming from vehicle emissions. The level of endotoxin in our highest treatment condition (30 mg/cm 2 ) was 0.433 EU/ml, as previously described (20).
Although treatment with high concentrations of PM 2.5 (! 50 mg/cm 2 ) induced apoptosis and cytotoxicity, as shown in our previous study, measures of lactate dehydrogenase release and cleaved poly-ADP ribosylation remained low at concentrations up to 30 mg/cm 2 in our experimental system (20). Levels of IL6, TNF, and TSLP robustly increased at this concentration, so treatment with 30 mg/cm 2 was chosen as the maximum concentration for our transcriptomic and methylomic studies.
We first treated BEAS-2B cells with either vehicle control or 30 mg/cm 2 of PM 2.5 for 24 h and performed RNA-Seq analysis to examine the transcriptomic changes that occur with PM 2.5 treatment. Compared to cells treated with vehicle control, cells treated with 30 mg/cm 2 of PM 2.5 demonstrated a statistically significant increase in 612 genes and decrease in 698 genes (Fig. 1). A list of all of the differentially expressed genes including foldchange and statistical significance is included in Supplementary Table S1. Examination of these differentially expressed genes included upregulation of CYP1A1, CYP1B1, and AHRR, which are commonly increased when cells have been exposed to polyaromatic hydrocarbons (7,35), such as that present in urban PM 2.5 . Genes such as HMOX1, NQO1, GCLM, SQSTM1, and FTL, which are often increased in response to oxidative stress (36,37), were also increased. Expression of MDM2 and CDKN1A were also elevated after treatment with PM 2.5 , which was also expected as a result of the ability of high concentrations of PM 2.5 to trigger DNA damage (5). The identification of these upregulated genes in our dataset, which have been shown individually in other studies, validates the reliability and accuracy of our global RNA-Seq analysis.
To take a more unbiased view of potential pathways enriched by these differentially expressed genes, we used Advaita to perform pathway analysis and identified pathways such as Cytokine-Cytokine Receptor Interaction, Glycine Serine Threonine Metabolism, mitogen-activated protein kinase (MAPK) Signaling Pathway, and PI3K-Akt Signaling Pathway that were enriched in our dataset (Table 1). Supplementary Table S2 lists the differentially expressed genes associated with these pathways. Although elevation of cytokines and enrichment of pathways such as MAPK and PI3K-Akt have been described in other studies of PM 2.5 (9,11,12,16), pathways related to amino acid synthesis and metabolism and cardiovascular disease (Table 1) were uniquely found in our dataset. Network analysis of genes within these pathways allowed us to identify key interactions, including genes that may serve as regulatory hubs (Fig. 1C, Supplementary Fig.  S1).
The Advaita iPathway software allowed us to additionally identify "upstream" hubs or mediators that regulate many of the genes in this dataset that may be downstream. We identified P53 and IL6 as statistically significant downregulated and upregulated hubs, respectively (Fig. 1D, Supplementary Fig. S1). These unbiased enrichment analyses support what has often been described in the literature as the ability of PM 2.5 to induce DNA damage (often regulated by P53) (5,7,9) and cause increases in inflammatory cytokines (e.g. IL6) (17,18). Finally, the Advaita software also allows one to examine potential upstream chemicals or toxicants that are often associated with these differentially expressed genes, and upstream chemicals include organophosphorus compounds, polyaromatic hydrocarbons, and metals such as copper, silver, nickel, and cadmium, which are all often found in PM 2.5 .
We next used the GEO database and PubMed to identify other studies that also performed transcriptomic analysis in bronchial epithelial cells after PM exposure (Table 2). These studies utilized PM of varying size (2.5-10 mm) and from different locations [Saudi Arabia (15), China (11), United States (10, 13), Italy (12)], with treatments ranging from a concentration of 10-50 mg/cm 2 for 24 h. Excluding studies where data were incomplete or not publicly available, we compared the differentially expressed genes found in our study with those identified in other datasets (Fig. 2). Several genes, including NQO1, TXNRD1, CYP1A1, and CYP1B1, were found to be differentially expressed after PM exposure in all studies. Our dataset, however, also identified a larger number of genes (1176 of them) that were not identified in any other published dataset (Supplementary Table  S1). Although other studies (9, 16) identified enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways such as Cytokine-Cytokine Receptor Interaction, MAPK Signaling Pathway, and Rheumatoid Arthritis, that we also observed were enriched in our dataset, we found additional pathways such as Glycine Serine and Threonine Metabolism, Arrhythmogenic Right Ventricular Cardiomyopathy, Fluid Shear Stress and Atherosclerosis, and Complement and Coagulation Cascades that were enriched only in our dataset.
In all prior transcriptomic studies, the lowest concentration of PM 2.5 used for transcriptomic analysis was 10 mg/cm 2 PM 2.5 ; we and others, however, have shown that concentrations lower than 10 mg/cm 2 were sufficient to alter the expression of some genes, and in some instances, altered the expression of genes that were not affected by higher concentrations of PM 2.5 (20). We thus performed RNA-Seq on cells treated 24 h with 1 mg/cm 2 of PM 2.5 . The number of genes altered after exposure to a lowconcentration of PM 2.5 was considerably less compared to the number of differentially expressed genes after treatment with higher concentrations of PM 2.5 . In fact, of the 1310 genes that were differentially expressed after exposure to 30 mg/cm 2 of PM 2.5 , only 154 (12%) showed differential expression (based on log fold-change ! 0.6) after low-concentration PM 2.5 . Treatment with lower concentrations of PM 2.5 often resulted in gene expression changes that were in the same direction as higher concentrations, but with a lower magnitude of effect (Supplementary Table S3). However, some genes such as TIE1, IL1RL1, DUSP5, CXCL8, ADTRP, LEF1, GPR3, FOS, PTX3, SERPINE1, FOSL1, AQP3, PADI2, and BMF exhibited opposing directions of expression after either low or high concentrations of PM 2.5 . Using a less stringent statistical model based on likelihood ratios (which has traditionally been used in many transcriptomic studies) rather than F-tests, we identified 40 genes (using the same false discovery rate < 0.05) that were differentially expressed to a statistically significant degree after lowconcentration exposure to PM 2.5 (Table 3). Certain genes such as FYB, RCAN2, WNT9A, and CTGF were affected by lowconcentration but not high concentration of PM 2.5 , suggesting that expression of several genes may be uniquely sensitive to these lower concentrations.

Transcriptomic Changes in BEAS-2B Cells after Repeated Exposures for Seven Days
Experimental models that simulate chronic daily exposure have increasing relevance to humans who are exposed to low levels of pollution on an everyday basis. We have observed that expression of some genes is not altered when cells are exposed to a single dose of PM 2.5 , but increased when exposure is repeated daily (20). To examine the potential for transcriptomic changes that occur with this model, cells were treated every day with a low-concentration of PM 2.5 (1 mg/cm 2 ) for one week. To prevent the accumulation of additive doses, cells were washed each day between treatments. We identified 40 genes (16 upregulated and 24 downregulated) that were differentially expressed in cells treated with chronic daily exposure to PM 2.5 compared to   control ( Fig. 3A and Table 4). Analysis of gene ontology identified biological processes such as blood vessel development, cell migration, and response to external stimulus to be enriched in this group. The relative expression of genes within each gene ontology group is shown in Fig. 3B. Interestingly, of these genes affected by chronic daily exposure, only 4 genes, CYP1A1, CYP1B1, AHRR, and LINC00886, were also shown to be affected after a single 24-h exposure with the same concentration of PM 2.5 (1 mg/cm 2 ). The lack of many more genes in common between single (Table 3) and repeated (Table 4) treatments with low-concentrations of PM 2.5 emphasize the importance of considering duration of treatment on gene expression. Because there were more genes that were affected by a single 24-h exposure using higher concentrations (30 mg/cm 2 ) of PM 2.5 , we also compared the genes that were differentially expressed after a single exposure with high concentration of PM 2.5 with those that were differentially expressed after repeated, lowconcentration of PM 2.5 treatment. Only 16 of the 40 (40%) genes that were differentially expressed after chronic, repeated exposure to low-concentrations of PM 2.5 were also affected when cells were treated with a single, 24-h exposure of highconcentration PM 2.5 (Table 4). These data highlight the differing responses of bronchial epithelial cells when exposed to either single or repeated doses of PM 2.5 .

Single Exposure to PM 2.5 Induce DNA Methylomic Changes in BEAS-2B Cells
Exposure to air pollution, and PM 2.5 specifically, has been wellrecognized to cause changes in DNA methylation as measured in vivo (26)(27)(28) and in in vitro models (38,39). These DNA methylation changes are capable of altering the gene expression and phenotype of cells and organisms. We performed eRRBS to first examine the DNA methylomic changes in BEAS-2B cells after 24 h treatment with PM 2.5 . Overall coverage was similar between control and treated samples and among all replicates. An average of 5 million CpG sites was covered by eRRBS per sample. The distribution of coverage, annotated to CpG islands, shores, shelves, or intergenic regions is shown in Supplementary Fig. S2. Cells treated with 30 mg/cm 2 of PM 2.5 demonstrated widespread changes in DNA methylation throughout the genome (Fig. 4). There were 73 648 differentially methylated CpGs (DMCs) and 12 810 differentially methylated CpG regions (DMRs) between treatment and control (Supplementary Table S4). Among the DMRs, there were slightly more regions that were hypomethylated than hypermethylated. Hypermethylated regions were found disproportionately more among CpG islands and shores, whereas hypomethylated changes were found in greater proportion among intergenic regions (Fig. 4A and B). A breakdown of methylation changes relative to regions within genes is shown in Fig. 4B. Although PM 2.5 caused a slight increase in overall hypomethylation across the genome, a higher proportion of hypermethylation was found in the 1-5 kb upstream regions, promoters, and exon/intron boundaries. Conversely, a greater proportion of hypomethylation was found in 5 0 -untranslated regions. Volcano plots of differentially methylated regions and their annotations to gene region are shown in Fig. 4C. Pathway analysis among the genes that were differentially methylated identified a variety of different enriched pathways (Supplementary Table S5). Some of these were the same pathways as that observed among differentially expressed genes, such as MAPK Signaling Pathway, Pathways in Cancer, and Arrhythmogenic Right Ventricular Cardiomyopathy. However, there were many other additional pathways that were enriched among the differentially methylated dataset.
We next correlated changes in DNA methylation with gene transcription. When all of DMCs were taken into account, there were 9011 genes that had least one differentially methylated CpG locus. Likewise, the 12 810 DMRs annotated to 5525 gene IDs. Of the 1310 differentially expressed genes, 381 were found to have DMRs and nearly half (620/1310) were found to have at least one differentially methylated cytosine (DMC) (Fig. 4D). Supplementary Table S6 lists those genes with differential expression and at least one differentially methylated region (DMR). Among the genes with both differential expression and differential methylation, pathway analysis identified enrichment in the specific pathways Hypertrophic Cardiomyopathy, MAPK Signaling Pathway, PI3K-Akt Signaling Pathway, and Transcriptional Misregulation in Cancer (Fig. 4E).
DNA methylomic analysis was also performed on BEAS-2B cells treated with a low-concentration of PM 2.5 (1 mg/cm 2 ) for 24 h. Although the number of genes that were differentially expressed after exposure to PM 2.5 was less in the lowconcentration treatment group compared to the highconcentration group, many of the genes that were differentially expressed also demonstrated differential methylation (Table 5).

DNA Methylomic Changes in Cells Exposed to Daily Repeated Exposures
Cells exposed to PM 2.5 with repeated doses may develop different DNA methylation changes over time. We thus performed eRRBS on our seven-day exposure protocol. There were 18 641 DMCs and 3213 DMRs in cells that were treated daily for seven days with low-concentration PM 2.5 compared to control. Similar to what we observed with single-dose short-term exposure conditions, treatment of epithelial cells with repeated doses of PM 2.5 resulted in a greater number of hypomethylated than hypermethylated DMRs, and this difference was even more   pronounced than compared to our short-term treatment protocol (Fig. 5A). Although repeated treatment with PM 2.5 caused a greater degree of hypomethylation across the genome, a greater proportion of the hypermethylated regions were found in CpG islands and shores, with hypomethylated regions found more predominantly in intergenic regions (Fig. 5A-C). DNA hypermethylation was also found disproportionately more in promoters, exons, and exon/intron boundaries, whereas DNA hypomethylation was found more in 3 0 -and 5 0 -untranslated regions. This pattern was similar to what we had observed in our 24 h experimental protocols.
Overall, the 3213 DMRs that resulted after chronic, repeated exposure to PM 2.5 were annotated to 1844 genes. Nearly half of these genes were different from those that were differentially methylated after single exposure. The genes that were differentially methylated after chronic, repeated exposure to PM 2.5 also not surprisingly demonstrated enrichment of pathways that were different from the pathways identified from differentially methylated genes after single PM 2.5 exposure.
To determine whether changes in DNA methylation influences gene expression, we examined the DNA methylation levels of the 40 genes that were differentially expressed by PM 2.5 after repeated exposure. Of the 40 genes, AHRR, COL5A1, TNS1, and LINC00886 were identified as both differentially expressed and methylated (Table 5). We next correlated the degree and direction of DNA methylation changes with differential gene expression among cells treated with low-concentration of PM 2.5 , either for 24 h or repeatedly for seven days. Overall, genes that were hypomethylated were more often associated with increased expression and genes that were hypermethylated were associated with diminished expression (Fig. 5D).
Since treatment with PM 2.5 resulted in an overall greater degree of DNA hypomethylation across the genome, we next examined expression of DNA methyltransferases (DNMTs) and ten-eleven translocation (TET) enzymes, which promote hydroxymethylation and eventual demethylation, in our experimental conditions. Single treatment with PM 2.5 for 24 h had no effect on DNA methyltransferase (DNMT) expression and caused a decrease in TET expression, but daily treatment with low-concentration of PM 2.5 for seven days resulted in an increase in expression of TET1, TET2, and TET3 (Fig. 6).

Discussion
In this study, we describe the transcriptomic and DNA methylomic changes that occur in bronchial epithelial cells after 24 h of either high (30 mg/cm 2 ) or low (1 mg/cm 2 ) concentration treatment with PM 2.5 from Beijing, China. We also examine the transcriptomic and DNA methylomic changes that occur in cells after repeated exposure to 1 mg/cm 2 of PM 2.5 every day for seven days. Overall, we found widespread changes in the transcriptomic and DNA methylomic patterns of bronchial epithelial cells after PM 2.5 treatment. These changes differed depending on the concentration and duration of exposure. Treatment with 30 mg/cm 2 of PM 2.5 altered the expression of a large number of genes, many of which included cytokines and which were enriched in MAPK, PI3K-Akt, IL6, and P53 pathways. Although a lower number of genes were differentially expressed after treatment with 1 mg/cm 2 of PM 2.5 compared to 30 mg/cm 2 of PM 2.5 , repeated exposure to low-concentrations of PM 2.5 over several days resulted in the differential expression of many genes that were interestingly, different from those observed after 24 h of PM 2.5 treatment. Treatment with PM 2.5 resulted in extensive changes in DNA methylation, and many of these DNA methylation changes were found in genes that also demonstrated differential expression. Overall, there were more hypomethylated changes compared to hypermethylated changes after PM 2.5 exposure, especially in cells treated repeated for seven days with PM 2.5 . More hypermethylated changes were noted in CpG islands and shores while hypomethylated changes were predominantly in intergenic regions. Together, these findings demonstrate the ability of PM 2.5 to induce widespread changes in gene expression and DNA methylation under different concentrations and duration of PM 2.5 exposure. The extent of these changes provides insight into how PM 2.5 affects airway epithelial cell biology, inflammatory responses, and ultimately, disease.
Treatment with high concentrations of PM 2.5 (30 mg/cm 2 ) induced upregulation of a number of genes that were associated with toxicological responses to PM 2.5 and that have been described in the literature. We observed upregulation of CYP1A1, CYP1B1, and AHRR, which is consistent with the polyaromatic hydrocarbon content of PM 2.5 (7,35). PM 2.5 is also known to induce oxidative stress, commonly through the transcription factor NRF2 (36,37). NRF2 binds to the anti-oxidant response and -3b after single, 24-h treatment with PM 2.5 . *P < 0.05, one-way ANOVA with Tukey's multiple comparisons post-test. element and can activate genes such as HMOX1, NQO1, GCLM, SQSTM1, and FTL, which were all increased in the dataset. Perhaps not surprising were also increased expression of MDM2 and CDKN1A, which were expected to occur with PM-associated DNA damage (5), or increased expression of cytokines such as IL6, IL1B, TSLP, CXCL2, and CXCL3, which have all been demonstrated to be upregulated in other studies (11,12,16). Metals such as silver, copper, cadmium, and nickel, which are present in PM 2.5 , were not surprisingly identified in our unbiased pathway analyses as enriched upstream chemical regulators. These findings emphasize the validity and reproducibility of our results in the context of expected changes reported in the literature.
Our RNA-Seq data, however, also identified a much larger number of genes, 1310, that were differentially expressed after treatment with 30 mg/cm 2 of PM 2.5 , which demonstrates the broad extent by which PM 2.5 is capable of altering gene transcription in airway epithelial cells. Even when employing a more stringent cutoff for fold-change, at least 709 of these genes had at least greater than 2-fold-change in expression. We identified a much larger number of differentially expressed genes than that reported in most other studies (9)(10)(11)(12)(13)(14)(15)(16), which may be attributed to the fact that we employed RNA-Seq, which is a more sensitive technique and provides a more comprehensive and broader analysis of transcriptomic changes than microarray analysis. By identifying a larger number of genes that were not previously noted in prior studies of PM 2.5 , we could discover additional pathways based on gene ontology that are altered by PM 2.5 , and identify novel genes that may provide additional insights into the effects of PM 2.5 on epithelial cell biology. Certain genes such as CYP1A1, CYP1B1, NQO1, and TXNRD1 were differentially expressed in all studies (9)(10)(11)(12)(13)(14)(15)(16) including ours, demonstrating common, universal effects of PM 2.5 , even when they are from different regions of the world. Likewise, pathway analysis identified many gene ontology pathways that were commonly activated among different studies of PM 2.5 , including Cytokine-Cytokine Receptor Interaction and PI3K-Akt and MAPK Signaling. However, we also identified Glycine Serine and Threonine Metabolism, Arrhythmogenic Right Ventricular Cardiomyopathy, Fluid Shear Stress and Atherosclerosis, and Complement and Coagulation Cascades as potentially important pathways that were enriched only in our dataset.
Unique to our study was the fact that we examined gene expression changes with concentrations of 1 mg/cm 2 of PM 2.5 , which was lower than that used in other transcriptomic studies. The rationale to study the effects of PM 2.5 at this lowconcentration was based on the fact that we had observed differential gene expression changes occurring with this concentration (sometimes even without changes at higher concentrations) (20), and low-concentrations of PM 2.5 have been shown to exert important effects in other studies (17)(18)(19). We found fewer genes that were differentially expressed after treatment with a low-concentration of PM 2.5 compared to higher concentrations, and the magnitude of effect was often less after lower concentrations of PM 2.5 compared to higher concentrations. There were, however, certain genes that were affected by low-concentration PM 2.5 treatment that were not observed after high-concentration treatment, and this illustrates how low levels of PM 2.5 exposure cause changes that are distinct from higher concentrations of PM 2.5 exposure. This is consistent with what we had observed in focused studies of select genes (20), though the mechanism for this is unclear. One potential hypothesis is that low-concentrations of PM 2.5 might activate transcription factors that are opposed by other transcription factors that are activated at only higher concentrations of PM 2.5 . Additionally, signaling pathways activated by PM 2.5 at low concentrations might be inhibited by pathways activated by only high concentrations, such as when anti-oxidant pathways become activated to counter the effects of oxidative stress. Finally, different concentrations of PM 2.5 may activate signaling pathways with different kinetics and some genes may be activated by only certain pathways and not others. Ultimately, there were not enough genes activated by only low-concentration but not high-concentration PM 2.5 that allowed us to perform rigorous pathway analyses, but future studies with increased sample size may enhance the statistical power to identify even more genes uniquely affected by low-concentration of PM 2.5 .
Studying the effects of repeated treatment with lowconcentrations of PM 2.5 over time has the potential to further model real-world exposure, since much of the world population are exposed to low levels of PM 2.5 on a chronic, daily basis. Here, we found that repeated exposures induced gene expression changes that are distinct from what was observed after a single, 24-h exposure protocol. This was true even when comparing the effects of chronic exposure to the 24-h changes associated with higher concentrations of PM 2.5 , indicating that the effects of repeated exposure cannot be explained merely by the arithmetic total of accumulated PM 2.5 . This differential effect of PM 2.5 between single and chronic, repeated treatment protocols has been described in other studies as well, in the context of epithelial mesenchymal transition and histone modifications (17)(18)(19). Along with these other studies, our data emphasize both the need and importance of utilizing models of repeated treatment to assess the effects of chronic exposure to PM 2.5 and highlights the profound changes, both transcriptomic and epigenetic, that occur when cells are exposed to PM 2.5 over time. The mechanisms for this may vary but could include adaptive changes in cells that are chronically exposed to PM 2.5 , as often observed when anti-oxidant signals become upregulated in response to oxidant damage. Epigenetic modifications may account for some of these adaptive changes, and many studies have examined how air pollutants affect DNA methylation in humans, not just acutely, but over the life course of an individual (26,28,40). One study showed that exposure of BEAS-2B cells to biomass from a power plant for five weeks induced widespread DNA methylation and gene expression differences (38). Exactly how epigenetic and gene expression changes continue to evolve over time remains to be determined.
DNA methylation changes represent an important epigenetic mechanism that affects gene expression, often in a persistent manner, and that contributes to the development of many diseases including cancer (21) and asthma (23). Both in vivo and in vitro studies have shown that air pollution affects DNA methylation patterns in a variety of cell types (26)(27)(28). Here, we utilized eRRBS to examine the DNA methylomic changes in bronchial epithelial cells after PM 2.5 treatment and found that PM 2.5 induced widespread DNA methylomic changes, many of which correlated with changes in gene expression. These changes occurred with both short-term and long-term PM 2.5 exposures, though the genes that were affected differed between exposure protocols. As compared to other studies in air pollution that utilized microarrays for DNA methylation analysis (38,39), which are biased toward pre-determined CpG sites, we employed eRRBS, which assays a more diverse region of genes and broader coverage. Our dataset thus has the potential to provide richer insights into the epigenetic changes induced by PM 2.5 and the mechanisms by which PM 2.5 affects individual gene expression. However, eRRBS has its own limitations, and in our study, only 47.6% of Illumina 450 K array and 37.3% of Illumina EPIC array (850 K) sites are covered by eRRBS. Thus, when compared to the study by Shi et al. (39) , only 16 genes  (ABCA3, ABCG1, ANPEP, LMNB2, GALNT2, HMGA1, HPCAL1,  KIF18B, PTK2B, MID1, CENPE, SEMA6B, TNS3, TPX2, EPHB4,  TNFRSF10B) were noted to be differentially methylated and expressed by both their study and ours. This might be due to variations in experimental protocol, source of PM 2.5 , and methods for DNA methylation analysis (Illumina array vs. eRRBS).
Although PM 2.5 induced an overall greater number of CpG sites that were hypomethylated compared to hypermethylated, there was increased hypermethylation found among CpG islands and CpG shores. This differential effect suggests that DNA methylation changes induced by PM 2.5 occur in a locusspecific manner rather than stochastically throughout the genome. Within genes, a greater proportion of hypermethylation was found among promoters and exons whereas hypomethylation was found more among intergenic regions and 3 0 -and 5 0untranslated regions. This overall greater degree of hypomethylation with certain regions biased toward hypermethylation has been observed in other in vitro studies as well (38,39). Together, these findings further emphasize the non-random, directed nature by which methylation changes occur throughout the genome.
Measurements of DNMT and TET expression after PM 2.5 exposure revealed that expression of TET1, TET2, and TET3 were elevated after treatment with PM 2.5 when given repeatedly for seven days. This is consistent with the observation that PM 2.5 induces an overall greater degree of hypomethylation across the genome, since TET enzymes are recognized to participate in demethylation by the addition of hydroxyl groups to methylcytosine followed by base-excision and repair. The effects of PM 2.5 on TET expression and DNA hypomethylation are consistent with certain other studies that demonstrate an overall increase in hypomethylation associated with air pollution (41,42), though the data in vivo tend to be mixed. In one, well-designed double cross-over interventional study performed in humans, PM 2.5 exposure resulted in a decrease in the predominant DNMT and TET isoforms, and this was associated with a decrease in both methylated and hydroxymethylated cytosine (43). Although our data with chronic repeated exposures demonstrated an increase in TET expression, a short-term single exposure was associated with a decrease in TET expression, which would thus be consistent with this recently published finding in humans. The differing effects of PM 2.5 on TET expression in our study once again highlight the importance of considering duration of exposure in different models. Although several studies demonstrate that PM 2.5 induces a global level of DNA hypomethylation (41,42), our study, in addition to these others, identified specific CpG sites that were hypermethylated as well. No changes in DNMT expression were noted after PM 2.5 exposure in our experiments, and future studies are needed to determine the exact mechanism by which specific loci are targeted for either hyper-or hypomethylation.
Several of the observed changes in DNA methylation correlated with changes in gene expression, though there were also many changes in DNA methylation that were annotated to genes whose expression was not significantly altered. The effects of changes in DNA methylation often depend on the number of CpG loci affected, the magnitude of change, and the location of these changes relative to any given gene, and our study may have identified many methylation changes that were insufficient to effect gene expression by itself. However, these DNA methylation changes may render genes poised for further modulation, and additional studies can test whether these methylation changes affect the sensitivity of cells to subsequent treatment with various mediators, cytokines, or other signals that activate transcription factors. The fact that DNA methylation patterns persist and can be passed during cell division indicates that these PM 2.5 -induced changes may have important implications in long-term health and disease, even after the cessation of PM 2.5 exposure. Additional studies would be needed to determine the durability and long-term stability of these DNA methylation changes over time.
Our study examined the gene expression and DNA methylation changes induced in BEAS-2B cells, but use of primary airway epithelial cells, from healthy subjects or individuals with airway disease, could produce different findings. We also treated cells directly with PM 2.5 dissolved in aqueous medium, and experiments performed using an air-liquid interface may produce different results, as suggested by other studies (44,45). Use of air-liquid interface may also obviate the need to account for how PM 2.5 is extracted, as varying levels of organic fractions has been shown to independently affect transcriptomic results (13). The lowest concentration we used was 1 mg/cm 2 , which was much lower than many other in vitro studies. This may be similar to the low-concentrations of PM 2.5 that humans are exposed to in developed countries, though calculating the equivalency between in vitro and in vivo concentrations has its limitations. We utilized PM 2.5 collected from air monitoring stations in Beijing, China, and similarities between our dataset and others suggest that PM 2.5 , even from different sources, exert common, universal effects. By contrast, transcriptomic and methylomic changes that were uniquely identified in our study could also indicate changes specific to PM 2.5 from Beijing. DNA methylation changes have been shown to occur in controlled, in vivo studies, performed in both rats (39) and humans (27). Although there are limitations to the applicability of our study to real-world, in vivo exposures, in vitro studies provide the advantage of systematically examining the effects of different concentration and exposure protocols of PM 2.5 on a single cell type in isolation. Repeated exposures also better approximate real-world conditions in which humans are often exposed to PM 2.5 on a daily basis. The fact that we observed very distinct gene expression and DNA methylation signatures from cells treated on a repeat vs. single basis highlights the need for more studies utilizing chronic models of exposure. Future studies that examine the upstream signals that are triggered by PM 2.5 , such as generation of reactive oxygen species or activation of nuclear factor kappa B or aryl hydrocarbon receptor, can help shed light on the mechanism(s) by which PM 2.5 causes these transcriptomic and DNA methylomic changes.
In conclusion, we identified transcriptomic and DNA methylomic changes that occur when bronchial epithelial cells are exposed to PM 2.5 from Beijing, China. Using RNA-Seq and eRRBS and comparing high vs. low-concentration of PM 2.5 at both single and repeated exposures, we identified different genomewide changes at the gene expression and epigenetic level that are unique and distinct from existing literature. These datasets demonstrate the extensive nature by which PM 2.5 influences the gene transcriptome and epigenome of airway epithelial cells, which may have important implications in long-term health and disease.