Disease relevant modifications of the methylome and transcriptome by particulate matter (PM_2.5) from biomass combustion

ABSTRACT

Exposure to particulate matter (PM) is recognized as a major health hazard, but molecular responses are still insufficiently described. We analyzed the epigenetic impact of ambient PM_2.5 from biomass combustion on the methylome of primary human bronchial epithelial BEAS-2B cells using the Illumina HumanMethylation450 BeadChip. The transcriptome was determined by the Affymetrix HG-U133 Plus 2.0 Array. PM_2.5 induced genome wide alterations of the DNA methylation pattern, including differentially methylated CpGs in the promoter region associated with CpG islands. Gene ontology analysis revealed that differentially methylated genes were significantly clustered in pathways associated with the extracellular matrix, cellular adhesion, function of GTPases, and responses to extracellular stimuli, or were involved in ion binding and shuttling. Differential methylations also affected tandem repeats. Additionally, 45 different miRNA CpG loci showed differential DNA methylation, most of them proximal to their promoter. These miRNAs are functionally relevant for lung cancer, inflammation, asthma, and other PM-associated diseases. Correlation of the methylome and transcriptome demonstrated a clear bias toward transcriptional activation by hypomethylation. Genes that exhibited both differential methylation and expression were functionally linked to cytokine and immune responses, cellular motility, angiogenesis, inflammation, wound healing, cell growth, differentiation and development, or responses to exogenous matter. Disease ontology of differentially methylated and expressed genes indicated their prominent role in lung cancer and their participation in dominant cancer related signaling pathways. Thus, in lung epithelial cells, PM_2.5 alters the methylome of genes and noncoding transcripts or elements that might be relevant for PM- and lung-associated diseases.

KEYWORDS:

• biomass combustion
• disease ontology
• methylome
• particulate matter
• transcriptome

Introduction

Exposure to ambient particulate matter (PM) is recognized as a major health hazard, associated with acute respiratory infections, lung cancer, and chronic respiratory or cardiovascular diseases.^1,2 Whereas industrial and traffic-related pollutants have been extensively investigated, the number of studies evaluating the impact of air pollution resulting from burning biomass is limited.

Worldwide, biomass fuel represents about 10% of direct human energy consumption and will remain the main source of energy for most of humanity.^3-6 According to the World Health Organization, more than 3 billion people are affected by emissions from biomass fuels, claiming globally over 4 million deaths per year.⁷ Epidemiologic and exposure studies have associated wild land fires, agricultural burning, and residential wood combustion with a variety of adverse respiratory health effects, including asthma, chronic obstructive lung disease, and acute lower respiratory tract infections, but also with cardiovascular effects and lung cancer.^3,6 Thus, emissions from residential energy use have the largest impact on premature mortality worldwide.⁵ The use of biomass fuels and wood burning is still expected to rise due to high fossil energy costs, the desire to use renewable, CO₂-neutral energy sources, the rise in population, need for cropland, and climate changes.³ Even in developed countries, residential biomass combustion is a major source of particle emissions.³

Emissions from industrial or traffic related sources are essentially under legislation, but residential combustion of biomass fuels, agricultural burning, and forest fires are largely uncontrolled and may result in particle concentration up to 3 mg/m³ or even higher.^3,8 Particles less than 2.5 μm (fine particles; PM_2.5), are considered to be particularly harmful, as they efficiently evade the mucociliary defense system and are deposited in the peripheral airways, where they can exert their pathophysiological mechanisms, linking exposure to ambient fine particles to diseases.⁹

Emissions from biomass fuels usually contain this harmful size range. Their chemical composition substantially depends on fuel composition and combustion characteristics and is different from industrial or traffic-related particles on which most health studies have focused.³ However, molecular studies have demonstrated that particulate matter from different sources may all induce oxidative stress and cellular inflammatory responses, potential mechanisms linking PM to respiratory or cardiovascular diseases and lung cancer.^10-13

Adaption to environmental changes occurs by epigenetic chromatin modifications resulting in durable alterations in gene expression.¹⁴ Indeed, for most of PM-associated diseases, epigenetic aberrations have been described.¹⁵ Correspondingly, changes in the DNA methylation pattern have been observed in the presence of ambient PM. Initially, studies determining methylation changes and its disease relevance have mainly focused on specific candidate target genes or repetitive elements in blood DNA, e.g., after exposure to black carbon, cigarette condensate, or traffic-related particles.^16-18 Only lately, global methylation changes were determined upon high-occupational PM_2.5 exposure.¹⁹ Despite the significant impact of PM from biomass combustion on health, comprehensive genome-wide methylation studies are still missing.

In our study, we generated profiles of the human methylome and transcriptome by microarray technologies in the presence of PM_2.5 from biomass combustion using human airway epithelial cells (BEAS-2B) that represent the first barrier to inhaled extraneous particles. Exposure to PM_2.5 from biomass combustion resulted in genome-wide changes of the methylome including repetitive elements and noncoding transcripts that may be relevant for lung cancer, inflammation, asthma, and other PM-associated diseases. Additionally, we filtered 152 differentially expressed genes, which were either hyper- or hypomethylated. From these genes, 66 were related to lung diseases, especially to lung cancer. Our results provide first insights on how changes in the methylation pattern of genes induced by chronic exposure to PM from biomass combustion may especially influence the development of various lung diseases.

Results

