Introduction

Asthma is one of the most common chronic disorders in childhood1, currently affecting an estimated 6.2 million children under 18 years, of which 3.1 million suffered from an asthma attack or episode in 20152. The annual cost of asthma is estimated to be $50.1 billion, representing $3,100 per person per year, more than half of which is related to prescription medication costs3,4,5. Although multiple asthma therapies exist, the heterogeneous nature of asthma and great variability among patients’ individual therapeutic responses highlight the need for alternative individualized therapies based on a better understanding of disease mechanisms6,7. A growing body of epidemiological studies, including ours8,9,10, have identified epigenetic variations associated with asthma and asthma severity in various tissues and suggested an important role for DNA methylation (DNAm) in asthma pathogenesis8,9,10,11,12,13. Studies in murine models demonstrated the role of DNA methylation in T cell polarization and cytokine production14,15, and that pharmacologic inhibition of DNA methylation in CD4+ T cells protected against the development of allergic airway inflammation15. Studies of lung tissues from a mouse model of experimental asthma identified DNA methylation changes related to smooth muscle functions in mice with allergic asthma16,17. However, further mechanistic studies are needed to understand the pathways through which DNA methylation contribute to asthma features, especially in the regulation of the airway epithelium, a major player in asthma due to its unique interface with the environment and interaction with immune cells18,19,20.

Our previous publications have associated TET1 promoter methylation with childhood asthma8,9. Tet1 is known to regulate gene expression in many cell types including embryonic stem cells (ESC) and HEK293 cells21,22,23,24,25. Tet1 belongs to the mammalian demethylase family (TET1/2/3) that hydrolyzes 5-methyl-cytosine (5mC) to generate 5hmC, usually resulting in activation of gene expression26,27,28,29,30. In addition to its catalytic activity, the TET1 protein recruits histone modifying protein complexes (e.g. OGT/SET1, SIN3A/HDAC and EZH2/SUZ12/EED) to alter histone marks and chromatin accessibility, leading to both activation and repression of gene expression31,32,33,34. TET1 is involved in many biological processes and diseases, including stem cell maintenance, T and B cell development, genomic imprinting, neural activity and cancer35,36,37,38,39,40,41,42,43. The expression of TET1 is regulated by asthma-related exposures such as diesel exhaust particles and house dust mite (HDM) in human bronchial epithelial cells8,9 and in whole lungs of mice16. However, the contributions of TET1 to asthma remain unknown.

In this paper, we examined the role of Tet1 in allergic airway inflammation using animal models and cell lines (Fig. 1, study design). We applied integrative transcriptomic analysis, pathway and protein interaction network analysis to identify genes and pathways regulated by Tet1 in HDM-exposed murine lungs and in human bronchial epithelial cells. We then examined the specific regulation of candidate genes in identified pathways by TET1 in bronchial epithelial cells at base line and following acute HDM challenges. Finally, we performed DNA methylation studies and functional genomic analyses to understand how Tet1 regulates gene expression in the airways. Collectively, our data from both mouse models and human bronchial epithelial cells strongly support that Tet1 suppresses allergic airway inflammation by transcriptional regulation of Interferon (IFN) and Aryl hydrocarbon Receptor (AhR) signaling pathways through changes in DNA methylation, and interactions with particular transcription factors and histone modifiers.

Figure 1
figure 1

Summary of study design.

Results

Loss of Tet1 exacerbates HDM-induced AHR and lung eosinophilia

To determine the role of Tet1 in asthma development, we utilized an experimental model of allergic airway inflammation that was previously established17,44 (Fig. 2a). Loss of Tet1 significantly exacerbated HDM-induced AHR at methacholine doses of 6, 12 and 25 mg/ml (Fig. 2b). No significant increase in AHR was observed in the HDM-treated Tet1+/− mice compared to HDM-treated Tet1+/+ mice (Fig. 2b). There were significantly more Bronchoalveolar lavage fluid (BALF) cells in both the HDM-challenged Tet1+/− and Tet1−/− mice (Fig. 2c). Significantly more eosinophils were observed in the BALF of Tet1−/− and Tet1+/− mice compared to Tet1+/+ mice (Fig. 2d), which may be driven by increased total cell number in BALF. In summary, these data suggest that loss of Tet1 increases HDM-induced AHR and lung inflammation. Sensitization was not significantly affected as exposure to HDM did not significantly alter total IgE, total IgG1, HDM-specific IgE and HDM-specific IgG1 levels in serum from Tet1−/− and Tet1+/− mice when compared to Tet1+/+ mice (Fig. 2e–h).

Figure 2
figure 2

Loss of Tet1 exacerbated airway hyperresponsiveness and lung inflammation in a mouse model of asthma. (a) Treatment protocol. (b) Airway hyperresponsiveness. (c) Total BALF cells. (d) Number of each cell type. (e) Total IgE. (f) Total IgG1. (g) HDM-specific IgE (1:5 dilution). (h) HDM-specific IgG1 (1:2000 dilution). Tet1+/+ Saline, n = 8–15; Tet1+/+ HDM, n = 7–24; Tet1+/− Saline, n = 9–14; Tet1+/− HDM, n = 9–24; Tet1−/− Saline, n = 6–7; Tet1−/− HDM, n = 9–16. mean ± SEM for each group is shown. Unpaired student t-test was applied with Bonferroni corrections for multiple testing. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Loss of Tet1 promotes the expression of pro-Th2 and Th17 cytokines and epithelial repair genes induced by HDM

