SARS-CoV-2 infection perturbs enhancer mediated transcriptional regulation of key pathways

Despite extensive studies on the effects of SARS-CoV-2 infection, there is still a lack of understanding of the downstream epigenetic and regulatory alterations in infected cells. In this study, we investigated changes in enhancer acetylation in epithelial lung cells infected with SARS-CoV-2 and their influence on transcriptional regulation and pathway activity. To achieve this, we integrated and reanalyzed data of enhancer acetylation, ex-vivo infection and single cell RNA-seq data from human patients. Our findings revealed coordinated changes in enhancers and transcriptional networks. We found that infected cells lose the WT1 transcription factor and demonstrate disruption of WT1-bound enhancers and of their associated target genes. Downstream targets of WT1 are involved in the regulation of the Wnt signaling and the mitogen-activated protein kinase cascade, which indeed exhibit increased activation levels. These findings may provide a potential explanation for the development of pulmonary fibrosis, a lethal complication of COVID-19. Moreover, we revealed over-acetylated enhancers associated with upregulated genes involved in cell adhesion, which could contribute to cell-cell infection of SARS-CoV-2. Furthermore, we demonstrated that enhancers may play a role in the activation of pro-inflammatory cytokines and contribute to excessive inflammation in the lungs, a typical complication of COVID-19. Overall, our analysis provided novel insights into the cell-autonomous dysregulation of enhancer regulation caused by SARS-CoV-2 infection, a step on the path to a deeper molecular understanding of the disease.


Introduction
The COVID-19 pandemic, which emerged in December 2019 has had a profound global impact, resulting in millions of deaths and affecting the lives of countless individuals. This devastating disease is caused by the SARS-CoV-2 virus, and displays a range of clinical manifestations, varying from asymptomatic cases to severe illness. Viruses evolved mechanisms to manipulate their host cells, to ensure their maintenance, replication, and transmission. These manipulations include regulation of the cell cycle, control of apoptosis, and manipulation of the immune system. Since epigenetic modifications play a pivotal role in regulating cellular gene expression, some viruses use epigenetic mechanisms to manipulate the regulation of the host's and their own genes during infection [1][2][3].
Tremendous recent efforts to unravel the molecular effects of SARS-CoV-2 infection have led us to an improved understanding of the changes in gene expression and biological processes that contribute to disease severity. Investigations into SARS-CoV-2 infection demonstrated large-scale chromatin structural changes which redistribute active and inactive chromatin. Moreover, infection causes a reduction in the overall levels of histone 3 lysine 27 acetylation (H3K27ac), a key marker of active enhancers [4,5]. Other studies have revealed the virus's epigenetic impact on monocytes, a crucial component of the immune response [6]. Previous studies on SARS viruses [7] highlighted the value of epigenetic and regulatory analyses to improve understanding of the disease and devise effective strategies to combat it. However, studies of how SARS-CoV-2 infection affects histone modifications and enhancer regulation in infected lung epithelial cells, which genomic regions are affected, are still lacking, and hold great promise to further uncover mechanisms of pathogenesis and suggest novel therapies.
In this study, we focused on elucidating the cell-autonomous impact of SARS-CoV-2 infection on enhancer regulation in lung epithelial cells. To achieve this, we performed a comprehensive reanalysis and integration of publicly available data on changes in gene expression and H3K27 acetylation following SARS-CoV-2 infection to understand the impact of the infection on enhancer activity and gene regulation. Our findings reveal the upregulation of genes involved in cell adhesion and pro-inflammatory cytokine production, alterations of transcription factor (TF) networks, and differential acetylation of corresponding enhancers. By shedding light on the specific mechanisms through which SARS-CoV-2 infection impacts enhancer regulation in lung epithelial cells, our study contributes to a deeper understanding of the disease pathogenesis and highlights potential targets for therapeutic interventions.

