TET1 contributes to allergic airway inflammation and regulates interferon and aryl hydrocarbon receptor signaling pathways in bronchial epithelial cells

Previous studies have suggested a role for Tet1 in the pathogenesis of childhood asthma. However, how Tet1 contributes to asthma remains unknown. Here we used mice deficient for Tet1 in a well-established model of allergic airway inflammation and demonstrated that loss of Tet1 increased disease severity including airway hyperresponsiveness and lung eosinophilia. Increased expression of Muc5ac, Il13, Il33, Il17a, Egfr, and Tff2 were observed in HDM-challenged Tet1-deficient mice compared to Tet1+/+ littermates. Further, transcriptomic analysis of lung RNA followed by pathway and protein network analysis showed that the IFN signaling pathway was significantly upregulated and the aryl hydrocarbon receptor (AhR) pathway was significantly downregulated in HDM-challenged Tet1−/− mice. This transcriptional regulation of the IFN and AhR pathways by Tet1 was also present in human bronchial epithelial cells at base line and following HDM challenges. Genes in these pathways were further associated with changes in DNA methylation, predicted binding of transcriptional factors with relevant functions in their promoters, and the presence of histone marks generated by histone enzymes that are known to interact with Tet1. Collectively, our data suggest that Tet1 inhibits HDM-induced allergic airway inflammation by direct regulation of the IFN and AhR pathways.

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.
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   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.  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).
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 pathogenesis 18,52,53 , airway epithelial cells are known to activate interferon and AhR signaling pathways following HDM challenges [54][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α isoforms 54 (Fig. 6c,d) and IFNβ (Fig. 6e). IFNλ1/2 expression was not detectable.
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 signaling 64,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 promoters [21][22][23][24] , or binding to transcriptional co-factors such as Suz12 or the SIN3A complex to change histone modifications [31][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, 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. www.nature.com/scientificreports www.nature.com/scientificreports/ Continued www.nature.com/scientificreports www.nature.com/scientificreports/ 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.
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 package 66 and motifs taken from the Cis-BP database 67 , 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 method 9 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 cells [21][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 structure [68][69][70][71][72][73][74][75] . "Readers" of 5hmC marks activate downstream gene expression networks [21][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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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  www.nature.com/scientificreports www.nature.com/scientificreports/ 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 expression [31][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 lines 76 . 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 lungs 16 . 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 AHR 77 , 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 responses 78 , 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 vivo 79,80 . In addition, loss of IRF7 did not significantly alter HDM-induced allergic immune responses in mice 81 . However, it has been well established that respiratory viral infection leads to asthma exacerbation in human studies and in mouse models 82 . 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 damage 84 . This free dsDNA binds to the TLR receptor and triggers signaling through IRF7 to promote downstream IFN responses 55 . 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 cells 54,55 , is known to regulate innate and adaptive type I and II immune responses (dendritic cells/Th1/Th2/Treg) following viral/bacterial/allergen challenges [85][86][87][88] and contribute to asthma development/exacerbation 89,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 function 91 . This is consistent with the presence of a T1-high group among severe asthmatics based on transcriptomic analysis 92 . Interestingly, 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.
www.nature.com/scientificreports www.nature.com/scientificreports/ Irf7 −/− and Ifnar2 −/− (deficient for type I IFN receptor) mice displayed normal dendritic cell development and allergic airway sensitization in response to HDM 81 , 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 multiciliogenesis 93 and oxidative stress responses 60 , and contributes to asthma development and exacerbation 46,62,94 . Specifically, AhR ligands attenuate allergic airway inflammation in OVA-challenged animal models 94,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 SOD1 96 . 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 lungs 59,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 particles 8,10 and cigarette smoke 98 . 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 miRNAs 34 , our findings may promote the development of new asthma therapies.

Methods
Murine model of allergic airway inflammation. Mice heterozygous for Tet1 tm1.1Jae /J (B6/129S4) were purchased from Jackson Labs (Bar Harbor, ME) and littermates of wild type (Tet1 +/+ ), heterozygous (Tet1 +/− ) and homozygous (Tet1 −/− ) for Tet1 tm1.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. www.nature.com/scientificreports www.nature.com/scientificreports/ 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). Total serum IgE and IgG1, HDM-specific IgE and IgG1 levels were measured as previously described 99,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 described 99,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. 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.

ELISA in serum, BALF and total lung homogenate.
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 RSEM 101 . Post-alignment, raw read counts were extracted and normalized to account for difference in read depth using DESeq 102 . 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). (2019) 9:7361 | https://doi.org/10.1038/s41598-019-43767-6 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 tool 103,104 and 1 st 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 described 10 .
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 sites 105 . P values were calculated by logistic regression for each individual CpG site. The sliding linear model (SLIM method 106 ) 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.

Data Availability
RNA-seq and RRBS data has be deposited in GEO (GSE124922).