Multi-OMICs landscape of SARS-CoV-2-induced host responses in human lung epithelial cells

Summary COVID-19 pandemic continues to remain a global health concern owing to the emergence of newer variants. Several multi-Omics studies have produced extensive evidence on host-pathogen interactions and potential therapeutic targets. Nonetheless, an increased understanding of host signaling networks regulated by post-translational modifications and their ensuing effect on the cellular dynamics is critical to expanding the current knowledge on SARS-CoV-2 infections. Through an unbiased transcriptomics, proteomics, acetylomics, phosphoproteomics, and exometabolome analysis of a lung-derived human cell line, we show that SARS-CoV-2 Norway/Trondheim-S15 strain induces time-dependent alterations in the induction of type I IFN response, activation of DNA damage response, dysregulated Hippo signaling, among others. We identified interplay of phosphorylation and acetylation dynamics on host proteins and its effect on the altered release of metabolites, especially organic acids and ketone bodies. Together, our findings serve as a resource of potential targets that can aid in designing novel host-directed therapeutic strategies.


INTRODUCTION
The rapid emergence of the COVID-19 pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) still shows no abatement. 1 Despite the availability of several vaccines, the emergence of newer variants with enhanced transmissibility and pathogenesis is a major cause of concern. Current efforts to develop antiviral therapeutics, especially host-directed therapies, are vital to significantly reduce the impact of both the current and future coronavirus epidemics. Several classes of FDA-approved drugs and drugs currently in clinical trials have been repurposed, showing promising results in hospitalized patients with severe diseases. [2][3][4][5] Development of host-directed therapeutics largely depends on the available information as the viral genome continuously evolves to adapt to the host environment and evade pre-or post-countermeasures. Several host factors that promote or restrict SARS-CoV-2 replication have been identified by genome-wide CRISPR knockout screens. [6][7][8][9][10][11][12][13][14] Furthermore, multi-OMICS studies have provided vital insights into understanding the viral profile, mechanisms of pathogenesis, and identifying host factors. [15][16][17][18][19][20][21][22][23] Notably, the role of protein post-translational modifications (PTMs), are increasingly recognized as key mechanisms viruses employ to target signaling pathways regulating host immune responses. [24][25][26] Nevertheless, strainspecific differences and their effects on the host cellular machinery, especially in evading the immune response and facilitating further transmissibility and clinical manifestations, are only recently being explored. 27 Given the evolving number of variants, a multi-OMICS approach focusing on the dynamic interplay of PTMs is vital aiding in therapeutics development and is helpful in pandemic preparedness. Most studies have focused on elucidating the phosphorylation dynamics and alterations in ubiquitin and glycosylation. 18,20,[27][28][29][30][31] Some are restricted to characterizing the glycosylation patterns of virus/membrane proteins. [32][33][34][35] Emerging evidence highlights protein lysine acetylation as a key regulatory mechanism originally thought to modulate epiproteome. Epigenetic responses are now increasingly identified as promising therapeutic targets in human viral infections. [36][37][38] However, alterations in global host acetylation dynamics  Figure S2B). Specifically, interferons-IFNB1 (Type I IFN) and IFNL1, -2, -3, -4 (Type III IFN) were induced as early as 12 hpi with the overexpression of antiviral factors, and interferon-stimulated genes (ISGs) such as OAS1-3, IFIT1-3, IFITM at both transcriptomes and proteome levels ( Figures 2B and 2C).
Notably, CMPK2, an immunomodulatory ISG associated with antiviral responses 56-59 in immune cells, and a kinase involved in mtDNA synthesis 60 and DNA repair, 61 was upregulated at 24 and 48 hpi both at transcriptome and proteome levels. Using immunoblot analysis, we confirmed the induction of OAS1, ISG15, TRIM5a, RNase-L, MX1, and CMPK2. ( Figure 2D). Comparing the ISG expression profile with that of the early lineage and Alpha variant 27 revealed a similar extent of expression in the Trondheim strain and Alpha variant compared to the early lineage viruses (IC19 and VIC), indicating strain-specific differences in inducing interferon response ( Figure 2E).
Although it is well established that Interferons activate Toll-like receptor (TLR) gene expression in macrophages during viral infections, 62 including in SARS-CoV-2, the expression of TLR genes in response to viral infections in epithelial cells is less known. Several members of the TLR family, including TLR1 (24 and 48 hpi), TLR3 (24 and 48 hpi), TLR4 (6 and 24 hpi), and TLR6 (24 and 48 hpi), were found to be upregulated in the Calu-3 transcriptome in response to infection. TLR2 has been reported to sense the SARS-CoV-2 envelope protein and thereby produce inflammatory cytokines, causing inflammation and damage in the lungs. 63 Although we did not identify differential expression of TLR2, increased expression of TLR1 and TLR6, known to form heterodimers with TLR2, 64 was observed. Notably, TLR3 activation has been reported in SARS-CoV-2-infected Calu-3 lung epithelial cells, 65 whereas there is indirect evidence indicating the involvement of TLR4 in epithelial cells during SARS-CoV-2 infection. 66,67 Concerning the expression of the pro-inflammatory cytokine signature (TNF, IL1B, IL6, CCL2, CXCL9, CXCL10, CXCL11) previously shown to be upregulated in plasma and/or BAL of severe COVID-19 patients, 68