Characterization of PM_2.5 from biomass combustion

PM_2.5 was obtained from a biomass power plant. Mineralogical analysis revealed that the particles mainly consisted of crystalline arcanite (K₂SO₄) or amorphous components (for detailed composition see Table S1 and Dornhof et al.).²⁰

Methylome of PM_2.5-exposed BEAS-2B cells

To analyze whether PM affects the human methylome, human bronchial epithelial BEAS-2B cells were persistently exposed to PM_2.5 from biomass combustion to accumulate DNA methylation changes. Cytotoxic effects were excluded by the absence of caspase-3/7 activity and DNA fragmentation (Fig. S1). Genome-wide DNA methylation was determined by using the Illumina Infinium HumanMethylation450 BeadChip (450K) Array. After 5 weeks, differential DNA methylation was observed in PM_2.5-exposed BEAS-2B cells at 1.29% of the annotated sites (5501 CpGs) (Fig. 1A). The majority of alterations in DNA methylation were related to hypomethylation (decrease in Δβ > 0.02; 4234 CpGs; Fig. 1B). Some individual differentially methylated CpGs were annotated to several different sites in relation to gene region or distance to CpG island (Fig. 1C/D).

Figure 1. Methylome of BEAS-2B cells exposed to PM_2.5 for 5 weeks. (A) Percentage of differentially methylated CpG sites in untreated BEAS-2B cells vs. PM_2.5-exposed cells (change in Δβ>0.02). (B) Percentage of hypermethylation and hypomethylation. (C) Functional genomic distribution of hypomethylated CpGs (upper pie chart) or in relation to CpG island regions (lower pie chart). (D) Functional genomic distribution of hypermethylated CpGs (upper pie chart) or in relation to CpG island regions (lower pie chart).

CpG hypo- and hypermethylation events in PM_2.5-exposed BEAS-2B cells occurred frequently in the body of genes or unrelated to any gene region. Because methylation changes at CpG island-associated promoter sites have high relevance for gene expression,¹⁴ these sites were further investigated. The promoter region, including CpGs located in the 5′UTR, the TSS200, and the TSS1500, was hypomethylated at 35.45% (1501 CpGs) and hypermethylated at 28.02% (355 CpGs) of the differentially methylated sites (Fig. 1C/D). However, most of the hypo- and hypermethylation events, observed in the proximal promoters of PM_2.5-exposed cells occurred in promoters not associated with CpG islands (Figs. S2/S3). Only 9.6% (144 CpGs) of the hypomethylated sites and 15.21% (54 CpGs) of the hypermethylated sites in the promoter region were associated with CpG islands.

Additionally, the number of the differentially methylated CpG sites in the regions island, N-shelf, N-shore, open sea, S-shelf, and S-shore were compared with the expected number of CpG sites, assuming that methylation changes appear equally distributed in all regions. However, it became evident that CpG sites located in islands appeared less frequently (negative enrichment, p<0.05) compared with a relatively high number of differentially methylated CpG sites (positive enrichment, p<0.05) located in the shelf or open sea region (Table S2).

Chromosomal distribution of DNA methylation of PM_2.5-exposed BEAS-2B cells

Next, the chromosomal distribution of PM_2.5 dependent DNA methylation was studied (Fig. 2, Table S3). DNA methylation changes were distributed on all analyzed chromosomes. Most methylation changes were detected on chromosome 1 (≥10% of either the hypo- or hypermethylated sites). In relation to the total amount of annotated CpGs per chromosome, hypomethylation was observed more frequently on chromosomes 1, 16, 17, and 22 (each ≥1.25%; enrichment tested by χ²-test, p < 0.05), whereas hypermethylation was enriched on chromosome 5 (≥0.4%; p < 0.05). Moreover, differential methylation was observed in both G- and R-bands. The fraction of hypomethylated CpGs compared with hypermethylated CpGs was significantly higher in R-bands (2868 hypomethylated CpGs; 746 hypermethylated CpGs) than in G-bands, which stain positive for Giemsa (1353 hypomethylated CpGs; 507 hypermethylated CpGs; χ²-test, p < 0.05). R-bands are known for their high gene density.

Figure 2. Ideogram of the human karyotype. R-bands are shown in white, G-bands in gray; the centromere and the coding region for rRNA are shown in red, regions of hypomethylation in blue, and regions of hypermethylation in yellow. The gender-specific XY-chromosomes were not further analyzed, because BEAS-2B cells originate from various individuals.

Functional clustering of differentially methylated genes

Differently methylated genes were functionally clustered by a gene ontology analysis, using the Enrichr web application.²¹ Genes with hypomethylated CpGs were most significantly enriched in pathways associated with GTPase activity, regulatory signal transduction to growth factor stimuli and acidic chemicals, or coded for proteins involved in the architecture of the extracellular matrix (p <1.0e-6; Fig. 3A). Genes with hypermethylated CpG sites were significantly enriched in pathways connected to cellular adhesion and pathways responsible for ion binding, transport, or channel activity (p < 1.0e-5; Fig. 3B). The affected genes are listed in Tables S4/S5.

Figure 3. Gene ontology of hypomethylated genes (A) or hypermethylated genes (B) following persistent exposure to 100 µg/ml PM_2.5 for 5 weeks. The 12 most significant functionally linked biologic processes, the number of clustered genes, the p-values and combined scores (in brackets) are shown.

