key: cord-0842867-zf096duh
authors: Stukalov, Alexey; Girault, Virginie; Grass, Vincent; Karayel, Ozge; Bergant, Valter; Urban, Christian; Haas, Darya A.; Huang, Yiqi; Oubraham, Lila; Wang, Anqi; Hamad, M. Sabri; Piras, Antonio; Hansen, Fynn M.; Tanzer, Maria C.; Paron, Igor; Zinzula, Luca; Enghleitner, Thomas; Reinecke, Maria; Lavacca, Teresa M.; Ehmann, Rosina; Wölfel, Roman; Jores, Jörg; Kuster, Bernhard; Protzer, Ulrike; Rad, Roland; Ziebuhr, John; Thiel, Volker; Scaturro, Pietro; Mann, Matthias; Pichlmair, Andreas
title: Multilevel proteomics reveals host perturbations by SARS-CoV-2 and SARS-CoV
date: 2021-03-15
journal: bioRxiv
DOI: 10.1101/2020.06.17.156455
sha: 006c69d5d26577104425656690c84d1c021dac20
doc_id: 842867
cord_uid: zf096duh

The global emergence of SARS-CoV-2 urgently requires an in-depth understanding of molecular functions of viral proteins and their interactions with the host proteome. Several individual omics studies have extended our knowledge of COVID-19 pathophysiology1–10. Integration of such datasets to obtain a holistic view of virus-host interactions and to define the pathogenic properties of SARS-CoV-2 is limited by the heterogeneity of the experimental systems. We therefore conducted a concurrent multi-omics study of SARS-CoV-2 and SARS-CoV. Using state-of-the-art proteomics, we profiled the interactome of both viruses, as well as their influence on transcriptome, proteome, ubiquitinome and phosphoproteome in a lung-derived human cell line. Projecting these data onto the global network of cellular interactions revealed crosstalk between the perturbations taking place upon SARS-CoV-2 and SARS-CoV infections at different layers and identified unique and common molecular mechanisms of these closely related coronaviruses. The TGF-β pathway, known for its involvement in tissue fibrosis, was specifically dysregulated by SARS-CoV-2 ORF8 and autophagy by SARS-CoV-2 ORF3. The extensive dataset (available at https://covinet.innatelab.org) highlights many hotspots that can be targeted by existing drugs and it can guide rational design of virus- and host-directed therapies, which we exemplify by identifying kinase and MMPs inhibitors with potent antiviral effects against SARS-CoV-2.

). In particular, we discovered that SARS-CoV-2 targets a number of key innate immunity regulators Table 2) . For instance, we recapitulated the known interaction between SARS-CoV NSP2 and prohibitins (PHB, PHB2) 12 but this was not conserved in SARS-CoV-2 NSP2, suggesting that the two viruses differ in their ability to modulate mitochondrial function and homeostasis through NSP2 (Extended data Fig. 2a) . The exclusive interaction of SARS-CoV-2 ORF8 with the TGFB1-LTBP1 complex is another interaction potentially explaining the differences in pathogenicity of the two viruses (Extended data Fig. 1f, 2b) . Notably, disbalanced TGF-β signaling has been linked to lung fibrosis and oedema, a common complication of severe pulmonary diseases including COVID-19 [13] [14] [15] [16] . NSP6. Intriguingly, cholesterol metabolism was recently shown to be implicated in SARS-CoV-2 replication and suggested as a promising target for drug development [19] [20] [21] . Beside perturbations at the pathway level, viral proteins specifically modulated single host proteins, possibly explaining more distinct molecular mechanisms involved in viral protein function.

Focusing on the 180 most affected host proteins, we identified RCOR3, a putative transcriptional corepressor, as strongly upregulated by NSP4 of both viruses (Extended data Fig. 2d, 3a) . Remarkably, the apolipoprotein B (APOB) was substantially regulated by ORF3 and NSP1 of SARS-CoV-2, suggesting its importance for SARS-CoV-2 biology (Extended data Fig. 3b ).

While interactome and effectome provide in-depth information on the activity of individual viral proteins, we wished to directly study their concerted activities in the context of viral infection. To this end, we infected ACE2-expressing A549 cells (Extended data Fig. 4a, b) with SARS-CoV-2 and SARS-CoV, and profiled the impact of viral infection on mRNA expression, protein abundance, ubiquitination and phosphorylation in a time-resolved manner ( Figure 2 a- b) .

In line with previous reports 9, 22 , both SARS-CoV-2 and SARS-CoV share the ability to down-regulate type-I interferon response and activate a pro-inflammatory signature at transcriptome and proteome levels (Figure 2a Tables 5, 6 ). Among these downregulated proteins, EDF1 has a pivotal role in the maintenance of endothelial integrity and may be a link to endothelial dysfunctions described for 25 . Profound regulation of cellular signaling pathways was also observed at the phosphoproteomic level: among 16 399 total quantified phosphorylation sites, 4 643 showed significant changes after SARS-CoV-2 or SARS-CoV infection (Extended data Fig. 5b , c; Supplementary Table 7) .

Highly regulated sites were identified for the proteins of the MAPK pathways (e.g. Tables 7 -8) . Notably, only SARS-CoV-2 but not SARS-CoV led to phosphorylation of the antiviral kinase EIF2AK2/PKR at the critical regulatory residue S33 28 . This differential activation of EIF2AK2/PKR could contribute to the difference in growth kinetics of the two SARS viruses (Supplementary Table 4 , 5).

Our data clearly point to an interplay of phosphorylation and ubiquitination patterns on individual host proteins. EGFR, for instance, showed increased ubiquitination on six lysine residues at 24 hours post-infection (h.p.i.) accompanied by increased phosphorylation of T693, S695 and S991 after 24 and 36 hours (Figure 2e , f). Ubiquitination of all six lysine residues on EGFR was more pronounced upon SARS-CoV-2 infection. Moreover, vimentin, a central co-factor for coronavirus entry 29 and pathogenicity 30,31 , displayed distinct phosphorylation and ubiquitination patterns on several sites early (e.g. S420) or late (e.g.

S56, S72, K334) in infection (Extended data Fig. 6a, b) . These discoveries underscore the value of testing different post-translational modifications simultaneously and suggest a concerted engagement of regulatory machineries to modify target protein's functions and abundance.

The majority of viral proteins were also post-translationally modified. Of the 27 detected SARS coronavirus proteins, 21 were ubiquitinated, among which N, S, NSP2, and NSP3