SARS-CoV-2 mediated modulation of PTM dynamics affects key host signaling pathways contributing to viral pathogenesis
To assess the impact of SARS-CoV-2 infection on host cellular signaling dynamics, serial enrichment analysis of phosphorylation and lysine acetylation was performed. Supervised clustering revealed 8 distinct clusters corresponding to early (clusters 5 and 6), intermediate (cluster 6), late (hyperphosphorylated: clusters 1, 7, and 8; hypophosphorylated: cluster 2), and sustained responses (clusters 3 and 4) ( Figure 3A). A minimal overlap was observed when comparing the hyperphosphorylated sites identified in the current study with other studies on the SARS-CoV-2 phosphoproteome 27,29,31 ( Figure S3A). Given the differences in the cell types, strains and platforms, it is not surprising that a large majority of sites identified vary across studies. Nevertheless, pathway analysis revealed enrichment of cytokine signaling, activation of DNA damage and repair pathways, Hippo signaling, cell cycle regulation, RNA splicing, and regulation of translation and transcription across different clusters indicative of viral-mediated alterations of host cellular machinery ( Figure 3B, Table S5).
Notably, differential phosphorylation of proteins related to the antiviral immune response (Cluster 1), including increased phosphorylation of STAT1 and STAT3 at S727 12 hpi, was observed with sustained levels seen even at 48 hpi ( Figure 3A). Immunoblot analysis further validated our findings ( Figure 3D). Phosphorylation and activation of STAT1/3 are required for the complete transcriptional activity of the ISGF  Figure S1).  70 Increased phosphorylation of IRF9 at S136 and IRF3 at S396 (Table S3, Figure 3D) was also observed. IRF9 interacts with phosphorylated STAT1 and IRF3 to form the ISGF3 complex, indicating induction of antiviral immune response that correlates with increased expression of ISGs. These findings strongly establish links to increased ISG levels observed both at the transcript and proteome level, implicating the establishment of a host antiviral state. [71][72][73] Cluster 2 revealed significant enrichment of time-dependent decrease in phosphorylation of proteins involved mainly in the regulation of cell cycle processes, DNA damage and repair pathways and chromatin modification such as CDK1, CDK2 (T14, Y15), CDKN1A (S130), DNA topoisomerases (TOP2A (S1106, S1213, S1504), TOP2B (T1431) and transcription factors/co-regulators-RB1, APC, TP53). Enrichment of DNA damage response pathway at the early time point corroborates with an earlier study that reported that infectious bronchitis virus (IBV) activated ATR-dependent cellular DNA damage response to induce cell-cycle arrest for the enhancement of viral replication and progeny production. They further demonstrated that the interaction between Coronavirus nsp13 and DNA polymerase d was essential for the induction of DNA damage response and cell-cycle arrest at the S phase. 74 Responses enriched as early as 3hpi indicate enrichment of protein autophosphorylation required for activation of kinases modulating downstream cellular processes.