Next, we assessed the impact of Tet1 genotype on genes involved in pulmonary HDM induced immune responses, including IL-13, a central mediator of allergic asthma45. The expression levels of Il13 in lung tissues were significantly elevated in Tet1−/− mice compared to Tet1+/+ mice treated by HDM (Fig. 3a). HDM-induced Il4 expression was similar between genotypes (Fig. 3a), consistent with the observed lack of differences in HDM-specific IgE levels (Fig. 2g). Among other pro-Th2 cytokines, HDM-treated Tet1−/− mice had significantly higher Il33 expression, and no significant changes in Tslp. In addition, we observed a significant increase in Il17a expression, but not Il17f. Expression of Il6 and Il1b, two cytokines involved in Th17 differentiation, were not changed. Since genes involved in epithelial repair are upstream of the activation of the NF-κB pathway and they promote the generation of pro-Th2 and Th17 cytokines46,47, we next examined the expression of Egfr, Tgfa, and Tff2 and found that the expression of Egfr and Tff2 were both significantly increased HDM-treated Tet1+/− mice (Fig. 3a). Consistent with changes in mRNA levels, we observed increased IL13 protein levels in BALF and significantly increased IL33 protein levels in lung homogenates of HDM-exposed Tet1−/− mice (Fig. 3b). IL5 protein levels in BALF remained unchanged. Collectively, these data suggest that loss of Tet1 promotes the expression of pro-Th2 and Th17 cytokines and genes involved in epithelial damage and repair in the lungs.

Figure 3
figure 3

Tet1 deficiency increased the expression of pro-Th2/Th17 and epithelial repair response genes. (a) Expression levels of indicated genes were measured by RT-qPCR 48 hrs after the last i.t. HDM challenges. Expression values were normalized to the expression of mRpl13a. Mean ± SEM for each group is shown. Tet1+/+ Saline, n = 4–10; Tet1+/+ HDM, n = 11–22; Tet1+/− Saline, n = 4–14; Tet1+/− HDM, n = 12–21; Tet1−/− Saline, n = 4–7; Tet1−/− HDM, n = 10–19. Student t-tests were applied to analyze normally distributed data for Il4, Il1b, Tslp, Il17f and Egfr with Bonferroni corrections for multiple testing. For comparisons involving data that are not normally distributed in all other genes, Mann Whitney tests were applied. (b) Protein levels of IL13, IL5 and IL33 measured by ELISA. IL13 and IL5 were measured in BALF, and IL33 was measured in lung homogenate. Tet1+/+ Saline, n = 10–11; Tet1+/+ HDM, n = 8–10; Tet1+/− Saline, n = 10; Tet1+/− HDM, n = 10; Tet1−/− Saline, n = 7; Tet1−/− HDM, n = 10–11. Mean ± SEM for each group is shown. Data are normally distributed and unpaired student t-tests were applied with Bonferroni corrections for multiple testing. *p < 0.05, **p < 0.01. ns represents not significant.

Loss of Tet1 increases the expression of Muc5ac

The strong Th2 response observed after HDM exposure induced significant mucus production; however, no significant differences in PAS stained large airways was observed between HDM challenged Tet1+/+, Tet1+/− and Tet1−/− mice (representative images shown in Fig. 4a). To determine if Tet1 deficiency was associated with minor alterations in mucin composition, we assessed lung expression of Muc5ac, which is involved in allergen-induced AHR48 and increased in asthmatics49,50, as well as Muc5b, whose expression is decreased in asthmatics49,50 and functions to maintain normal mucus clearance51. The expression of Muc5ac was significantly increased in HDM-challenged Tet1+/− mice and Tet1−/− mice compared to Tet1+/+ mice (Fig. 4b), consistent with IL-13 BALF levels (Fig. 3b). No significant difference was observed for Muc5b expression (Fig. 4c). Taken together with increased lung IL-13 levels in HDM-treated Tet1−/− mice, these data suggest that IL-13 mediated Muc5ac expression may contribute to the increased AHR in the HDM-treated Tet1−/− mice.

Figure 4
figure 4

Tet1 deficiency increased the expression of Muc5ac in the lungs. (a) Periodic acid–Schiff (PAS) staining of lung sections. Representative images (4 slides per animal, 2–3 animals per group) are shown. (b) Expression of Muc5ac (b) and Muc5b (c) were measured by RT-qPCR. Expression values were normalized to Rpl13a. Mean ± SEM for each group is shown. Tet1+/+ Saline, n = 8–10; Tet1+/+ HDM, n = 14–16; Tet1+/− Saline, n = 12–14; Tet1+/− HDM, n = 16–21; Tet1−/− Saline, n = 6–7; Tet1−/− HDM, n = 15–17. Unpaired student t-tests (Muc5ac) or Mann Whitney tests (Muc5b) were applied. *p < 0.05. ns represents not significant. Bar represents 100 μm.

Lung transcriptomic analysis reveals interferon signaling as a major network upregulated in Tet1−/− mice challenged by HDM

To gain mechanistic insight into the role of Tet1 in allergic airway responses, we performed an exploratory RNA-seq experiments on total lung RNA samples from saline- and HDM-challenged Tet1+/+ and Tet1−/− mice (two animals in each group). Compared to Tet1+/+ mice treated with saline, saline-challenged Tet1−/− mice showed differential expression in 104 genes, among which 22 genes were downregulated and 82 genes were upregulated (Supplementary Table S1). Pathways affected by the genes with elevated expression include Calcium signaling, Agranulocyte Adhesion and Diapedesis, Actin Cytoskeleton Signaling, TR/RXR Activation and Role of IL-17A in Psoriasis (Supplementary Table S2). In the HDM-treated mice, 72 genes showed reduced expression and 61 genes showed increased expression in Tet1−/− compared to Tet1+/+ mice (Supplementary Table S3). Nineteen genes out of 133 genes overlap with genes dysregulated in the saline-treated mice, suggesting that these 19 genes are responsive to Tet1 but may not be specific to Tet1 deficiency in the context of HDM challenge. As our focus is to understand the impact of Tet1 on allergen-induced genes, these 19 genes were excluded from further pathway and gene network analysis. Among the downregulated genes in HDM-treated Tet1−/− mice, the most significantly enriched pathways include LPS/IL-1 Mediated Inhibition of RXR Function, Nicotine Degradation II, Aryl Hydrocarbon Receptor Signaling, Xenobiotic Metabolism Signaling and Glutathione-mediated Detoxification (Fig. 5a and Supplementary Table S4A), all of which share genes such as Cyp1a1, Gsta3/Gstm3/Gsto1, Fmo2/Fmo3 and Aldh1a1. Among the upregulated genes in HDM-treated Tet1−/− mice, top significantly enriched pathways include Activation of IRF by Cytosolic Pattern Recognition Receptors, Interferon Signaling, Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses and Role of RIG1-like Receptors in Antiviral Innate Immunity (Fig. 5b and Supplementary Table S4A), all of which include components of the interferon (IFN) signaling pathway such as upstream regulator Irf7 and downstream interferon-stimulated genes (ISGs) such as Isg15, Ifit3 and Oas1/Oas2/Oas3. This was further supported by the upstream regulator analysis, where many components of IFN signaling, including type I/II/III interferon signaling, were identified as activated (Supplementary Table S4B). Accordingly, functions related to infectious diseases (viral infection and replication) were predicted to be repressed whereas functions related to inflammatory responses including T cell, B cell, macrophage and inflammatory diseases (rheumatic disease and multiple sclerosis) were predicted to be activated (Table 1 and Supplementary Table S4C).

