1 SARS-CoV-2. In this study, we analyzed multiple transcriptomes obtained from cell lines infected with MERS-CoV, SARS-CoV and SARS-CoV-2, and from COVID-19 patient-derived samples. Using integrative analyses of gene co-expression networks and de-novo pathway enrichment, we characterize different gene modules and protein pathways enriched with Transcription Factors or Epifactors relevant for SARS-CoV-2 infection. We identified EP300, MOV10, RELA and TRIM25 as top candidates, and more than 60 additional proteins involved in the epigenetic response during viral infection that have therapeutic potential. Our results show that targeting the epigenetic machinery could be a feasible alternative to treat COVID-19.
The coronavirus family (CoV) are non-segmented, positive-sense and enveloped RNA viruses that have been identified as the cause of multiple enteric and respiratory diseases in both animals and humans 1 . By being intracellular pathogens, viruses' infection strategy requires the continuous subordination and exploitation of cellular transcriptional machinery and metabolism in order to ensure its own expansion. To do so, the host genome expression must be used and, to be successful, this will depend on chromatin dynamics and transcription regulation, which are principally ruled by epigenetic mechanisms, such as DNA methylation, histone post translational modifications (HPTM), and transcription factors (TF) 5 . During a viral infection, it has been reported that epigenetic and transcriptional changes occur for both sides: the infected cell promotes an antiviral environmental response, leading to the induction of pathways to survive, while the virus switches off the expression of critical anti-viral host cell genes 6,7 . 2 Several studies have reported the importance of epigenetic modifications in viral infections. In influenza virus, specific gene promoter methylation 8 , decreased H3K4me3 (a hallmark of active chromatin) 9 , histone acetylation in H3 and H4 histones, and increased levels of H4K20me2 and unmodified H3K36 and H4K79 have been reported 10 . Interestingly, these HPTM do not always trigger the same mechanisms and lead to similar phenotypes; for example, depletion of H3K79me2, an epigenetic mark that is usually increased upon viral infections due to an upregulation of DOT1L, results in impaired viral growth in human cytomegalovirus infection 11 , while enhancing the replication in influenza virus 10 . However, these mechanisms usually lead to host transcriptional inactivation, which contribute to the altered cellular transcription produced by viral infections.
Regarding CoVs, few experimental studies have been conducted to unravel the epigenetic proteins and marks involved in their infection and pathogenesis in MERS-CoV and SARS-CoV, being specially scarce in SARS-CoV-2 due to its recent appearance. For MERS-CoV and SARS-CoV, different outcomes have been reported, such as the mechanisms used to control the interferon-stimulated genes, which involves H3K27 methylation in MERS-CoV but not in SARS-CoV 12 , and the ones used to down-regulate antigen-presenting molecules, which involves DNA methylation in MERS-CoV and not in SARS-CoV 8 .
These studies show that epigenetic mechanisms are highly important in the host gene expression control carried out by the virus and that, despite the phylogenetic closeness, these mechanisms can be very different between strains, highlighting the need to understand the epigenetic processes that play a role in SARS-CoV-2 infection.
Integrative computational methods are promising approaches used to generate research hypotheses, generate consensus regulatory networks and to describe deregulated processes in SARS-CoV-2 infection 13, 14 . Nevertheless, they have overlooked key epigenetic and transcription factors that underlie the infected phenotype. Since drugs that target the epigenetic landscape of diseased cells have shown a great potential and have proved to be game-changing as complementary treatments of complex diseases, such as cancer 15 , the identification of these key epigenetic and transcription factors becomes highly important 3 in our current context, where popular regimen candidates for treating COVID-19, such as Remdesivir, Hydroxychloroquine, Lopinavir and Interferon have shown to have little or no effect on reducing mortality of hospitalized COVID-19 individuals 16 .
In this work, we gathered publicly available RNA-seq data from SARS-CoV-2, SARS-CoV and MERS-CoV infected cell lines and patient samples, and performeddifferential expression analyses together with weighted gene co-expression network analysis to identify unique and shared central epigenetic players in SARS-CoV-2, SARS-CoV and MERS-CoV. Candidate genes were further prioritized by integrating DEGs enrichment tests, gene-coexpression network and viral-host protein-protein interaction network analysis to propose potential key epigenetic proteins involved in SARS-CoV-2 infection. Finally, we identified currently approved drugs that target key epigenetic drivers of SARS-CoV-2 infection, and thus they are potential new therapeutic approaches for COVID-19.
Raw sequencing data was trimmed with Trimmomatic version 0.39 17 using the mapping parameters suggested in Jin et al. 20 : (--outFilterMultimapNmax 100 --winAnchorMultimapNmax 100). To estimate the abundance of the transcripts accounting for coding and non-coding genes as well as repetitive elements, we used TETranscripts version 2.1.4 20 with the multi mode. Raw count tables were used for differential expression analysis using DESeq2 21 . Differentially expressed genes (DEGs) were identified with a p adj. < 0.05 and abs(log 2 fold change) > log 2 (1.5).
4 Viral transcriptome was constructed with the 11 gene sequences reported in SARS-CoV-2 genome (NCBI Reference Sequence NC_045512.2). Viral transcript expression was quantified in each trimmed RNA-seq file of SARS-CoV-2 infected samples with Salmon v 1.3.0 22 .
Virus-associated gene sets were obtained with the intersection of DEGs identified in all the cell lines infected with the corresponding virus, except for SARS-CoV-2. For SARS-CoV, the intersection between the cell lines infected consisted of 182 genes (SARS-CoV-DEGs); for MERS-CoV the intersection was 1139 genes (MERS-CoV-DEGs); and for SARS-CoV-2, the intersection between at least 3 out of the 4 cell lines was used instead and consisted in 909 genes (SARS-CoV-2-DEGs) (Supplementary Table 2 ).
Patient-associated gene set was obtained with the DEGs in at least 2 out of the 3 conditions (543 genes, patient-DEGs) (Supplementary Table 2 ).
To build the Epigenes catalogue, 4 different databases were used: EpiFactors 23 , Histome 24 , dbEM 25 and the manually curated TF list from Lambert et al. 26 . TFs' functional annotation was taken from Lambert et al. 26 . The final list consisted of 2161 genes (776 epifactors, 1348 TFs and 41 categorized as both TF and epifactor).
Count matrices of the analyzed cell lines were filtered to remove low-expressed genes using the function filterByExpr from edgeR 27 while accounting for the treatment (i.e. virus infection) and cell type in the filtering design. Following, normalization of gene counts was performed with vst function from DESeq2 21 (treatment and cell type of each sample were included into the design matrix and accounted for these effects with the blind argument).The gene co-expression network was built with the log 2 fold changes of each biological sample compared with the controls of the same biological condition by applying the formula (1).
(1) log 2 FC i = log 2 ( SC i / ACC i ) 5
Where SC and ACC correspond to the normalized counts of gene i in the infected and controls samples respectively. The resulting matrix containing the log2FoldChanges per sample was used to construct the weighted gene co-expression network with the WGCNA package 28 . A soft threshold of 9 was used to construct the network and modules were identified with a minimum size of 20. Modules whose expression was similar were merged using a dissimilarity threshold of 0.25, resulting in a total of 24 modules. Finally, the module-eigengene pearson correlation of each module with the viruses was tested.
Gene Ontology enrichment analyses were performed using clusterProfiler 29 in virus associated and patient gene sets. For the differential expression analyses of infected cell lines, the enrichment of GO terms in DEGs was tested using the expressed genes on each particular comparison as background. For the co-expression network, the enrichment of GO terms was tested in each module using the genes of the full network as background.
Epigenes, virus-associated DEGs and TF-target enrichment analyses were performed with gProfiler2 30 using a custom gmt file or the TRANSFAC database included in the package for TF-target enrichment.
The correction method used was g:SCS and an adjusted p-value significance threshold of 0.05. As background, all the genes annotated in the co-expression network were used for epigenes and TF-target enrichment and the expressed genes in each virus for virus-associated DEG enrichment.
Virus-human interactions were obtained from Gordon et al. 31 and Stukalov et al. 32 . The human protein-protein interaction network (PPI) was obtained from IID version 2018-11 33 using only the experimentally validated interactions ("exp", "exp;ortho", "exp;ortho;pred" or "exp;pred"). After homogenizing the viral protein nomenclature, the three sources of interactions were merged to create the entire virus-human PPI, followed by removal of duplicated edges and self loops. (1) The shortest path length was calculated between all pairs of viral proteins and human proteins in the PPI network with the igraph package version 1.0.0 35 . The retained epigenes were the ones whose shortest path length with at least one viral protein was less than 3. K=2 for module 6, and K=3 for modules 8, 10 and 11.
All approved and non-approved drugs targeting the genes/proteins contained in each network were obtained with CoVex 37 by mapping the gene names to uniprot IDs, using the closeness centrality algorithm and the following parameters: result size= 50000, disabled hub penalty, disabled max degree, include indirect drugs=FALSE and include non-approved drugs=TRUE.
Raw RNA-seq data was obtained from the Sequence Read Archive ( www.ncbi.nlm.nih.gov/sra ) of the National Center for Biotechnology Table 1) .
In order to identify the genes that change their expression in pulmonary cell lines (Calu-3, MRC-5, A549
and NHBE) due to infection of Coronaviruses such as MERS-CoV, SARS-CoV or SARS-CoV-2, differential expression analysis was performed in RNA-seq data.
As a first approach, we evaluated the overlapping DEGs identified for each virus regardless of the cell Furthermore, the analysis showed that the differentially expressed epifactors belong to a wide range of functional categories, such as histone writers, histone readers, histone erases, Polycomb group proteins, chromatin remodeling, DNA modifications, among others ( Figure 1D ). In addition,regarding differentially expressed TFs, cell lines infected with SARS-CoV-2 show differential expression of TFs of the STAT (mediators of the cellular response to cytokine) and IRF (interferon-regulatory factor) family, which are not differentially expressed in MERS-CoV and SARS-CoV ( Figure 1E ). We noted that most of the TFs that are differentially expressed and shared between two or more of the Coronaviruses infected cells are 
Afterwards, we evaluated the transcriptional response in patient samples infected with SARS-CoV-2 to assess their resemblance to the previously observed results in cell lines. For this purpose, we obtained datasets from samples of bronchoalveolar lavage fluid (BALF), lung and peripheral blood mononuclear cells (PBMC). From the differential expression analysis, we found 16 DEGs shared among all three sample types (Figure 2A ), the majority of them with a different fold change direction in PBMC compared to LUNG and BALF, and 6 genes with a consistent fold change direction ( Figure 2B ). This is most likely explained by the higher cell-lineage similarity between the Lung and BALF samples than with PBMC, 10 which could lead to a more similar transcriptional response to the infection. In addition, GO enrichment analysis shows that most of the DEGs were related to the immune response to viral infection such as leukocyte migration and humoral immune response ( Figure 2C ). Following, we evaluated the similarity of these results with the data observed for SARS-CoV-2 in infected cell lines by comparing the overlap between the virus-associated genes with the DEGs present in at least 2 out of the 3 sample types in patients (hereinafter referred as patient-DEGs, 543 genes). We found 61, 15 and 33 genes in common with SARS-CoV-2, MERS-CoV and SARS-CoV, respectively. In particular for SARS-CoV-2 and patients, 6 TFs ( STAT5A, MAFF, IRF9, MXD1, FOSB and STAT4) and no epifactors were identified. We also found, at a lesser extent, shared DEG between samples of SARS-CoV-2 infected patients and MERS-CoV and SARS-CoV infected cell lines, which are likely to be non-specific viral-responding immune genes ( Figure 2D ). Finally, we found that 804, 1 and 20 repeat elements are differentially expressed in BALF, PBMC and LUNG samples, respectively, being LTR elements the most differentially expressed in all samples (Supplementary Table 3 ). Collectively, these results show that the gene expression changes promoted by SARS-CoV-2 infection in patients are diverse among tissue types;
however, despite these differences, immune response processes are the main processes affected across all sample types.
the same gene co-expression modules, which recapitulate the expression profiles in COVID-19
So far, our analyses showed that cell lines and patient samples infected with SARS-CoV-2 exhibited DEGs related to immunological processes, which has been previously described by Blanco-Melo et al. 42 which is congruent with our results. However, differential expression analysis often overlooks the subtle differences in several genes that altogether can be responsible for major changes in global transcriptional regulation. Weighted gene co-expression network analysis overcomes this limitation by studying the expression of thousands of genes in the same analysis 28 . Thus, we expanded our previous results by including a co-expression analysis to identify gene modules associated with each viral infection and genes that play central roles within them.
We constructed the co-expression network with the log2 fold changes of each sample compared to its controls. After identifying the modules, we calculated the correlation between each module and the different traits, where we found that out of the 24 total modules identified, 13 Regarding SARS-CoV-2, module 9 was uniquely associated with it, but also other modules standed out. such as the epifactors MOV10 and EP300, the TF RELA and the enzyme TRIM25, since they are exception connectors (genes that do not belong to the specific module, but are important in the protein 15 pathway found) in more than one module (modules 6, 10 and 11), and belong to modules enriched in genes that participate in histone H3-K4 methylation and in the response to interferon gamma. EP300 is a histone acetyltransferase that was also identified in SARS-CoV-2 infected cell lines 13 . Additionally, we found that MOV10, a putative helicase, also participates in SARS-CoV-2 infection. The TF RELA has been increasingly recognized as a crucial modulator of the response to SARS-CoV-2 infection 13, 45 and is part of the NF-κB complex, along with RELB 46 , which is differentially expressed in MERS-CoV and SARS-CoV-2 infected cell lines (Supplementary Figure 1D ). TRIM25 is a ubiquitin ligase required for production of INF-1 and is inhibited by Nucleocapsid of SARS-CoV 47 . Finally TRIM28 (also known as KAP1) has been shown to interfere with viral integration into host genome 48 and represses the expression of repeat elements of the LINE family, in particular L1NA4 41 , which was previously identified as differentially expressed in cell lines infected with the three Coronaviruses.
Afterwards, we evaluated whether the proteins in the networks had annotated drugs targeting them. We Table 1 .
Acid, whose derivatives show anti-HIV activity 49, 50 . EP300 is targeted by curcumin, a molecule with anti-inflammatory properties 51 . The latter proposes RELA and EP300 as new potential drug target candidates for SARS-CoV-2 infection, not only because they participate in immune-related processes, but also because they belong to the cellular epigenetic machinery used by the virus during infection. We further identify at the protein interaction level, that several TFs take part mainly in the immunological response to viral infection. One example is NF-κB, whose p65 subunit (also known as RELA ) is a central part in the protein interaction network for SARS-CoV-2. NF-κB induces the expression of several pro-inflammatory cytokines, including IL-6, CCL2 and CCL3 53 , which had been found in high levels in COVID19 patients 54 . On the other hand, TRIM25, an ubiquitinase, is essential for the activation of NF-κB and the production of IL-6 55 . TRIM25 is over-expressed in cell lines infected with SARS-CoV-2 but not in those with MERS-CoV, which furthermore suggest that NF-κB could be a medullary part of the host immune response against SARS-CoV-2. The previous observation is reinforced by the fact that it was observed that RELA directly interacts with histone acetyltransferase EP300, and both proteins interact with various components of the AP-1 complex such as FOS, JUND, and FOSL1. AP-1, EP300 and NF-κB regulate chromatin accessibility in the proximal promoter region of IL-6 and CCL2, both pro-inflammatory cytokines 56, 57 . The p300/CBP complex is one of the best characterized cofactors of NF-κB and specifically binds RELA and acetylates it along with the surrounding histones 46 . It is known that adults 18 older than 65 years have higher NF-κB levels compared to younger adults 58 and some authors had suggested that this may be one reason older adults are more susceptible to develop the severe form of COVID-19 59 .
According to our results, SARS-CoV-2 infection modifies the expression of several TFs of the interferon regulatory factor (IRF) and STAT families, which are primarily involved in the immune response against pathogens. STAT1 and STAT2 are key elements of the signaling induced by type I interferons, this proteins form a dimer upon interferon mediated phosphorylation and, together with IRF9, form the complex ISGF3 that activates the transcription of interferon stimulated genes 60 . Our results also showed that IRF9 is upregulated in cell lines infected with SARS-CoV-2; however, module 12's interactome showed that STAT2 and STAT1 interact with IRF2. IRF2 is a negative regulator of IFNα and its inhibition causes an increase in the antiviral response induced by IFNα 61 . This fact further suggests an impairment of interferon type I stimulated genes activation, as previously described as a hallmark of SAR-CoV-2 infection 62 . On the other hand, IRF1 and IRF7 were also upregulated in SARS-CoV-2 infected cell lines.
IRF7 is a key TF for IFNα expression, and it has been previously identified as a hub gene for SARS-CoV-2 infection together with IFR9 and STAT1 63 . It is also interesting that IRF7 lost of function mutations were associated with severe COVID-19 patients 64 and with the development of life-threatening influenza in children 65 which suggest that inhibition of IRF7 activity is crucial for SARS-CoV-2 pathology.
Viruses have been reported to use epigenetic machinery to take advantage of the cell and hijack its regulatory capacity for their own benefit 12 . The epigenetic machinery can be affected by coronaviruses in this same sense, and this can happen either by promoting alterations in the epigenetic code, such as DNA methylation and post-translational modifications of histones, or directly by promoting the dysregulation of enzymes and other proteins associated with the epigenome.
We found that among the deregulated epifactors with histone acetylation function are HDAC9 and SIRT1 enzymes. In this sense, it has recently been reported that the SIRT1 protein (a class 3 HDAC) was positively regulated in the lung of patients with severe COVID-19 comorbidities 66 . Likewise, another work demonstrated that under conditions of cellular energy stress, SIRT1 can epigenetically regulate the ACE2 receptor 67 . Also, it has been observed that treatment with non-steroidal anti-inflammatory drugs can inhibit SIRT1 activity, which in turn could affect ACE2 expression 68 . Accordingly, it has been postulated that some diseases such as lupus where the epigenetic dysregulation is implicit in this disease could facilitate the entrance of SARS-CoV-2 into the host cells 69 .
Interestingly, the enzymes HAT1, HDAC2 and KDM5B have been reported to also potentially regulate ACE2 in human lungs. KDM5B has gained interest, because it is associated with other viral infections such as the hepatitis B virus 70 , and potentially with SARS-CoV-2 66 . Remarkably, in breast cancer cells, it has been shown that inhibition of this enzyme triggers a robust interferon response that results in resistance to infection by DNA and RNA viruses 71 . In this regard, we observed several deregulated KDMs in the different coronavirus infections, in which KDM6B stands out by being deregulated in both MERS-CoV and SARS-CoV2 infection. KDM6B is a specific demethylase of H3K27me3, which acts as a repressive histone mark. Although it remains to be fully studied, it is associated with the regulation of a wide range of genes involved in inflammatory agents, development, cancer, viral infection response, senescence and is an important host response against environmental, cellular stress 72 . Therefore, adding to the above, it is suggested that demethylases, such as KDM6B, are potential epigenes that are affected during SARS-CoV-2 infection and can be presented as potential targets for the treatment of COVID-19.
However, this should be further studied.
Several epigenes previously involved in response to viral infections stood out in our protein interaction analysis, such as BRD4, TOP2A, and TRIM28. Bromodomain protein 4 (BRD4) is a histone acetylation reader and writer that plays an important role in DNA replication, transcription, and repair 73 . This epigene is critical for the maintenance of the higher-order chromatin structure, since its inhibition leads to chromatin decondensation and fragmentation, and it also can stimulate innate antiviral immunity 73 . BRD4
complexes with RelA and CDK9 and is functionally required for effective activation of NF-kB-dependent immediate early cytokine genes in response to viral patterns. In this sense, our results show a protein interaction with EP300, which involves the p300 / CBP complex, which is one of the best characterized cofactors of NF-kB and binds specifically to RelA 74 . Validating the possible importance of this system in 20 infection with SARS-CoV2. Examples like this suggest that the virus, through these epigenetic remodelers, is promoting that the cell induces a chromatin remodeling that in many cases will lead to opening, both at the local level of specific genes and at the global level. Accordingly, an indicator of global changes is the increased expression of transcripts from repeated sequences such as LINE1. If this is so, then the virus is manipulating the chromatin aperture to promote the expression of genes that support its invasion. In this regard, other work has suggested the importance of LINE1 elements. Where this type of repetitive elements are very relevant in gene regulation, especially when these elements are in proximity to neighboring genes, since they could alter their expression. Therefore, the dysregulation of repeated elements such as LINE1 could indirectly change the cellular transcriptome 75 .
Furthermore, we find epigenes that interact with the viral proteins directly or very closely. This connection suggests a virus-promoted modulation to affect the epigenome of the host cell's interactome. Which reinforces the idea that the virus strategy is partly to take advantage of the epigenetic machinery. In general, our data suggest that the SARS-CoV-2 infection deregulates the epigenetic master machinery.
One of the points that should be taken into consideration in the future is that if this epigenetic machinery is not re-established after disease courses it could generate other diseases such as cancer in the long term. This is based on the fact that many of the genes that we found in our study have been proposed as epigenetic hallmarks in various neoplasms.
Our last key finding is the identification of driver epigenetic proteins and TFs involved in SARS-CoV-2 infection that can be targeted by existing drugs. We identified S-adenosyl-L-homocysteine (SAH)
targeting several epigenetic components of the host response to SARS-CoV-2 infection. SAH is the product of the chemical reaction performed by methyltransferases using nucleic acids or proteins as substrates, and has been previously suggested as a potential treatment for viral infections such as ZIKA, [84] [85] [86] . It is worth mentioning that these drugs are being tested in clinical trials based on their described inhibitory activity of enzymes related to the activation of immune response and inflammation, such as growth receptors 87 . The latter, together with our results, suggests that drugs targeting epigenetic mechanisms could be also effective to treat SARS-CoV-2 by modulating their kinase activity.
Finally, we also identified Bortezomib and betulinic acid associated with RELA. Bortezomib is a proteasome inhibitor that has been proposed as COVID-19 therapy given its capacity to inhibit (although only marginally) the papain-like protease (NSP3) of SARS-CoV, which also has deubiquitinase activity [88] [89] [90] . Likewise, betulinic acid has been proposed as a target of NSP3 in SARS-CoV-2 91 .
Together, we have supporting evidence that current drug-based therapies to treat COVID-19 also target the transcriptional response to infection by the modulation of the epigenetic proteins identified in this study. Furthermore, we provide additional new potential drug targets and drug candidates which could be effective and whose potential use has not been exploited yet ( Figure 5 ). These results provide comprehensive evidence that epigenetic therapy could aid in restoring the transcriptional changes 22 observed during infection. By using epigenetic drugs, a therapeutic effect can be achieved due to their systemic effects, which can be advantageous to treat a disease that targets different tissues and cellular mechanisms, as observed in COVID-19.
In this study, we used a blend of bioinformatic approaches to comparatively analyze transcriptomic data from SARS-CoV-2, SARS-CoV and MERS-CoV pulmonar infected cell lines and COVID-19 patient-derived samples.In particular, we focused on the epigenetic processes and transcriptional factors, since these have been widely proposed as the master regulators of the expression of most genes. We found that the transcriptional response elicited by MERS-CoV and SARS-CoV-2 is more similar to that observed for SARS-CoV, which is observed both in infected cell lines and in samples from patients with COVID-19. At the same time, we identified specific altered modules in the response to infection with SARS-CoV2 that could serve as an aliasing for the proposal of different therapeutic strategies based on epigenetic therapy. Thus, our results add a piece to the puzzle of the strategies used by the different coronaviruses to manipulate the gene regulation capacity of the cell. Although the pathways are differential between them, the virus objective is to take advantage of the TFs and various chromatin remodelers to avoid being detected and prevail in the invasion. This is a very fine strategy that the virus uses and it has been poorly studied in both its biological importance and its future therapeutic application.
This could open a new window of opportunities for treatment and thus close the chapter on this pandemic disease. 