ll
It has been previously shown that sustained RAF/MAPK signaling results in the downregulation of ROCKI and Rho-kinase, two-Rho effectors required for stress fiber formation and promote cytoskeleton reorganization. 75 Inhibition of MEK functionally restored the activity of ROCK1/Rho-kinase in promoting cytoskeleton reorganization in NRK/RAS cells. We observed sustained MAPK1 activation on SARS-CoV-2 infection with decreased levels of pROCK1 24 hpi, thereby suggesting potential alterations in cytoskeleton organization (cluster 4). Increased phosphoERK1/2 (T202, Y204) was confirmed by immunoblotting ( Figure 3D). Clusters C2 and C6 demonstrate a significant time-dependent decrease in phosphorylation, specifically at 24 hpi. Proteins involved in Hippo/YAP signaling, Signaling by Rho GTPases, vesicle-mediated transport, tight junction proteins, including TJP1, TJP2, CTTN, SYMPK, among others, were significantly enriched. Hippo/YAP1 Signaling was enriched specifically in cluster 6, corroborating with hypophosphorylation of YAP1, an essential component of the Hippo/YAP1 signaling pathway, at S61. We also observed hypophosphorylation of other Hippo signaling proteins -TJP1, TJP2, DLG1, SCIB, among others. A time-dependent decrease in pYAP1 (S61) was independently confirmed by immunoblotting ( Figure 3D). Hippo signaling is further discussed in subsequent sections.
We also observed temporal changes in the SARS-CoV-2-responsive kinome and phosphatome ( Figure 3C), with several protein kinases demonstrating decreased phosphorylation levels over time. On the contrary, the phosphorylation levels of protein phosphatases such as PTPN3, CTDSPL, and CTDSPL2 increased with time, with 9 phosphatases hyperphosphorylated at 24 hpi. Of the 107 kinases found to be differentially regulated, 18 have been described as substrates of viral proteins, 23,29 including kinases involved in PI3K/ AKT signaling (GSK3B, RAF1, PAK4), cell cycle (PRKDC, CDK2, and TTK) and MARK kinase signaling (MARK2 and MARK3). Further PTM signature enrichment analysis using PTMSigDB revealed enrichment of signatures of MAPK signaling pathways including MAPK14, MAPKAP2, leptin pathway, PRKACA, cell cycle kinases CDK1, CDK2, and Aurora B at early time points of infection correlating with activation of DNA damage and repair pathways and altered cell cycle regulation ( Figure S3B). Furthermore, a substantial number of transcription factors were differentially phosphorylated in response to SARS-CoV-2 infection (E) A graph comparing transcript levels of interferon-stimulated genes (ISGs) between the Trondheim strain (S15), alpha, IC19, and VIC strains from Thorne et al. 27 (F) Cytokine array data showing differential levels of Ang-1, Dkk-1, uPAR, CXCL10, VEGF, MIC-1, and CD147 in Calu-3 cells post-infection. (See also Figure S2, Tables S1 and S2).  Figures S3C and S3D). Notably, TAF6 (S672) and TAF7 (S264), components of the DNA-binding general transcription factor complex TFIID, were hypophosphorylated on SARS-CoV-2 infection, especially TAF7, as early as 6 hpi. Phosphorylation of TAF7 (S264) mediated by TAF1 has been shown to influence the levels of cyclin D1 and cyclin A gene transcription by increasing TAF1 histone acetyltransferase (HAT) activity and histone H3 acetylation levels. 76 Emerging evidence strongly indicates the role of protein lysine acetylation as a regulatory mechanism in viral infection 77,78 and, therefore, can likely serve as potential therapeutic targets. 79,80 However, the dynamic alterations in protein acetylation on SARS-CoV-2 infection have not been explored. Analysis of protein lysine acetylation dynamics revealed 112 acetylated peptides differentially regulated on SARS-CoV-2 infection ( Figure 3E). Of these, site-specific alterations in acetylation of 53 proteins were largely independent of protein abundance changes except for CDK1 (K33), which was found to be hypoacetylated at 24 hpi. In contrast, the protein expression was downregulated only at 48 hpi. On the contrary, the acetylation levels of TMSB4X at K26 and K39 decreased at 6 hpi with the proteome expression downregulated 24 hpi. k-means clustering of the acetylome measurements enabled the segregation of the regulated genes into early, late, and mid-late responders to viral infection. We largely observed significant downregulation of acetyl sites on proteins as early as 6 hpi with respect to Clusters 2 and 3. In contrast, Cluster 1 included acetylation sites on proteins that were predominantly downregulated at a late time-point (48 hpi). Clusters 4 and 5 revealed hypoacetylation at 6 hpi in comparison with 3 hpi followed by increased acetylation (log2FC>0.5) post 12 hpi.
Pathway enrichment analysis largely revealed the enrichment of proteins involved in various biological processes ( Figure 3F). Notably, a time-dependent decrease in the acetylation levels across proteins belonging to several classes but primarily histones, epigenetic modifiers, and proteins involved in metabolic regulation. Cluster 1, which largely showed a delayed response, was enriched in proteins involved in glucose metabolism and hexose biosynthetic process, as well as histone subunits (H2B and H4C), acetylated on multiple sites. In addition to proteins enriched in glucose/hexose metabolism, Cluster 3 consisted of hypoacetylated proteins involved in the regulation of the mitochondrial organization. Acetylation of GAPDH at K251 was observed as early as 3 hpi, followed by a persistent decline at later time points of infection. Acetylation at this site is known to be mediated by PCAF and is required for nuclear translocation of GAPDH during apoptotic stress. 81 In the case of other metabolic enzymes, a progressive decline in the acetylation status was observed across time points, with a significant reduction post 24 h. Cluster 2 and Cluster 4 were enriched in proteins involved in viral transcription and translation termination, further corroborating the role of viruses inducing host protein translation shutoff activity. In addition, proteins localized to the membrane were enriched in Cluster 4, whereas proteins involved in macromolecule complex disassembly and translational elongation were enriched in Cluster 2. Acetylation of EIF5A2, a translation initiation factor essential for cell growth at K47, is responsible for regulating its subcellular localization with the hypoacetylated form predominantly localized to the cytoplasm. 82 Cluster 5 includes histones-H2AZ1 acetylated at K8, K12, and K14, H3C1 (K24), as well as keratins KRT17 (K219) and KRT19 (K215) were found to be hyperacetylated 12 hpi. H3C1 is a core component of the nucleosome and is involved in the process of post-transcriptional and translational regulation of genes. Furthermore, studies have revealed that the site H3C1 (K24) is amenable to acyl modifications and is highly responsive and reversibly regulated by nutrient availability. 83 Overall, our data is indicative of SARS-CoV-2 mediated triggering of deacetylases and likely restricted availability of acetate pool that results in significant decrease in protein acetylation at later time points of infection. These results correlate with an earlier report suggesting that active deacetylase activity is required to induce ISG expression and antiviral immune responses, 84 thereby creating an environment conducive to increased viral replication. iScience Article We next assessed if there was an overlap of dysregulated phosphorylation and acetylation datasets. Our data clearly shows interplay of phosphorylation and acetylation on host proteins. In all, we identified 37 proteins with differentially regulated multiple PTM sites. These include proteins involved in innate immune signaling (PPIA, HSP90AA1, ILF3, ANXA2, ANXA4, and CFL1), regulation of cellular metabolic processes as well as proteins involved in the regulation of mRNA stability, including ANP32A, a multifunctional transcriptional regulator and a component of the inhibitor of histone acetyltransferases complex. 85,86 Significant hypoacetylation was observed at 24-48 hpi, with the phosphorylation levels on T15 and S17 showing an opposite trend ( Figure 3G). Of interest, the levels of ubiquitination at K99, as demonstrated by Stukalov et al., were downregulated as early as 6 hpi with a minimal increase over time. Host proteins, especially ANP32A, interact with influenza virus RNA polymerase components, forming a replication platform essential to promoting vRNA synthesis. 87,88 Furthermore, ANP32A is known to be a component of the inhibitor of histone acetyltransferases complex and modulates the HAT activity of EP300/CREBBP (CREB-binding protein) and EP300/CREBBP-associated factor 89, 90 and is involved in the positive regulation of ISG transcription 85 which we observe at both transcript and protein levels. It remains to be determined if the differential acetylation at this site located in the LRR domain plays a role in CoV-2 replication and infectivity. As described previously, we also identified distinct phosphorylation and acetylation patterns on several sites on Vimentin (VIM), an essential host factor responsible for the entry and pathogenesis of SARS-CoV-2. 91 Decreased acetylation levels as early as 6 hpi (K294, K402) whereas distinct phosphorylation on several sites early (S73) or late (S83, S325, S420, and S430) were observed ( Figure S3E). Of interest, Stukalov et al. report decreased ubiquitination levels following infection with SARS-CoV-2 and SARS-CoV.
A comparison of the changes in the acetylation data with the ubiquitinome performed in ACE2-A549 cells 29 revealed 26 proteins common across the datasets. Of these, ribosomal proteins-RPL13A (K159), RPL34 (K36), RPL7 (K29), RPS25 (K57); LRRC59 (K135), and UBE2N (K82) were differentially modified at the same residue. However, the PTM dynamics observed were opposite in trend, especially at 24 hpi, wherein the acetylation levels were downregulated as opposed to increased ubiquitination. Notably, LRRC59, a leucine-rich repeat-containing ER membrane protein, plays an important role in innate immune signaling by modulating DDX58-mediated type I IFN signaling 92 and regulating trafficking of nucleic-acid sensing TLRs. 93 It was found to be hypoacetylated at 6 hpi with acetylation reaching basal level at 12 hpi. On the contrary, the ubiquitination levels at K135 and others (K71, K207) remain unchanged at 6 hpi but gradually increase across time points, with significantly elevated levels observed at 24 hpi ( Figure S3F). Overall, our findings emphasize the need to assess the effects of different post-translational modifications simultaneously, as the functional outcome of the cellular response is dependent on the concerted action of regulatory machinery.