were the most frequently modified proteins in both viruses (Extended data Fig. 6c , Supplementary Table 6 ). Many of these ubiquitination sites were shared between the two viruses. Around half of the sites specifically regulated in either of the two viruses were conserved but differentially ubiquitinated, while the other half was encoded by either of the two pathogens, indicating that such acquired adaptations are also post-translationally modified and could recruit cellular proteins with appropriate functions (Figure 3a ). Our interactome data identified several host E3 ligases (e.g. SARS-CoV-2 ORF3 with TRIM47, WWP1/2, STUB1; M and TRIM7; NSP13 and RING1) and deubiquitinating enzymes (e.g. 

Our unified experimental design in a syngeneic system permitted direct time-resolved comparison of SARS-CoV-2 and SARS-CoV infection across different levels. Integrative pathway enrichment analysis demonstrated that both viruses largely perturb the same cellular processes at multiple levels albeit with varying temporal patterns (Extended data Fig. 7a ). Transcriptional downregulation of proteins involved in tau-protein kinase activity and iron ions sequestration at 6 h.p.i., for instance, was followed by a decrease in protein abundance after 12 h.p.i. (Supplementary Table 8) . RHO GTPase activation, mRNA processing and role of ABL in ROBO-SLIT signaling appeared to be regulated mostly through phosphorylation (Extended data Fig. 7a ). In contrast, processes connected to cellular integrity such as the formation of senescence-associated heterochromatin foci, apoptosisinduced DNA fragmentation and amino acid transport across the plasma membrane were modulated through concomitant phosphorylation and ubiquitination events, providing insights into the molecular relationships of these post-translational modifications. Ion transporters, especially the SLC12 family (cation-coupled chloride cotransporters), previously identified as cellular factors in pulmonary inflammation 34 , were also regulated at multiple levels, evidenced by reduced protein abundance as well as differential posttranslational modifications (Extended data Fig. 7a ).

The pathway enrichment analysis provided a global and comprehensive picture of how SARS-CoV-2 and SARS-CoV affect the host. We next applied an automated approach to systematically explore the underlying molecular mechanisms contained in the viral interactome and effectome data. We mapped the measured interactions and effects of each viral protein onto the global network of cellular interactions 35 and applied a network diffusion approach 36 ( Figure 4a ). Such analysis utilizes known cellular protein-protein interactions, signaling and regulation events to identify connection points between the interactors of the viral protein and the proteins affected by its expression (Extended data Fig.   1b, 2d, Supplementary Tables 2, 3 ). The connections inferred from the real data were significantly shorter than for randomized data, confirming both the relevance of the approach and the data quality (Extended data Fig. 8a, b) . Amongst many other findings, this approach pointed towards the potential mechanisms of autophagy regulation by ORF3 and NSP6; the modulation of innate immunity by M, ORF3 and ORF7b; and the Integrin-TGF- This inhibition may be due to the interaction of the ORF3 protein with the HOPS complex (VPS11, , which is essential for autophagosome-lysosome fusion, as well as by the differential phosphorylation of regulatory sites (e.g. on TSC2, mTORC1 complex, ULK1, RPS6, SQSTM1) and ubiquitination of key components (MAP1LC3A, GABARAP(L2), VPS33A, VAMP8) ( Figure 4c , Extended data Fig. 8g ). This inhibition of autophagosome function may have direct consequences for protein degradation. The abundance of APOB, a protein degraded via autophagy 38 , was selectively increased after SARS-CoV-2 infection or expression of the SARS-CoV-2 ORF3 (Extended data Fig. 3b , 8h). Accumulating APOB levels could exacerbate the risk of arterial thrombosis 39 , one of the main complications contributing to lung, heart and kidney failure in COVID-19 patients 40 . The inhibition of the IFN-α/β response observed at transcriptional and proteome levels was similarly explained by the network diffusion analysis (Extended data Fig. 8i ), which implicated multiple proteins of SARS-CoV-2 in the disruption of antiviral immunity.

Additional experiments functionally corroborated the inhibition of IFN-α/β induction or signaling by ORF3, ORF6, ORF7a, ORF7b, ORF9b (Extended data Fig. 8j) . Upon virus infection, we observed the regulation of TGF-β and EGFR pathways modulating cell survival, motility and innate immune responses (Extended data Fig. 9a -d) . Specifically, our network diffusion analysis revealed a connection between the binding of the ORF8 and ORF3 proteins to TGF-β-associated factors (TGFB1, TGFB2, LTBP1, TGFBR2, FURIN, BAMBI), the differential expression of ECM regulators (FERMT2, CDH1) and the virusinduced upregulation of fibrinogens (FGA, FGB), fibronectin (FN1) and SERPINE1

(Extended data Fig. 9a , b) 41 . The increased phosphorylation of proteins involved in MAPK (e.g. SHC1-S139, SOS1-S1134/1229, JUN-S63/S73, MAPKAPK2-T334, p38-T180/Y182) and receptor tyrosine kinase signaling (e.g. phosphorylation of PI3K complex members, PDPK1 (S241) and RPS6KA1 (S380)) as well as a higher expression of JUN, FOS and EGR1 are further indicative of TGF-β and EGFR pathways regulation (Extended data Fig.   9a , c, d). In turn, TGF-β and EGFR signaling are known to be potentiated by integrin signaling and activation of YAP-dependent transcription 42 , which we observed to be regulated in a time-dependent manner upon SARS-CoV-2 infection (Extended data Fig. 9a ).

Besides promoting virus replication, activation of these pathways has been implicated in fibrosis 13-15 , one of the hallmarks of COVID-19 16 .

Taken together, the viral-host protein-protein interactions and pathway regulations observed at multiple levels identify potential vulnerability points of SARS-CoV and SARS-CoV-2 that we decided to target by well-characterized selective drugs for antiviral therapies. To test antiviral efficacy, we established time-lapse fluorescent microscopy of SARS-CoV-2 GFPreporter virus infection 43 . Inhibition of virus replication by IFN-α/β treatment corroborated previous conclusions that efficient SARS-CoV-2 replication involves an inactivation of this pathway at an early step and confirmed the reliability of this screening approach (Extended data Fig. 10a) This drug screen demonstrates the value of our combined dataset that profiles SARS-CoV-2 infection at multiple levels. We hope that further exploration of these rich data by the scientific community and additional studies of the interplay between different omics levels will substantially advance our molecular understanding of coronaviruses biology, including the pathogenicity associated with specific human coronaviruses, such as SARS-CoV-2 and SARS-CoV. Moreover, this resource, together with complementary approaches by the community 26,52-54 , will streamline the search for antiviral compounds and serve as a base for rational design of combination therapies that target the virus from multiple synergistic angles, thus potentiating the effect of individual drugs while minimizing potential sideeffects on healthy tissues. 2  p  r  o  t  e  i  n  i  n  t  e  r  a  c  t  i  o  n  m  a  p  r  e  v  e  a  l  s  t  a  r  g  e  t  s  f  o  r  d  r  u  g  r  e  p  u  r  p  o  s  i  n e t a l . 

1 . G o r d o n , D . E . e t a l . A S A R S - C o V -A S A R S - C o V - 2 - h o s t p r o x i m i t y i n t e r a c t o m e . b i o R x i v 2 0 2 0 . 0 9 . 0 3 . 2 8 2 1 0 3 ( 2 0 2 0 ) d o i : 1 0 . 1 1 0 1 / 2 0 2 0 . 0 9 . 0 3 . 2 8 2 1 0 3 . 6 . L a u r e n t , E . M . N . e t a l . G l o b a l B i o I D - b a s e d S A R S - C o V - 2 p r o t e i n s p r o x i m a l i n t e

HEK293T, A549, Vero E6 and HEK293-R1 cells and their respective culturing conditions were described previously 55 . All cell lines were tested to be mycoplasma-free. Expression constructs for C-terminal HA tagged viral ORFs were synthesized (Twist Bioscience and BioCat) and cloned into pWPI vector as described previously 56 with the following modifications: starting ATG codon was added, internal canonical splicing sites were replaced with synonymous mutations and C-terminal HA-tag, followed by amber stop codon, was added to individual viral open reading frames. C-terminally hemagglutinin(HA)tagged ACE2 sequence was amplified from an ACE2 expression vector (kindly provided by Stefan Pöhlmann) 57 into the lentiviral vector pWPI-puro. A549 cells were transduced twice, and ACE2-expressing A549 (A549-ACE2) cells were selected with puromycin. Lentiviruses production, transduction of cells and antibiotic selection were performed as described previously 52 . RNA-isolation (Macherey-Nagel NucleoSpin RNA plus), reverse transcription (TaKaRa Bio PrimeScript RT with gDNA eraser) and RT-qPCR (Thermo-Fisher Scientific PowerUp SYBR green) were performed as described previously 54 . RNA-isolation for NGS applications was performed according to manufacturer's protocol (Qiagen RNeasy mini kit, RNase free DNase set). For detection of protein abundance by western blotting, HA-HRP (Sigma-Aldrich; H6533; 1:2500 dilution), ACTB-HRP (Santa Cruz; sc-47778; 1:5000 dilution), MAP1LC3B (Cell Signaling; 3868; 1:1000 dilution), MAVS (Cell Signaling; 3993; 1:1000 dilution), HSPA1A (Cell Signaling; 4873; 1:1000 dilution), TGFβ (Cell Signaling; 3711; 1:1000 dilution), phospho-p38 (T180/Y182) (Cell Signaling; 4511; 1:1000 dilution), p38 (Cell Signaling; 8690; 1:1000 dilution) and SARS-CoV-2/SARS-CoV N protein (Sino Biological; 40143-MM05; 1:1000 dilution) antibodies were used Secondary antibodies detecting mouse (Cell Signaling; 7076; 1:5000 dilution/Jackson ImmunoResearch; 115-035-003; 1:5000 dilution), rat (Invitrogen; 31470; 1:5000 dilution), and rabbit IgG (Cell Signaling; 7074; 1:5000 dilution) were horseradish peroxidase (HRP)coupled. For AP-MS and AP-WB applications, HA-beads (Sigma-Aldrich and Thermo Fisher Scientific) and Streptactin II beads (IBA Lifesciences) were used. WB imaging was performed as described previously 58 . For the stimulation of cells in the reporter assay, recombinant human IFN-α was a kind gift from Peter Stäheli, recombinant human IFN-γ were purchased from PeproTech and IVT4 was produced as described before 59 . All compounds tested during the viral inhibitor assay are listed in Supplementary Table 9 .

SARS-CoV-Frankfurt-1, SARS-CoV-2-MUC-IMB-1 and SARS-CoV-2-GFP strains 43 were produced by infecting Vero E6 cells cultured in DMEM medium (10% FCS, 100 ug/ml Streptomycin, 100 IU/ml Penicillin) for 2 days (MOI 0.01). Viral stock was harvested and spun twice (1000g/10min) before storage at -80°C. Titer of viral stock was determined by plaque assay. Confluent monolayers of VeroE6 cells were infected with serial five-fold dilutions of virus supernatants for 1 hour at 37 °C. The inoculum was removed and replaced with serum-free MEM (Gibco, Life Technologies) containing 0.5% carboxymethylcellulose (Sigma-Aldrich). Two days post-infection, cells were fixed for 20 minutes at room temperature with formaldehyde directly added to the medium to a final concentration of 5%. Fixed cells were washed extensively with PBS before staining with H2O containing 1% crystal violet and 10% ethanol for 20 minutes. After rinsing with PBS, the number of plaques was counted and the virus titer was calculated.

A549-ACE2 cells were infected with either SARS-CoV-Frankfurt-1 or SARS-CoV-2-MUC-IMB-1 strains (MOI 2) for the subsequent experiments. At each time point, the samples were washed once with 1x TBS buffer and harvested in SDC lysis buffer (100 mM Tris HCl pH 8.5; 4% SDC) or 1x SSB lysis buffer (62.5 mM Tris HCl pH 6.8; 2% SDS; 10% glycerol; 50 mM DTT; 0.01% bromophenol blue) or RLT (Qiagen) for proteome-phosphoproteomeubiquitinome, western blot, and transcriptome analyses, respectively. The samples were heat-inactivated and frozen at -80°C until further processing, as described in the following sections.

To determine the interactomes of SARS-CoV-2 and SARS-CoV and the interactomes of an accessory protein (encoded by ORF4/ORF4a of HCoV-229E or ORF3 of HCoV-NL63) that presumably represents a homolog of the ORF3 and ORF3a proteins of SARS-CoV-2 and SARS-CoV, respectively, four replicate affinity purifications were performed for each HAtagged viral protein. A549 cells (6×10 6 cells per 15-cm dish) were transduced with lentiviral vectors encoding HA-tagged SARS-CoV-2, SARS-CoV or HCoV-229E/NL63 proteins and protein lysates were prepared from cells harvested three days post-transduction. Cell pellets of two 15-cm dishes were lysed in lysis buffer (50 mM Tris-HCl pH 7.5, 100 mM NaCl, 1.5 mM MgCl 2 , 0.2% (v/v) NP-40, 5% (v/v) glycerol, cOmplete protease inhibitor cocktail (Roche), 0.5% (v/v) 750 U/µl Sm DNAse) and sonicated (5 min, 4°C, 30 sec on, 30 sec off, low settings; Bioruptor, Diagenode SA). Following normalization of protein concentrations of cleared lysates, virus protein-bound host proteins were enriched by adding 50 µl anti-HAagarose slurry (Sigma-Aldrich, A2095) with constant agitation for 3 hours at 4°C. Nonspecifically bound proteins were removed by four subsequent washes with lysis buffer followed by three detergent-removal steps with washing buffer (50 mM Tris-HCl pH 7.5, 100 mM NaCl, 1.5 mM MgCl 2 , 5% (v/v) glycerol). Enriched proteins were denatured, reduced, alkylated and digested by addition of 200 µl digestion buffer (0.6 M guanidinium chloride, 1 mM TCEP, 4 mM CAA, 100 mM Tris-HCl pH 8, 0.5 µg LysC (WAKO Chemicals), 0.5 µg trypsin (Promega) at 30°C overnight. Peptide purification on StageTips with three layers of C18 Empore filter discs (3M) and subsequent mass spectrometry analysis was performed as described previously 55,56 . Briefly, purified peptides were loaded onto a 20 cm reverse-phase analytical column (75 µm diameter; ReproSil-Pur C18-AQ 1.9 µm resin; Dr. Maisch) and separated using an EASY-nLC 1200 system (Thermo Fisher Scientific). A binary buffer system consisting of buffer A (0.1% formic acid in H 2 O) and buffer B (80% acetonitrile, 0.1% formic acid in H 2 O) with a 90 min gradient (5-30% buffer B (65 min), 30-95% buffer B (10 min), wash out at 95% buffer B (5 min), decreased to 5% buffer B (5 min), and 5% buffer B (5 min)) was used at a flow rate of 300 nl per min. Eluting peptides were directly analysed on a Q-Exactive HF mass spectrometer (Thermo Fisher Scientific). Data-dependent acquisition included repeating cycles of one MS1 full scan (300-1 650 m/z, R = 60 000 at 200 m/z) at an ion target of 3×10 6 , followed by 15 MS2 scans of the highest abundant isolated and higher-energy collisional dissociation (HCD) fragmented peptide precursors (R = 15 000 at 200 m/z). For MS2 scans, collection of isolated peptide precursors was limited by an ion target of 1×10 5 and a maximum injection time of 25 ms. Isolation and fragmentation of the same peptide precursor was eliminated by dynamic exclusion for 20 s. The isolation window of the quadrupole was set to 1.4 m/z and HCD was set to a normalized collision energy of 27%.

For the determination of proteome changes in A549 cells expressing SARS-CoV-2, SARS-CoV or HCoV-229E/NL63 proteins, a fraction of 1×10 6 lentivirus-transduced cells from the affinity purification samples were lysed in guanidinium chloride buffer (6 M GdmCl, 10 mM TCEP, 40 mM CAA, 100 mM Tris-HCl pH 8), boiled at 95°C for 8 min and sonicated (10 min, 4°C, 30 sec on, 30 sec off, high settings). Protein concentrations of cleared lysates were normalized to 50 µg and proteins were pre-digested with 1 µg LysC at 37°C for 1 hour followed by a 1:10 dilution (100 mM Tris-HCl pH 8) and overnight digestion with 1 µg trypsin at 30°C. Peptide purification on StageTips with three layers of C18 Empore filter discs (3M) and subsequent mass spectrometry analysis was performed as described previously 55, 56 . Briefly, 300 ng of purified peptides were loaded onto a 50 cm reversed phase column (75 μ m inner diameter, packed in house with ReproSil-Pur C18-AQ 1.9 μ m resin [Dr. Maisch GmbH]). The column temperature was maintained at 60°C using a homemade column oven. A binary buffer system, consisting of buffer A (0.1% formic acid (FA)) and buffer B (80% ACN, 0.1% FA), was used for peptide separation, at a flow rate of 300 nl/min. An EASY-nLC 1200 system (Thermo Fisher Scientific), directly coupled online with the mass spectrometer (Q Exactive HF-X, Thermo Fisher Scientific) via a nanoelectrospray source, was employed for nano-flow liquid chromatography. Peptides were eluted by a linear 80 min gradient from 5% to 30% buffer B (0.1% v/v formic acid, 80% v/v acetonitrile), followed by a 4 min increase to 60% B, a further 4 min increase to 95% B, a 4 min plateau phase at 95% B, a 4 min decrease to 5% B and a 4 min wash phase of 5% B. To acquire MS data, the data-independent acquisition (DIA) scan mode operated by the XCalibur software (Thermo Fisher) was used. DIA was performed with one full MS event followed by 33 MS/MS windows in one cycle resulting in a cycle time of 2.7 seconds. The full MS settings included an ion target value of 3×10 6 charges in the 300 -1 650 m/z range with a maximum injection time of 60 ms and a resolution of 120 000 at m/z 200. DIA precursor windows ranged from 300.5 m/z (lower boundary of first window) to 1 649.5 m/z (upper boundary of 33rd window). MS/MS settings included an ion target value of 3×10 6 charges for the precursor window with an Xcalibur-automated maximum injection time and a resolution of 30 000 at m/z 200.

To generate the proteome library for DIA measurements purified peptides from the first and the fourth replicates of all samples were pooled separately and 25 µg of peptides from each pool were fractionated into 24 fractions by high pH reversed-phase chromatography as described earlier 60 . During each separation, fractions were concatenated automatically by shifting the collection tube every 120 seconds. In total 48 fractions were dried in a vacuum centrifuge, resuspended in buffer A* (0.2% TFA, 2% ACN) and subsequently analyzed by a top12 data-dependent acquisition (DDA) scan mode using the same LC gradient and settings. The mass spectrometer was operated by the XCalibur software (Thermo Fisher). DDA scan settings on full MS level included an ion target value of 3×10 6 charges in the 300 -1 650 m/z range with a maximum injection time of 20 ms and a resolution of 60 000 at m/z 200. At the MS/MS level the target value was 10 5 charges with a maximum injection time of 60 ms and a resolution of 15 000 at m/z 200. For MS/MS events only, precursor ions with 2-5 charges that were not on the 20 s dynamic exclusion list were isolated in a 1.4 m/z window. Fragmentation was performed by higher-energy C-trap dissociation (HCD) with a normalized collision energy of 27eV.

Frozen lysates of infected A549-ACE2 cells harvested at 6, 12 and 24 hours (also 36 hours only in phosphoproteomics study) post-infection were thawed on ice, boiled for 5 min at 95°C and sonicated for 15 min (Branson Sonifierer). Protein concentrations were estimated by tryptophan assay 61 . To reduce and alkylate proteins, samples were incubated for 5 min at 45°C with TCEP (10 mM) and CAA (40 mM). Samples were digested overnight at 37°C using trypsin (1:100 w/w, enzyme/protein, Sigma-Aldrich) and LysC (1:100 w/w, enzyme/protein, Wako). For proteome analysis, 10 µg of peptide material were desalted using SDB-RPS StageTips (Empore) 61 . Briefly, samples were diluted with 1% TFA in isopropanol to a final volume of 200 µl and loaded onto StageTips, subsequently washed with 200 µl of 1% TFA in isopropanol and 200 µl 0.2% TFA/ 2% ACN. Peptides were eluted with 75 µl of 1.25% Ammonium hydroxide (NH4OH) in 80% ACN and dried using a SpeedVac centrifuge (Eppendorf, Concentrator plus). They were resuspended in buffer A* (0.2% TFA/ 2% ACN) prior to LC-MS/MS analysis. Peptide concentrations were measured optically at 280 nm (Nanodrop 2000, Thermo Scientific) and subsequently equalized using buffer A*. 1µg peptide was subjected to LC-MS/MS analysis.

The rest of the samples was four-fold diluted with 1% TFA in isopropanol and loaded onto SDB-RPS cartridges (Strata™-X-C, 30 mg/ 3 ml, Phenomenex Inc), pre-equilibrated with 4 ml 30% MeOH/1% TFA and washed with 4 ml 0.2% TFA. Samples were washed twice with 4 ml 1% TFA in isopropanol, once with 0.2% TFA/ 2% ACN and eluted twice with 2 ml 1.25% NH4OH/ 80% ACN. Eluted peptides were diluted with ddH 2 O to a final ACN concentration of 35%, snap frozen and lyophilized.

For phosphopeptide enrichment, lyophilized peptides were resuspended in 105 µl of equilibration buffer (1% TFA/ 80% ACN) and the peptide concentration was measured optically at 280nm (Nanodrop 2000, Thermo Scientific) and subsequently equalized using equilibration buffer. The AssayMAP Bravo robot (Agilent) performed the enrichment for phosphopeptides (150µg) by priming AssayMAP cartridges (packed with 5 µl Fe(III)-NTA) with 0.1% TFA in 99% ACN followed by equilibration in equilibration buffer and loading of peptides. Enriched phosphopeptides were eluted with 1 % Ammonium hydroxide, which was evaporated by Speedvac'ing samples for 20 minutes. Dried peptides were resuspended in 6 µl buffer A* and 5 µl was subjected to LC-MS/MS analysis.

For diGly peptide enrichment, lyophilized peptides were reconstituted in IAP buffer (50 mM MOPS, pH 7.2, 10 mM Na 2 HPO 4 , 50 mM NaCl) and the peptide concentration was estimated by tryptophan assay. K--GG remnant containing peptides were enriched using the PTMScan® Ubiquitin Remnant Motif (K--GG) Kit (Cell Signaling Technology). Crosslinking of antibodies to beads and subsequent immunopurification was performed with slight modifications as previously described 62 . Briefly, two vials of crosslinked beads were combined and equally split into 16 tubes (~31 µg of antibody per tube). Equal peptide amounts (600 µg) were added to crosslinked beads and the volume was adjusted with IAP buffer to 1 ml. After 1 hour of incubation at 4°C and gentle agitation, beads were washed twice with cold IAP and 5 times with cold ddH 2 O. Thereafter, peptides were eluted twice with 50 µl 0.15% TFA. Eluted peptides were desalted and dried as described for proteome analysis with the difference that 0.2% TFA instead of 1%TFA in isopropanol was used for the first wash. Eluted peptides were resuspended in 9 µl buffer A* and 4 µl was subjected to LC-MS/MS analysis.

Samples were loaded onto a 50 cm reversed phase column (75 μ m inner diameter, packed in house with ReproSil-Pur C18-AQ 1.9 μ m resin [Dr. Maisch GmbH]). The column temperature was maintained at 60°C using a homemade column oven. A binary buffer system, consisting of buffer A (0.1% formic acid (FA)) and buffer B (80% ACN plus 0.1% FA) was used for peptide separation, at a flow rate of 300 nl/min. An EASY-nLC 1 200 system (Thermo Fisher Scientific), directly coupled online with the mass spectrometer (Orbitrap Exploris 480, Thermo Fisher Scientific) via a nano-electrospray source, was employed for nano-flow liquid chromatography. The FAIMS device was placed between the nanoelectrospray source and the mass spectrometer and was used for measurements of the proteome and the PTM-library samples. Spray voltage was set to 2 650 V, RF level to 40 and heated capillary temperature to 275°C.

For proteome measurements we used a 100 min gradient starting at 5% buffer B followed by a stepwise increase to 30% in 80 min, 60% in 4 min and 95% in 4 min. The buffer B concentration stayed at 95% for 4 min, decreased to 5% in 4 min and stayed there for 4 min. The mass spectrometer was operated in data-independent mode (DIA) with a full scan range of 350-1 650 m/z at 120 000 resolution at 200 m/z, normalized automatic gain control (AGC) target of 300% and a maximum fill time of 28 ms. One full scan was followed by 22 windows with a resolution of 15 000, normalized automatic gain control (AGC) target of 1 000% and a maximum fill time of 25 ms in profile mode using positive polarity. Precursor ions were fragmented by higher-energy collisional dissociation (HCD) (NCE 30%). Each of the selected CVs (-40, -55 and -70) was applied to sequential survey scans and MS/MS scans; the MS/MS CV was always paired with the appropriate CV from the corresponding survey scan.

For phosphopeptide samples, 5 µl were loaded and eluted with a 70 min gradient starting at 3% buffer B followed by a stepwise increase to 19% in 40 min, 41% in 20 min, 90% in 5 min and 95% in 5 min. The mass spectrometer was operated in data-independent mode (DIA) with a full scan range of 300-1 400 m/z at 120 000 resolution at 200 m/z and a maximum fill time of 60 ms. One full scan was followed by 32 windows with a resolution of 30 000. Normalized automatic gain control (AGC) target and maximum fill time were set to 1 000% and 54 ms, respectively, in profile mode using positive polarity. Precursor ions were fragmented by higher-energy collisional dissociation (HCD) (NCE stepped 25-27.5-30%). For the library generation, we enriched A549 cell lysates for phosphopeptides and measured them with 7 different CV settings (-30, -40, -50, -60, -70, -80 or -90 V) using the same DIA method. The noted CVs were applied to the FAIMS electrodes throughout the analysis.

For the analysis of K--GG peptide samples, half of the samples were loaded. We used a 120 min gradient starting at 3% buffer B followed by a stepwise increase to 7% in 6 min, 20% in 49 min, 36% in 39 min, 45% in 10 min and 95% in 4 min. The buffer B concentration stayed at 95% for 4 min, decreased to 5% in 4 min and stayed there for 4 min. The mass spectrometer was operated in data-independent mode (DIA) with a full scan range of 300-1 350 m/z at 120 000 resolution at m/z 200, normalized automatic gain control (AGC) target of 300% and a maximum fill time of 20 ms. One full scan was followed by 46 windows with a resolution of 30 000. Normalized automatic gain control (AGC) target and maximum fill time were set to 1 000% and 54 ms, respectively, in profile mode using positive polarity. Precursor ions were fragmented by higher-energy collisional dissociation (HCD) (NCE 28%). For K--GG peptide library, we mixed the first replicate of each

Raw MS data files of AP-MS experiments conducted in DDA mode were processed with MaxQuant (version 1.6.14) using the standard settings and label-free quantification enabled (LFQ min ratio count 1, normalization type none, stabilize large LFQ ratios disabled). Spectra were searched against forward and reverse sequences of the reviewed human proteome including isoforms (UniprotKB, release 2019.10) and C-terminally HA-tagged SARS-CoV-2, SARS-CoV and HCoV proteins by the built-in Andromeda search engine 63 .

In-house Julia scripts 64 were used to define alternative protein groups: only the peptides identified in AP-MS samples were considered for being protein group-specific, protein groups that differed by the single specific peptide or had less than 25% different specific peptides were merged to extend the set of peptides used for protein group quantitation and reduce the number of protein isoform-specific interactions.

Spectronaut version 13 (Biognosys) with the default settings was used to generate the proteome libraries from DDA runs by combining files of respective fractionations using the human fasta file (Uniprot, 2019.10, 42 431 entries) and viral bait sequences. Proteome DIA files were analyzed using the proteome library with the default settings and disabled cross run normalization.

Spectronaut version 14 (Biognosys) 65 was used to generate the libraries and analyze all DIA files using the human fasta file (UniprotKB, release 2019.10) and sequences of SARS-CoV-2/SARS-CoV proteins (UniProt, release 2020.08). Orf1a polyprotein sequences were split into separate protein chains according to the cleavage positions specified in the UniProt. For the generation of the PTM-specific libraries, the DIA single CV runs were combined with the actual DIA runs and either phosphorylation at Serine/Threonine/Tyrosine or GlyGly at Lysine was added as variable modification to default settings. Maximum number of fragment ions per peptide was increased to 25. The proteome DIA files were analyzed using direct DIA approach with default settings and disabled cross run normalization. All PTM DIA files were analyzed using their respective hybrid library and either phosphorylation at Serine/Threonine/Tyrosine or GlyGly at Lysine was added as an additional variable modification to default settings with LOESS normalization and disabled PTM localization filter.

A collection of in-house Julia scripts 64 were used to process the elution group (EG) -level Spectronaut reports, identify PTMs and assign EG-level measurements to PTMs. The PTM was considered if at least once it was detected with ≥ 0.75 localization probability in EG with q-value ≤ 10 -3 . For further analysis of given PTM, only the measurements with ≥ 0.5 localization probability and EG q-value ≤ 10 -2 were used.

Unless otherwise specified, the bioinformatic analysis was done in R (version 3.6), Julia (version 1.5) and Python (version 3.8) using a collection of in-house scripts 64,66 .

MaxQuant and Spectronaut output files were imported into R using in-house maxquantUtils R package 67 . For all MS datasets, the Bayesian linear random effects models were used to define how the abundances of proteins change between the conditions. To specify and fit the models we employed msglm R package 68 , which utilizes rstan package (version 2.19) 69 for inferring the posterior distribution of the model parameters. In all the models, the effects corresponding to the experimental conditions have regularized horseshoe+ priors 70 , while the batch effects have normally distributed priors. Laplacian distribution was used to model the instrumental error of MS intensities. For each MS instrument used, the heteroscedastic intensities noise model was calibrated with the technical replicate MS data of the instrument. These data were also used to calibrate the logit-based model of missing MS data (the probability that the MS instrument will fail to identify the protein given its expected abundance in the sample). The model was fit using unnormalized MS intensities data. Instead of transforming the data by normalization, the inferred protein abundances were scaled by the normalization multiplier of each individual MS sample to match the expected MS intensity of that sample. This allows taking the signal-to-noise variation between the samples into account when fitting the model. Due to high computational intensity, the model was applied to each protein group separately. For all the models, 4 000 iterations (2 000 warmup + 2 000 sampling) of the No-U-Turn Markov Chain Monte Carlo were performed in 7 or 8 independent chains, every 4th sample was collected for posterior distribution of the model parameters. For estimating the statistical significance of protein abundance changes between the two experimental conditions, the p-value was defined as the probability that a random sample from the posterior distribution of the first condition would be smaller (or larger) than a random sample drawn from the second condition. No multiple hypothesis testing corrections were applied, since this is handled by the choice of the model priors.

The statistical model was applied directly to the MS1 intensities of protein group-specific LC peaks (evidence.txt table of MaxQuant output). In R GLM formula language, the model could be specified as

where APMS effect models the average shift of intensities in AP-MS data in comparison to full proteome samples, Bait is the average enrichment of a protein in AP-MS experiments of homologous proteins of both SARS-CoV and SARS-CoV-2, and Bait:Virus corresponds to the virus-specific changes in protein enrichment. MS1peak is the log-ratio between the intensity of a given peak and the total protein abundance (the peak is defined by its peptide sequence, PTMs and the charge; it is assumed that the peak ratios do not depend on experimental conditions 71 ), and MSbatch accounts for batch-specific variations of protein intensity. APMS, Bait and Bait:Virus effects were used to reconstruct the batch effect-free abundance of the protein in AP-MS samples. The modeling provided the enrichment estimates for each protein in each AP experiment. Specific AP-MS interactions had to pass the two tests. In the first test, the enrichment of the candidate protein in a given bait AP was compared against the background, which was dynamically defined for each interaction to contain the data from all other baits, where the abundance of the candidate was within 50%-90% percentile range (excluding top 10% baits from the background allowed the protein to be shared by a few baits in the resulting AP-MS network). The non-targeting control and Gaussia luciferase baits were always preserved in the background. Similarly, to filter out any potential side-effects of very high bait protein expression, the ORF3 homologs were always present in the background of M interactors and vice versa. To rule out the influence of the batch effects, the second test was applied. It was defined similarly to the first one, but the background was constrained to the baits of the same batch, and 40%-80% percentile range was used. In both tests, the protein has to be 4 times enriched against the background (16 times for highly expressed baits: ORF3, M, NSP13, NSP5, NSP6, ORF3a, ORF7b, ORF8b, HCoV-229E ORF4a) with the p-value ≤ 10 -3 .

Additionally, we excluded the proteins that, in the viral protein expression data, have shown upregulation, and their enrichment in AP-MS data was less than 16 times stronger than observed upregulation effects. Finally, to exclude the carryover of material between the samples sequentially analyzed by MS, we removed the putative interactors, which were also enriched at higher levels in the samples of the preceding bait, or the one before it. For the analysis of interaction specificity between the homologous viral proteins, we estimated the significance of interaction enrichment difference (corrected by the average difference between the enrichment of the shared interactors to adjust for the bait expression variation). Specific interactions have to be 4 times enriched in comparison to the homolog with p-value ≤ 10 -3 .

The statistical model of the viral protein overexpression data set was similar to AP-MS data, except that protein-level intensities provided by Spectronaut were used. The PCA analysis of the protein intensities has identified that the 2nd principal component is associated with the batch-dependent variations between the samples. To exclude their influence, this principal component was added to the experimental design matrix as an additional batch effect.

As with AP-MS data, the two statistical tests were used to identify the significantly regulated proteins (column "is_change" in Supplementary Table 3) . First, the absolute value of median log 2 -fold change of the protein abundance upon overexpression of a given viral protein in comparison to the background had to be above 1.0 with p-value ≤ 10 -3 . The background was individually defined for each analyzed protein. It was composed of experiments, where the abundance of given protein was within the 20%-80% percentile range of all measured samples. Second, the protein had to be significantly regulated (same median log 2 -fold change and p-value thresholds applied) against the batch-specific background (defined similarly to the global background, but using only the samples of the same batch).

An additional stringent criterion was applied to select the most significant changes (column "is_top_change" in Supplementary Table 3 ; Extended data Fig. 1i ).

For each protein we classified bait-induced changes as:

• "high" when |median log 2 fold-change| ≥ 1 and p-value ≤ 10 -10 both in background and batch comparisons • "medium" if 10 -10 < p-value ≤ 10 -4 with same fold-change requirement and • "low" if 10 -4 < p-value ≤ 10 -2 with same fold-change requirement, all other changes were considered non-significant.

We then required that "shared" top-regulated proteins should have exactly one pair of SARS-CoV-2/SARS-CoV "high"-or "medium"-significant homologous baits among the baits with either up-or downregulated changes and no other baits with significant changes of the same-type.

We further defined "SARS-CoV-2-specific" or "SARS-CoV-specific" top-regulated proteins to be the ones with exactly one "high"-significant change, and no other significant changes of the same sign. For "specific" hits we additionally required that in comparison of "high"significant bait to its homolog |median log 2 fold-change| ≥ 1 and p-value ≤ 10 -3 . When the homologous bait was missing (SARS-CoV-2 NSP1, SARS-CoV ORF8a and SARS-CoV ORF8b), we instead required that in the comparison of the "high"-significant change to the background |median log 2 fold-change| ≥ 1.5.

The resulting network of most affected proteins was imported and prepared for publication in Cytoscape v.3.8.1 72 .

Similarly to the AP-MS DDA data, the linear Bayesian model was applied to the elution group (EG) level intensities. To model the protein intensity, the following linear model (in R notation) was used:

where • after(t i ) effect corresponds to the protein abundance changes in mock-infected samples that happened between t i-1 and t i h.p.i. and it is applied to the modeled intensity at all time points starting from t i ;

• infection:after(t i ) (t i =6, 12, 24) is the common effect of SARS-CoV-2 & SARS-CoV infections occurred between t i-1 and t i ;

• CoV2:after(t i ) is the virus-specific effect within t i-1 and t i h.p.i. that is added to the log intensity for SARS-CoV-2-infected samples and subtracted from the intensity for SARS-CoV ones;

• EG is the elution group-specific shift in the measured log-intensities. The absolute value of median log 2 fold change between the conditions above 0.25 and the corresponding unadjusted p-value ≤ 10 -3 were used to define the significant changes at a given time point in comparison to mock infection. We also required that the protein group is quantified in at least two replicates of at least one of the compared conditions. Additionally, if for one of the viruses (e.g. SARS-CoV-2) only the less stringent condition (|median log 2 fold-change| ≥ 0.125, p-value ≤ 10 -2 ) was fulfilled, but the change was significant in the infection of the other virus (SARS-CoV), and the difference between the viruses was not significant, the observed changes were considered significant for both viruses.

The data from single-double-and triple-modified peptides were analyzed separately and, for a given PTM, the most significant result was reported. The data was analyzed with the same Bayesian linear model as proteome SARS-CoV/-2 infection data. In addition to the intensities normalization, for each replicate sample the scale of the effects in the experimental design matrix was adjusted, so that on average the correlation between log fold-changes of the replicates was 1:1. The same logic as for the proteome analysis, was applied to identify significant changes, but the median log 2 fold change had to be larger than 0.5, or 0.25 for the less stringent test. We additionally required that the PTM peptides are quantified in at least two replicates of at least one of the compared conditions. To ignore the changes in PTM site intensities that are due to proteome-level regulation, we excluded PTM sites on significantly regulated proteins if the direction of protein and PTM site changes was the same and the difference between their median log 2 fold changes was less than 2. Phosphoproteomics data were further analyzed with Ingenuity Pathway Analysis software (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis) Transcriptomic analysis of SARS-CoV-2 and SARS-CoV infected A549-ACE2 cells As for the analysis of the transcriptome data, Gencode gene annotations v28 and the human reference genome GRCh38 were derived from the Gencode homepage (EMBL-EBI). Viral genomes were derived from GenBank (SARS-CoV-2 -LR824570.1, and SARS-CoV -AY291315.1). Dropseq tool v1.12 was used for mapping raw sequencing data to the reference genome. The resulting UMI filtered count matrix was imported into R v3.4.4. CPM (counts per million) values were calculated for the raw data and genes having a mean cpm value less than 1 were removed from the dataset. A dummy variable combining the covariates infection status (mock, SARS-CoV, SARS-CoV-2) and time point was used for modeling the data within Limma (v3.46.0) 73 .

Data was transformed with the Voom method 73 followed by quantile normalization. Differential testing was performed between infection states at individual timepoints by calculating moderated t-statistics and p-values for each host gene. A gene was considered to be significantly regulated if the FDR adjusted p-value was below 0.05. The data for this study have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB38744.

We have used Gene Ontology, Reactome and other EnrichmentMap gene sets of human proteins (version 2020.10) 74 as well as protein complexes annotations from IntAct Complex Portal (version 2019.11) 75 and CORUM (version 2019) 76 . PhosphoSitePlus (version 2020.08) was used for known kinase-substrate and regulatory sites annotations, Perseus (version 1.6.14.0) 77 was used for annotation of known kinase motifs. For transcription factor enrichment analysis (Extended data Fig. 2e ) the significantly regulated transcripts were submitted to ChEA3 web-based application 78 and ENCODE data on transcription factortarget gene associations were used 79 . To find the nonredundant collection of annotations describing the unique and shared features of multiple experiments in a dataset (Figure 1d , Extended data Fig. 2l, m) , we have used inhouse Julia package OptEnrichedSetCover.jl 80 , which employs evolutionary multi-objective optimization technique to find a collection of annotation terms that have both significant enrichments in the individual experiments and minimal pairwise overlaps. The resulting set of terms was further filtered by requiring that the annotation term has to be significant with the specified unadjusted Fisher's Exact Test p-value cutoff at least in one of the experiments or comparisons (the specific cutoff value is indicated in the figure legend of the corresponding enrichment analysis). The generation of diagonally-split heatmaps was done with VegaLite.jl package (https://github.com/queryverse/VegaLite.jl).

For matching the PTMs of SARS-CoV-2 and SARS-CoV the protein sequences were aligned using the BioAlignments.jl Julia package (v.2.0, https://github.com/BioJulia/BioAlignments.jl).

Needleman-Wunsch algorithm with BLOSUM80 substitution matrix, -5 and -3 penalties for the gap and extension, respectively. As for the cellular proteins, we required that the viral phosphorylation or ubiquitination site is observed with q-value ≤ 10 -3 and localization probability ≥ 0.75. For the PTMs with lower confidence (q-value ≤ 10 -2 and localization probability ≥ 0.5) we required that the same site is observed with high confidence at the matching position of the orthologous protein of the other virus.

To systematically detect functional interactions, which may connect the cellular targets of each viral protein (interactome dataset) with the downstream changes it induces on proteome level (effectome dataset), we have used the network diffusion-based HierarchicalHotNet method 36 as implemented in Julia package HierarchicalHotNet.jl 81 . Specifically, for network diffusion with restart, we used the ReactomeFI network (version 2019) 35 of cellular functional interactions, reversing the direction of functional interaction (e.g. replacing kinase→substrate interaction with substrate→kinase). The proteins with significant abundance changes upon bait overexpression (|median(log 2 fold change)| The restart probability was set to 0.4, as suggested in the original publication, so that the probability of the random walk to stay in the direct neighborhood of the node is the same as the probability to visit more distant nodes. To find the optimal cutting threshold of the resulting hierarchical tree of strongly connected components (SCCs) of the weighted graph corresponding to the stationary distribution of signal diffusion and to confirm the relevance of predicted functional connections, the same procedure was applied to 1 000 random permutations of vertex weights as described in Reyna et al. 36 (vertex weights are randomly shuffled between the vertices with similar in-and out-degrees). Since cutting the tree of SCCs at any threshold t (keeping only the edges with weights above t) and collapsing each resulting SCC into a single node produces the directed acyclic graph of connections between SCCs, it allowed efficient enumeration of the paths from the "source" nodes (proteins strongly perturbed by viral protein expression with vertex weight w, w ≥ 1.5) to the "sink" nodes (interactors of the viral protein). At each threshold t, the average inverse of the path length from source to sink nodes was calculated as:

where N src is the number of "sources", N sink is the number of "sinks", L SCC (p) is the number of SCCs that the given path p from source to sink goes through, and the sum is for all paths from sources to sinks. The metric changes from 1 (all sources and sinks in the same SCC) to 0 (no or infinitely long paths between sources and sinks). For the generation of the diffusion networks we were using the t opt threshold that maximized the difference between In the generated SCC networks, the direction of the edges was reverted back, and the results were exported as GraphML files using in-house Julia scripts 64 . The catalogue of the networks for each viral bait is available as Supplementary Data 1.

To assess the significance of edges in the resulting network, we calculated the p-value of the edge g i g j as the probability that the permuted data-based transition probability between the given pair of genes is higher than the real data-based one:

This p-value was stored as the "prob_perm_walkweight_greater" edge attribute of GraphML output. The specific subnetworks predicted by the network diffusion (Figure 4b -d) were filtered for edges with p-value ≤ 0.05.

When the g i g j connection was not present in the ReactomeFI network, to recover the potential short pathways connecting g i and g j , ReactomeFI was searched for intermediate g k nodes, such that the edges g i g k and g k g j are present in ReactomeFI. The list of these short pathways is provided as "flowpaths" edge attribute in GraphML output.

The GraphML output of network diffusion was prepared for publication using yEd (v.3.20, www.yworks.com).

The intersection between the data generated by this study and other publicly available datasets was done using the information from respective supplementary tables. When multiple viruses were used in a study, only the comparisons with SARS-CoV and SARS-CoV-2 were included. For time-resolved data, all time points up to 24 h.p.i. were considered. The dataset coverage was defined as the number of reported distinct protein groups for proteomic studies and genes for transcriptomic studies. Confident interactions/significant regulations were filtered according to the criteria specified in the original study. A hit was considered as "confirmed" when it was significant both in this and external data and showed the same trend.

RNA isolation from SARS-CoV and SARS-CoV-2 infected A549-ACE2 cells was performed as described above (Qiagen). 500 ng total RNA was used for reverse transcription with PrimeScript RT with gDNA eraser (Takara). For relative transcript quantification PowerUp SYBR Green (Applied Biosystems) was used. Primer sequences can be provided upon request.

HEK293T cells were transfected with pWPI plasmid encoding single HA-tagged viral proteins, alone or together with pTO-SII-HA expressing host factor of interest. 48 hours after transfection, cells were washed in PBS, flash frozen in liquid nitrogen and kept at -80°C until further processing. Co-immunoprecipitation experiments were performed as described previously 55, 56 . Briefly, cells were lysed in lysis buffer (50 mM Tris-HCl pH 7.5, 100 mM NaCl, 1.5 mM MgCl 2 , 0.2% (v/v) NP-40, 5% (v/v) glycerol, cOmplete protease inhibitor cocktail (Roche), 0.5% (v/v) 750 U/µl Sm DNAse) and sonicated (5 min, 4°C, 30 sec on, 30 sec off, low settings; Bioruptor, Diagenode SA). HA or Streptactin beads were added to cleared lysates and samples were incubated for 3 hours at 4°C under constant rotation. Beads were washed six times in the lysis buffer and resuspended in 1x SDS sample buffer 62,5 mM Tris-HCl pH 6.8, 2% SDS, 10% glycerol, 50 mM DTT, 0.01% bromophenol blue). After boiling for 5 minutes at 95°C, a fraction of the input lysate and elution were loaded on NuPAGE™ Novex™ 4-12% Bis-Tris (Invitrogen), and further submitted to western blotting using Amersham Protran nitrocellulose membranes. Imaging was performed by HRP luminescence (ECL, Perkin Elmer).

SARS-CoV-2 infected A549-ACE2 cell lysates were sonicated (10 min, 4°C, 30 sec on, 30 sec off, low settings; Bioruptor, Diagenode SA). Protein concentration was adjusted based on Pierce660 assay supplemented with ionic detergent compatibility reagent. After boiling for 5 min at 95°C and brief max g centrifugation, the samples were loaded on NuPAGE™ Novex™ 4-12% Bis-Tris (Invitrogen), and blotted onto 0,22 µm Amersham ™ Protran ® nitrocellulose membranes (Merck). Primary and secondary antibody stainings were performed according to the manufacturer's recommendations. Imaging was performed by HRP luminescence using Femto kit (ThermoFischer Scientific) or Western Lightning PlusECL kit (Perkin Elmer).

N CTD dimers of SARS-CoV-2 (PDB: 6YUN) and SARS-CoV (PDB: 2CJR) were superimposed by aligning the α -carbons backbone over 111 residues (from position 253/254 to position 364/365 following SARS-CoV-2/SARS-CoV numbering) by using the tool MatchMaker 82 as implemented in the Chimera software 83 . Ubiquitination sites were visually inspected and mapped by using the PyMOL software (https://pymol.org). Phosphorylation on Ser310/311 was simulated in silico by using the PyTMs plugin as implemented in PyMOL 84 . Inter-chain residue contacts, dimer interface area, free energy and complex stability were comparatively analyzed between non-phosphorylated and phosphorylated SARS-CoV-2 and SARS-CoV N CTD by using the PDBePISA server 85 . Poisson-Boltzmann electrostatic surface potential of native and post-translationally modified N CTD was calculated by using the PBEQ Solver tool on the CHARMM-GUI server by preserving existing hydrogen bonds 86 . Molecular graphics depictions were produced with the PyMOL software.

The following reporter constructs were used in this study: pISRE-luc was purchased from Stratagene, EF1-α-ren from Engin Gürlevik, pCAGGS-Flag-RIG-I from Chris Basler, pIRF1-GAS-ff-luc, pWPI-SMN1-flag and pWPI-NS5 (ZIKV)-HA was described previously 56,87 .

For the reporter assay, HEK293-R1 cells were plated in 24-well plates 24 hours prior to transfection. Firefly reporter and Renilla transfection control were transfected together with plasmids expressing viral proteins using polyethylenimine (PEI, Polysciences) for untreated and treated conditions. In 18 hours cells were stimulated for 8 hours with a corresponding inducer and harvested in the passive lysis buffer (Promega). Luminescence of Firefly and Renilla luciferases was measured using dual-luciferase-reporter assay (Promega) according to the manufacturer's instructions in a microplate reader (Tecan).

Total amounts of IFN-α/β in cell supernatants were measured by using 293T cells stably expressing the firefly luciferase gene under the control of the mouse Mx1 promoter (Mx1luc reporter cells) 88 . Briefly, HEK293-R1 cells were seeded, transfected with pCAGGS-flag-RIG-I plus viral protein constructs and stimulated as described above. Cell supernatants were harvested in 8 hours. Mx1-luc reporter cells were seeded into 96-well plates in triplicates and were treated 24 hours later with supernatants. At 16 hours post-incubation, cells were lysed in the passive lysis buffer (Promega), and luminescence was measured with a microplate reader (Tecan). The assay sensitivity was determined by a standard curve.

A549-ACE2 cells were seeded into 96-well plates in DMEM medium (10% FCS, 100 ug/ml Streptomycin, 100 IU/ml Penicillin) one day before infection. Six hours before infection, or at the time of infection, the medium was replaced with 100μl of DMEM medium containing either the compounds of interest or DMSO as a control. Infection was performed by adding 10μl of SARS-CoV-2-GFP (MOI 3) per well and plates were placed in the IncuCyte S3 Live-Cell Analysis System where whole well real-time images of mock (Phase channel) and infected (GFP and Phase channel) cells were captured every 4 hours for 48 hours. Cell viability (mock) and virus growth (mock and infected) were assessed as the cell confluence per well (Phase area) and GFP area normalized on cell confluence per well (GFP area/Phase area) respectively using IncuCyte S3 Software (Essen Bioscience; version 2019B Rev2).

For comparative analysis of antiviral treatment activity against SARS-CoV and SARS-CoV-2, A549-ACE2 cells were seeded in 24-well plates, as previously described. Treatment was performed for 6 hours with 0.5ml of DMEM medium containing either the compounds of interest or DMSO as a control, and infected with SARS-CoV-Frankfurt-1 or SARS-CoV-2-MUC-IMB-1 (MOI 1) for 24 hours. Total cellular RNA was harvested and analyzed by qRT-PCR, as previously described.

The raw sequencing data for this study have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB38744 (https://www.ebi.ac.uk/ena/data/view/PRJEB38744). The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE 89 partner repository with the dataset identifier PXD022282, PXD020461 and PXD020222.

The protein interactions from this publication have been submitted to the IMEx (http://www.imexconsor-tium.org) consortium through IntAct 90 with the identifier IM-28109.

The data and analysis results are accessible online via the interactive web interface at https://covinet.innatelab.org. t  r  u  c  t  u  r  e  o  f  t  h  e  S  A  R  S  C  o  r  o  n  a  v  i  r  u  s  N  u  c  l  e  o  c  a  p  s  i  d  P  r  o  t  e  i  n  R  N  A  -b  i  n  d  i  n  g   D  i  m  e  r  i  z  a  t  i  o  n  D  o  m  a  i  n  S  u  g  g  e  s  t  s  a  M  e  c  h  a  n  i  s  m  f  o  r  H  e  l  i  c  a  l  P  a  c  k  a  g  i  n  g  o  f  V  i  r  a  l  R  N  A  .   J  o  u  r  n  a  l  o  f   M  o  l  e  c  u  l  a  r  B  i  o  l  o  g  y   3  6  8  ,  1  0  7  5  -1  0  8  6  (  2  0  0 a  n  t  l  o  ,  E  .  ,  B  u  k  r  e  y  e  v  a  ,  N  .  ,  M  a  r  u  y  a  m  a  ,  J  .  ,  P  a  e  s  s  l  e  r  ,  S  .  &  H  u  a  n  g  ,  C  .  A  n  t  i  v  i  r  a  l  a  c  t  i  v  i  t  i  e  s  o  f  t  y  p  e  I   i  n  t  e  r  f  e  r  o  n  s  t  o  S  A  R  S  -C  o  V  -2  i  n  f  e  c  t  i  o  n  .   A  n  t  i  v  i  r  a  l  R  e  s   1  7  9  ,  1  0  4  8  1  1  (  2  0  2 the relationship between the minimal allowed edge weight of the random walk graph (xaxis) and the mean inverse length of the path from the regulated proteins to the host targets of the viral protein along the edges of the resulting filtered subnetwork (y-axis). The red curve represents the metric for the network diffusion analysis of the actual data. The grey band shows 50% confidence interval, and dashed lines correspond to 95% confidence interval for the average inverse path length distribution for 1000 randomised datasets.

Optimal edge weight threshold that maximizes the difference between the metric based on real data and its 3rd quartile based on randomized data is highlighted by the red vertical line. 

Expression of HA-tagged viral proteins, in stably transduced A549 cells, used in AP-MS and proteome expression measurements. When several bands are present in a single lane, * marks the band with expected molecular weight (n = 4 independent experiments). For gel source data, see Supplementary Figure 1. (b) Extended version of the virus-host protein

SARS-CoV proteins, as well as ORF3 of HCoV-NL63 and ORF4 and ORF4a of HCoV

Host targets regulated upon viral protein overexpression are highlighted (see the in-plot legend). (c-f) Co-precipitation experiments in HEK293T cells showing a specific enrichment of (c) endogenous MAVS co-precipitated with C-terminal HA-tagged ORF7b of SARS-CoV-2 and SARS-CoV (negative controls: SARS-CoV-2 ORF6-HA, ORF7a-HA), (d) ORF7b-HA of SARS-CoV-2 and SARS-CoV

For gel source data, see Supplementary Figure 1. AP-MS: affinity-purification coupled to mass spectrometry; MD: Macro domain; NSP: Non-structural protein

Extended data Figure 2 | SARS-CoV-2 and SARS-CoV proteins trigger shared and specific interactions with host factors, and induce changes to the host proteome

Differential enrichment of proteins in (a) NSP2 and (b) ORF8 of SARS-CoV-2 (x-axis) vs SARS-CoV (y-axis) AP-MS experiments (n=4 independent experiments)

Ontology Biological Processes enriched among the cellular proteins that are up-(red arrow) or down-(blue arrow) regulated upon overexpression of individual viral proteins. (d) The A549 cells (see materials and methods for the exact protein selection criteria). Homologous viral proteins are displayed as a single node

Extended data Figure 3 | RCOR3 and APOB regulation upon SARS-CoV-2 and SARS-CoV protein over-expression. (a-b) Normalized intensities of selected candidates specifically perturbed by individual viral proteins: (a) RCOR3 was upregulated both by SARS-CoV-2 and SARS-CoV NSP4 proteins, (b) APOB was upregulated by ORF3 and downregulated by NSP1 specifically to SARS-CoV-2. The box and the whiskers represent 50% and 95% confidence intervals

Extended data Figure 4 | Tracking of virus-specific changes in infected A549-ACE2 cells by transcriptomics and proteomics. (a) Western blot showing ACE2-HA expression levels in A549 cells untransduced (wild-type) or transduced with ACE2-HA-encoding lentivirus (n = 2 independent experiments). For gel source data, see

mRNA expression levels of SARS-CoV-2 N relative to RPLP0 as measured by qRT-PCR upon infection of wild-type A549 and A549-ACE2 cells at the indicated MOIs. Error bars represent mean and standard deviation (n=3 independent experiments). (c) Volcano plot of mRNA expression changes of A549-ACE2 cells, infected with SARS-CoV-2 at an MOI of 2 in comparison to mock infection at 12 h.p.i. Significant hits are highlighted in gray

Diamonds indicate that the actual log 2 fold change or p-value were truncated to fit into the plot. (d) Expression levels, as measured by qRT-PCR, of SARS-CoV-2/SARS-CoV N and host transcripts relative to RPLP0 in infected (MOI 2) A549-ACE2 cells with SARS-CoV-2 (orange) and SARS-CoV (brown) at indicated time points. Error bars correspond to mean cells (n=3 independent experiments). For gel source data, see

as well as (h) an increase in total protein abundance of APOB with indicated median, 50% and 95% confidence intervals (n=4 independent experiments). (i) Overview of perturbations to hostcell innate immunity-related pathways, induced by distinct proteins of SARS-CoV-2, derived from the network diffusion model and overlaid with transcriptional, ubiquitination and phosphorylation changes upon SARS-CoV-2 infection. (j) Heatmap showing the effects of the indicated SARS-CoV-2 proteins on type-I IFN expression levels, ISRE and GAS promoter activation in HEK293-R1. Accumulation of type-I IFN in the supernatant was evaluated by testing supernatants of PPP-RNA (IVT4) stimulated cells on MX1-luciferase reporter cells, ISRE promoter activation -by luciferase assay after IFN-α stimulation, and GAS promoter activation -by luciferase assay after IFN-γ stimulation in cells expressing SARS

Overview of perturbations to host-cell Integrin-TGF-β-EGFR-RTK signaling, induced by distinct proteins of SARS-CoV-2, derived from the network diffusion model and overlaid with transcriptional, ubiquitination and phosphorylation changes upon SARS-CoV-2 infection. (b) Profile plots of total protein levels of SERPINE1 and FN1 in SARS-CoV-2 or SARS-CoV-infected A549-ACE2 cells at 6, 12, and 24 h.p.i., with indicated median

We thank Stefan Pöhlmann for sharing ACE2 plasmids and Robert Baier for technical assistance, Monika Trojan for advice on drugs, Juan Pancorbo and Johannes Albert-von der Gönna from the Leibniz Supercomputing Centre (www.lrz.de) for technical assistance. We further thank the Karl Max von Bauernfeind -Verein for support for the screening microscope. Work in the authors' laboratories was supported by an ERC consolidator grant (ERC-CoG ProDAP, 817798), the German Research Foundation (PI 1084/3, PI 1084/4, PI