Figure 5
figure 5

Pathway and protein network analysis of genes differentially expressed between HDM-challenged Tet1+/+ mice and Tet1−/− mice. (a) Top 11 significantly enriched canonical pathways among the genes that were downregulated in HDM-treated Tet1−/− mice. (b) Significantly enriched canonical pathways among the genes that were upregulated in HDM-treated Tet1−/− mice. (c) Minimal network containing only direct protein-protein interactions between seed proteins from the input dataset. 133 genes differentially expressed genes were used as input, and the links between them indicate known protein-protein interactions. Green color marks genes with downregulated expression, red color marks genes with upregulated expression, grey color marks genes that are not included in the input dataset. (d) Significantly enriched pathways among genes included in c. (e) RT-qPCR validation of selected differentially expressed genes from RNA-seq. Tet1+/+ Saline, n = 6–7; Tet1+/+ HDM, n = 6–7; Tet1−/− Saline, n = 5–6; Tet1−/− HDM, n = 7. Expression levels were measured by RT-qPCR and were normalized to the expression of mRpl13a. Mean ± SEM for each group is shown. Data are normally distributed and unpaired student t-tests were applied with Bonferroni corrections for multiple testing. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. ns represents not significant.

Table 1 IPA predicted direction of change in functions by genes differentially expressed comparing HDM-challenged Tet1−/− mice with Tet1+/+ mice.

We further performed network analysis to identify significantly enriched protein-protein interactions among the genes differentially expressed between Tet1+/+ and Tet1−/− mice challenged by HDM. Among proteins encoded by differentially expressed genes (DEGs) in our dataset (seed nodes) and proteins directly interacting with them, there are several hub nodes centered around Irf7, Foxo3, Lepr, Pml, Cola1, and Fos. After trimming the network to only keep nodes that are necessary to connect the seed nodes (input genes), a smaller dense network including Irf7 and multiple ISGs (such as Isg15, Mx2, Ifit1, Oas2/3) was identified (Fig. 5c). Functional analyses (KEGG, Reactome and Molecular function) revealed that this sub-network is enriched for genes in Toll-like receptor signaling and IFN signaling (Fig. 5d and Supplementary Table S5). The expression differences in Irf7, Isg15 and Cyp1a1 were further validated by RT-qPCR (Fig. 5e). Collectively, these analyses identify IFN signaling as a major network that is significantly upregulated by loss of Tet1 following HDM challenges.

TET1 regulates IFN and AhR signaling pathways in bronchial epithelial cells

As an important player in asthma pathogenesis18,52,53, airway epithelial cells are known to activate interferon and AhR signaling pathways following HDM challenges54,55,56,57,58,59,60,61,62. Because our transcriptomic analysis of lung RNA identified Irf7 as a major component of IFN signaling regulated by Tet1 in mice challenged by HDM, we next determined whether such regulation of Irf7 by Tet1 occurs in bronchial epithelial cells treated by HDM. As shown in Fig. 6a,b, partial knockdown of TET1 (~60%) in human bronchial epithelial cells (HBECs)63 significantly increased IRF7 expression compared to siRNA treated controls at base line and following 24 hr HDM challenge (25 ug/ml). Similar increases were observed in the expression of IFNα isoforms54 (Fig. 6c,d) and IFNβ (Fig. 6e). IFNλ1/2 expression was not detectable.

Figure 6
figure 6

Knockdown of TET1 in human bronchial epithelial cells regulates the expression of genes in IFN and AhR pathways. Expression levels of TET1, IRF7, IFNα1, IFNα2, IFNβ, CYP1A1, and ALDH1A1 in HBECs with or without TET1 knockdown and treatment of HDM were measured at indicated time points. Expression values of six technical replicates were normalized to the expression of GADPH. Mean ± SEM for each group is shown. Data are normally distributed and student t-tests were applied with Bonferroni corrections for multiple testing. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. ns represents not significant.

In addition, we also observed downregulation of genes involved in AhR signaling, including Cyp1a1 and Aldh1a1 in lungs from HDM-treated Tet1−/− mice comparing to Tet1+/+ mice. Since the expression of CYP1A1 and ALDHA1 (Phase I enzymes) represents an early response following activation of aryl hydrocarbon receptor signaling64,65, we examined the regulation of these genes by TET1 in HBECs at both 1 hr and 24 hr time points. Downregulation of TET1 in HBECs significantly reduced the expression of CYP1A1 in the presence of HDM at 1 hr (Fig. 6f,h). For ALDHA1, we observed significant reduction in expression at 1 hr without HDM exposure and a trend of decrease with HDM exposure (Fig. 6h). Interestingly, both genes showed a significant increase following 24 hr exposure to HDM when TET1 was knocked down (Fig. 6i,j), underscoring the temporally dynamic changes of these genes in response to HDM. In summary, our data showed that TET1 transcriptionally regulates the IFN and AhR signaling pathways in human bronchial epithelial cells following HDM challenges, consistent with the findings from the mouse model.