Modulation of secretory effectors involved in SARS-CoV-2 infection
To identify molecular/metabolic pathways perturbed following SARS-CoV-2 exposure, we measured the extent of release of 100 metabolites in the culture supernatants of Calu-3 cells treated with SARS-CoV-2 for the time points described above ( Figure 1A). Of these, 85 were identified and quantitated across five time points after SARS-CoV-2 infection and were considered for further analysis. A total of 163 differential events were identified across the five time points, of which 146 were increased, and 17 were decreased (Figure 4A). The metabolites that were significantly changed in response to SARS-CoV-2 infection have been provided in Table S6.
Changes in the metabolite release and metabolic pathway perturbations were analyzed using targeted metabolome-wide association. Metabolite set enrichment analysis (MSEA) indicated significant enrichment of several metabolic pathways, including the TCA cycle, Glycolysis/Gluconeogenesis, butanoate metabolism, and metabolism pathways of several amino acids ( Figure 4B). Metabolites involved in amino acid metabolism, including branched-chain amino acids, and aromatic amino acids, were among the significantly enriched pathways at early and mid-time points of infection. At later time points, arginine metabolism and metabolites of the TCA cycle pathways were found to be significantly enriched. However, the extent of changes in the levels of arginine was not significant across the time points. We further compared and correlated our findings from the exometabolome data with intracellular metabolomics studies of SARS-CoV2 mediated dysregulation in cellular metabolism, 94,95 primarily focusing on glycolysis, TCA cycle and amino acid metabolism as it is known that cellular metabolism specifically energy metabolism supports the demand for viral replication. 96 iScience Article partial partial correlation with respect to TCA cycle metabolites and no correlation with metabolites of amino acid metabolism was observed ( Figure S4). Elevated levels of glycolytic metabolites, fumarate and malate, but not the other TCA intermediates were observed in the intracellular metabolome datasets.
In the exometabolome dataset, elevated levels of glycolytic metabolites were observed with the exception of pyruvic acid (decreased levels 12-48 hpi). Of interest, intermediates of the TCA cycle were found to be increased in the extracellular mileu in response to SARS-CoV-2 infection ( Figure 4C). An increase in succinate levels was observed at 24-48 hpi, whereas fumaric acid levels increased at 6-48 hpi, and alpha-ketoglutarate levels were increased across all the studied time points. However, cis-aconitate was only increased at 48 hpi. Further, we observed a slight decrease in the glutamine levels 24-48 hpi indicative of glutaminolysis and significant decrease in Serine, and Cysteine levels ( Figure S4) which can be converted to pyruvate and further to lactate. Notably, increased levels of secreted lactate and alanine were observed indicative of the host epithelial cells demonstrating a Warburg-like effect. The Warburg effect has been shown to support oncogenic viruses' infection and speculated to likely supports replication of SARS-CoV-2 replication in airway cells. 98 Furthermore, lactic acid has been shown to influence immune cell function (reviewed in 99 ). From these findings, it is imperatve that, SARS-CoV-2 infection results in metabolomic reprogramming of Calu-3 cells toward an aerobic glycolysis phenotype and is in concordant with earlier reports demonstrating the role of energy metabolism on viral replication 100 and rewiring of antiviral immune responses. 101 In addition, metabolites including the vitamin Pyridoxal, phenylpyruvate, organic acids such as, methylmalonic acid, and branched-chain amino acid metabolism intermediates: 3-methyl-2-oxovaleric acid and Ketoleucine and 3-Hydroxybutyrate showed increased levels across the studied time points (Figures 4D-4I). 3-Hydroxybutyrate, also known as Beta-hydroxybutyrate (BHB) is an endogenous ketone body acts as a highly efficient oxidative fuel. Recent studies implicate its vital role in immunomodulation. 102 Of interest, we observed increased levels of sialic acid N-acetylneuraminate as early as 3 hpi, and it has been recently demonstrated that N-acetylneuraminic acid (Neu5Ac) serves as a plausible alternative receptor for SARS-CoV as a key domain in CoV-2 spike protein binds to Neu5Ac, a process essential for viral entry into cells. 103,104 The increased secretion of several acids and ketone bodies in response to SARS-CoV-2 infection at all time points suggests metabolic ketoacidosis.