Many TFs are significantly enriched at the differential acetylation peaks
We reanalyzed data from an H3K27ac ChIP-seq experiment performed on A549-ACE2 cells 8 hours and 24 hours after SARS-CoV-2 infection (see Methods). We found that 8 hours postinfection there were 744 putative enhancers that gained H3K27 acetylation, and 514 putative enhancers whose acetylation was decreased (FDR<0.05, fold-change > 150%). At 24h postinfection we found 3773 enhancers gaining acetylation and acetylation loss in 8741 enhancers ( Fig 1A and S1 Table).
To understand which transcription factors may be involved in the differentially acetylated enhancers, we performed motif enrichment analysis of regions displaying differential H3K27ac activity for each time point (FDR < 0.05). In enhancers that gained H3K27ac, we found an enrichment of motifs recognized by nuclear factor κB (NF-κB), GATA family, STAT family, IRF family, and EBF family, which are related to immune response. In addition, we found enrichment of motifs recognized by AP-1, p53, E2F, ETS, bHLH, and MYB ( Fig 1B and Table). In enhancers that lost their acetylation, we found enrichment of motifs from the bHLH family, ETS family, and MYB family (Fig 1C and S2 Table).
To identify the actual TFs that may bind the differentially acetylated enhancers, we obtained the binding sites of 27 DNA binding factors experimentally measured in A549 cells from ENCODE. Binding sites of 17 out of these factors were found to be significantly enriched in enhancers that lost acetylation after 24h, and a few of them were also found to be enriched in enhancers that lost or gained acetylation after 8h (Fisher exact test, FDR < 0.05, Table 1). The 17 enriched factors include factors related to transcriptional machinery and epigenetic regulation (EP300, EHMT2, SMC3, KDM1A, RAD21, TEAD4, CHD4, FOSL2, JUNB, JUN) and immune-related processes (CEBPB, GATA3, MAFK, MAX, NR3C1). Notably, for MAFK (acetylation loss after 24h), JUNB (acetylation gain after 8h), GATA3 (acetylation gain after 8h), and MAX (acetylation loss after 24h), the binding sites enrichment is also supported by a matching motif enrichment.

Differentially acetylated enhancers are adjacent to differentially expressed genes
To test the correlation of enhancer changes with gene expression and focus on clinically relevant changes, we obtained gene expression data of ex-vivo infected human lung tissues and human COVID-19 patients' lungs (see Methods). We looked for significantly differentially expressed genes near enhancers with significant acetylation change. For lungs infected ex-vivo, we found that 87 genes that were differentially expressed 24h after SARS-CoV-2 infection compared to mock control (FDR < 0.1, fold-change > 150%, Fig 2A) and that were located near enhancers with significant acetylation change (FDR < 0.05, fold-change > 150%). 80 were more highly expressed, and 7 genes had lower expression. We found that enhancers that lost acetylation tend to be near genes with significantly lower expression (Fisher exact test, p < 0.04, Fig 2B) and enhancers that gained acetylation tend to be near (< 100kb) genes that had a significantly higher expression (Fisher exact test, p < 0.004, Fig 2B).
When we looked at the distribution of acetylation levels around the transcription start site (TSS) of genes with significant expression change, we found that acetylation was indeed observed mostly at the promoter region upstream of the TSS, and importantly gain of acetylation at promoters was significantly higher in genes that gained expression (t-test, p < 0.03, Fig  2C and 2D).

Acetylation changes linked to cell adhesion and inflammatory response pathways
To understand which biological processes are affected by changes in enhancer activity, we performed gene set enrichment analysis for differentially expressed genes near differentially acetylated enhancers. We analyzed differentially expressed genes separately for ex-vivo infected lung tissues and COVID-19 patients' lung data. Several processes were consistently enriched between the COVID-19 patient data and the lung tissues data, including genes involved in the regulation of cell adhesion, positive regulation of ERK1 and ERK2 cascade, positive regulation of NF-kappaB signaling, and many pathways related to immune processes, including interferon signaling (Fig 3 and S3 and S4 Tables).