PM_2.5 treatment affects tandem repeats and the miRNA-specific methylome

Differential methylation of repetitive DNA by environmental pollutants might have serious consequences for the structural integrity of the genome and on transcriptional activity, consequently affecting human health.²² Therefore, the methylation profile of these tandem repeats was also analyzed. Interspersed repeats (long interspersed nuclear elements, LINEs; short interspersed nuclear elements, SINEs; long-terminal repeats, LTRs) satellite DNA (simple sequence repeats, SSRs; satellite) and low complexity regions (LCRs) were included in the study (Table 1). Reduced methylation of tandem repeats was prominent, but also some DNA loci were hypermethylated, indicating potentially impaired functionality of these sites.

Table 1. PM-dependent methylation changes of tandem repeats.

Tandem repeat	LINE	SINE	LTR	SSR	LCR	Satellite DNA
Hypomethylated repeats	202 (1.47%)	199 (1.35%)	92 (1.10%)	8 (0.39%)	3 (0.10%)	0 (0.00%)
Hypermethylated repeats	74 (0.54%)	54 (0.37%)	46 (0.55%)	3 (0.15%)	3 (0.10%)	4 (1.50%)
Analyzed repeats in total	13771	14793	8338	2043	2898	267

Furthermore, the differential methylation of noncoding transcripts (miRNAs, miRs) was analyzed using positional data from the miRBase database (http://mirbase.org). Persistent exposure to PM_2.5 significantly modified the methylation profile at 45 CpG loci annotated to 43 miRs. A large part of these CpGs were located proximal to the promoter and were hypomethylated, indicating a potential mechanism for their transcriptional activation (Table 2).

Table 2. Disease relevance of miR with a PM-dependent, differential DNA methylation pattern.

miRNA	CpG locus	Δβ-value^d	p-value	Lung cancer/neoplasms	Asthma	Interstitial lung disease	Pulmonary embolism	Cancer/neoplasms	Inflammation	Infection, viral response	Heart disease/heart failure	Cardiomyopathy, muscular disorders	Lipid disorders	Atherosclerosis	Vascular disease	Hypertension	Diabetes	Stroke	Neurological disorders
miR-23a	TSS200/S_Shelf	−0.029	0.0008100	x				x	x	x	x	x			x
miR-24-2	TSS1500/S_Shelf	−0.029	0.0008100	x				x	x	x	x	x			x		x		x
miR-27a	TSS1500/S_Shelf	−0.029	0.0008100	x				x		x	x	x			x
miR-122	TSS1500/OpenSea	−0.077	0.0003100					x		x	x		x			x		x	x
miR-126	TSS1500/N_Shore	−0.051	0.0005800	x	x			x	x		x	x		x	x		x		x
miR-134	TSS1500/OpenSea	0.029	0.0004700	x			x	x			x								x
miR-299	Body/OpenSea	0.049	0.0000600					x			x	x
miR-323b	TSS1500/OpenSea	0.039	0.0002700					x
miR-338	TSS200/N_Shelf	−0.131	0.0000010	x				x											x
	TSS200/N_Shelf	−0.102	0.0000300	x				x											x
miR-365a	TSS200/OpenSea	−0.072	0.0006900	x				x	x	x	x
miR-382	TSS1500/OpenSea	0.029	0.0004700					x		x	x
miR-425	TSS1500/N_Shore	−0.033	0.0002500					x
miR-485	TSS1500/OpenSea	0.039	0.0002700	x	x			x											x
miR-570	TSS1500/OpenSea	−0.062	0.0000010					x
miR-591	TSS200/OpenSea	−0.072	0.0006400					x
miR-657	TSS1500/N_Shelf	−0.076	0.0007600					x
	TSS1500/N_Shelf	−0.131	0.0000010					x
	TSS1500/N_Shelf	−0.102	0.0000300					x
miR-661	TSS200/S_Shore	−0.034	0.0004000					x
	TSS200/S_Shore	−0.059	0.0004000					x
	TSS200/S_Shore	−0.056	0.0003600					x
miR-668	TSS200/OpenSea	0.039	0.0002700			x
miR-675	TSS200/Island	−0.057	0.0001300					x	x
miR-1249	TSS200/N_Shore	−0.062	0.0005800					x
	TSS1500/N_Shore	−0.067	0.0000090					x
	TSS1500/N_Shore	−0.062	0.0006700					x

Table

miRNA	CpG locus	Δβ-value^d	p-value	miRNA	CpG-locus	Δβ-value^d	p-value
miR-380	TSS1500/OpenSea	0.049	0.0000600	miR-4512	TSS1500/S_Shore	−0.022	0.0000100
miR-584	TSS1500/OpenSea	−0.052	0.0003300		TSS1500/S_Shore	−0.031	0.0000400
miR-662	TSS1500/OpenSea	−0.047	0.0004200	miR-4686	TSS1500/OpenSea	0.062	0.0004800
miR-1208	Body/OpenSea	−0.071	0.0000600	miR-5091	TSS1500/Island	0.053	0.0005300
miR-1250	TSS200/N_Shelf	−0.047	0.0006000	miR-6753	TSS1500/OpenSea	−0.070	0.0000001
miR-3065	TSS1500/N_Shelf	−0.088	0.0001800		TSS1500/OpenSea	−0.166	0.0000020
	TSS1500/N_Shelf	−0.076	0.0007600	miR-6775	TSS1500/OpenSea	−0.042	0.0001300
miR-3189	TSS1500/N-Shore	−0.035	0.0001500	miR-6787	TSS1500/N_Shore	−0.054	0.0000700
miR-3195	TSS1500/Island	−0.050	0.0000300		TSS200/Island	−0.047	0.0002700
miR-3912	TSS1500/N_Shore	0.026	0.0007200	miR-6887	TSS1500/N_Shelf	−0.069	0.0000900
miR-4297	TSS1500/OpenSea	0.380	0.0008400	miR-7107	TSS200/Island	−0.021	0.0000500
miR-4425	TSS1500/OpenSea	−0.027	0.0005200		TSS1500/Island	−0.032	0.0003900
miR-4496	TSS1500/OpenSea	−0.030	0.0004600	miR-7846	TSS1500/N_Shore	−0.068	0.0001300
miR-4496	TSS1500/OpenSea	−0.071	0.0000003

^a miRs identified by the miRBase database (http://mirbase.org);

^b disease association according to the human microRNA Disease Database v.2.0 (http://cuilab.cn/hmdd);

^c miRs with known contribution to PM-associated diseases are listed in the upper part of Table 2, miRs with yet unidentified disease relevance are listed below;

^d negative values: hypomethylated; positive values: hypermethylated

The direct disease association of miR was examined by the human microRNA Disease Database v.2.0 (http://cuilab.cn/hmdd). We found that miRs with a distinct methylation pattern were prominently involved in lung and PM-associated diseases, such as lung cancer and neoplasms, asthma, interstitial lung disease, or pulmonary embolism (Table 2). Moreover, miRs were linked to different types of cancer, regulation of inflammation and viral responses, which may also adversely affect lung physiology. A high proportion of miRs could be allocated to heart diseases. Several miRs contribute to the metabolic syndrome, and affect hypertension, diabetes, or lipid metabolism. All these diseases can be related to high ambient PM exposure.^1,2,23

Effect of gene methylation on gene transcription

To determine the biologic significance of differential DNA methylation, gene expression profiling of PM_2.5-exposed BEAS-2B cells was performed using the Human Genome U133 Plus 2.0 Array. Subsequent correlation of the transcriptome with changes in methylation resulted in 152 individual genes that show significant congruence (p-value threshold ≤ 0.05; Tables S6/S7). Mapping of expression vs. CpG methylation of these genes indicated a trend toward transcriptional activation by hypomethylation (Fig. 4). Accordingly, 78.42% of the hypomethylated CpGs were positioned on upregulated genes and located in the upper left quadrant. However, a correlation of hypermethylation and gene repression was not evident (lower right quadrant, 51.28% of hypermethylated CpGs). This trend did not change significantly, when only CpG island or promoter-associated alterations in the methylation pattern were included in the analysis (data not shown).

Figure 4. Mapping of gene expression vs. DNA methylation (p-value threshold ≤ 0.05). Each dot represents a single CpG, located on a gene that is both differentially methylated and expressed upon persistent exposure to 100 µg/ml PM_2.5 for 5 weeks. Hypomethylated CpGs located on genes that cluster in the upper left quadrant were upregulated, whereas hypermethylated CpGs located on genes that were downregulated clustered in the lower right quadrant. The number of CpGs in each quadrant is indicated (UL, upper left quadrant; UR, upper right quadrant; LL, lower left quadrant; LR, lower right quadrant).

Gene ontology analysis using the Enrichr web application allocated genes that show both differential CpG methylation and transcription to functional properties (Fig. 5, Tables S8/S9). Functional enrichment of downregulated genes revealed that PM_2.5 largely affected common pathways related to cytokine signaling, cell motility, differentiation, and development. Upregulated genes were largely involved in biologic responses to extracellular stimuli, e.g., chemotaxis, cell growth, wound healing, inflammation, angiogenesis, and responses to foreign particles or molecules.

Figure 5. Functional enrichment of genes that show both differential methylation and transcription following persistent exposure to 100 µg/ml PM_2.5 for 5 weeks identified by the Enrichr web server. The top 14 of the most significant functionally linked pathways, p-values, combined score, and number of clustered genes are indicated.

PM_2.5-associated changes in the methylome and transcriptome are linked to human disease markers

Exposure to PM_2.5 has been epidemiologically related to pulmonary and cardiovascular disease, metabolic disorders, and lung cancer.^1,2,23 Therefore, a disease ontology approach was used to find a causal relationship between alterations in the methylome and transcriptome of PM_2.5-exposed BEAS-2B cells and diseases known to be exacerbated by PM. The 152 genes that showed both exposure-associated alterations in DNA methylation and gene expression were included in the analysis. Biocomputational profiling using the human disease methylation database (http://bioinfo.hrbmu.edu.cn/diseasemeth), the Dragon Database for Methylated Genes in Diseases (http://cbrc.kaust.edu.sa/ddmgd/), and the Human Gene Database GeneCards® (http://genecards.org) revealed 66 molcular regulators mainly associated with lung diseases, in particular with lung cancer (Table 3). Extracted genes involved in lung independent diseases are listed in Table S9.

Table 3. Differentially methylated and expressed genes and allocated lung-specific and lung-associated diseases.

^# for references see supplemental Table S8;

^a DiseaseMeth (http://bioinfo.hrbmu.edu.cn/diseasemeth);

^b GeneCard^® (http://genecards.org);

^c Dragon Database for Methylated Genes and Diseases (http://cbrc.kaust.edu.sa/ddmgd/);

^d according to data base search;

+promoted/overexpressed;

−suppressed/downregulated;

* when absent or inhibited; x involved but not specified

Next, the cellular and molecular relevance of the genes associated with lung cancer was investigated by a Medline search in more detail (data are included in Table 3, for references see Table S10). A significant number of these genes are upregulated in tumors and modulated by methylation. Some of these genes are implicated in tumor suppression and improved prognosis, whereas other genes were found to be largely involved in cellular transformation, proliferation and cell growth, tumor progression, metastasis, angiogenesis, inhibition of apoptosis, and therapeutic resistance. These tumor-promoting genes are prominently associated with Wnt, Akt, growth factor or hypoxia-related signaling pathways, epithelial-mesenchymal transition (EMT), and the circadian system. All of these pathways are central in promoting cellular survival and tumorigenesis.^24-30

Discussion

This study demonstrates that long-term exposure to PM_2.5 results in genome-wide alterations of the methylome and significantly alters gene expression. Predominantly, we observed hypomethylation of CpG sites, which is consistent with other studies, analyzing tandem repeats or individual candidate genes involved in cancer or cardiovascular disease upon exposure to polluted air.^18,31 Extensive DNA hypomethylation can be induced by oxidative stress.³² Correspondingly, we detected increased intracellular levels of reactive oxygen species in chronically PM_2.5-exposed lung epithelial cells (Regina Dornhof, unpublished observation).

Although methylation changes were also found in the promoter region, only a minor part of those were located within CpG islands. The methylation pattern of CpG islands within the promoter region seems to be most critical for gene expression,¹⁴ but methylation outside this region may also be functionally relevant, because these CpG sites are evolutionarily conserved and not subject to negative selection. Accordingly, it was shown that changes in the methylation pattern of non-CpG island promoters have almost an identical impact on gene expression and disease relevance as changes in CpG island promoters.³³

PM-mediated variations in DNA methylation were distributed among most chromosomes, but accumulated at regions with high gene density. The fraction of hypomethylated to hypermethylated CpG sites was enriched in the transcriptionally competent R-bands compared with G-bands and hypomethylation was most prominent on the chromosomes 17 and 19, which contain the highest gene density in the genome. The observed intra-chromosomal clustering of DNA methylation changes may be a specific biologic response to PM_2.5. However, it may also demonstrate that changes in methylation is enriched at loci of high gene density or simply is a matter of technical bias by annotation of variable CpG sites to these regions.

Gene ontology enrichment analysis revealed that hypomethylated genes prominently affected pathways involved in G-protein signaling, cellular structure and extracellular interaction. Correspondingly, we observed PM-dependent cytoskeletal changes, such as the formation of stress fibers,²⁰ which might progress to epithelial barrier disruption in the respiratory tract, systemic spread of inhaled particles and lung-independent diseases. Hypermethylated genes were also associated with pathways involving cellular adhesion, indicating that these pathways are a dominant field of epigenetic regulation by PM. Loss of cellular adhesion by repressive hypermethylation might promote cancer cell progression, permeability of epithelial barriers, and metastasis. Furthermore, hypermethylation also clustered for pathways responsible for ion binding and shuttling, which might represent a biologic response to the high salt content in PM_2.5.

Hypomethylation was evident at repetitive element-associated DNA. Likewise, exposure to ambient air pollution or cigarette smoke induces hypomethylation of interspersed elements (LINE-1) or tandem repeats.^18,31,34 Reduced methylation of interspersed repeats is associated with genomic instability, overexpression of oncogenes, and a higher risk of developing and dying from cancer.^35,36 It has been proposed that hypomethylation of repetitive elements may also be used as a marker for many other human diseases, including metabolic disorders, lung injury, cardiovascular diseases, and associated risk factors.^16,37-39 Thus, genome-wide hypomethylation of repetitive elements might represent a major mechanism to promote PM-associated diseases.

miRs posttranscriptionally affect gene expression and their deregulation may be involved in the pathogenesis of human disease.⁴⁰ Transcriptional activation of miRs can be regulated epigenetically⁴¹ and exposure to environmental pollutants or cigarette smoke may result in differential expression of miRs.^42,43 We found methylation changes in some of these miR loci but also identified others for which relationship to PM exposure has not yet been reported (miRs listed in Table 2). Methylation changes were located near to the transcription start point suggesting functional relevance. A miR-specific disease database predicted a connection to PM-associated diseases. Almost all hits were related to the development or progression of cancer. miR-loci with a tumor suppressive allocation were frequently hypermethylated,^44-48 whereas some oncogenic miR-loci were hypomethylated.^49-55 However, the biologic functions of miRs are complex and most often tumor-type specific. Several miRs were also involved in the regulation of inflammation, response to pathogens, cardiovascular diseases, and associated risk factors, or various lung specific diseases (see Table 2). Some of these miRs have been described as useful biomarkers, such as for early detection of non-small cell lung cancer (NSCLC),⁵⁶ COPD,⁵⁷ or type 2 diabetes,⁵⁸ indicating that their dysregulation might represent an early event in disease development. Furthermore, dysregulation of some of these miRs frequently correlates with outcome of cancer and is predictive for refractory therapy,^59-61 development or relapse of disease before clinical manifestation,⁶² and poor prognosis.^50,51,63 Additionally, it may lead to hypertension⁶⁴ and is also indicative of increased risk of death upon chronic heart disease⁶⁵ or post infarction.⁶⁶ However, further experiments are needed to prove that methylation changes affect miR and target gene expression, and thus represent a prominent mechanism for PM-associated diseases.

Hypomethylation of genes is associated with transcriptional activation.¹⁴ Accordingly, most genes that were hypomethylated also showed increased gene expression. In contrast, only less than half of the hypermethylated genes showed transcriptional repression. Because hypermethylation of CpG islands within promoter regions has a high impact on transcriptional repession,¹⁴ these genes were separately analyzed, but also showed no statistical bias for gene repression (data not shown). Controversially, increased gene body methylation or high methylation of shore and shelf regions, flanking demethylated CpG islands are associated with high gene expression, which might partly explain lack of correlation.^67,68 Nevertheless, methylated CpG dinucleotides are considered to be mutational hot spots that cause genetic diseases and cancer through spontaneous deamination of 5-methylcytosine (5mC) and its transition to thymine.¹⁴

We functionally clustered differentially methylated genes according to their up- or downregulated expression. Interestingly, most of the selected genes were repeatedly identified by different GO terms with potential impact on PM-associated diseases (genes listed in Tables S8/S9 and referenced in Table S10). Some antiviral genes were downregulated, whereas other genes, encoding viral docking sites, were upregulated. This may facilitate respiratory infections, frequently associated to ambient PM or cigarette smoke exposure.^69,70 Furthermore, genes that promote respiratory or cardiovascular diseases were primarily upregulated, whereas suppressors of these diseases were transcriptionally repressed. Encoded proteins are known to be critically involved in the pathogenesis of acute and chronic inflammation, asthma, COPD, atherosclerosis, or lung injury, but also often have pleiotropic functions and thus may additionally contribute to epithelial-mesenchymal transition, chronic inflammation and, consequently, to a pro-tumorigenic microenvironment.^71,72 Other allocated genes may directly promote tumorigenesis by facilitating proliferation, invasiveness, metastasis, angiogenesis, immune evasion or resistance and are associated with poor patient survival when upregulated in tumor cells, tumor stroma, or surrounding tissue (Tables S8/S9). Transcriptional and functional activation by cigarette smoke or airborne particulate matter and regulation by methylation have been independently confirmed for some of these genes (Tables S9/S10). In contrast, tumor suppressive genes were transcriptionally repressed in the presence of PM (Tables S8/S10). Interestingly, some upregulated genes contradictorily showed hypermethylation whereas some repressed genes were hypomethylated (Tables S6-S9). This indicates that these genes may not be exclusively epigenetically regulated, but that additional, unknown mechanisms might predominate for their expression. However, most of these genes were repeatedly identified by multiple gene ontology terms with highest significance highlighting their potential relevance. Because they represent central mediators of severe diseases, their reciprocal interaction and pathophysiological impact on adverse health effects by PM should be further examined.

In addition, the disease relevance of differentially methylated and expressed genes was confirmed by disease ontology. Some of the previously functionally clustered genes were allocated again and numerous additional genes were identified (Table 3). Genes, significant for cancer formation and progression were most highly enriched, indicating that exposure to PM may facilitate pro-tumorigenic conditions and transform lung epithelial cells to tumor cells. The cancer-associated molecular functions of these genes are referenced in Table S10. Again, tumor-promoting genes and genes associated with poor outcome in cancer patients were generally upregulated, whereas tumor-suppressive genes were often downregulated. However, the observed upregulation of some tumor-suppressive genes might represent an adaptive response to protect from cellular transformation. Independent studies have shown that many of these genes can be regulated by DNA methylation (see Table 3 and Table S10) and, interestingly, most of the allocated genes were associated with the same dominant pro-tumorigenic pathways, such as EGF/TGF-β-signaling, EMT including downregulation of E-cadherin, b-catenin activation, and response to hypoxia. Since these pathways are substantially linked to malignant progression, metastasis, angiogenesis, therapy resistance, and poor clinical outcome,^25-29 it is crucial to further verify their activation and pathologic function during PM exposure.

In addition to cancer, dysregulation of some of these disease relevant genes may also exacerbate pulmonary inflammation, lung injury, severity of respiratory infections, and progression of COPD and ARDS (Table S10), but also lung-independent diseases, such as cardiovascular disease, hypertension, metabolic syndrome, neurologic disorders and chronic inflammation, or cancer in different tissue types (Table S11) and, thus, may contribute to other common pathological features symptomatic for PM exposure. However, here we used bronchial epithelial cells and therefore our discussion focuses on lung diseases. Whether PM_2.5 induces similar biologic responses in different cell types needs further investigation.

Our results may open the door for preventive or therapeutic considerations, e.g., by using dietary compounds to reduce or reverse detrimental methylation changes when exposed to high PM_2.5 concentrations. It has already been shown that folic acid can enhance methylation and that DNA methylation inhibitors can significantly reduce carcinogen-induced lung cancer development in animal models.^73,74

However, the presented data may have some limitations by the use of a single fraction of PM from biomass combustion. Particle composition and cellular responses can vary depending on the type of biomass and combustion conditions.⁷⁵ Additionally, the concentration of PM and duration of exposure may have a direct impact on biologic responses. The given concentration (100 μg/ml), which had no significant impact on cell viability upon prolonged treatment, was also used in other studies investigating biocombustion.¹¹ Moreover, our approach for efficient identification of putative target sites for methylation changes is additionally supported by the fact that ambient PM concentrations from biomass combustion may significantly exceed industrial or traffic-related emissions,^3,8 and that inhaled PM <2.5 μm may accumulate in the lung parenchyma.⁹ Consistently, health effects are especially evident in regions with recurrent and intense pollution, but rarely after a single and short exposure.¹ Nevertheless, further studies will be necessary for targeted validation of the most important top hits. We also acknowledge that the use of an immortalized cell line and cell culture conditions may not reflect in vivo processes. Moreover, both 5mC and 5-hydroxymethylcytosine (5hmC) are interrogated by the Illumina 450K array. The 5hmC modification, which is an intermediate of 5mC demethylation by TET enzymes,⁷⁶ may be responsible for some of the activation, correlated with higher promoter methylation. Finally, the Illumina 450K array does not completely cover the entire epigenome.

In conclusion, using BEAS-2B cells as a model, we could observe genome wide DNA methylation changes at CpG sites residing in disease relevant genes, repetitive elements, and miRs. Consequently, DNA methylation analysis should be included in future efforts to understand the interactions between environmental exposures and PM-associated diseases. By identifying putative targets and their allocation to selective prominent disease relevant pathways, there is now an urgent need to independently confirm the biologic effects of methylation changes on gene expression, to clarify the interplay of the affected mediators, and to demonstrate how they contribute to adverse health effects. This might finally result in the establishment of new sensitive and reliable diagnostic biomarkers and therapeutic targets for the prevention and/or treatment of PM-associated diseases.

Materials and methods

Preparation of PM_2.5

Preparation of PM_2.5 has been described previously.²⁰ Briefly, bulk fly ash was collected from an electrostatic precipitator of a medium-scale biomass power plant with a nominal thermal output of 1.7 MW (Bürger Energie St. Peter eG, St. Peter, Schwarzwald, Germany), which exclusively combusts chips of soft wood (mainly spruce) derived from the local forests. The wood chips included debarked stem wood and branches with a minimum of 7 cm in diameter (merchantable wood) without leaves and twigs. The collected ash was subsequently size-fractionated by a cyclone with an aerodynamic cut-off diameter of 2.5 µm (Labor für Partikeltechnologie/Mechanische Verfahrenstechnik, Hochschule Konstanz, Technik, Wirtschaft und Gestaltung). The resultant size fraction (PM_2.5) was then used for all biologic assays.

Cell culture and PM_2.5 treatment

Immortalized bronchial epithelial cells (BEAS-2B; LGC Standards GmbH, #ATCC-CRL-9609) were maintained in Dulbecco's Modified Eagle Medium/F12 with stable Glutamine (GE Healthcare, #E15–889) containing 5% fetal bovine serum (Life Technologies, #10270–106) 100 IU streptomycin and 100 IU penicillin (Sigma-Aldrich, #11074440001) at 37°C in a humidified incubator. BEAS-2B cells were plated at a density of 0.7 × 10⁶ cells per 75 cm² cell culture flasks (Greiner bio-one, #658175) and were split twice a week after rinsing 2 times with phosphate buffered saline and subsequently using trypsin, 0.25% EDTA (Life Technologies, #25200056) for detachment. For the experiments, adherent cells were repeatedly treated with PM_2.5 every first day after each passage. Particles were immersed in BEAS-2B growth medium at a concentration of 200 μg/ml and sonicated for 20 min (SONOREX RK510H, BANDELIN electronic GmbH & Co. KG, Berlin, Germany) immediately before dilution to the cells at a final concentration of 100 µg/ml PM_2.5. The cells were maintained in the same PM_2.5-containing medium until the next passage. For the last passage, cells were seeded in 6-well plates and were exposed on the next day with 100 µg/ml PM_2.5 for a final time period of 48 h. Control cells were similarly treated in the absence of PM_2.5. After 5 weeks of exposure, cells were harvested for DNA (DNeasy Blood & Tissue Kit, Qiagen, #69504) and RNA (RNeasy® Plus Mini Kit, Qiagen, #74136) extraction following the manufacturer's instructions. To avoid RNA contamination, 400 μg RNase A (Carl Roth, #7156.1) were added during cellular lysis for genomic DNA isolation. For each condition (control and PM_2.5-treated cells) 3 independent replicates have been performed for DNA and for RNA isolation.

DNA methylation chip preparation and data analysis

The Illumina HumanMethylation450 BeadChip (450K) array was performed by AROS Applied Biotechnology A/S (Aarhus, Denmark) resulting in IDAT files containing the raw data. Analysis of Illumina 450K v1.2 data was performed with the R minfi package.⁷⁷ For the preprocessing of the raw data, the preprocessQuantile function was used. After calculation of β-values, 56,235 cross hybridizing probes were removed.⁷⁸ Differentially methylated positions were identified using the dmpFinder function and a False Discovery Rate (FDR) of less than 5% was obtained (q-value cutoff of 0.05). The differences of averaged β-values of the 4 replicates of each control and particular matter sample were defined with a cutoff for minimum difference of 0.02 (decrease in Δβ > 0.02, hypomethylated; increase in Δβ>0.02, hypermethylated). Probes were annotated using the getAnnotation function implemented in the minfi package. Enrichment calculations (Table S2) were based on the χ²-test (p < 0.05).

Gene expression and chip preparation

Affymetrix HG-U133 Plus 2.0 oligonucleotide microarrays (Affymetrix, #900467) were used to determine transcription levels of more than 54,000 transcripts and variants, representing over 20,000 genes. RNA was quality assessed (Agilent RNA 6,000 Nanokit, Agilent Technologies, #AGI 5067–1511; Agilent 2100 Bioanalyzer, Agilent Technologies, Waldbronn, Germany) and quantified (NanoDrop ND-1,000 spectrophotometer, Thermo Fisher Scientific, Berlin, Germany). The RNA quality revealed values of RIN ≥ 9.9 for all samples. After reverse transcription of 300 ng total RNA (Invitrogen™ Ambion™ RETROscript™ Reverse Transcription Kit, Thermo Fisher Scientific, #AM1710) and generation of cRNA probes (3′ IVT Plus Amplification Kit, Affymetrix, #902395), 50 μg/mL of biotinylated and fragmented cRNA was hybridized to HG-U133 Plus 2.0 arrays. After washing and fluorescence labeling (GeneChip Hybridization, Wash & Stain Kit, Affymetrix, #900720), arrays were scanned (Affymetrix GeneChip Scanner 3000, Affymetrix).

Gene expression analysis

Analysis of the raw expression data was performed using the simpleaffy⁷⁹ and the affy R packages⁸⁰ for uploading and RMA normalization of the data. The lmFit, contrasts.fit, as well as the eBayes functions of the limma package⁸¹ were used to preprocess the data, and, finally, the topTable function was applied to identify differentially expressed genes, with a false discovery rate defined by the Benjamini-Hochberg method of smaller than 0.05 and each p-value smaller than 0.05. Annotation data of probes were obtained from the supplier's homepage (http:/www.affymetrix.com/support/technical/byproduct.affx?product = hg-u133-plus, file: “HG-u133_Plus_2.na35.annot.csv,” creation_date = 2014–10–06). If several probe sets were assigned to the same LocusID either all sets were kept individually since variances between those values might be a result of alternative splicing of the related transcript (Fig. 4)⁸² or the average of the transcripts were calculated (Tables 3, S6, and S7) for a clearer presentation.

Data plotting

Mean β values of the control probes were plotted against the β values of the PM-treated samples of all differentially methylated CpG sites using the ggplot2 package of R.⁸³

Transcription start site (TSS) positions were obtained from the UCSC website (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/; file: “knownGene.txt.gz”) taking into account the strandedness of annotated genes. Distances of all differentially methylated CpG positions to their nearest TSSs were calculated with the help of the DistanceToNearest-function of the “GenomicRanges” package,⁸⁴ and plotted against their corresponding Dbvalue.

The phenogram shown in the manuscript is a collage of 3 different phenograms (for hypomethylated CpGs, hypermethylated CpGs, and Giemsa bands). Each phenogram was generated using the web application PhenoGram.⁸⁵ The input data table for the Giemsa-stained cytobands was obtained from the PhenoGram website (cytoBand-human.txt, http://visualization.ritchielab.psu.edu/phenograms/document).

The names of all differentially methylated or differentially expressed genes were converted with the help of the web application GeneCards (http://www.genecards.org/) and the db2db database (https://biodbnet-abcc.ncifcrf.gov/db/db2db.php), and subsequently analyzed for enrichment using the default parameters of the Enrichr platform.²¹ Plotting of the tabular output of Enrichr was performed using Python's matplotlib module (http://matplotlib.org/).

Overlaps of all differentially methylated CpG positions with different kinds of repetitive elements (“LINE,” “SINE,” “LTR,” “Simple_repeat,” “Low_complexity,” “Satellite”) were counted using the GenomicRanges package in R.⁸⁴ The repetitive element positions were downloaded from the UCSC website (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/, file: rmsk.txt.gz). Overlaps with miRs were counted and incorporated with the input data from the miRBase website (http://www.mirbase.org/ftp.shtml, file: hsa.gff3; gff-version 3 - miRBase v21, build GRCh38, hg38). The CpG positions were converted from hg19 to hg38 by using the liftOver tool (http://genome.ucsc.edu/cgi-bin/hgLiftOver).

Circos plots were generated using the default settings of related web application (http://mkweb.bcgsc.ca/tableviewer/).

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Authors contribution

Katharina Heßelbach provided biologic samples, was involved in data processing, their interpretation and assisted with manuscript preparation. Gwang-Jin Kim and Stefan Günther were responsible for bioinformatics data processing, statistical analysis and manuscript preparation. Stephan Flemming established the bioinformatics data processing and was involved in the interpretation of results. Thomas Häupl and Marc Bonin conducted the Affymetrix HG-U133 Plus 2.0 transcriptome analysis. Regina Dornhof provided biologic samples. Irmgard Merfort and Matjaz Humar were responsible for study design, sample preparation, data analysis, interpretation and presentation of results and manuscript preparation.

Funding

This work was supported by the EU Interreg IV Programm “Oberrhein” (project C35 BIOCOMBUST) and the German Research Foundation (DFG, CRC992, Medical Epigenetics).