Global alterations in the viral proteome profile
In addition to the changes in host cellular proteins, analysis of upregulated DEGs revealed a significant increase in the expression of viral genes observed mainly at 24 hpi (Table S7). This was demonstrated by the low alignment rate observed against the human reference database, indicating that the RNA fraction at these time points was taken over by the virus. At the proteome level, nine canonical SARS-CoV-2 proteins were consistently identified and quantified across six time points of infection. In concordance with previous reports, over 3-fold induction of viral proteins with similar expression trends across all nine viral proteins was observed ( Figure 1C). The dynamics of protein expression across time points indicated that the viral protein synthesis increased continuously post-infection, with a peak observed at 24 hpi except for Orf8, with a slight decrease in expression observed at 48 hpi ( Figure 1C, Table S2). SARS-CoV-2 N, M, S, and Orf9b were among the most abundant proteins consistently detected.
Analysis of the PTM dynamics revealed 41 phosphorylation sites on eight viral proteins with maximum changes in the phosphorylation dynamics observed post 12 hpi ( Figure 5A). The trend observed was in accordance with the increase in protein abundance. Comparison with previously reported SARS-CoV-2 phosphoproteome enabled us to identify ten novel phosphorylation sites on 6 viral proteins, including ORF1a (S142, S2517) and ORF3a (T24, S272), ORF7a (S44), ORF8 (S103) ( Figure 5B). Replicase protein 1a (P0DTC1) is a polyprotein that is proteolytically cleaved to generate several non-structural proteins. We found the polyprotein to be modified on 3 sites; one corresponding to NSP1 (S142) and 2 sites corresponding to NSP3 (S2517, S2644). In concordance with earlier studies, multiple sites of phosphorylation on the  Figure S4 and iScience Article nucleocapsid protein, including 3 novel sites corresponding to T245, S327, and S379, clustered in the linker region between the RNA-binding (RBD) and dimerization domains as well as located in the dimerization domain (S327) were observed ( Figure 5C). The nucleocapsid protein plays an essential role in the replication, transcription, and assembly of the SARS-CoV-2 genome 105 and is also known to control host cell cycle machinery. 23, 106 The identified phosphorylation motifs were bioinformatically assessed using NetPhos 3.1 to predict potential host kinases phosphorylating the viral proteins ( Figure 5D). Consistent with previous findings, we too observed several motifs in nucleoprotein to be modified by CMGC kinases, including members of MAPK family CDKs and GSK3. Of interest, we observed several novel phosphorylation sites (N, ORF1a, ORF8) mapping to recognition motifs of the AGC family-PKA, PKB, and PKC and DNAPK, a critical player involved in the DNA damage response pathways. It is known that DNA-PK acts as a DNA sensor that activates innate immunity. However, its role in phosphorylating RNA virus N protein, especially at the novel site T379 located C-terminal of the dimerization domain, remains to be explored. 107 In concordance with previous results, we also observed an increased phosphorylation level of nucleocapsid protein. A recent study observed the highest percentage of disorder compared to other viral proteins as well as a high number of variable molecular recognition features (MoRFs) and classified as a highly disordered protein with a central role in viral pathogenesis. Furthermore, the viral acetylome profile revealed hyperacetylation of 6 sites on nucleoprotein, including at K65, K248, K249, K256, K341, and K362. Of these, two sites on the nucleoprotein (K248, and K249) have been recently demonstrated to be acetylated by P300/CBP-associated factor (PCAF) and general control nonderepressible 5 (GCN5). 108 The hyperacetylation of N protein coincides with the hypoacetylation of ANP32A and the induction of ISGs, which strongly demonstrates the critical role of ANP32A in modulating the host PTM dynamics, thereby contributing to viral pathogenesis.
We next looked at the positions of phosphorylated residues in the available and predicted structural models of SARS-CoV-2 proteins ( Figure S5). For those residues that could be mapped on the available structures, all the identified phosphorylation sites were plausible because they were located at the surface of the proteins and, therefore, should be accessible to kinases. Similar to the previous reports, 20 most of these residues were located in loops or at the edge of secondary structure elements. The majority of the side chains were not engaged in any intramolecular interactions (for example, S63 of ORF9b or T141 and S327 of N protein), though some residues formed hydrogen bonds that could support the positioning of some flexible loops (for example, S103 of ORF8 or S44 of ORF7a) or contribute to an interaction with another protein (T72 of ORF9b bound to the backbone of V556 of Tom70). In addition, many phosphorylation sites were found in regions that were missing from crystal or cryo-electron microscopy structures and therefore presumed to be likely in unstructured regions (for example, T24, S272, and T248 of ORF3a). Indeed, one such region was S176-S180 of Protein N, which was modeled as a flexible linker in the NMR ensemble 109 and is part of the serine/arginine-rich (SR) domain at the start of one of the predicted intrinsically disordered regions (IDR) (residues 1-44, 175-254, 365-$400). The SR region may participate in RNA binding 109 and is a hotspot for mutations in Protein N. 110 Of interest, S327 of N Protein, a phosphorylation site unique to our study, has also been found to be frequently mutated, typically to Leu. 110 For some of the most interesting cases for which experimental models were missing, we used computational modeling with TrRosetta. 111 Although computational models should be approached as only approximate, especially in the positioning of flexible regions and the orientation of independent domains, they all showed that the phosphorylation sites were located on protein surfaces and in flexible regions, including T172, S212, and S213 of Protein M, T24, S272 and T248 of ORF3a, S50 of ORF6, S26, T245 and T379 of Protein N, S142 of Rep1a (Nsp1), S2517 and S2644 of Rep1a (Nsp3). In summary, the phosphorylation sites observed in this study were predicted to be unlikely to play an important structural role in supporting the architecture of SARS-CoV-2 proteins, which corresponded well to their exposed positions and the lack of conservation noted in some cases. Instead, phosphorylation might play a role in regulating interactions mediated by these residues -the change in electrostatic potential or in the pattern of hydrogen bonds could alter the binding of other proteins or RNA. Whether these could be beneficial or detrimental to the SARS-CoV-2 remains to be validated by functional studies.  Figure S5 and Table S7).