Association with Wnt signaling, MAPK signaling, NF-kappaB and adhesion
Since the gene set enrichment analysis indicated enrichment in Wnt signaling, MAPK signaling, NF-kappaB signaling, and cell adhesion, we further examined the differential genes (FDR < 0.1, fold-change > 150%) located near (< 100kb) differentially acetylated enhancers and involved in these pathways (S5 Table, Fig 4A). Due to a lower statistical power to detect differential expression in the ex-vivo data, we relied on the COVID-19 patient data for this analysis. Overall, we saw a loss of negative regulators and gained expression of positive regulators of adhesion, NF-kappaB, and Wnt signaling, while for MAPK signaling changes in expression were more complex and not in a clear direction (Table 2). For the Wnt signaling pathway we found upregulation of positive regulators, such as BCL9L, CDH3, CSNK1E, and SDC1, as well as downregulation of negative regulators, such as CAV1, PLPP3, and SFRP4 (S5 Table). For NF-kappaB signaling, we found upregulation of many positive regulators and of the TF subunits: APP, CASP10, EEF1D, EIF2AK2, IRAK3, MAVS, NPM1, RPS3, and TRIM22 as well as downregulation of negative regulators, such as CTNNB1, PELI1, and RORA (S5 Table). Again, we also found downregulation of positive regulators and subunits: AKAP13, CAV1, CD36, MAP3K14, PELI2, and REL and upregulation of negative regulators: TRIM38, TNIP1, and HSPB1. For MAPK signaling pathway we found upregulation of positive regulators and genes that are involved in the MAPK cascade: ADAM9, APP, EIF2AK2, GAS6, KITLG, GDF15, NOTCH1, and TRIM5, and downregulation of negative regulators: CAV1, CNKSR3, DUSP1, DUSP6, PDE8A, SPRED2, and PELI2. In a few cases changes in expression did not support MAPK activation when positive activators, such as AKAP13, PLCE1, RIPK1, and AREG, were downregulated and negative regulators like STK38 and TNIP1 were upregulated (S5 Table). Together this suggested that SARS-CoV-2

WT1 is enriched in H3K27ac peaks that lose acetylation after covid infection
To focus on SARS-CoV-2-specific effects, we repeated the analysis but removed genes that are also modulated by influenza infection (fold-change > 150%, FDR < 0.1), relying on the exvivo data, which provided matched SARS-CoV-2 and influenza control at 24 h. Only four genes passed this criterion: IER2, DLGAP1, WT1, and BLACAT1. All these genes were located near (< 100kb) enhancers that lose acetylation. WT1, DLGAP1, and BLACAT1 were also significantly downregulated at this time point (Fig 4B). Since WT1 (Wilms Tumor Suppressor) is a transcription factor that is downregulated during SARS-CoV-2 infection but not during influenza infection, we focused on the transcriptional effects of its downregulation. The WT1 gene encodes a zinc finger transcription factor and RNA-binding protein that directs the development of several organs and tissues. WT1 regulates numerous target genes that are involved in growth, differentiation, and apoptosis. It can serve both as a transcriptional activator and a suppressor. As mentioned above (Fig 1B and  1C), we found that the motif recognized by WT1 is enriched both at enhancers that lose acetylation and at enhancers that gain acetylation, suggesting WT1 indeed serves both as an activator and a repressor in this context as well.
We found 5 downregulated genes near WT1 enhancers that lost acetylation: DLGAP1, MAVS, EPHX1, PXMP4, and SPN. Finally, we found that WT1 enhancers that lost their acetylation tend to be found next to significantly downregulated genes (Fisher exact test, p < 0.0002, Fig 4E) and that WT1 enhancers that gained acetylation tend to be found next to upregulated genes (Fisher exact test, p < 10 −13 , Fig 4F).
To better understand the function of WT1 targets in this context, we performed gene set enrichment analysis for differentially expressed genes (SARS-CoV-2 vs. mock, fold-change 150%, FDR<0.1) near (< 100kb) differentially acetylated WT1 enhancers. We found enrichment of genes related to regulation of MAPK cascade, phosphorylation, cell adhesion, cytokine production, and other immune processes (Fig 5, S6 Table).
Together, this data suggested that WT1 was down-regulated after SARS-CoV-2 infections, which leads to major changes in the acetylation of WT1 bound enhancers, and the transcriptional regulation of numerous targets. This process was specific for SARS-CoV-2 and is not shared with other respiratory infections such as influenza, may be responsible for some of the unique features of COVID-19 such as excessive acute inflammation, and may support cell-tocell infection.

