Loss of TET1 decreases chromatin accessibility in human bronchial epithelial cells.
In order to understand how TET1 regulates gene expression in HBECs, we knocked down TET1 using siRNA (48 hour incubation led to ~ 72% downregulation, Supplementary Fig. 1) and examined the impact on chromatin accessibility using ATAC-seq. Compared to control WT saline cells, cells with deficient TET1 expression (KD saline) showed a general decrease in chromatin accessibility (Fig. 1). Overall, we identified 6975 differentially accessible peaks (FDR ≤ 0.05, fold-change ≥ 2) in the KD saline vs. WT saline comparison, which were annotated to 4669 unique genes (Supplementary Table 1, selected candidates shown in Fig. 2 and Supplementary Fig. 2). The majority of these peaks (84.7%) were less accessible in KD saline cells. IPA analysis on the genes associated with significantly differentially accessible ATAC-seq peaks revealed 111 significantly enriched canonical pathways (FDR ≤ 0.05; Supplementary Table 2, top pathways by cumulative absolute activation z-score shown in Fig. 3A). Some enriched pathways of interest included Leukocyte extravasation signaling, IL-15 production, and Pulmonary Fibrosis Idiopathic Signaling Pathway. Genes involved in oxidative stress response (e.g. TLR3, TLR5, TLR10, ALDH1A1, ALDH1B1, ALDH2, and ALDH7A1), type I immune response (e.g. STAT3, STAT4, STAT5B, IFNGR1, and OAS1), type II immune response (e.g. STAT3, STAT4, and STAT5B) and type 17 immune response (e.g. AHR, AHRR, and IRF4) all showed changes in accessibility (mostly reduced) when TET1 was low (Supplementary Table 1). Together, these data show that TET1 promotes chromatin accessibility in HBECs.
Loss of TET1 leads to changes in methylation, primarily in regions where chromatin accessibility did not change.
It is well established that TET1 regulates DNA methylation (56), and that DNA methylation has a direct impact on transcription factor binding (57, 58). Some evidence suggests that DNA methylation may directly influence chromatin accessibility as well (59). Therefore, we performed genome-wide bisulfite sequencing to evaluate the DNA methylation landscape in TET1-deficient HBECs. Given that TET1 is a demethylase, we expected that TET1 knock down samples would have generally higher methylation levels than samples where TET1 was intact. However, this was not the case. When comparing TET1 KD saline samples to WT saline samples, we identified 4359 significant DMRs (empirical p-value ≤ 0.01, ≥10% difference) annotated to 3224 unique genes (Supplementary Table 3). Of these DMRs following TET1 knockdown, 49.1% were hypermethylated (Fig. 1). These DMR-associated genes were enriched for 94 IPA pathways (Supplementary Table 4, top pathways in Fig. 3B). Enriched pathways of interest include Pulmonary Fibrosis Idiopathic Signaling Pathway, IL-8 Signaling, and IL-15 Production. Out of these 4359 DMRs between KD saline and WT saline, only 14 had overlapping coordinates with one of the 5479 differentially accessible regions between these two sample groups (Fig. 4A). Collectively, although WGBS measures the sum of 5mC and 5hmC and may miss regions that undergo simultaneous and similar opposite changes in both marks (60), these data suggest that TET1 independently influences chromatin accessibility and DNA methylation in HBECs.
TET1 knockdown-induced chromatin accessibility and methylation changes influence gene expression.
We next performed RNA-sequencing to assess the effect of TET1 knockdown on gene expression and to attempt to correlate changes in gene expression with changes in chromatin accessibility and/or DNA methylation. We identified 5242 genes that were differentially expressed in the KD saline vs. WT saline comparison (Supplementary Table 5). There was a fairly even split of upregulation and downregulation in these DE genes (51.9% upregulated, Fig. 1). These differentially expressed genes were enriched for IPA pathways of interest such as IL-8 Signaling, HIF1α Signaling, Aryl Hydrocarbon Receptor Signaling, and NRF2-mediated Oxidative Stress Response (Supplementary Table 6, top pathways in Fig. 3C). Importantly, these changes in gene expression correlated with changes in DNA methylation and chromatin accessibility following TET1 knockdown. Genes near regions with changes in accessibility tended to have reduced expression in KD saline (Fig. 5A). Genes with more accessible promoters in KD saline showed a trend for increased expression (not significant, Fig. 5B), while genes with less accessible promoters were significantly more likely to have decreased gene expression (Fig. 5C). Genes near DMRs were significantly more likely to be downregulated in KD saline (Fig. 5D). Genes with hypermethylated promoters in KD saline showed a trend of having decreased expression (Fig. 5E), while genes with hypomethylated promoters showed a trend of having increased expression (Fig. 5F).
To further compare the effects of methylation and accessibility on gene expression, we also performed IPA pathway analyses on genes that were both associated with a DMR and differentially expressed following TET1 knockdown (but not differentially accessible, DMR + DE only). These genes (7.2% of all DE genes, Fig. 4B) were enriched for 9 IPA canonical pathways (Supplementary Table 7). We compared these pathways to enriched pathways in genes that were simultaneously associated with a differentially accessible peak and differentially expressed, but not differentially methylated (ATAC + DE only). These genes (14.3% of all DE genes, Fig. 4B) were enriched for 132 IPA canonical pathways (Supplementary Table 7). Out of the nine enriched pathways from the DMR + DE only set, seven were also enriched in the ATAC + DE only set. On the other hand, some of the 126 pathways enriched in ATAC + DE only but not DMR + DE only include EIF2 Signaling, IL-8 Signaling, NRF2-mediated Oxidative Stress Response, and Xenobiotic Metabolism Signaling. There were also 228 genes (4.4% of all DE genes, Fig. 4B) that were simultaneously DE, differentially methylated, and differentially accessible (ATAC + DMR + DE). There was a physical overlap between the DMR and the differentially accessible region in only 14 of these genes. The ATAC + DMR + DE genes were only enriched for one IPA canonical pathway (Protein Kinase A Signaling), and this pathway was not enriched in either the ATAC + DE only or DMR + DE only groups. Overall, these results imply that changes in chromatin accessibility and DNA methylation following TET1 knockdown both impact gene expression, but they often affect different genes, pathways, and regions of the genome.
HDM challenge closes chromatin in a similar fashion to TET1 knockdown in HBECs.
HDM promotes allergy and asthma in humans (61, 62). We previously found that loss of TET1 exaggerated HDM-induced allergic airway inflammation in mice (9). In our current experiments, HDM significantly downregulated Tet1 expression following 24hrs of challenge (Supplementary Fig. 1). Therefore, we hypothesized that TET1 loss and HDM treatment have a similar impact on chromatin accessibility. We observed 2894 ATAC-seq peaks that were significantly different between WT HDM and WT saline cells (FDR ≤ 0.05, fold-change ≥ 2; Fig. 1, Supplementary Table 1, annotated to 2328 genes). Similar to the changes we observed to chromatin accessibility following TET1 knockdown, HDM treatment primarily led to decreased chromatin accessibility (88.3%) in regions where accessibility was altered. 65.8% of these peaks overlapped with a significant peak from the KD saline vs. WT saline comparison, significantly more than expected by chance (p ~ = 0, Fig. 6A). Even more striking is that for each overlap, the direction of change in the HDM-treated group was the same as the direction of change in the KD saline group (i.e. if accessibility increased in KD saline vs. WT saline, it also increased in WT HDM vs. WT saline). A principal component analysis on the chromatin accessibility data revealed that WT HDM samples clustered closest with KD saline samples, with these two groups split from WT saline and KD HDM along PC1 (Supplementary Fig. 3). Genes associated with differentially accessible peaks were significantly enriched for 88 IPA pathways (Fig. 3A, Supplementary Table 2). The top five most enriched pathways in genes associated with differentially accessible peaks following HDM treatment were all in the top 32 most enriched pathways following TET1 knockdown. Additionally, all five of these pathways had the same inferred direction of change. 65 (73.0%) enriched pathways following HDM treatment were also enriched following TET1 knockdown, including TGF-β Signaling and Pulmonary Fibrosis Idiopathic Signaling Pathway (Supplementary Table 2). Of these 65, 55 had similar z-scores (i.e. the inferred direction of change in activation was the same), 8 had no activation z-scores calculated (due to pathway limitations), and 2 had opposing activation z-scores. Collectively, our data show that HDM challenge and TET1 deficiency have similar impacts on chromatin accessibility in airway epithelial cells, with both likely impacting proinflammatory pathways and lung function.
HDM treatment also alters DNA methylation patterns in a similar fashion to TET1 knockdown in HBECs.
We previously showed that HDM challenges led to changes in 5mC and 5hmC in HBECs using microarrays, and changes in 5mC responsive to HDM and diesel particles significantly correlated with changes in gene expression (14). When we compared methylation patterns between WT HDM HBECs and WT saline HBECs using WGBS, we identified 4718 significant DMRs (empirical p-value ≤ 0.01, at least 10% difference in methylation) annotated to 3428 unique genes (Supplementary Table 3). Of these DMRs, 48.8% were hypermethylated in WT HDM compared to WT saline (Fig. 1). There were significantly more positional overlaps than expected between DMRs in the KD saline vs. WT saline comparison and DMRs in the WT HDM vs. WT saline comparison (p = 1.2 x 10− 92, Fig. 6B). For these positional overlaps, methylation changed in the same way (i.e. increased or decreased) after either TET1 knockdown or HDM treatment 94.8% of the time. These positional overlaps were annotated to genes such as CLU, TLR7, and IL13RA1 (Supplementary Table 8). At the gene level, out of the genes that were associated with DMRs following TET1 knockdown, 36.7% were also associated with at least one DMR following HDM treatment. The DMR-associated genes following HDM treatment were enriched for 70 IPA canonical pathways (Supplementary Table 4). These pathways were generally similar to enriched pathways from DMR-associated genes in the KD saline vs. WT saline comparison; there were six overlapping pathways in the top ten enriched pathways from each comparison. Out of the enriched pathways from DMR-associated genes following TET1 knockdown, 59.6% were also enriched in DMR-associated genes following HDM treatment, including Pulmonary Fibrosis Idiopathic Signaling Pathway, Nitric Oxide Signaling in the Cardiovascular System, Leukocyte Extravasation Signaling, and IL-8 Signaling. Overall, HDM treatment and TET1 knockdown had similar effects on the methylome, potentially affecting genes and pathways contributing to inflammation and lung function.
HDM treatment led to few changes in expression, but these changes were similar to changes following TET1 knockdown.
We next performed RNA-sequencing to assess the effect of HDM treatment on gene expression, and to compare these results with how TET1 knockdown affected gene expression. We identified 104 DE genes in the WT HDM vs. WT saline comparison (Supplementary Table 5). There were many fewer DE genes following HDM treatment alone compared to after TET1 knockdown alone (104 vs. 5242 genes). Of these DE genes, 64.4% were upregulated following HDM treatment. Out of these DE genes, 77.9% were also DE in the KD saline vs. WT saline comparison, indicating that most genes with changed expression following HDM treatment were also affected by TET1 knockdown (Fig. 6C). Further supporting this notion was the finding that out of these shared DE genes, 91.4% changed in the same direction (i.e. HDM treatment and TET1 knockdown both led to downregulation or both led to upregulation). These genes whose expression was altered similarly following HDM treatment or TET1 knockdown included IL1A (increased), IL1B (increased), IFNGR2 (increased), and AHR (decreased) (Supplementary Fig. 4A). We also performed IPA pathway analysis on the DE genes from the WT HDM vs. WT saline comparison to see if expression levels in similar pathways were affected by both HDM treatment and TET1 knockdown (Fig. 3C, Supplementary Table 6). Of 22 enriched IPA pathways, 11 (50%) were also enriched in DE genes from the KD saline vs. WT saline comparison. Out of these 11 shared pathways, 7 had estimated activation z-scores in the WT HDM vs. WT saline comparison (Supplementary Fig. 4B); in all 7 cases, the activation z-scores for both comparisons indicated consistent pathway level changes. For example, the activation z-scores for Aryl Hydrocarbon Receptor Signaling were negative in both cases, indicating that expression changes following HDM treatment alone and TET1 knockdown alone both likely led to some deactivation of that pathway. Overall, despite HDM treatment leading to many fewer DE genes than knocking down TET1 in HBECs, genes and pathways that were affected by HDM treatment were frequently also affected by TET1 knockdown as well.
Loss of TET1 results in a transcriptomic signature enriched in stress and immune response following HDM challenges in HBECs.
As our data demonstrated that TET1 loss tends to have a similar impact to HDM challenge on chromatin accessibility, DNA methylation and gene expression, we next compared expression in KD HDM HBECs to WT HDM HBECs through RNA-seq analysis. In a pairwise analysis, we identified 4457 genes that were differentially expressed (FDR ≤ 0.05, FC ≥ 1.2, Fig. 1, Supplementary Table 5). Of these DE genes, 49.9% were significantly upregulated in KD HDM compared to WT HDM HEBCs (Fig. 1). The majority of these genes (67.6%) were also differentially expressed when comparing KD saline vs. WT saline. Out of the overlapping DE genes, 99% showed changes in the same direction when TET1 was knocked down, suggesting that these gene expression changes were primarily driven by TET1 loss regardless of the presence of HDM (Supplementary Fig. 3). Pathway analysis via IPA on all DE genes from the KD HDM vs. WT HDM comparison (Fig. 3C, Supplementary Table 6) revealed the enrichment of 269 IPA pathways. These genes were enriched for pathways of interest such as Oxidative Phosphorylation, Mitochondrial Dysfunction, HIF1α signaling, NRF2-mediated Oxidative Stress Response, and Aryl Hydrocarbon Receptor Signaling. Additionally, the significant enrichment and strong activation (activation z-score = 3.97) of EIF2 Signaling in DE genes in KD HDM vs. WT HDM HBECs suggests an upregulation of protein translation. Interestingly, we also observed significant enrichment for genes involved in the regulation of cell cycle and division (Cell Cycle: G2/M DNA Damage Checkpoint Regulation, Cyclins and Cell Cycle Regulation, Cell Cycle: G1/S Checkpoint Regulation, etc.). Significant downregulation was observed in several cyclin and cyclin-dependent kinases (CCND1, CCND3, CDC25A and CDK4/6) and transcriptional regulators required for E2F-mediated transcription (E2F6), which may block G1/S cycle progression, and this is similar to what we observed in KD saline compared to WT saline.
To further elucidate the role of TET1 in regulating responses to HDM, we separately analyzed the 1444 genes that were DE in KD HDM vs. WT HDM but not in KD saline vs. WT saline. There were several DE genes in this set that are subunits of the main complexes in the electron transport chain, including ATP5F1C, ATP5MC1, COX5A, NDUFA4, NDUFA7, UQCR11, and UQCRQ. There were also immune and detoxification genes such as ARNT, FLOT2, GSTO2, JUN, JUNB, IL17RE, IRF2BP2, SMAD4, TLR5, and TRAF4 (examples shown in Supplementary Fig. 5). These genes were enriched for seven IPA pathways, including Oxidative Phosphorylation, Mitochondrial Dysfunction, and Autophagy (Supplementary Table 9). All of these pathways were also enriched in the KD saline vs. WT saline DE genes (Supplementary Table 6) even though overlapping genes were specifically filtered out for this analysis, indicating that TET1 broadly regulates these pathways but that the specific genes that are impacted by the loss of TET1 are at least partially context-dependent.
HDM influences the regulation of chromatin accessibility and DNA methylation by TET1.
In contrast to the mostly reduced accessibility observed in KD saline cells compare to WT saline cells, we observed a strong pattern of increased accessibility in KD HDM cells compared to WT HDM (Fig. 1). There were 6404 significantly different ATAC-seq peaks (Supplementary Table 1, annotated to 4449 unique genes), and 92.4% of the peaks showed increased peak intensity in KD HDM (Fig. 1). Of these genes associated with differential chromatin accessibility in KD HDM vs. WT HDM, 51.8% were also associated with differential accessibility in the KD saline vs. WT saline comparison. Only 10.5% of these genes overlapping genes, however, showed consistent changes in direction in the two comparisons (i.e. a gene increasing in accessibility in KD HDM vs. WT HDM and also increasing in KD saline vs. WT saline). Out of the differentially accessible peaks in KD HDM vs. WT HDM, 30.8% overlapped positionally with a differentially accessible peak in KD saline vs. WT saline (p ~ = 0). For these positional overlaps, in only 1.4% did TET1 knockdown have a consistent effect of either increasing or decreasing accessibility regardless of whether HDM was present or not. These results indicate a likely interaction between TET1 knockdown and the presence of HDM since the effect of TET1 knockdown was different depending on the environment. This was further supported by a principal component analysis where PC1 (which explained 65.8% of the variance) separated WT saline and KD HDM samples from WT HDM and KD saline samples (Supplementary Fig. 3). Pathway analysis via IPA on genes associated with changes in accessibility revealed significant enrichment for 112 pathways (Supplementary Table 2), 61.6% of which also were enriched in significant ATAC peaks from KD saline vs. WT Saline. Of the pathways where an activation z-score could be assigned (assuming a positive correlation between chromatin accessibility and gene expression), 96.8% showed opposite patterns of activation in KD HDM vs. WT HDM and KD saline vs. WT saline (Supplementary Fig. 4C). This indicates that knocking down TET1 consistently leads to changes in chromatin accessibility in certain pathways, but that precisely how accessibility is affected depends on whether HDM is present or not.
When we compared genome-wide DNA methylation levels in KD HDM samples to WT HDM samples, we identified 4086 DMRs (Fig. 1) annotated to 3037 unique genes (Supplementary Table 3). The 3037 DMR-associated genes were enriched for 68 IPA pathways (Fig. 3B, Supplementary Table 4). Of the genes that were associated with a DMR in this comparison (KD HDM vs. WT HDM), 34.4% were also DMR-associated in the KD saline vs. WT saline comparison. Of these overlapping genes, the average methylation changes were in the same direction following TET1 knockdown 49.4% of the time. There were 44 positional overlaps between KD HDM vs. WT HDM DMRs and KD saline vs. WT saline DMRs (more than expected by chance, p = 6.5 x 10− 15). There were consistent directional changes in methylation regardless of whether HDM was present in 43.2% of these positional overlaps. These results imply that, similar to chromatin accessibility, methylation changes following TET1 loss were influenced by the presence of HDM. There were only 10 positional overlaps between KD HDM vs. WT HDM DMRs and differentially accessible regions, adding another piece of evidence that the changes we observed in chromatin accessibility following TET1 knockdown generally happened in different regions of the genome than changes in DNA methylation.
Loss of TET1 and HDM treatment lead to similar changes in H3K27ac, and these changes are enriched for asthma-associated changes in H3K27ac.
Since H3K27ac is highly enriched at active enhancers (63), we performed H3K27ac ChIP-seq and searched for regions with changes in H3K27ac signal intensity. We identified 2848 differential regions in which KD saline and WT saline groups had similar signal strength, with only 16 regions showing increased H3K27ac in KD saline and 2832 regions showing decreased H3K27ac in KD saline (Supplementary Table 10). Out of these 2848 regions, 34.2% were annotated to a differentially expressed gene, of which 57.1% showed the generally expected positive correlation between H3K27ac levels and expression (enriched pathways in Supplementary Table 11). Interestingly, 33 of these differences overlapped positionally with changes in chromatin accessibility (p = 0.03, Fig. 7A), including an overlap at MAP3K14 (Fig. 2B), which stimulates NFKB signaling (64).
As enhancers may cluster together and form super-enhancers to drive transcription of genes implicated in cell identity (65, 66) and disease (51, 67), we next compared the presence of super-enhancers in KD saline vs. WT saline samples. We identified 109 consensus “gained” super-enhancers in KD saline samples that did not overlap with any consensus super-enhancers in WT saline samples, and 329 “lost” consensus super-enhancers in WT saline samples with no overlaps in the KD saline group (Supplementary Table 12). These 438 “gained” or “lost” super-enhancer regions overlapped with 140 DE genes (i.e. the closest gene to the super-enhancer was DE) and 146 overlapped physically with regions with changes in chromatin accessibility (p = 1.75x10− 62, Fig. 7A). A positive correlation between “gain” or “loss” of a nearby super-enhancer and changes in expression might be generally expected, but we observed this expected correlation only 41.4% of the time. Overall, we generally observed decreased H3K27ac signal in KD saline samples compared to WT saline samples, but these changes were not well correlated with changes in gene expression, possibly because other regulatory changes were underlying gene expression changes, changes in expression came after the time point that we sampled, or our “nearest gene” annotation strategy might not always connect a regulatory region to its affected gene.
Similar to what we observed following TET1 knockdown, a majority of regions with significant differences in H3K27ac showed decreased levels in WT HDM compared to WT saline (374/504, 74.2%) (Supplementary Table 10). These changes in H3K27ac overlapped with 9 DE genes, and there was a positive correlation between changes in H3K27ac and gene expression in 7 of these instances (77.8%). These regions largely overlapped with differential H3K27ac regions from the KDSal vs. WTSal comparison (Fig. 6D) Also similar to what we observed following TET1 knockdown, there were many more “lost” super-enhancers in WT HDM (329) than gained super-enhancers (86) in WT HDM (Supplementary Table 12). These 415 “gained” or “lost” super-enhancer regions overlapped with only 4 DE genes (there were only 104 DE genes in WT HDM vs. WT saline in total), and in only 1 out of 4 of these DE genes did we observe a positive correlation between super-enhancer presence and expression changes. When comparing KD HDM to WT HDM, out of 1479 regions showing significant changes in H3K27ac in regions where both groups had H3K27ac signal, 99.4% showed decreased levels of H3K27ac in KD HDM (Supplementary Table 10). There were 482 of these regions with significant differences in H3K27ac that overlapped with a significant DE gene, and 67% showed a positive correlation with changes in expression. However, unlike what we observed following TET1 knockdown and HDM treatment alone (both led to many more “lost” super-enhancers based on analysis of H3K27ac ChIP-seq), the number of “lost” super-enhancers (194) was similar to the number of “gained” super-enhancers (226) in KD HDM (Supplementary Table 12). These 420 “gained” or “lost” super-enhancer regions overlapped with 116 DE genes (there were only 104 DE genes in WT HDM vs. WT saline in total), but only 47.4% of the time did we observe a positive correlation between super-enhancer presence and expression changes. These super-enhancer results provide another example of the combination of TET1 knockdown and HDM treatment having an interactive effect (rather than additive) on epigenomic regulation. Similar to the results from KD saline vs. WT saline, there was only limited correlation between super-enhancer gain or loss and expression, but there was some correlation between changes in H3K27ac levels and expression. Overall, this indicates that TET1 knockdown and HDM treatment both lead to similar decreases in H3K27ac presence, again highlighting that these two conditions have a similar impact on the overall regulatory landscape.
Besides enhancer-specific markers, we also observed significant overlap between TET1-regulated differentially accessible regions and all A549 ENCODE (25) predicted cCREs (candidate cis-regulatory elements, defined by chromatin accessibility, CTCF binding, H3K27ac and H3K4me3) and CTCF-bound cCREs (Fig. 7A). As CTCF is important in setting up 3-D chromatin interactions and topologically associating domains (TADs), we also quantified the overlap between differentially accessible regions and TADs in A549 cells (GSE92819) and found significant overlap (Fig. 7A). In addition, TET1 regulates the expression of genes involved in asthma to protect against allergic asthma (9, 68). Therefore, we searched for the presence of asthma-associated variations in DNA methylation and histone acetylation found in airway epithelial cells among TET1-loss induced differentially accessible regions and DMRs. We observed significant enrichment for asthma-associated H3K27ac changes (51), but not asthma-associated DMRs or DMPs (52–54) (Fig. 7A and 7B). Therefore, our data support that TET1-mediated chromatin accessibility might influence H3K27ac marks and cis-regulatory elements and CTCF binding to regulate gene expression, which may explain its role in protection from inflammation and asthma.
Altered chromatin accessibility and DNA methylation following loss of TET1 and HDM treatment may affect CTCF and CEBP binding.
As altered chromatin accessibility may impact the binding of TFs and therefore lead to changes in gene expression, we next performed unbiased HOMER and RELI transcription factor enrichment analyses. Using a custom HOMER database of all human transcription factor binding motifs obtained from the CisBP database (31) (see Methods), we compared TF binding motif enrichment in peaks that had changed accessibility following TET1 KD with peaks that did not show any significant change in accessibility (Fig. 8). We also performed a similar analysis with RELI, which utilizes results from publicly available ChIP-seq experiments instead of TF binding motifs (Fig. 8). In regions where KD saline had reduced accessibility, both analyses indicated lower enrichment of CTCF and CEBP binding sites compared to regions that were not changing in accessibility (Fig. 8). In regions where KD saline had increased accessibility, however, both analyses showed increased enrichment of CEBP binding sites and decreased enrichment of CTCF binding sites. These data suggest the TET1-regulated chromatin accessibility may specifically impact CTCF and CEBP binding. In contrast, we observed no enrichment for either CTCF or CEBP binding sites within DMRs (Supplementary Fig. 6).
In order to further test the association between CTCF and CEBP binding and changes in chromatin accessibility, we compared binding sites for CTCF (GSM1003606) and CEBPB (GSM935630) in A549 cells to our differentially accessible peaks. The differentially accessible peaks overlapped much more frequently than expected by chance with the CTCF peaks (Fig. 7A). DMRs between KD saline HBECs and WT saline HBECs, however, did not overlap more frequently than expected with the CTCF binding sites from A549 cells (Fig. 7B), emphasizing how the loss of TET1 likely affects gene regulation in different ways in different types of genomic environment. There was also significant overlap between both differentially accessible peaks (Fig. 7A) and DMRs (Fig. 7B) and CEBPB binding sites in A549 cells. The much lower p-value and much greater odds ratio for the differentially accessible peaks, however, imply a more significant overlap between differentially accessible peaks and CEBPB sites than between DMRs and CEBPB sites. Overall, TET1 loss-associated changes in accessibility were clearly enriched for both CTCF and CEBPB motifs and ChiP-seq peaks, while changes in methylation showed some evidence for enrichment of CEBPB ChIP-seq peaks.
Similarly, TF motif enrichment analysis (HOMER) revealed that in peaks where HDM treatment led to increased accessibility, there was more enrichment for CEBP binding sites and fewer CTCF binding sites than in peaks where there were no changes in accessibility (Fig. 8). In peaks where HDM treatment led to reduced accessibility, both CTCF and CEBP binding sites were less than expected (Fig. 8). These HOMER results were very similar to the HOMER results from the KD saline vs. WT saline comparison, implying that chromatin changes caused by HDM treatment and TET1 knockdown have similar impacts on TF binding. The RELI analysis (which compares significant results to publicly available ChIP-seq studies rather than searching for TF binding motifs) showed similar enrichment patterns as the HOMER analysis (Fig. 8). Despite observing opposing changes in accessibility, similar to the results when TET1 was knocked down in saline, we observed less overlap with CTCF binding sites than expected in regions with increased accessibility in KD HDM (Fig. 8). Together, these results support that altered chromatin accessibility and DNA methylation following loss of TET1 and HDM treatment may alter the binding of CTCF and CEBP.
A TET1-mediated regulatory element at the IL1B locus enhances gene expression to the same extent as HDM challenge.
HDM and TET1 loss both activate the acute phase response pathway and downregulate the AHR signaling pathway (Supplementary Fig. 4B). IL1B, a major stimulator of the acute phase response pathway (69) that has been linked to asthma (70), showed increased gene expression following both TET1 loss and HDM treatment (Supplementary Fig. 4A). Interestingly, a region downstream of IL1B became significantly less accessible following HDM treatment and became significantly more accessible in KD HDM compared to WT HDM. This region overlapped with a predicted cCRE from A549 ENCODE data (Fig. 9A). Reporter assays confirm that this region enhances the expression of a florescence gene reporter gene at baseline, to the same extent of the effect of HDM treatment alone (Fig. 9B). These data further support that TET1 may regulate chromatin accessibility and enhancer activity to control gene expression and response to allergen in HBECs.