OPEN ACCESS
iScience 26, 105895, January 20, 2023 13 iScience Article Integrated analysis of multi-omics datasets identified aberrant hippo signaling, DNA damage/repair, and ubiquitination machinery of host cells Integrated analysis of key signaling pathways and processes identified in our datasets using pre-defined gene sets from MSigDB 112 revealed significant changes at multiple levels on proteins involved in Hippo signaling, DNA/damage/response, protein ubiquitination, alternate splicing pathways ( Figures 6A-6G and S6A-S6H) and metabolic pathways involved in energy metabolism ( Figure S7). Although the phosphoproteome data showed a time-dependent decrease in the phosphorylation status of Hippo signaling, other pathways, including DNA damage/response, protein ubiquitination, and spliceosome pathways, demonstrated increased signaling. The proteome changes showed reduced expression, whereas the transcriptome showed increased transcription of genes belonging to these pathways.
Although in cancers, the Hippo signaling pathway is a tumor suppressor by nature, 113 various studies have indicated its crucial role in viral infections such as those caused by Hepatitis B virus (HBV), Hepatitis C virus (HCV), Human papillomavirus (HPV), Epstein-Barr virus (EBV), Zika virus (ZIKV), amongst others. 114 During the course of Hippo signaling, MST1/2 kinases (STK3/STK4) phosphorylate and activate Lats1/2 kinases, which in turn phosphorylate and inhibit essential transcriptional coactivators-YAP/TAZ (YAP1/ WWTR1). 113,115 When upstream kinases are repressed, dephosphorylated YAP/TAZ accumulate in the cell nucleus and associate with several transcriptional factors that stimulate genes involved in cell survival and proliferation, including TEADs, ERBB4, EGR1, TBX5, and members of the SMAD and RUNX family. 116 Notably, YAP was found to negatively regulate the antiviral immune responses in Sendai virus (SeV), vesicular stomatitis virus (VSV), and herpes simplex virus type 1 (HSV-1) infections. 117 Wang et al. further demonstrated that depletion of YAP in macrophages increased interferon beta expression. Also, YAP overexpression led to repression of IRF3 dimerization through YAP association and decreased nuclear localization of IRF3. Furthermore, viral-mediated activation of IKKε phosphorylated YAP led to its degradation, thereby relieving cells of YAP-mediated antiviral response. We observed altered phosphorylation of Hippo signaling-related proteins in response to SARS-CoV-2 infection. Recently Garcia et al. showed that activation of the Hippo signaling pathway during SARS-CoV-2 infection contributed to host antiviral response. 118 Furthermore, using a pharmacological inhibitor of YAP, the infection mediated by SARS-CoV-2 was significantly reduced. This study corroborates our findings on the importance of Hippo signaling in the context of SARS-CoV-2 infection and its regulation of the innate immune response. 119 In the current study, we observed hypophosphorylation of upstream kinases involved in Hippo signaling, including STK3 (S386, T354, S370/S371), and STK4 (S40, T41, S43). Furthermore, we observed YAP1 (S61) and WWTR1 (S62, T67) phosphorylation levels progressively decreasing and the lowest at 48 hpi. In addition, validation using immunoblot corroborated the hypophosphorylation of YAP1 at S61. Furthermore, we observed hypophosphorylation of PRKAA1 or AMPK (S486, T490), an upstream kinase that phosphorylates YAP1 at S61. 120 We also found other regulators of the Hippo signaling pathway such as KIBRA or WWC1 (S947), FRMD6 (S375), NF2 (S518), and ZO-2 or TJP2 (S966, S1027) to be hypophosphorylated. Signaling crosstalk between the Hippo and TGFb pathway is known. 121 In our study, TGFb receptor TGFBR2 was progressively hyperphosphorylated, reaching maximal hyperphosphorylation at 48 hpi. LATS-mediated phosphorylation of YAP results in its interaction and binding to 14-3-3, leading to cytoplasmic retention. 122 We identified hypoacetylation of both YWHAB (14-3-3b) at K13 and YWHAZ (14-3-3z) at K11. Whether altered acetylation directly affects YAP localization and activity 123 remains to be determined.
Although it is well known that metabolic cues can control YAP/TAZ regulation, 124,125 increasing evidence suggests YAP influences cellular metabolism. YAP has been demonstrated as a central regulator of glucose metabolism, 126 mediating the regulation of glucose transporter 1 (glut1) in zebrafish. Yap mutant zebrafish were reported to have impaired glucose uptake, nucleotide synthesis and glucose tolerance in adults. 127 Our study indicates high levels of glycolysis intermediates, including glucose 6-phosphate, fructose 6-phosphate, dihydroxyacetone phosphate, and lactate at 48 hpi. These findings likely suggest that YAP and, thereby, hippo signaling could profoundly impact glucose metabolism during SARS-CoV-2 infection. We further observed increased alpha-ketoglutarate, proline, and uracil levels in response to SARS-CoV-2 infection and likely associate this observation with dysregulated YAP/Hippo signaling as YAP has been found previously to regulate fatty acid oxidation and amino acid metabolism. 128 YAP inhibition in the skeletal muscles decreased levels of various metabolites, including undecanoic acid, capric acid, 2-octanoic acid, 2-oxoglutaric acid (alpha-ketoglutarate), amino acids-lysine, serine, proline, aspartate and uracil intermediate, 3-ureidopropionic acid. YAP has also been shown to regulate glutamine metabolism to stimulate nucleotide biosynthesis through increased expression of glutamine synthetase in zebrafish. 129  Furthermore, the DNA damage response pathway was enriched across multiple omics datasets in response to SARS-CoV-2 infection. We identified p53 S314/315, not S15, and the phosphorylation level was downregulated (log 2 FC À0.5). Immunoblot analysis revealed decreased phosphorylation of p53 S15 and corresponded to the reduced abundance of total p53 levels indicating CoV-2 infection potentially mediates p53 degradation ( Figure 6H). At later time points of infection, we observed hyperphosphorylation of Ser/Thr kinase PRKDC (also known as DNAPK) at multiple sites. DNAPK acts as a molecular sensor for DNA damage and is also known to regulate DNA virus-mediated innate immune response. We identified hyperphosphorylation of several DNAPK sites including S3205, an ATM target phosphorylated either through autophosphorylation 130 or through PLK1 131 and S2672, which is part of the ABCDE cluster, which promotes DNA end processing, and autophosphorylation of this cluster is required to initiate DNA damage repair. 130 Regulation of DNA damage response mechanism by DNAPK is mediated by phosphorylation of S139/S140 of histone variant H2AX/H2AFX. 132,133 Concurrently, we observed increased phosphorylation of gH2AX at S140 from 12 hpi with no changes in the total protein levels. We further confirmed our results through immunoblot analysis ( Figure 6H).YAP1/Hippo signaling has also been known to cross-talk with the DNA damage response pathways, [134][135][136][137][138] and this could help explain why DNA damage/repair pathways were enriched in Calu-3 on SARS-CoV-2 infection. Cross talk between YAP and P53 via a SIRT-1 mediated deacetylation mechanism is involved in the control of cell cycle in A549 lung cells. 139 In addition, TAZ (or WWTR1) has been shown to reduce P53 activity through p300-mediated acetylation negatively. Notably, ANP32E has previously been reported to remove H2A.Z from DNA double-strand breaks, promoting DNA repair and nucleosome organization. 140 All these findings suggest a deeper connection between DNA damage and YAP/TAZ signaling.
Ubiquitination is essential in regulating the innate immune response. E3 ligases, such as TRIM25 and RIPLET, have been shown to mediate RIG-I ubiquitination and type I IFN induction. 138 We observed several Ubiquitin (Ub) E2 and E3 ligases such as UBE2C, UBE2J1, UBE2S, PJA2, RNF5, RNF14, and RNF115, among others, to be significantly downregulated as early as 12 hpi ( Figure S6I). To assess if this affected the status of protein ubiquitination, the pan-Ub profile was explored using Anti-Ub FK2, which detects K29-, K48-, and K63-linked mono-and polyubiquitinylated proteins. Similar to the expression profiles of E3 ligases, a significant decrease in ubiquitination was observed, especially at 24 and 48 hpi ( Figure 6I), suggesting CoV-2mediated hijacking of the host ubiquitin system to maximize their own survival likely through mechanisms adopted by DNA tumor viruses and Adenoviruses. A recent study aimed to identify substrates of the Mpro from three coronaviruses revealed many E3 ubiquitin ligases were cleaved by the SARS-CoV-2 main protease Mpros including RNF20, ITCH, UBE3A, confirming our findings that albeit increased expression at the transcript level, profound downregulation of E3 ligases owing to degradation was essential to counteract the host innate immune response. 138,141 In addition, Coronaviruses encode papain-like protease (PLP) that acts as a cysteine protease as well as possesses intrinsic deubiquitinating and deISGylating activities required for viral replication and the evasion of host responses.
We observed time-dependent dephosphorylation of several proteins involved in the cell cycle at multiple sites, including Rb, which acts as a key driver of cell cycle regulation. A time-dependent decrease in phosphorylation at S249, T821, and T826 is known to be mediated by members of type 1 serine/threonine protein phosphatases (PP1) 142 was observed. Earlier studies indicate that endogenous Rb is dephosphorylated on T821 when cells undergo apoptosis. 143 Changes in the phosphorylation dynamics were accompanied by downregulation at the proteome level for most proteins involved in regulating cell cycle processes. (Figure S6J). iScience Article In addition to the various cellular processes regulated at multiple levels, changes in the expression dynamics as well as alterations in the phosphorylation/acetylation status of genes involved in energy metabolism specifically glycolysis, TCA cycle and amino acid metabolism was evaluated ( Figure S7). We observed extensive transcriptional changes of genes belonging to these pathways, however the impact on the protein expression seemed to be minimal. Notably, altered PTM status indicating hyperphosphorylation and hypoacetylation especially at later time points of infection of several glycolysis and TCA cycle genes was observed. The altered PTM dynamics correlate with increased viral gene expression signifying their role in the metabolic reprogramming of human airway epithelial cells in response to SARS-CoV-2 infection. 78,94,100 However, further studies are warranted to understand the precise role of these post-translational modifications in the regulation of these metabolic pathways.