Next, as we observed relatively small but significant changes induced by Tet1 knockdown in HBECs without HDM at 24 hr time point (Fig. 6b,c,i,j), we examined whether activation of TET1 regulates the IFN signaling and AhR signaling pathways in HBECs without HDM challenges. Using a dCas9-SAM system to epigenetically activate the endogenous locus of TET1, we were able to increase TET1 expression by ~14 fold (Supplementary Fig. S1a). The activation of TET1 expression led to significant downregulation of IRF7 (Supplementary Fig. S1b) and other genes in IFN signaling pathway (Supplementary Fig. S1c, Supplementary Table S6 and S7). Genes involved in the AhR pathway were also regulated, including CYP1A1 (Supplementary Table S6 and S7). A protein-protein interaction network including IRF7 and ISGs was identified by network analysis (Supplementary Fig. S1d, highlighted by blue circles). These data provide further support for the transcriptional regulation of IFN and AhR signaling pathways by TET1 in bronchial epithelial cells.

Tet1 regulates IFN and AhR signaling in airway epithelial cells possibly through 5mC/5hmC, transcriptional factor binding and histone modification

Tet1 regulates gene expression by altering 5mC/5hmC at enhancers and promoters21,22,23,24, or binding to transcriptional co-factors such as Suz12 or the SIN3A complex to change histone modifications31,32,33,34 (Fig. 7a). We therefore performed genome-wide methylation analysis comparing lung DNA from HDM-treated Tet1−/− and Tet1+/+ mice (two animals in each group). Our genome-wide analysis identified 10,087 differentially methylated CpG sites, which are located near 2,738 genes (Supplementary Table S8). The majority of these sites (10, 047 out of 10,087 CpG sites) showed increased DNAm in HDM-treated Tet1−/− mice compared to Tet1+/+ mice, consistent with the established role of Tet1 in DNA demethylation. Although these pathways were not significant enriched (p > 0.05), we did find genes involved in Activation of IRF by Cytosolic Pattern Recognition Receptors/IFN signaling (Tlr9, Irf8, Ifna4 and Ifnar1, etc) and LPS/IL-1 Mediated inhibition of RXR function and Aryl hydrocarbon receptor signaling (Sod3, Nqo1, Gsto2, Aldhl1/Aldhl2/Aldha2/Aldh3a1, etc) that were differentially methylated. However, we did not observe significant expression changes in these genes between HDM-treated Tet1−/− and Tet1+/+ mice. In addition, Irf7, Dhx58 and Oas2 in IFN signaling pathway and Aldh1a1, Fos, Gsto1 in AhR signaling pathway showed significant changes in DNAm at the level of p < 0.05. As shown in Fig. 7b, 5 CpG sites spanning the Irf7 gene locus (chromosome coordinates in Supplementary Table S9) were found to be differentially methylated in RRBS and two of them were validated in additional samples (Fig. 7c). Collectively, our analysis support that the Tet1 regulates genes in IFN and AhR signaling pathways through DNAm, even though at certain loci such regulation may not correlate with detectable gene expression changes in mouse lungs.

Figure 7
figure 7

Possible mechanisms through which Tet1 regulates IFN and AhR signaling pathways. (a) Schematic showing two possible pathways that Tet1 utilizes to silence gene expression. (b) DNAm at 5 CpG sites across Irf7 gene (mm10) from RRBS data at p < 0.05. Open bar: Tet1+/+ HDM; Grey bar: Tet1−/− HDM. (c) Measurement of DNAm levels at CG4 and CG5 in b by bisulfite pyrosequencing. N = 10 for both groups and unpaired student t-tests were applied on normally distributed data. *p < 0.05, **p < 0.01. Open bar: Tet1+/+ HDM; Grey bar: Tet1−/− HDM. (d) Transcription factor binding motif enrichment in promoter regions of Tet1-regulated genes in human bronchial epithelial cells. The top 30 significantly enriched binding motifs are grouped by their corresponding TFs (p < 10−2). The size of each square or circle indicates fold-enrichment over background. (e) Overlap significance between promoter regions (2 kb upstream transcription start site) of Tet1-regulated genes in bronchial epithelial cells and ChIP-seq datasets. The top 30 significant results are shown regardless of cell type, along with results from lung related datasets passing a p < 10−6 significant threshold. Lung-related cell types are shown in multi-color, and non-lung cell types are shown in dark blue. (f) Overlap significance between gene promoters and histone marks. The significance of the degree of overlap between promoters and each member of a large library of histone mark ChIP-seq datasets was estimated. Histone marks with at least one significant result (p < 10−2) are shown. The Y-axis indicates the histone mark, in decreasing order of significance. The X-axis indicates the significance (−log P-value) of the overlap of the given dataset. The size of each circle indicates the fold-enrichment relative to background.

We next explored other potential mechanisms that may account for the regulation of gene expression by TET1. To this end, we applied two complementary approaches to identify regulatory molecules present in the promoters of genes with expression levels dependent on TET1. First, we performed TF binding motif enrichment analysis using the HOMER software package66 and motifs taken from the Cis-BP database67, revealing an enrichment for interferon signaling (IRF and STAT) and AhR signaling (ARNT) TF binding motifs (Fig. 7d and Supplementary Table S10A). We next applied a previously described computational method9 to search for enriched ChIP-seq peaks for transcription factors, co-factors and histone marks at the promoter regions of genes transcriptionally regulated by TET1 in human bronchial epithelial cells. Consistent with our pathway, protein network and motif enrichment analyses, the top enriched TFs include significant players in IFN signaling (Fig. 7e and Supplementary Table S10B, STAT1/2/3/4 and IRF1/2/3/5/7/8/9). Many of these enriched motifs and ChIP-seq datasets are present in the promoter of IRF7 (including IRF4/5, STAT1/3/4/5 and ARNT). Moreover, enriched binding was observed for EZH2 (Fig. 7e), a component of the PRC2 complex that generates H3K27me3 marks and is known to interact with TET1, was found. Consistent with the predicted binding of EZH2, H3K27me3 is by far the most significantly enriched histone marks around the promoter regions of TET1-regulated genes in HBECs (Fig. 7f and Supplementary Table S10C). For the IRF7 and CYP1A1 promoters, the binding of SIN3A, a known binding partner of TET1 to silence gene expression, was found (GSM803525). Histone modifications of H3K4, H3K9 and H3K27 that are indicative of promoter and enhancers were also found at the IRF7 promoter in lung-related cells and tissues. Collectively, our functional genomics analyses suggest that TET1 may regulate genes in the IFN signaling pathway through interactions with TFs and histone modifiers in bronchial epithelial cells.