Discussion and conclusion
The COVID-19 pandemic affected the lives of most of humanity in recent years. Despite record breaking efforts to reveal the molecular impact of SARS-CoV-2 infection, our understanding of the impact on enhancer mediated transcriptional regulation remains poor. Given the crucial role of enhancer regulation and chromatin modifications in the regulation of transcriptional networks and cellular gene expression, our study aimed to elucidate the impact of SARS-CoV-2 infection on H3K27ac, a histone modification associated with active enhancers, and downstream gene expression changes.
Our findings revealed substantial changes in enhancer acetylation at 8 and 24 hours postinfection. Both gain and loss of acetylation were common and associated with differential expression of nearby genes. While H3K27ac data is currently limited to in vitro data due to technical difficulties, we demonstrated here its relevance to human biology by establishing a significant agreement with gene expression data in COVID-19 patients and ex-vivo infected lung samples. Furthermore, our analysis of motif and TF binding data at differential enhancers identified several transcription factors whose activity underwent significant alterations. While some of these effects are related to the general immune response of the infected cells, we identified several changes that are unique to SARS-CoV-2 and were not observed after influenza infection, most prominently downregulation of WT1 and its effect on WT1 bound enhancers and their target genes.
We discovered that differentially acetylated enhancers are adjacent to differentially expressed genes. Gene set enrichment analysis revealed changes in the inflammatory response,

PLOS COMPUTATIONAL BIOLOGY
SARS-CoV-2 infection perturbs regulation of key pathways cell adhesion, Wnt signaling, MAPK cascade, and NF-kappaB signaling. In particular, careful examination revealed the upregulation of positive regulators and the downregulation of negative regulators of Wnt signaling, NF-kappaB signaling, and cell adhesion, suggesting coordinated regulation to activate these pathways. While previous studies have already uncovered SARS-CoV-2 infection activates Wnt [8][9][10], MAPK [11][12][13], and NF-kappaB [11,14,15], our research unveils the involvement of enhancer mediated transcriptional networks in these processes for the first time. Furthermore, the activation of cell adhesion suggests that SARS-CoV-2 infection orchestrates transcriptional networks to support cell-to-cell transmission, a critical infection route for some SARS-CoV-2 strains, including those investigated in our study [16].
Importantly, our investigation demonstrated significant downregulation of WT1 during SARS-CoV-2 infection, and that this is not part of the general response observed in respiratory infections as it does not occur in influenza infections. WT1 is known as a negative regulator of MAPK signaling pathway [17], acting through the activation of SPRY1 [18] and DUSP6 (MKP3) [19]. Indeed, differentially expressed genes near WT1 enhancers were enriched in these processes, implying that the loss of WT1 contributes, at least in part, to the activation of the MAPK pathway. The novel finding of WT1's impact on increased cell adhesion warrants further investigation and may be a significant factor in facilitating cell-to-cell infection of SARS-CoV-2.
Our analysis provides valuable novel insights into the cell-autonomous transcriptional effects of SARS-CoV-2 infection, an important step toward a more comprehensive molecular understanding of the disease. Furthermore, our findings generate several hypotheses that merit further exploration and validation. This study showcases the potential of enhancer analysis in uncovering the regulatory effects of SARS-CoV-2 infection, as well as its applicability to other viral infections. Specifically, our analysis revealed the role of enhancer mediated transcriptional regulation after SARS-CoV-2 infection in manipulating WT1 activity, cell-to-cell infection, MAPK signaling, Wnt signaling and pro-inflammatory cytokines. Overall, our work demonstrated how integrative analysis of enhancer acetylation and gene expression after infection can reveal new pathways and transcriptional cascades that are perturbed by SARS-CoV-2 infection, moving us forward toward a deeper understanding of the molecular underpinning of COVID-19.
A union set of peaks for each comparison (ctrl vs. 8h and ctrl vs. 24h), was generated using bedtools merge [25]. Reads in each peak for each experimental condition were counted using bedtools multicov function. The obtained read counts were utilized as input for DESeq2 (default parameters) [26] to identify differentially acetylated peaks. A cutoff of FDR adjusted p-value < 0.05 was applied to identify significantly different H3K27ac peaks for further analysis. HOMER (version 4.11) was used to identify enriched motifs in the sequences of the differentially acetylated enhancers. The merged set of enhancer peaks was used as background to control for common motifs in A549 cells and focus on differential motifs.
analyses. In order to compare the enrichment level of the pathways we calculated the relative enrichment score by: (number of target genes in term * total number of genes)/(total number of target genes * number of genes in term), as suggested in [31].

WT1 enhancers identification
FIMO tool (version 5.5.0) [32] of MEME suite [33] was used to find enhancers that contain WT1 known motif (MCTCCCMCRCAB). Enhancers that exhibited the presence of the WT1 known motif with false discovery rate (FDR) < 0.1 were classified as WT1 enhancers.
Supporting information S1