Conclusions
We describe a systems-level approach to discover signaling pathways modulated on infection of Calu-3 cells with SARS-CoV-2 S15 Trondheim strain belonging to the wave 1 isolates. In-depth systems-wide and time-resolved characterization of the host and viral changes throughout productive infection revealed temporal changes in the transcriptome, proteome, phosphoproteome, acetylome, and exometabolome. Although other studies have extensively explored SARS-CoV-2 mediated alterations in the proteome and phosphoproteome of infected cells, 27,29,31,53,54 to our knowledge, this study is the first to perform comprehensive profiling, including exploring the acetylome and exometablome dynamics in response to SARS-CoV-2 infection. Multiple genes and functional pathways identified in our data were previously reported to promote SARS-CoV-2 mediated pathogenesis, 18,20,22,31,44,46,53,144 validating the rigor of our approach and providing further support for the role of these specific host factors. In concordance with earlier studies, we too observed a robust induction of the innate immune response accompanied by increased interferon signaling, increased expression of ISGs, TLRs, and cytokines and chemokines. 63,145,146 Furthermore, our results are concordant with studies showing that the pulmonary viral load positively correlates with the expression of ISGs, with reduced risk and disease severity. [147][148][149] Integrated analysis revealed changes in Hippo signaling, DNA damage response, ubiquitination, and spliceosome machinery pathways across omics datasets. In response to infection, several novel phosphorylations and acetylation sites were identified on vital host proteins. We also provide evidence of host-mediated PTM modification of viral proteins, including time-dependent alterations in the acetylation levels of nucleoprotein for the first time. Altered PTM levels have been demonstrated to influence virus-mediated control of host cellular dynamics 78,150,151 and the role of novel acetylation sites warrants further investigations. Altered release of TCA cycle intermediates with increased secretion of organic acids especially lactate strongly indicates enhanced aerobic glycolysis or Warburg-like effect that likely supports SARS-CoV-2 replication in airway epithelial cells. 98 Furthermore, secretion of ketone bodies such as BHB could likely be a counteractive measure to restrict SARS-CoV2 infection. 152 Finally, the data from the study shows evidence of potential crosstalk between Hippo signaling, DNA damage response, glucose metabolism and post-translational modifications. However, the site-specific alterations in the PTM dynamics and it effect on the signaling crosstalks in response to SARS-CoV2 infection will need to be deciphered in greater detail by future investigations.
Although we observed a small overlap with other in vitro Calu-3 infections studies, a comparison of both the transcriptome and proteome of uninfected Calu-3 cells revealed a high degree of correlation across datasets. Furthermore, differences were also observed with several published multi-omics studies on SARS-COV-2 infection of other cell types. 20,27,29,31,53 Overall, our findings demonstrate that the variance observed is unlikely because of pronounced expression program differences of uninfected Calu-3 cells rather influenced by protocols used for infection, sample preparation, and data acquisition strategies. Nevertheless, the processes and functional pathways identified in our data have been previously reported to promote SARS-CoV-2-mediated pathogenesis. A detailed compilation of the information would provide a holistic overview on the discreet molecular level regulation. Specifically, the measurement of PTM dynamics across the time course of infection can be leveraged to investigate immune evasion mechanisms and enhanced transmission. Furthermore, the findings from this study and other studies could serve as a resource for future investigations on newer SARS-COV-2 strains and a point of reference during the evolution of existing strains. Such an approach could aid in future pandemic preparedness in terms of finding novel drug targets.