Discussion

In this paper, we found that loss of Tet1 increased the severity of allergic airway disease, most notably AHR, in a mouse model exposed to the common allergen HDM. Accordingly, there was increased expression of Muc5ac, pro-Th2 and Th17 cytokine (e.g., Il13, Il33 and Il17a), increased BALF eosinophilia, as well as increased genes involved in epithelial repair responses in the lungs (e.g. Egfr and Tff2). Integrative transcriptomic analysis, pathway and network analyses revealed the upregulation of genes in the IFN signaling pathway and downregulation of genes in the AhR signaling pathway in Tet1-deficient lungs challenged by HDM. A gene network centered around IRF7 was identified as the major hub among genes affected by Tet1. Consistent with data from the mouse lungs, knockdown of TET1 in bronchial epithelial cells (HBECs), a major cell type that contributes to AHR and initiates lung inflammation, significantly upregulated expression of genes in the IFN signaling (IRF7 and type I interferons) and altered genes in the AhR pathway (CYP1A1 and ALDHA1) at base line and/or following HDM challenges. This transcriptional regulation of IFN and AhR pathways by Tet1 was further supported by transcriptomic analysis following activation of endogenous TET1 expression in HBECs. To investigate possible mechanisms that Tet1 may regulate Irf7 and other genes in vivo, we performed genome-wide methylation studies in mouse lung DNA and found changes in DNAm cross the Irf7 gene and many other genes in IFN and AhR pathways. In addition, promoters of genes regulated by TET1 in HBECs are enriched for predicted binding of TFs that interact with TET1 or function in IFN and AhR signaling, as well as the presence of histone marks indicative of promoters and other regulatory elements. Collectively, our data suggest that TET1 prevents allergic airway inflammation through inhibition of the IFN signaling and activation of the AhR pathway in airway epithelial cells.

TET1 is known to regulate gene expression in many cell types including embryonic stem cells (ESC) and HEK293 cells21,22,23,24,25. This regulation could be achieved through 5hmC because TET1 can directly bind DNA to convert 5mC to 5hmC and 5hmC is known to promote DNA demethylation at enhancers and gene bodies by altering the local chromatin structure68,69,70,71,72,73,74,75. “Readers” of 5hmC marks activate downstream gene expression networks21,22,23,24. Consistent with this mechanism, we observed small but significant changes in DNAm within the Irf7 gene in HDM-challenged-Tet1−/− mice (Fig. 7b,c). Similar changes were observed for several genes in the IFN and AhR pathways that showed gene expression changes between HDM-treated Tet1−/− and Tet1+/+ mice, although further validation studies are needed due to the limited number of samples included in RRBS. In addition, we also observed substantial DNAm changes (up to 60%) at other genes involved in IFN and AhR signaling whose expression changes were not detected following HDM challenges in the absence of Tet1. The lack of changes in gene expression may be due to the study of RNA expression in a mixture of lung cells instead of a single cell type, resulting in loss of signal. This could also explain why we did not observe expression changes in IFNα and IFNβ in total lung as in HBECs (differences in timing of measurement could be another reason). Alternatively, but not mutual exclusively, TET1 recruits histone modification enzymes (including OGT/SET1, SIN3A and EZH2 complexes) to regulate both active and repressive histone marks at promoters and enhancers in ESC and HEK293 cells and alter gene expression31,32,33,34. In support of this mechanism, we observed enriched binding of SIN3A and EZH2 within the promoters of Tet1-regulated genes in human bronchial epithelial cells and histone marks associated with active regulatory regions. Further studies investigating the spatial and temporal recruitment of TFs, histone marks, chromatin accessibility and 5mC/5hmC marks at genes regulated by Tet1 in airway epithelial cells, in combination with gene expression at the single-cell level, will significantly expand our understanding of Tet1’s gene regulatory roles in airway-related diseases.

In addition to Tet1, Tet2 and Tet3 may also contribute to the regulation of target genes due to their overlapping modulation of 5mC and 5hmC in other cell lines76. In addition, the DNA methyltransferases (DNMT1, 3A and 3B) are also enzymes that are important for maintenance and generation of 5mC marks. We found that Tet1 and Tet2 were significantly downregulated following HDM challenges in total lung RNA (a trend of downregulation was observed for Tet3), while Dnmt1, Dnmt3a and Dnmt3b remained unchanged (Supplementary Fig. S2). Therefore, the transcriptional responses to HDM in wildtype mouse lungs is likely mediated by TET, as opposed to DNMT proteins. This is different from a previous report using a relatively chronic HDM induced mouse model (1 i.p. sensitization by 100 ug HDM followed by 15 i.t. administration of 100 ug HDM over 5 weeks), where Tet1 was significantly upregulated and Dnmt3a was significantly downregulated in HDM-treated mouse lungs16. This discrepancy is probably due to the differences between mouse models and the dynamic changes in these enzymes involved in DNAm in response to HDM. Importantly, loss of Tet1 did not significantly influence the expression of Tet2 and Tet3 in mouse lungs at base line or following HDM challenges (p = 0.21 and 0.35), suggesting that the exacerbated allergic airway inflammation we observed in HDM-challenged Tet1−/− mice compared to Tet1+/+ mice were not mediated by Tet2 or Tet3. Interestingly, in the absence of Tet1, HDM treatment significantly reduced the expression of Dnmt1, Dnmt3a and Dnmt3b in addition to Tet2 and Tet3, suggesting complex interactions between Tet proteins and DNMTs in mediating the responses to HDM. Although loss of Tet1 significantly increased the expression of Dnmt1 and Dnmt3b in saline-treated mice, Dnmt1 and Dnmt3b were downregulated following HDM challenges (comparing HDM-challenged Tet1−/− mice with Tet1+/+ mice). Because inhibition of DNMT activity using 5-aza before sensitization in a mouse model of experimental asthma alleviated allergic asthmatic features including lung inflammation, mucus production and AHR77, the exacerbated asthma features we observed in HDM-challenged Tet1−/− mice compared to the HDM-challenged Tet1+/+ mice are unlikely due to the downregulation of Dnmt1, Dnmt3a, and Dnmt3b.

