SARS-CoV-2 Nsp14 mediates the effects of viral infection on the host cell transcriptome

Viral infection involves complex set of events orchestrated by multiple viral proteins. To identify functions of SARS-CoV-2 proteins, we performed transcriptomic analyses of cells expressing individual viral proteins. Expression of Nsp14, a protein involved in viral RNA replication, provoked a dramatic remodeling of the transcriptome that strongly resembled that observed following SARS-CoV-2 infection. Moreover, Nsp14 expression altered the splicing of more than 1,000 genes and resulted in a dramatic increase in the number of circRNAs, which are linked to innate immunity. These effects were independent of the Nsp14 exonuclease activity and required the N7-guanine-methyltransferase domain of the protein. Activation of the NFkB pathway and increased expression of CXCL8 occurred early upon Nsp14 expression. We identified IMPDH2, which catalyzes the rate-limiting step of guanine nucleotides biosynthesis, as a key mediator of these effects. Nsp14 expression caused an increase in GTP cellular levels, and the effect of Nsp14 was strongly decreased in presence of IMPDH2 inhibitors. Together, our data demonstrate an unknown role for Nsp14 with implications for therapy.


INTRODUCTION
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the virus responsible for the COVID-19 pandemic that began in 2019. As of early February, 2022, COVID-19 has caused 5.7 million deaths worldwide. Coronaviruses are enveloped, relatively small (60-140 nm diameter), positive-stranded RNA viruses belonging to the Coronaviridae family. They derive their name from the crown-like appearance (corona means crown in Latin) that results from the spike glycoproteins in their envelope (V'kovski et al., 2021a). The SARS-CoV-2 RNA genome is 30-kb long and has 14 open reading frames (ORFs) that encode 29 proteins (16 non-structural proteins, 4 structural proteins, and 9 accessory factors, Kim et al., 2020). During the first step of the viral infection, the spike glycoprotein mediates attachment and fusion with the cellular membrane. Once the virus is in the cytoplasm of the host cells, the host ribosome machinery is recruited for the synthesis of viral proteins. Non-structural proteins are required for viral genome replication and are generated by proteolytic cleavage of the polyprotein encoded by ORF1a and ORF1ab. Once the viral genome is replicated, virions are assembled in the host endoplasmic reticulum-Golgi intermediate complex. Finally, new viral particles are incorporated into vesicles and secreted by the host cells (V'kovski et al., 2021a). Viral infection triggers a variety of pathways in the host cells that ultimately lead to the hijacking of the cellular machineries and escape from immune surveillance. SARS-CoV-2 infection elicits a peculiar gene expression response that first involves activation of interferon pathway (Blanco-Melo et al., 2020;Vanderheiden et al., 2020;Wyler et al., 2021), and then the NFkB pathway (Kircheis et al., 2020;Hariharan et al., 2021;Wyler et al., 2021), as well as expression of specific cytokines such as IL6 and IL8 (Wang et al., 2007;Blanco-Melo et al., 2020;Coperchini et al., 2020;Park and Lee, 2020). 4 Recent and intense research on SARS-CoV-2 characterized the role of viral proteins during viral replication, and showed that these functions are often conserved across coronaviruses (V'kovski et al., 2021a). Less is known about the roles of the individual proteins in modulating host cell pathways (Gordon, Hiatt, et al., 2020;Gordon, Jang, et al., 2020). For instance, recent studies have proposed that Nsp16 is a splicing modulator (Banerjee et al., 2020) and that Nsp1 and Nsp14 are translational repressors (Schubert et al., 2020;Hsu et al., 2021).
Nsp14 is a 60-kDa protein conserved among coronaviruses that is involved both in viral replication and immune surveillance escape (Ogando et al., 2020). The N-terminal region of Nsp14 contains an exonuclease domain (ExoN) that excises mismatched nucleotides to ensure accurate replication of the viral genome (Ogando et al., 2020). As a result of this proofreading mechanism, coronaviruses have a lower mutation rate than other RNA viruses (error rate of 10 6 -10 7 versus 10 3 -10 5 ) (Sanjuán et al., 2010;Robson et al., 2020). Loss of function of the ExoN domain results in increased sensitivity to the RNA mutagen 5-fluorouracil (Eckerle et al., 2010) and attenuated virulence (Graham et al., 2012). Furthermore, the interaction with Nsp10 augments Nsp14 ExoN activity up to 35 fold, and inhibition of the interaction between Nsp10 and Nsp14 leads to reduced replication fidelity (Ma et al., 2015a;Smith et al., 2015). Nsp14 is also involved in assembly the cap at the 5' end of the viral RNA genome, which is crucial for evading immune surveillance. The C-terminal region of Nsp14 functions as an S-adenosyl methionine-dependent guanine-N7 methyl transferase that is independent of the ExoN activity (Chen et al., 2009). Both enzymatic domains are essential for successful viral replication, making Nsp14 an appealing drug target (Otava et al., 2021;Saramago et al., 2021).
Nsp14 is part of the Replication Complex, therefore it interacts with other SARS-CoV-2 proteins. As for other coronaviruses, replication of SARS-CoV-2 genome takes place in replication 5 organelles that provide a protective environment for the newly synthesized viral genome (V'kovski et al., 2021b). Notably, these organelles are formed in the cytoplasm and present convoluted double layered membranes that likely exchange material with the cytoplasm through pores (Wolff et al., 2020).
Guanosine-5'-triphosphate (GTP) is necessary for DNA replication and transcription and is used as energy source for translation and as mediator of signal transduction (Hesketh and Oliver, 2019).
How the physical interaction between Nsp14 and IMPDH2 impacts the host pathways is not completely understood (Li et al., 2021).
Interestingly, a recent study has shown that expression of Nsp14 results in global translation inhibition, but it is not clear if this is a direct effect or a downstream consequence of the changes that the expression of this protein provokes to the cellular environment (Hsu et al., 2021). Indeed, no direct interaction between Nsp14 and ribosomes or a known translational modulator has been reported, suggesting that the potential translational inhibition might be a downstream effect. 6 Here we undertook transcriptome analyses in cells that express each SARS-CoV-2 protein individually. Nsp14 altered the expression of about 4,000 genes, mostly involved in splicing, RNA metabolism, and cell-cycle control. Importantly, the effect of Nsp14 on cellular gene expression resembled the transcriptional changes that occur upon SARS-CoV-2 infection and included the activation of the NFkB pathway and the expression of CXCL8 (encoding IL8), a marker of acute severe respiratory distress syndrome in COVID-19 patients (Adcock et al., 2015;Blanco-Melo et al., 2020;Kircheis et al., 2020). Intriguingly, we also detected an increase in circRNAs expression upon Nsp14 expression; recent studies indicate that circRNAs can act as modulators of the innate immune response during viral infections (Li et al., 2017;Liu et al., 2019;Chen et al., 2020;Yan and Chen, 2020). Moreover, we showed that the cellular enzyme IMPDH2 mediates the gene expression response induced by Nsp14. We found that IMPDH2 mRNA is downregulated upon Nsp14 expression and that the cellular GTP concentration strongly increases, indicating that Nsp14 might activate IMPDH2 enzymatic activity. In accordance with our hypothesis, we showed that treatment with IMPDH2 inhibitors (Mycophenolic acid and Mizoribine, Seungheon  partially rescued the changes in gene expression induced by Nsp14.

Expression of individual SARS-CoV-2 proteins specifically remodels the transcriptome
Infection of cells with SARS-CoV-2 induces strong and specific changes in the transcriptome of host cells and tissues (Blanco-Melo et al., 2020;Wyler et al., 2021). It is assumed that this response results from the hijacking of the cellular systems by the virus as well as from the defense by the host. To identify unknown functions of the individual SARS-CoV-2 proteins and to determine how much each protein contributes to the takeover of cellular systems, we determined how the transcriptome changed when we individually express each viral protein in a human cell line. Specifically, we expressed individual SARS-CoV-2 proteins (Gordon, Jang, et al., 2020) in HEK293T cells, and after 48 hours we purified RNA, generated and sequenced 3' RNA-seq libraries, and identified differentially expressed genes (DEGs; Figure 1A). Extensive cell death or changes in morphology did not occur upon expression of individual proteins, as assessed visually.
Interestingly, expression of most proteins resulted in modest or no changes in the transcriptome of the HEK293T cells. Indeed, we detected less than 300 DEGs upon expression of 17 of the 25 tested proteins (Table 1 and Figure 1B). Expression of seven viral proteins, M, Nsp9, E, ORF9b, ORF3a, Nsp13, and Nsp1, modestly altered the transcriptome (between 300 and 700 DEGs). Interestingly, these DEGs tended to be upregulated rather than downregulated ( Figure 1B).
Striking, Nsp14 altered the expression of more than 4,000 RNAs (1,862 upregulated and 2,161 downregulated; Figure 1B). The profound impact of Nsp14 expression on the transcriptome of HEK293T cells suggests that this protein has roles beyond its known functions in viral genome proofreading and immune system escape.
To identify cellular pathways affected by the expression of the individual SARS-CoV-2 proteins, we performed Gene Ontology (GO) analysis of the DEGs upon expression of the different viral proteins (Table 2 and Figure 1C). Expression of the viral proteins impacted mRNAs encoding 8 proteins related to different aspects of gene expression (regulation of transcription, translation, and RNA metabolism), cell metabolism, cell division, and innate immunity. For example, expression of Nsp7, an accessory protein of the RNA-dependent RNA polymerase (Kirchdoerfer and Ward, 2019;Hillen et al., 2020), deregulated expression of genes encoding proteins involved in DNA and RNA metabolism ( Figure 1C). Expression of Nsp9 and Nsp13, known to bind RNA and to regulate RNA metabolism, respectively (Egloff et al., 2004;Shu et al., 2020), resulted in the misregulation of mRNAs encoding proteins involved in RNA metabolism and translation ( Figure 1C). Furthermore, expression of the protein ORF3b, known to be an interferon inhibitor (Konno et al., 2020), resulted in alteration of genes involved in cell metabolism and immunity. Indeed, the interferon response is known to impact several metabolic pathways, including carbohydrate catabolism (reviewed in Konno et al., 2020), and genes involved in these pathways are downregulated upon ORF3b expression ( Figure 1C). In general, expression of the viral proteins mainly impacted pathways related to RNA metabolism, translation, and cell metabolism that were previously reported to be affected upon viral infection (Blanco-Melo et al., 2020;Gardinassi et al., 2020;Wyler et al., 2021). Nsp14 affected the expression of genes involved in RNA splicing, metabolism, and processing, translation, cell-cycle control, and the cytoskeleton. Nsp14 also impacted the expression of genes involved in general metabolism, especially on genes implicated in nucleotide metabolism ( Figure 1C).

Nsp14 expression induces transcriptional changes that resemble SARS-CoV-2 infection
We decided to focus on Nsp14, as expression of this protein resulted in significantly more DEGs than any other tested protein. To obtain a more comprehensive view of the transcriptome changes provoked by Nsp14, we complemented the 3' RNA-seq data with total RNA-seq data 9 generated from cells transfected with Nsp14. The new dataset strongly resembled the one obtained by 3' RNA-seq ( Figure 2A, Table 3, and Figure 2 -figure supplement 1A) but also contained information about the full-length RNAs, non-poly-adenylated RNAs, and pre-mRNAs.
To gain further insights into the mechanism and pathways altered by Nsp14, we performed Gene Set Enrichment Analysis (GSEA) of the total RNA-seq data from Nsp14-expressing cells.
Strikingly, we found that for genes upregulated upon expression of Nsp14, three of the four most highly enriched GSEA datasets were those of cells infected with SARS-CoV-2 (Blanco-Melo et al., 2020) ( Figure 2B, 2C, and Figure 2 -figure supplement 1B). Moreover, we observed a smaller but significant enrichment of genes upregulated upon infection with MERS, a related virus (Blanco-Melo et al., 2020). Further, the mRNAs downregulated upon expression of Nsp14 tended to be strongly enriched for those downregulated following SARS-CoV-2 viral infection ( Figure   2B). We did not detect enrichment when our dataset were compared to data from cells infected with influenza A virus (IAV) (Blanco-Melo et al., 2020) ( Figure 2B). These results demonstrate that expression of a single viral protein can recapitulate a considerable portion of response of a host cell to SARS-CoV-2 infection and highlights the potential importance of Nsp14 in hijacking host gene expression.

Nsp14 alters gene expression mostly at the transcriptional level
To evaluate whether the transcriptome changes were due to alterations at the transcriptional or post-transcriptional level or a combination of both, we utilized intronic signals from the total RNA-seq experiment as a proxy for transcriptional activity (Stuart . For each DEG in the 3' RNA-seq dataset, we determined whether the intronic signal in the total RNA-seq dataset was altered in the same direction. We observed a strong correlation (R = 0.73 and p-value = 2.2e −16 ) 10 between these two measures ( Figure 2D), indicating that a large part of the response is transcriptional. For example, of the 1,862 genes that were upregulated upon Nsp14 expression, 1,006 also had higher intronic signal (at least 20% signal increase), whereas only 57 of the showed decreased intronic signal. Similarly, of the 2,161 genes downregulated upon Nsp14 expression, 1,242 displayed lower intronic signal (at least 20% signal decrease), whereas only 109 showed higher intronic signal (Table 4). We then performed a similar analysis using the exonic signal of DEGs in the same total RNA-seq dataset. We obtained an even stronger correlation between the changes in exonic signal and total RNA signal for each gene (R=0.87 and p-value= 2.2e −16 , Figure   2E). In sum, we observed that more than 50% of the changes in the transcriptome upon Nsp14 expression are transcriptional. Notably, this might be an underestimation, as changes in intronic signal could be lower than the steady-state RNA levels or some introns could be very efficiently and quickly removed, so that splicing intermediates are not detected. To further confirm the transcriptional effect, we determined the levels of a subset of DE genes from chromatin-bound nascent RNA from cells transfected with a control or a Nsp14-expressing plasmid. As expected a pre-mRNA (pre-TBP), was enriched in the chromatin bound fraction, U6 was abundant in the nuclear compartment (nucleoplasm and chromatin bound), whereas 18S rRNA and a circRNA and CXCL8) were also upregulated in the chromatin bound fraction, while the downregulated ones (SH2D2A and COL13A) were downregulated indicating that at least those mRNAs are regulated at the transcriptional level by Nsp14. These results verify the genomic observations and strongly suggest that the gene expression changes observed upon Nsp14 are mainly transcriptional ( Figure   11 2F). Interestingly, we also observed that there are more than 1,000 genes that display higher intronic signals with no changes in gene expression, indicating that, in addition to the transcriptional effects, Nsp14 might also affect splicing (Table 5).

Nsp14 expression provokes changes in alternative splicing and circRNAs production
Genes upregulated upon Nsp14 expression are enriched in genes encoding proteins with GO terms related to RNA metabolism and, more specifically, splicing ( Figure 1C). Moreover, for more than 1,000 genes, there were significant increases in intron signal upon expression of Nsp14 without changes in mRNA levels, suggesting that Nsp14 alters the splicing of these pre-mRNAs (Table 5). Indeed, we found that expression of Nsp14 strongly altered the inclusion patterns of almost 2,300 exons, with more than 2,000 exons displaying lower inclusion in the mature mRNA and 238 showing higher levels of inclusion when we expressed Nsp14 ( Figure 3A, Table 5, and . Furthermore, we also identified genes which used alternative acceptor or donor splice sites upon expression of Nsp14 ( Figure 3A and Table 5). Moreover, we observed an increase in the retention of more than 2,000 introns following Nsp14 expression ( Figure 3A and Figure 3C). Although the effects on exon inclusion and use of alternative splice sites clearly indicate that Nsp14 influences splicing, the increase in intron retention could be secondary to changes in transcription. However, as most of the introns with increased retention are within genes that were not differentially expressed ( Figure 3B and Figure 3 -figure supplement 1B), we reasoned that expression of Nsp14 leads to changes in splicing in this subset of genes as well. Moreover, the effect of Nsp14 on alternative splicing appears to be specific to particular introns and exons in each gene, as the majority (~62%) of identified genes had a single spicing event altered ( Figure 3D and Table 5 To further characterize the effects of Nsp14 expression on alternative splicing, we looked for genomic features associated with the altered patterns of splicing in the presence of Nsp14. Interestingly, we found that the affected introns tended to be at the 5' end of the transcript and were embedded in genomic regions with higher-than-average GC content (Figure 3 -figure supplement 1C and Table 6). The GC content may be a confounding factor of the location of these introns, as the GC content is higher at the 5' end of genes and around the transcription start sites (Zhang et al., 2004). Therefore, we concluded that Nsp14 expression has a strong and specific effect on splicing efficiency for a subset of genes.
Furthermore, most of the affected exons are shorter and have higher GC content around their splice sites than a randomized subset of exons not affected by Nsp14 (Figure 3 -figure supplement 1D and Table 6). The higher GC content suggests that stable RNA structures form around the splice sites in these exon, which correlates with alternative splicing propensity (Zhang, Kuo and Chen, 2011;Lin, Taggart and Fairbrother, 2016).
Previous research showed that SARS-CoV-2 infection disrupts splicing, mostly by inducing intron retention (Banerjee et al., 2020). To check whether the effects on splicing induced by Nsp14 were comparable to the ones that happen during SARS-CoV-2 infection, we reanalyzed a published dataset of total RNA sequencing from HEK293T cells expressing human ACE2 and infected with SARS-CoV-2 (Sun et al., 2021). We found that there is a significant overlap (p-value <10 -15 , see Material and Methods) between the alternative splicing events during the infection and in our model. Specifically, 10% of the alternative splicing events changing with SARS-CoV-2 infection (Table 7) also change (and in the same direction) when we express Nsp14 13 (Table 5). In sum, we showed that expression of Nsp14 partially recapitulates the gene expression changes, as well as some of the alternative splicing events that occur upon SARS-CoV-2 infection.
As circRNAs are generated by back-splicing, a process that competes with linear pre-mRNA splicing (Ashwal-Fluss et al., 2014), we checked whether Nsp14 also alters circRNAs expression. Using the total RNA-seq data, we found that there is a strong increase (more than 2 fold) in the total number of circRNA reads upon Nsp14 expression ( Figure 3E). Most of the 246 circRNAs that were differentially expressed upon expression of Nsp14 were upregulated (put the number of upregulated; Figure 3F and Table 8). These deregulated circRNAs were contained within 190 loci (some loci host multiple circRNAs). Interestingly, the levels of most of the mRNAs that host these circRNAs were unchanged ( Figure 3G and Table 4). Indeed, Nsp14 did not increase the transcription of loci hosting the upregulated circRNAs (we observed increased levels of intronic sequences in only 20 out of the 190 loci with upregulated circRNAs; Figure 3H and Table   4). Together, these results strongly suggest that the increased circRNA levels are caused by increased biosynthesis or stability of the circRNAs. In sum, our data shows that expression of Nsp14 influences alternative splicing and back-splicing on the host cells.

The effect of Nsp14 on gene expression is independent of the exonuclease activity but requires the N7-guanine-methyltransferase domain
Nsp14 has two separated enzymatic activities: (1) it works as a proofreading exonuclease, and (2) it is a N7-guanine-methyltransferase required for the modification of the viral RNA cap.
To determine whether the exonuclease activity is required for the strong effects of Nsp14 on gene expression, we performed two different experiments. First, we tested whether co-expression of Nsp10, which dramatically increases the exonuclease activity of Nsp14 (Bouvet et al., 2012(Bouvet et al., , 2014 14 Ma et al., 2015a), resulted in an enhanced effect of Nsp14. Briefly, we transfected HEK293T cells with plasmids for expressing Nsp14 and/or Nsp10. As a control, we used cells transfected with an empty vector. Co-expression of Nsp14 and Nsp10 did not impact levels of Nsp14 protein ( Figure   4A). We then prepared RNA from these samples and quantified the expression of some mRNAs that were altered upon expression of Nsp14 alone. For instance, the levels of FGF-18 increased upon expression of Nsp14 alone, whereas the levels of SH2D2A were significantly downregulated Moreover, co-expression of Nsp10 did not result in additional differentially expressed mRNAs 15 ( Figure 4B and Figure 4C) and did not alter of the number of circRNAs reads ( Figure 4D) compared to expression of Nsp14 alone. These data strongly suggests that the Nsp14 exonuclease activity is not responsible for the dramatic effects of Nsp14 expression on the transcriptome.
To definitely test whether the exonuclease activity is related to the remodeling of gene inhibit the exonuclease activity of Nsp14 (Ma et al., 2015a;Ogando et al., 2020). We generated and sequenced 3' RNA-seq libraries from cells that expressed wild-type or ExoN-mutant Nsp14 proteins or eGFP, as a control. PCA indicated that the two mutants have overlapping transcriptomes that differed from the control (  Figure 4E). Moreover, we did not detect significant differences in the expression of the three tested circRNAs between samples that expressed the wild-type Nsp14 (Nsp14 WT) and the ExoN mutants (Figure 4 -figure supplement 2D). Taken together, these data indicate that the exonuclease activity of Nsp14 is not required for the observed transcriptome effects.
We also generated a N7-guanine-methyltransferase mutant (Nsp14 D331A), by altering a conserved site previously shown to be essential for the N7-guanine-methyltransferase activity (Jin et al., 2013) (schematic Figure 4 -figure supplement 2A). Then, we checked by RT-qPCR how Nsp14 D331A affects the expression of some mRNAs and circRNAs that are deregulated with Nsp14 WT. Interestingly, overexpression of Nsp14 D331A does not alter the expression of the tested targets ( Figure 4G and 4H). The lack of effect is not due to low expression or stability of the Nsp14 D331A mutant, as it was expressed at similar levels than Nsp14 WT as assessed by Western Blot (Figure 4 -figure supplement 2E). In summary, our results indicate that the N7guanine-methyltransferase domain of Nsp14 is required for driving the gene expression changes described in this study.

Dissecting the transcriptome changes upon Nsp14 expression
To gain further insights into the gene expression program orchestrated by Nsp14, we performed a time-course experiment. Briefly, we prepared and sequenced 3' RNA-seq libraries using RNA isolated from cells transfected with the Nsp-14 expressing plasmid 6, 8, 12, 24, and 48 hours after transfection. As a control, we collected data from cells transfected with a control plasmid (eGFP) harvested at the same time points. Importantly, we observed detectable levels of  Table 9), as well as an increase in the number of significantly enriched GO terms ( Figure 5C and Table 10). For example, expression of genes related to cell metabolism, RNA metabolism, and cell cycle were altered relative to the control only at 24-and 48-hours post transfection ( Figure 5C).
We reasoned that the mRNAs that were differentially expressed at early time points could provide cues regarding the mechanism of action of Nsp14. We observed only 19 DEGs at the 8-hour time point and 200 DEGs at the 12-h time point ( Figure 5B). Interestingly, CXCL8, which encodes IL8, was the most highly induced mRNA at the 8-hour time point (fold change =11.5, pvalue =1.2e -5 , Table 8). Furthermore, we found that CXCL8 was induced by more than 12-fold consistently between 8 and 48 hours post transfection ( Figure 5D). It was previously reported that CXCL8 mRNA is significantly upregulated following infection of cells with SARS-CoV-2 and that the levels of IL8 are increased in COVID-19 patients (Blanco-Melo et al., 2020;Coperchini et al., 2020;Del Valle et al., 2020;Gardinassi et al., 2020;Li et al., 2021).
To investigate whether Nsp14 expression induces the transcription of CXCL8, we generated a firefly luciferase reporter under the control of the CXCL8 promoter (-500 bp to +80 bp around the transcription start site, Figure  CXCL8 expression is controlled by several transcription factors, including NFkB, a master regulator of inflammation (Kunsch and Rosen, 1993). Recent reports have shown that NFkB is activated during infections with SARS-CoV and SARS-CoV-2, although the mechanism and the importance of this activation are not well-established (Liao et al., 2005;Kircheis et al., 2020;Park and Lee, 2020;Hariharan et al., 2021;Hsu et al., 2021;Li et al., 2021). Strikingly, some of the top upregulated RNAs at 8-hours post transfection of Nsp14 were NFKBIA, JUN, and ATF3 (Figure 5 -figure supplement 1D), which encode proteins involved in the activation and modulation of the NFkB pathway (Stein et al., 1993;Baldwin, 1996;Pahl, 1999).
To determine if NFkB activation contributes to the early response of cells following Nsp14 expression, we performed an unbiased analysis for transcription factor binding site enrichment on the promoter regions of the genes with strong upregulation at the 8-hour time point. We found a 18 significant enrichment in those promoters for binding sites recognized by components of the NFkB transcription complex and the NFkB regulatory network ( Figure 5F, Table 11, and Table 12, Stein et al., 1993;Baldwin, 1996;Pahl, 1999). We next tested whether expression of Nsp14 activates a luciferase reporter containing a minimal promoter and NFkB sites (Wilson et al., 2013). Indeed, we observed that expression of Nsp14 lead to an approximately 1.6-fold increase in the luciferase signal ( Figure 5 -figure supplement 1F), demonstrating that this viral protein activates the NFkB pathway.
Then, for the mRNAs strongly upregulated early after Nsp14 transfection, we performed an analysis of their extended promoter regions (1,000 bp upstream of the transcription start site).We found that these regions were enriched for binding sites for factors involved in the interferon response such as EGR2, STAT2, and SP3 (Table 13). Previous work showed that SARS-CoV-2 infection modulates both the interferon and the NFkB pathways (Blanco-Melo et al., 2020;Vanderheiden et al., 2020;Wyler et al., 2021) and that Nsp14 affects these pathways (Hsu et al., 2021;Li et al., 2021). Our data indicates that the initial response to Nsp14 expression involves modulation of both the NFkB and the interferon pathways.

IMPDH2 partially mediates the transcriptional changes induced by Nsp14 expression
Given the strong transcriptional changes we observed upon Nsp14 expression, it is possible that this viral protein acts directly as a transcription factor. However, previous research showed that Nsp14 localizes in the cytoplasm when expressed in HEK293T cells (Gordon, Hiatt, et al., 2020). We then performed subcellular fractionation and chromatin purification, followed by Western Blot at different timepoints after transfection, and we showed that Nsp14 localizes in the cytoplasm, and it does not associate with chromatin demonstrating that the transcriptional effect is indirect ( Figure 6 -figure supplement 1A).
Recent proteomics studies found that Nsp14 interacts with the host protein IMPDH2 (Gordon, Hiatt, et al., 2020;Gordon, Jang, et al., 2020), which catalyzes the rate-limiting step of guanine biosynthesis, but did not reveal whether this interaction results in any metabolic or cellular change (Li et al., 2021). One possibility is that Nsp14 activates IMPDH2 which then modulates gene expression directly or indirectly. In support of this possibility, expression of Nsp14 changed the levels of genes involved in purine metabolism ( Fig 1C). Furthermore, we found IMPDH2 mRNA expression significantly downregulated upon Nsp14 expression ( Figure 6A). As GTP levels inversely regulate IMPDH2 expression (Glesne, Collart and Huberman, 1991), we hypothesized that Nsp14 promotes the activation of IMPDH2, leading to higher GTP levels that, in the end, negatively regulate IMPDH2 mRNA. To test whether expression of Nsp14 alters IMPDH2 activity, we performed metabolomic analysis of lysates of cells transfected with a plasmid for expression of Nsp14 or a control plasmid. Consistent with modulation of IMPDH2 activity by Nsp14, we detected an approximately 3-fold increase in the levels of xanthine diphosphate and GTP ( Figure 6B and Table 14), the downstream products of XMP. Furthermore, we observed increased concentrations of other nucleotides that directly or indirectly derive from IMP ( Figure 6B and Table 14), probably due to activated compensatory pathways.
To determine if IMPDH2 activity mediates the Nsp14-driven changes in gene expression, we compared the gene expression changes upon Nsp14 expression in the presence of an IMPDH2 inhibitor. Briefly, we transfected cells with a Nsp14 or with a control plasmid (eGFP) and 8 hours later treated the cells with Mycophenolic acid (MPA), a non-isoform-selective pan-IMPDH inhibitor, or DMSO. We harvested the cells 40 hours after the treatment. We used a concentration 20 of MPA of 0.5 µM to avoid cytotoxic effects (Qasim et al., 2011). Although we treated the cells with a low concentration of MPA, MPA treatment could reduce the concentration of cellular GTP, which is the primary energy source for mRNA translation. We confirmed that MPA treatment did not alter Nsp14 or IMPDH2 levels by Western Blot ( Figure 6C). We next performed 3' RNA-seq on the same samples. Interestingly, we observed that MPA treatment diminishes the effects induced by expression of Nsp14 ( Figure 6D, Figure  Finally, we investigated whether IMPDH2 also mediates the effects induced by Nsp14 on alternative splicing. To do so, we treated cells with MPA or MZR and we checked by RT-qPCR the abundance of some of transcripts that displayed intron retention in the total RNA seq data ( Figure 3B and Table 5). As expected, we detected increased abundance of retained introns for In summary, these results indicate that IMPDH2 mediates both the gene expression changes, and the alternative splicing effects induced by Nsp14.

DISCUSSION
Viral infection leads to a complex set of events orchestrated by multiple viral proteins. In the last two years, new studies shed light in some unexpected functions of diverse SARS-CoV-2 proteins. For instance, protein interactome experiments defined the interaction map of SARS-CoV-2 proteins with host cellular proteins (Gordon, Hiatt, et al., 2020;Gordon, Jang, et al., 2020).
Here we used an alternative approach; we performed an expression screen to determine how the individual expression of SARS-CoV-2 proteins influenced host cell gene expression.
Nsp14 displayed the most striking effects, deregulating the expression of approximately 4,000 mRNAs and altering splicing. While we acknowledge that the conditions under which we performed our experiments might lead to higher levels of Nsp14 than those observed during infection, we are convinced that our findings are highly relevant. Indeed, both the effects on gene expression and alternative splicing significantly overlap with the ones occurring during SARS-CoV-2 infection. Furthermore, during the infection Nsp14 localizes in the replication organelles that are formed in the cytoplasm (V'kovski et al., 2021b), and in our model we also see that Nsp14 localizes in the cytoplasm. Overall, our data indicate that our model, despite some limitations, recapitulates some gene expression changes that occur during SARS-CoV-2 infection.
Regarding the other tested proteins, when expressed individually M, Nsp9, E, ORF9b, ORF3a, Nsp13, and Nsp1 also altered the host cells transcriptome in ways that will require further investigation. Although most of the tested proteins did not dramatically alter gene expression, this does not rule out their roles in remodeling of host cell gene expression. For instance, some may 23 require other viral proteins as co-factors. In the future, it will be interesting to test how different combinations of viral proteins that are known to work in complex may affect host cell gene expression.
Our findings unrevealed an unexpected function of the SARS-CoV-2 Nsp14 in regulation of host cell transcription, splicing, and circRNAs expression. Notably, we also showed that IMPDH2 mediates these effects, as pharmacological inhibition of this cellular enzyme partially reverts them.
Our analysis of the transcriptomes of HEK293T cells as a function of time after Nsp14 expression revealed that the transcriptional response is temporarily subdivided into two phases.
The early response involves NFkB activation, whereas the later phase includes effects on cellular metabolism, RNA metabolism, and cell-cycle control. Notably, NFkB is one of the first pathways activated upon SARS-CoV-2 infection (Li et al., 2021;Wyler et al., 2021). Furthermore, early after transfection of Nsp14, we detected increased CXCL8 expression, which encodes IL8. IL8 is elevated in plasma of COVID-19 patients (Blanco-Melo et al., 2020;Hariharan et al., 2021;Li et al., 2021) and is a marker of acute severe lung distress (Adcock et al., 2015), an inflammatory condition that can lead to severe complications or death, further confirming the relevance of our model in mimicking the molecular events triggered by SARS-CoV-2 infection. Besides, in accordance with our data, a recent study showed that IMPDH2 mediates the activation of CXCL8 expression and of NFkB transcriptional reporters during Nsp14 expression (Li et al., 2021).
Nevertheless, for a subgroup of genes the response appears to involve alterations of splicing patterns and other post-transcriptional events. Remarkably, splicing is strongly disrupted upon SARS-CoV-2 infection (Banerjee et al., 2020) and we see a significant overlap between the alternative splicing events induced by Nsp14 and the ones occurring during the infection. Furthermore, Nsp14 affects the expression of genes encoding proteins with functions in RNA processing and splicing. However, unlike the splicing modulation by Nsp16 which has been proposed to be general (Banerjee et al., 2020), the effect of Nsp14 seems to be specific to certain exons and introns as, in most cases, only one splicing event is modulated per gene. It is likely that Nsp16 and Nsp14 influence alternative splicing through different mechanisms, although we could not compare the effects of the two proteins as the plasmid for Nsp16 expression was not in the collection we utilized in our screening (Gordon, Jang, et al., 2020). In any case, while Nsp16 directly binds to pre-mRNA (Banerjee et al., 2020), Nsp14 may trigger a cascade of events that culminate in altered splicing patterns.
Surprisingly, we also detected increased levels of circRNAs upon Nsp14 expression.
Specifically, Nsp14 altered the expression of about 10% of the cellular circRNAs. The magnitude of this effect is striking, considering that deletion or knock-down of factors involved in circRNAs production generally affects the expression of a smaller percentage of circRNAs (Conn et al., 2015;Errichelli et al., 2017;Knupp et al., 2021). We showed that this increase is not due to enhanced transcription of the hosting RNAs, and our results indicate that Nsp14 influences biosynthesis and/or stability of those circRNAs. Published work showed that circRNAs accumulate in non-proliferating cells (Bachmayr-Heyda et al., 2015). Nsp14 alters the expression of genes regulating cell cycle and proliferation, and this may result in inhibition of pathways related to circRNA degradation or in inhibition of cell division, ultimately leading to circRNAs accumulation. However, this possible explanation will require future investigation.
Interestingly, recent studies indicate that circRNAs modulate some immune factors by repressing their pro-inflammatory activity in normal conditions (Liu et al., 2019). On the contrary, upon viral infection, circRNAs are degraded and these immune factors can orchestrate the 25 appropriate response in the cell (Liu et al., 2019). Indeed, circRNA titer dramatically decreases upon viral infection (Liu et al., 2019;Chen et al., 2020). The fact that we detected an increase in levels of certain circRNAs upon Nsp14 expression suggests a possible mechanism for the immune evasion orchestrated by Nsp14.
Nsp14 is conserved among coronaviruses and it participates to viral replication and immune surveillance escape (Ogando et al., 2020). Indeed, on the one end it works as an exonuclease (ExoN), important for the correct replication of the viral genome, on the other hand it functions as N7-guanine-methyltransferase for modifying the 5'-cap of the viral genome.
Notably, both domains are important for successful viral replication and are considered potential drug targets (Otava et al., 2021;Saramago et al., 2021). Our data indicate that N7-guaninemethyltransferase domain, but not the ExoN one, is necessary for Nsp14 effects on host gene expression.
Nsp14 interacts with the Replication Complex and in proximity with other enzymes, such as Nsp16, involved in the modification of the 5' cap (Romano et al., 2020). Although we excluded the importance of the co-factor Nsp10 (Ma et al., 2015b), we did not test how the combination of Nsp14 with other SARS-CoV-2 proteins modulate gene expression in the host cells. As mentioned above, in the future it would be interesting to co-express Nsp14 and Nsp16, since both are involved in the 5'cap modification (and, as mentioned before, the N7-guanine-methyltransferase domain of Nsp14 is important for the effects on gene expression), and both proteins disrupt splicing of cellular transcripts (Banerjee et al., 2020;Romano et al., 2020).
Finally, we showed that IMPDH2 meditates the effects on cell transcriptome, altered splicing, and circRNAs deregulation induced by Nsp14 expression, as pharmacological inhibition of IMPDH2 activity with MPA and Mizoribine mitigated the effects of Nsp14. Intriguingly, MPA 26 is a known immunosuppressant and antiviral agent (Pan et al., 2012;Dang et al., 2017), and it inhibits SARS-CoV-2 propagation (Kato et al., 2020).

Ideas and Speculations
Although our data clearly indicate that IMPDH2 mediates the effects induced by Nsp14, we still do not know the detailed mechanism. We reported altered expression of genes involved in purine synthesis, including a decrease in IMPDH2 mRNA levels, and we also detected elevated cellular GTP concentrations upon Nsp14 expression. Taken together, these data suggest that Nsp14 activates IMPDH2. Besides, previous studies showed that Nsp14 interacts with IMPDH2 (Gordon, Hiatt, et al., 2020;Gordon, Jang, et al., 2020), but we do not know whether the direct interaction between Nsp14 and IMPDH2 is sufficient for triggering the massive gene expression chenges or if this system requires a co-factor.
Notably, IMPDH2 binds nucleic acids (McLean et al., 2004;Mortimer and Hedstrom, 2005) and acts as a transcription factor in Drosophila (Kozhevnikova et al., 2012). However, although IMPDH2 is conserved (Hedstrom, 2009) and can localize in the nucleus of mammalian cells (Ahangari et al., 2021), the transcription factor activity of IMPDH2 has not been demonstrated in mammalian cells. Whether Nsp14 triggers cytoplasmic-nuclear shuffle of IMPDH2 in mammalian cells will require further investigation. Intriguingly, a recent study reported that IMPDH2 can bind a circRNA and that this interaction activates IMPDH2 enzymatic activity (Wang et al., 2021). In the future, it will be interesting to check the circRNA binding capacity of IMPDH2 in our model, where circRNAs expression is highly deregulated. Besides, IMPDH2 can also associate with polyribosomes and modulates the translation of specific mRNAs (Hsu et al., 2021).

27
In addition, another recent report indicates that Nsp14 blocks the interferon response, which is upstream of the NFkB pathway, by shutting down global translation 24 hours after Nsp14 transfection (Hsu et al., 2021). However, our data indicate that the transcriptional activation of the NFkB pathway occurs very early after transfection. Furthermore, we confirmed that some effects are transcriptional, and we detected alterations on splicing and circRNAs levels, events which are not necessarily downstream to translational inhibition. Moreover, we showed that the transcriptional effects occur soon after Nsp14 expression, as treatment with MPA or Mizoribine 8 hours after transfection partially recovered the gene expression profile of the cells. This indicates that some of the gene expression changes might precede the global translation shutoff (Hsu et al., 2021). Therefore, we speculate that Nsp14 can alter gene expression by disrupting both the transcriptome and translation. Another possibility is that the transcriptional and the translational effects depend on the levels of Nsp14, that cannot be compared between the two studies (Hsu et al., 2021 and this study).
Although we showed that the N7-guanine-methyltransferase domain of Nsp14 is crucial for mediating the gene expression changes in the host cell, previous studies indicated that this domain is important for mediating Nsp14 translational repression (Jin et al., 2013;Hsu et al., 2021). This leads us to further speculate that Nsp14 may disrupt transcription, alternative splicing, and translation. Interestingly, Nsp14 can methylate free GTP in the cell and convert it to m7GTP, that is known to block cap-dependent translation by competing with mRNA for eIF4E (Jin et al., 2013). Our results indicate that Nsp14 increases cellular GTP concentration by activating IMPDH2. Therefore, a compelling model would be that Nsp14 activates IMPDH2, causing not only to altered transcription and splicing, but also increased cellular GTP, that is then converted to m7GTP by Nsp14, ultimately leading to general translation repression. However, this hypothesis 28 will require further investigation and could only explain only a small fraction of the changes observed when expressing Nsp14, as expression of a viral factor that inhibits protein synthesis directly (Nsp1) altered the levels of a much smaller set of genes.
To sum up, we here provided novel insights regarding the function of Nsp14 in altering gene expression of the host cell through the interaction with the cellular enzyme IMPDH2 in the context of SARS-CoV-2 infection.

IMPDH2 inhibitors treatment
Mycophenolic acid (MPA) (Sigma Aldrich) was a kind gift from Hedstrom lab at Brandeis University and was reconstituted in DMSO (Sigma-Aldrich) and aliquots were stored at -20 o C. Random hexamers were used for the nascent RNA experiment showed in Figure 2F and  Fluorescence intensities were plotted versus the number of cycles by using an algorithm provided by the manufacturer. All the primers used in this assay are listed in Table 16. Regarding the RT-qPCR for detecting retained introns ( Figure 6G and Figure 6 -figure supplement 2G, we positioned one primer in the exonic sequence before or after the retained intron, and the other primer in the retained intronic sequence (indicated as exon-intron). In parallel, we also run a reaction using both primers in the exonic regions surrounding the retained intron. We normalized the abundance of the exon-intron and the exon-exon product on a housekeeping gene (TBP mRNA). Finally, we plotted the intronic-exonic ratio, defined as the ratio between the relative abundance of the exon-intron product on the relative abundance of the exon-exon one.

Luciferase Assay
The day before transfection 0.18*106 cells per well were plated into 24 well plates. Cells were transfected as aforementioned using 750ng of total DNA. Cells were transfected with 75ng of 42 Firefly reporter, 75ng of Renilla reporter, and 600ng of pLVX-EF1alpha-SARS-CoV-2-Nsp14-2xStrep-IRES-Puro or 600ng of empty vector. We used the following plasmids as reporters: pGL_CXCL8_promoter_Firefly, generated as indicated above. Coelenterazine (Goldbio.com) was added to the sample and luminescence was measured again.
Non-transfected cell lysate was used to measure the background signal, that was subtracted before calculating the firefly/renilla ratios.

Linear RNA alignment and quantification
Raw reads were aligned to the human genome (hg38) using STAR (Dobin et al., 2013). Mapped reads in the 3' UTR were quantified using End Sequence Analysis Toolkit (Gohr and Irimia, 2019)for 3' RNA libraries. FeatureCounts (Gohr and Irimia, 2019) was used to quantify mapped transcripts from total RNA-seq libraries. DEGs were found using DESeq2 as previously described (Sarrion-Perdigones, Gonzalez and Venken, 2020). DEGs with |L2FC| > 1, p-adjusted value <0.05 were considered significant and used as input in downstream Gene Ontology (GO) analysis (DAVID,v 6.8). For the analysis of the timepoint 3' RNA-seq experiment, we used |L2FC| > 0.5, p-adjusted value <0.05. HEK293T cell transcriptome was used as a reference to query for Enriched GO terms from up-and down-regulated DEG lists. Gene set enrichment analysis (GSEA) was performed using the fgsea package in R (Sergushichev et al., 2016). Gene rank was determined prior to GSEA by calculating -log(p-value)* sign(log2 fold change) (Sergushichev et al., 2016).
Gene sets used for GSEA were downloaded from Molecular Signatures Database (http://www.msigdb.org). To compare our RNA-seq results with other studies we used a gene set from MSigDB containing DEGs from published RNA-seq data. A curated list of terms from GSEA were plotted to assess the transcriptional response to viral protein expression. Data visualization was carried out using ggplot2 in R. Non-omics statistical analysis and data visualization was completed using Prism.

circRNA detection and differential expression analysis
Computational analysis of circRNAs was performed on total RNA-seq data. circRNAs were detected by searching unaligned reads for head-to-tail splice junctions as previously described (Memczak et al., 2013). Differentially expressed circRNAs were found by DESeq2 (Love, Huber and Anders, 2014). circRNA reads were normalized using size factors computed by DESeq2 together with all mapped reads.

Intronic reads quantification and analysis
44 Reads were aligned with STAR-aligner (Wilson et al., 2013) to the human genome and transcriptome version Hg38. Intronic reads were extracted and counted using Featurecounts function in R with an intronic region reference.

Alternative splicing analysis
Percentage of inclusion (PSI) quantification was done with Vast-Tools (Tapial et al., 2017). Delta PSI was then calculated from the mean of each condition. To determine events changing between Nsp14 and transfection control we chose a minimum of 15% difference in mean PSI and no overlapping between replicas. Introns and exons feature analysis was done using Matt (Gohr and Irimia, 2019).
With the same approach, we re-analyzed total RNA sequencing from HEK293T-hACE2cells infected with SARS-CoV-2 GSE169158 (Sun et al., 2021) . Then, we tested the significance of the overlap between the Alternative Splicing Analysis on our dataset and this dataset by performing a randomization of sets with equal size, evaluating their overlap to generate a null distribution, and calculating the probability of our overlap (10.000 replicas p-value <10 -15 ).For this analysis we used the function enrichment_test from the RVenn package in R.

Motif enrichment analysis
We extracted the promoters' sequences of upregulated genes with a fold change of at least 1.75 (log2FC > 0.8 and adjusted p-value > 0.05) using windows of 100, 500 and 1000 bases upstream of each annotated transcription start site for each gene. We then used the MEMEsuite (Bailey et al., 2009) using the tool AME(McLeay and Bailey, 2010) for motif enrichment analysis. We used the motif data base HOmo sapiens COmprehensive MOdel COllection (HOCOMOCO) v11 transcription factor (TF) binding models (binding profiles or binding motifs) for human transcription factors. As control for the input sequences, we used the tool provided for 45 scrambled sequences option. Finally, we used MEMEsuit tool FIMO (Grant, Bailey and Noble, 2011) to identify the position of these motifs in the analyzed 100bp promoters. We extracted the promoters' sequences of upregulated genes with a fold change of at least 1.75 (log2FC > 0.8 and adjusted p-value > 0.05) using windows of 100, 500 and 1000 bases upstream of each annotated transcription start site for each gene. We then used the MEMEsuite tool AME(McLeay and Bailey, 2010) for motif enrichment analysis using. We used the motif data base HOmo sapiens COmprehensive MOdel COllection (HOCOMOCO) v11 transcription factor (TF) binding models (binding profiles or binding motifs) for human transcription factors. As control for the input sequences, we used the tool provided scrambled sequences option. Finally. we used MEMEsuit tool FIMO (Grant, Bailey and Noble, 2011) to identify the position of these motifs in the analyzed 100bp promoters.

46
List of the tables Table 1: List of DEGs in the screening described in Figure 1. 3' RNA-seq. Table 2: List of GO terms described in Figure 1. 3' RNA-seq. Table 3: GSEA analysis described in Figure 2. Total RNA-seq. Table 4: Differential gene expression results described in Figure 2. Table 5: Alternative Splicing analysis results described in Figure 3. Table 6: Feature analysis with MATT for upregulated intron retentions and downregulated exon inclusion events described in Figure 3 -figure supplement 1.   Figure 3. Table 9: DEG at different time points described in Figure 5. 3' RNA-seq.  Figure 5. 3' RNA-seq. Table 11: Ame analysis with 100bp of promoter described in Figure 5. Table 12: Fimo analysis with 100bp promoter described in Figure 5. Table 13: Ame analysis with 1000bp promoter described in Figure 5.  : list of plasmids used for SARS-CoV-2 protein expression in the initial screening described in Figure 1.              GAPDH is used as cytoplasmatic marker, and H3K27me3 as a chromatin bound marker. See

Control Nsp14
A.