Limitations of the study
We performed a temporal multi-OMICs study of host response to SARS-CoV-2 infection using a human lung epithelial cellline model. iScience Article valuable findings in terms host-directed therapeutic targets, one of the key limitations of this study is that the viral strain thata we used is one of vast number of SARS-CoV-2 strains being studied, and is from the early days of the pandemic. The current SARS-CoV-2 strain in circulation could be very different from the strain that we studied and hence more multi-OMICS studies are required to find common targets for host-directed therapy.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:    158 The resulting gene lists were annotated and filtered for significantly differentially upand down-regulated genes.

Sample preparation for mass spectrometry analysis
The cell lysates were sonicated using a probe sonicator (Branson Digital Sonifier, Branson Ultrasonics Corporation, USA) on ice for 5-10 min (20% amplitude, 10 cycles). The lysates were then heated at 95 C and centrifuged at 12,000 rpm each for 10 min, respectively. The concentration of proteins was determined by the Bicinchoninic acid (BCA) assay (Pierce, Thermo Fisher Scientific, Waltham, MA). A total of 300 mg proteins per sample were used for downstream processing using the filter-aided sample preparation (FASP) method. 159  The peaks for metabolite were confirmed using Mass spectrometry metabolite library kit MSMLS-1EA (Ms Sigma Aldrich supplied by IROA Technologies). Peak integration was done with the TraceFinder 4.1 software (M/s Thermo Fischer Scientific, Waltham, MA, USA). The peak area data was exported as an Excel file for further analysis. Data quality was monitored throughout the run using pooled healthy human serum as Quality Control (QC) processed and extracted in the same way as unknown samples and interspersed throughout the run as every 10 th sample. After integrating QC data with TraceFinder 4.1, each detected metabolite was checked and %RSD was calculated, and the acceptance limit was set %20%. Blank samples for carryover were injected after every fifth unknown samples to monitor the metabolites' carryover effect and calculated against the mean QC area, and the acceptance % carryover limit was set %20% for each ll OPEN ACCESS iScience 26, 105895, January 20, 2023 iScience Article metabolite. Background noise blank (first solvent blank of the run) was injected, and % background noise was calculated against the mean QC area, and the acceptance % background noise limit was set %20% for each metabolite.

Cytokine array analysis
The cell culture supernatants from mock-infected and hCoV-19-S15 strain infected CALU3 cells were collected at indicated time points (0, 3, 6-, 12-, 24-, and 48-h post-infection (hpi) and subjected to centrifugation for 5 min at 14,000 rpm. Cytokines were analyzed using Proteome Profiler Human Cytokine Array Kit (R&D Systems) according to the manufacturer's instructions. The membranes were exposed to X-ray films, scanned, and analyzed using ImageJ software (NIH). Thefold change was calculated in comparison to the mock-infected sample.

Western blot analysis
Protein samples were run on pre-cast NuPAGE Bis-Tris gels (Invitrogen) with 1x MOPS buffer (Invitrogen) and transferred on nitrocellulose membranes using the iBlot2 Gel Transfer Device (Invitrogen). Membranes were washed in Tris-Buffered Saline with 0.1% Tween-X100 (TBS-T) and blocked with TBS-T containing 5% BSA (BSA, Sigma-Aldrich). Membranes were incubated with primary antibodies at 4 C overnight. The following primary antibodies were used:,anti-b-Actin (1:

QUANTIFICATION AND STATISTICAL ANALYSIS
The transcriptome FPKM data was subjected to differential analysis. Differentials at each time point were determined using a pvalue cut-off of <-0.01 and log2(fold-change) cut-off of +-1.5. Protein abundances and phosphosite abundances across multiple replicates were subjected to quantile normalization and differential expression using limma (v3.44.3) in R (v4.0.2 https://www.r-project.org/). Log2 fold changes and pvalues were calculated using the proDA (v1.2.0) package for R. Differentials at each time point were determined using a pvalue cut-off of <-0.05and log2(fold-change) cut-off of +-1.5. For the Acetylome data, abundance ratios obtained from Proteome Discoverer were considered for further analysis. The proteomics and phosphoproteomics datasets generated in the current study were compared with corresponding datasets from Thorne et al., 27 Hekman et al., 31 Stukalov et al. 29 , and Grossegesse et al. 53 Gene ontology and Pathway analysis for transcriptome total proteome, phosphoproteome and acetylome data were carried out using enrichR (v3.0) R package. The enrichment databases consisted of "GO_Biological_Pro-cess_2015","GO_Cellular_Component_2015", and "Reactome_2015" and significant enrichment used a pvalue cut-off of 0.05. A list of ISGs was compiled based on previous literature. 27,161 The gene sets pertaining to cytokines, chemokines, Hippo signaling, regulation of Hippo signaling, DNA damage response, DNA repair, protein ubiquitination, Regulation of protein mono and polyubiquitination, Alternative splicing by spliceosome, and Regulation of mRNA splicing via spliceosome were obtained from 112 MSigDB (v7.5.1), and comparisons were carried out. The differentials from the phosphoproteome data were subjected to gene set enrichment analysis (GSEA) against the PTM signatures database (PTMsigDB), 162 which consists of modification site-specific signatures of perturbations kinase activities, and signaling pathways curated from the literature.