In the HDM-challenged mouse model, we observed an increase in the Th2 response in the absence of Tet1, while the IFN signaling pathway negatively regulated by Tet1 is mostly involved in viral responses and has been linked to asthma severity and exacerbation. As a master regulator of type I interferon-dependent immune responses78, IRF7 mediates airway epithelial responses to respiratory viral infection and was identified as a major hub gene connecting the IFN responses with virus-induced asthma exacerbations in vivo79,80. In addition, loss of IRF7 did not significantly alter HDM-induced allergic immune responses in mice81. However, it has been well established that respiratory viral infection leads to asthma exacerbation in human studies and in mouse models82. This is possibly due to the release of double-stranded DNA (dsDNA) triggered by viral infection that comes from NETosis (the formation of Neutrophil extracellular traps)83 in addition to HDM-induced dsDNA release due to tissue damage84. This free dsDNA binds to the TLR receptor and triggers signaling through IRF7 to promote downstream IFN responses55. We showed that loss of Tet1 transcriptionally mimics viral infection responses in mouse lungs and in human bronchial epithelial cells, which may promote the release of dsDNA and therefore contribute to an exacerbated IFN response. Therefore, we measured extracellular free dsDNA in BALF from our mouse models. Consistent with previous literature, there was significantly more extracellular dsDNA in BALF of mice challenged by HDM (Supplementary Fig. S3). Loss of Tet1 did not significantly change the amount of dsDNA in BALF. These data suggest that the enhanced asthma-related features in HDM-treated Tet1−/− mice are not due to increased dsDNA ligand binding to toll receptors triggering downstream IFN signaling, but rather due to the loss of transcriptional modulation of components of the IFN signaling pathway by Tet1.

IFN signaling, including type I and III IFNs that can be produced by airway epithelial cells54,55, is known to regulate innate and adaptive type I and II immune responses (dendritic cells/Th1/Th2/Treg) following viral/bacterial/allergen challenges85,86,87,88 and contribute to asthma development/exacerbation89,90. However, we observed no obvious effects of Tet1 deletion on T effector/memory cell differentiation (data not shown), suggesting that the enhanced Th2 responses in HDM-challenged Tet1−/− mice may not be due to the influences of IFN signaling on antigen presentation and T cell differentiation. In support of this, a very recent transcriptomic analysis of nasal epithelial brushing samples identified prominent activation of ISGs in adult asthmatics that is independent of viral infection, unrelated to type II inflammation, and associated with reduced lung function91. This is consistent with the presence of a T1-high group among severe asthmatics based on transcriptomic analysis92. Interestingly, Irf7−/− and Ifnar2−/− (deficient for type I IFN receptor) mice displayed normal dendritic cell development and allergic airway sensitization in response to HDM81, suggesting independence between type I IFN signaling and HDM-induced airway allergy when Tet1 is present. On the other hand, the AhR pathway, the second most significantly enriched pathway in our analysis, regulates airway epithelial multiciliogenesis93 and oxidative stress responses60, and contributes to asthma development and exacerbation46,62,94. Specifically, AhR ligands attenuate allergic airway inflammation in OVA-challenged animal models94,95, suggesting that down-regulation of the AhR pathway by Tet1 would increase the severity of HDM-induced allergic airway inflammation. Additionally, interactions between the IFN and AhR signaling pathways exist: increased IFN signaling by respiratory syncytial virus infection in airway epithelial cells induces protein degradation of transcription factor NRF2, which downregulates Nrf2-mediated expression of ARE-containing genes catalase and SOD196. As the CYP, ALDH and GST genes identified in our RNA-seq analysis are also targets of Nrf2, the downregulation of the AhR pathway might be the consequence of upregulated IFN signaling in Tet1−/− mice treated by HDM. Due to limited animal numbers included in the RNA-seq analysis, future studies will be performed to validate additional candidate genes in these pathways other than Irf7, Isg15 and Cyp1a1 (Fig. 5e). Whether the upregulation of IFN signaling (including upstream PRR signaling) and/or downregulation of AhR signaling pathways by Tet1 in airway epithelial cells directly leads to exacerbated HDM-induced allergic airway inflammation, and identification of the associated mechanisms, require further studies. Moreover, whether Tet1-mediated gene regulation would explain the interplay between viral infection and HDM-induced responses in airway epithelial cells and mouse lungs59,97 will be investigated in future studies.

In conclusion, we identified a novel role for Tet1 in asthma development. Our data suggest that the regulation of the IFN and AhR signaling pathways by Tet1 may contributes to the lung and airway epithelial responses to HDM challenges and the establishment of allergic airway inflammation. Similar to HDM, Tet1 expression is also responsive to other environmental exposures that contribute to asthma development and exacerbation, including diesel exhaust particles8,10 and cigarette smoke98. Therefore, our studies suggest that these exposures may contribute to asthma through the function of Tet1 on gene regulation. As TET1 can be modulated by small molecules (vitamin C, L-cysteine, 2-HG, etc) and miRNAs34, our findings may promote the development of new asthma therapies.

Methods

Murine model of allergic airway inflammation

Mice heterozygous for Tet1tm1.1Jae/J (B6/129S4) were purchased from Jackson Labs (Bar Harbor, ME) and littermates of wild type (Tet1+/+), heterozygous (Tet1+/−) and homozygous (Tet1−/−) for Tet1tm1.1Jae/J aged 8–12 weeks were utilized in the experiments. House dust mite (HDM) extract (Dermatophagoides pteronyssinus) was purchased from Greer Laboratories (Lenoir, NC). For an acute model of allergic airway inflammation, on days 0 and 7, mice were sensitized intraperitoneally with 100 μg of HDM (representing 33 µg of protein; 11 µg of Der p1; 5 EU of endotoxin) in 100 μL of PBS plus alum or an equivalent amount of PBS alone. On days 14, 19 and 21, mice were challenged intratracheally with 100 μg of HDM in 50 μL of PBS or PBS alone. On day 23, airway hyperresponsiveness (AHR) was analyzed, bronchoalveolar lavage fluid (BALF) were collected, blood was collected for measurement of serum proteins, and lung tissues were harvested for histology, immunohistochemistry and DNA/RNA extraction. All mice were maintained in a barrier facility at Cincinnati Children’s Hospital Medical Center and handled under Institutional Animal Care and Use Committee-approved procedures (protocol approved by CCHMC IACUC committee). All methods were performed in accordance with the relevant guidelines and regulations.

Measurement of airway hyperresponsiveness (AHR)

Invasive measurements of AHR were made on a flexiVent apparatus (Scireq, Montreal, Canada). Mice were anesthetized with Ketamaine, Xylazine, and Acepromazine (100, 20 and 10 mg/ml mixed at a ratio of 4:1:1). Mouse tracheas were cannulated with a 20-gauge blunt needle and the mice were ventilated at 150 breaths/min, 3.0 cm water positive end expiratory pressure. Two total lung capacity perturbations were then performed for airway recruitment before baseline measurement and subsequent methacholine challenges were performed. Dynamic resistance (R) and compliance (C) were determined by fitting the data to a single compartment model of airway mechanics where Ptr = RV + EV + Po, and Ptr = tracheal pressure, V = volume, E = elastance, Po is a constant and C = 1/E. Measurements were made using a 1.25 s, 2.5 Hz volume-driven oscillation applied to the airways by a computer-controlled piston (SnapShot perturbation). PBS or methacholine was aerosolized for 15 s (Aeroneb ultrasonic nebulizer) followed by 15 s of ventilation and a SnapShot perturbation was performed followed by 10 s of ventilation. Twelve SnapShot/ventilation cycle measurements were made. The procedure was repeated for 0, 6.25, 12.5, 25, 50, and 100 mg/ml concentrations of methacholine. The maximum R value and minimum C value with a coefficient of determination of 0.9 or greater (as determined by the flexiVent software) were used to determine the dose-response curve.

Bronchoalveolar lavage fluid (BALF) collection and analysis

Bronchoalveolar lavage was performed by cannulation of the trachea. The lungs were lavaged with 1 ml PBS +0.5% BSA +2 mM EDTA. The collected BALF was centrifuged and the total cell numbers counted with a hemacytometer. Cells were spun onto slides and stained with the HEMA3 stain kit (Fisher Scientific, Kalamazoo, MI). 400 cells were counted and the total number and % of each cell type were calculated.

Histology and immunohistochemistry

The left lobe of the lung was fixed in formalin, paraffin embedded and cut into 5 µm sections. Sections were stained with hematoxylin and eosin (H&E) or PAS according to the manufacturer’s recommendations (PolyScientific).

Lung homogenate

Frozen lung tissue was homogenised mechanically (OmniPrep Rotor Stator Generator, Omni International, USA) and enzymatically by addition of RIPA lysis buffer (MilliporeSigma, Burlington, MA) with protease inhibitors (MilliporeSigma) added. Samples were centrifuged and supernatants were saved. Total protein concentration was determined using BCA assay (ThermoFisher Scientific).

ELISA in serum, BALF and total lung homogenate

Total serum IgE and IgG1, HDM-specific IgE and IgG1 levels were measured as previously described99,100. For measurement of HDM-specific IgE and IgG1 levels, wells were coated with 0.01% HDM (Greer Laboratories) overnight. Serum was diluted 1:5 for IgE and 1:2000 for IgG1. After 2 hours of incubation, plates were washed, and either horseradish peroxidase–conjugated anti-mouse IgG1 (X56; 1:1000; BD Biosciences PharMingen, San Jose, Calif) or biotin–anti-mouse IgE (R35-118; 1:250; PharMingen) was added for 1 hour, followed by an incubation with streptavidin–horseradish peroxidase (R&D DY998; 1:200) in the case of IgE. BALF cytokine were also quantified as previously described99,100. IL13 protein levels were measured using Mouse IL-13 DuoSet ELISA kit (R&D Systems, 62.5–4,000 pg/mL) per the manufacturer’s instructions. IL5 was measured using Mouse IL-5 ELISA MAX™ Standard kit (Biolegend, 7.8 pg/mL to 500 pg/mL) per the manufacturer’s instructions. IL33 protein levels in lung homogenates were measured using Mouse IL-33 DuoSet ELISA kit (R&D Systems, 15.6 pg/mL–1000 pg/mL) per the manufacturer’s instructions. The results were then presented in relation to total protein concentration for each sample.

Activation and knockdown of TET1 expression in human bronchial epithelial cells

Human bronchial epithelial cells (HBECs) is a cell line obtained from Dr. John Minna’s lab at UT Southwestern Medical Center63. This cell line was created from bronchial specimens from areas of the lung histologically not involved with lung cancer by ectopic expression of Cdk4 and hTERT. It is not carcinogenic, has epithelial morphology and the expression of epithelial markers, and has an intact p53 checkpoint pathway. Previous cytogenetic analysis and array comparative genomic hybridization profiling identified duplication of parts of chromosomes 5 and 20. Microarray gene expression analysis demonstrated that this cell line clustered with non-immortalized bronchial cells, distinct from lung cancer cell lines.

HBECs were grown to 70–90% confluence in keratinocyte serum free medium supplemented with human recombinant Epidermal Growth Factor and Bovine Pituitary Extract (Life Technologies, Carlsbad, CA). For the CRISPR plasmid Activation of TET1, control plasmid and TET1 CRISPR activation plasmid (Santa Cruz Biotechnology, Dallas, TX) was transfected into HBECs using Lipofectamine® 2000 (Thermo Fisher Scientific, Florence, KY) and Opti-MEM medium according to manufacturer’s protocol. Briefly, 1 ug plasmid DNA was diluted with Opti-MEM. Lipofectamine® 2000 Reagent was diluted, in a separate tube, with Opti-MEM. The DNA plasmid dilution was then mixed with the Lipofectamine® 2000 dilution (1:1 ratio) for 5 minutes. The DNA-lipid complex was added to cells in a 12-well plate and incubated for 48 hours before harvest. For the siRNA knockdown of TET1, two siRNA targeting TET1 products were combined to increase knockdown efficiency (Thermo Fisher scientific and Santa Cruz). Lipofectamine® RNAiMAX™ Transfection Reagent (Thermo Fisher Scientific, Florence, KY) was used.

RNA extraction and RT-qPCR

Mouse lungs were homogenized using ceramic beads (Omni Inc, Kennesaw, Georgia) and the Omni Bead Ruptor 24 (Omni Inc, Kennesaw Georgia). Total RNA was isolated from homogenized mouse lung using Trizol (Invitrogen, Carlsbad, CA) and RNA was purified using the RNeasy Mini Kit (Qiagen, Valencia, CA) according to manufacturer’s instructions. DNase treatment was performed (Qiagen, Valencia, CA) before being reverse transcribed with SuperScript™ IV Vilo Master Mix™ (Invitrogen, Carlsbad, CA). Reverse-Transcriptase quantitative PCR (RT-qPCR) analysis was done using LightCycler FastStart DNA master SYBR green I as a ready-to-use reaction mixture (Roche, Indianapolis, Indiana). Targeted genes were amplified from cDNA using the primers listed in Supplementary Table S11 and gene expression was normalized to Rpl13a (mouse) or GAPDH (human). Melting curves have been generated for all assays and one product for each gene has been identified.

RNA-seq library preparation and sequencing

1 ug RNA per sample was used for each RNA-seq experiment. Samples from two animals in each group (Tet1+/+ Saline, Tet1+/+ HDM, Tet1−/− Saline, Tet1−/− HDM) were included. Two technical replicates from the Tet1 activation experiment in HBECs were analyzed. The purity and concentration of the isolated RNA was quantified using the Nanodrop 1000 (Thermo Fisher Scientific, Wilmington, DE) and 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA). A total of 1 μg RNA for each sample with RNA integrity number (RIN) values greater than 8 was used for library construction. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Following purification, the mRNA was fragmented and reverse-transcribed to create the final cDNA library in accordance with the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, USA). cDNA libraries were sequenced on an Illumina HiSeq™ 2000 platform, and reads were generated in 100-bp single-end format.

RNA-seq data analysis and functional genomics analyses

All fastq files were adapter trimmed prior to alignment. Files were aligned against Human (GRch37/Ensembl) or Mouse reference genome (GRCm38/mm10). Sequence alignment was performed using Bowtie2 and RSEM101. Post-alignment, raw read counts were extracted and normalized to account for difference in read depth using DESeq102. Following normalization, gene read counts were filtered to remove low read values: maximum Reads Per Kilobase of transcript, per Million mapped reads (RPKM) across all samples must be ≥1. Differentially expressed genes (DEGs) were identified using DESeq, and genes with an FDR adjusted p-value (q-value) <0.01 and fold of change ≥1.2, were considered significant. P values were calculated by student t test for each individual gene, and Benjamini and Hochberg correction was applied to generate FDR (q-value).

Gene ontology and pathway analysis were performed in Ingenuity Pathway Analysis (IPA, Ingenuity Systems, Redwood City, CA) and a cutoff of 0.05 was used for statistical significance in IPA analysis. Protein-protein interaction networks and modules were exacted using a web-based tool103,104 and 1st order subnetwork or minimal networks were shown. Enriched transcription factor (TF) interaction and histone marks at promoter regions (2 kb upstream TSS) of differentially expressed genes over a background data set (containing genes that are not differentially expressed from the same RNA-seq analysis) were identified as previously described10.

Reduced representation bisulfite sequencing (RRBS) sample processing, alignment, methylation calling and differential methylation analysis

After mouse lungs were homogenized (see above in RNA extraction), DNA was extracted using DNeasy Blood and Tissue Kits (Qiagen, Valencia, CA) per the manufacturer’s instructions. Samples from two animals in Tet1+/+ HDM and Tet1−/− HDM groups were included. For RRBS, 100 ng gDNA and 10 units MspI restriction enzyme were added in 20 ul reaction, and then incubated at 37 °C for 2 hours. MspI digests the CCGG sites and generated CG overhang at 5′ end, which were end-repaired and A-tailed using the T4 DNA polymerase and Klenow fragments, and subsequently ligated with Illumina adapters. The DNA samples were then treated with sodium bisulfite which converts unmethylated cytosine to uracil while methylated cytosine remains unaffected. Size selection was performed to obtain the optimal fragments for genome coverage and remove restriction fragments that failed to ligate with the adapter. Purified DNA then underwent the minimum number of cycles to produce an evenly represented library. 20 million SE75 reads were generated from Illumina HiSeq 2500.

FASTQ files were obtained from the DNA Sequencing and genotyping Core facility of CCHMC. Quality control steps were performed to determine overall quality of the reads from the FASTQ files. Upon passing basic quality matrices, the reads were trimmed to remove adapters and low-quality reads using Trim galore. The trimmed reads were then mapped to the bisulfite converted mouse reference genome mm10 using bismark, which created a summary of alignment, overall methylation profile and output BAM files. Next, the methylation levels of all the covered CpG sites were obtained using methylation extractor option of Bismark tool. Their association with a CpG-Island and/or a promoter and read coverage of these sites were calculated using R. An R package called MethylKit was used to identify differentially methylated sites105. P values were calculated by logistic regression for each individual CpG site. The sliding linear model (SLIM method106) was applied to generate FDR (q-value). The CpG sites were annotated to genes using a custom PERL script.

Bisulfite Pyrosequencing

A total of 100 ng genomic DNA was subjected to sodium bisulfite treatment. Standard PCR amplification reactions were performed to amplify targeted gene fragments at an annealing temperature of 50 °C before being subjected to pyrosequencing. The generated pyrograms were automatically analyzed using PyroMark analysis software (Qiagen, Valencia, CA, USA). Pyrosequencing assay design and genomic coordinates are documented in Supplementary Table S9.

Statistical analysis

All statistical analyses were performed using GraphPad Prism 7.0 with a p value cutoff of 0.05 considered significant. Normality of data were checked using Shapiro-Wilk normality test and α was set to be 0.05. When data are normally distributed, groups were compared using student’s t-tests. For data that are not normally distributed, Mann Whitney tests were performed. For multiple comparisons, Bonferroni post hoc corrections were used.