Ubiquitin-Modified Proteome of SARS-CoV-2-Infected Host Cells Reveals Insights into Virus-Host Interaction and Pathogenesis.

The outbreak of coronavirus disease 2019 (COVID-19), which is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has posed a serious threat to global public health. The mechanism of pathogenesis and the host immune response to SARS-CoV-2 infection are largely unknown. In the present study, we applied a quantitative proteomic technology to identify and quantify the ubiquitination changes that occur in both the virus and the Vero E6 cells during SARS-CoV-2 infection. By applying label-free, quantitative liquid chromatography with tandem mass spectrometry proteomics, 8943 lysine ubiquitination sites on 3086 proteins were identified, of which 138 sites on 104 proteins were quantified as significantly upregulated, while 828 sites on 447 proteins were downregulated at 72 h post-infection. Bioinformatics analysis suggested that SARS-CoV-2 infection might modulate host immune responses through the ubiquitination of important proteins, including USP5, IQGAP1, TRIM28, and Hsp90. Ubiquitination modification was also observed on 11 SAR-CoV-2 proteins, including proteins involved in virus replication and inhibition of the host innate immune response. Our study provides new insights into the interaction between SARS-CoV-2 and the host as well as potential targets for the prevention and treatment of COVID-19.

Introduction

The outbreak of coronavirus disease 2019 (COVID-19) is a serious threat to global public health, leading to numerous deaths and massive economic impacts. COVID-19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is the third coronavirus to cause a human epidemic in the last 2 decades.[1] Coronaviruses (family Coronaviridae) are positive-sense, single-stranded RNA viruses that can infect a wide range of avian and mammalian hosts to cause diseases of varying severity, ranging from mild respiratory illnesses to severe fatal illnesses. The proteins encoded by coronaviruses include 4 main structural proteins, spike (S), envelope (E), membrane (M), and nucleoprotein (N), 16 major proteins. SARS-CoV-2 shares 79.5% genetic sequence identity with SARS-CoV and 96.2% homology with a bat coronavirus.[2,3] Excessive cytokine release after SARS-CoV-2 infection results in an uncontrolled “cytokine storm” and multi-organ damage that lead to death in 28% of fatal COVID-19 cases.[4] However, the underlying mechanism and the altered targets remain unclear.

Post-translational modifications (PTMs) are chemical modifications that are essential for precise regulation of the properties and functions of proteins and are key mechanisms that increase proteomic diversity in cells. Lysine ubiquitination (Kub) is a PTM in eukaryotic cells that is involved in diverse biological activities. In addition to mediating degradation via the ubiquitin–proteasome system (UPS), Kub also regulates the cell cycle, subcellular localization, enzyme activation, cell signal transduction, and other physiological cellular functions. Research on virus infection showed that ubiquitination modification plays a key role in the innate immune response to viral infections by inducing antiviral defenses and facilitating virus replication.[5−12] It has been shown that an active UPS is required for the replication of some viruses, including human immunodeficiency virus, rotavirus, and coxsackievirus.[5−7] Ubiquitination of the tombusvirus p33 replication protein was shown to be critical for virus replication during tombusvirus infection via its interaction with the host protein Vps23p.[9] Ubiquitination and deubiquitination of NP proteins have also been shown to regulate influenza A virus RNA replication.[8]

Despite limited evidence, ubiquitination has been implicated in the mechanisms underlying the process of coronavirus infection.[13−16] In one study, porcine coronavirus was found to evade host immune surveillance through reduction of host IFN-β production by blocking K63-linked polyubiquitination.[13] Siu et al.[14] reported activation of the NLRP3 inflammasome by SARS-CoV open reading frame 3a (ORF3a) protein after infection, which was mediated by TRAF3-dependent ubiquitination, resulting in a proinflammatory cytokine storm. In addition, SARS-CoV negatively regulates host antiviral defenses by disrupting STING-mediated signaling and IFN induction via downregulation of ubiquitinated STING, retinoic acid-inducible gene I (RIG-I), TBK1, and IRF-3.[16] However, little is known about the alterations in host ubiquitination caused by SARS-CoV-2 infection.

In the present study, we applied a quantitative proteomic technology to investigate the ubiquitination changes that occur in Vero E6 cells during SARS-CoV-2 infection. Using bioinformatics methods, we profiled the lysine ubiquitome in SARS-CoV-2-infected Vero E6 cells and identified 966 significant differential ubiquitination sites mapped to 551 proteins. Our findings improve our understanding of the mechanism of COVID-19 pathogenesis as it relates to protein ubiquitination as well as the identification of potential therapeutic targets.

Experimental Procedures

We performed biological replicates by label-free quantification (LFQ), ubiquitination enrichment technology, and high-resolution liquid chromatography–mass spectrometry combined with quantitative proteomics research strategy to analyze the intracellular protein and viral protein after SAR-CoV-2 infects the hosts. In this study, the quantitative values of each sample in three replicates were obtained by LFQ intensity. The first step is to calculate the differential expression of the protein between the two groups of samples and then calculate the ratio of the average values between the two groups of samples. The ratio is used as the final quantitation. The second step is to calculate the significant p value of differential expression between two groups of samples. The relative quantitative values of each sample were taken as log2 transform (so that the data conforms to the normal distribution), and p value was calculated by the two-sample two-tailed T-test method. p value <0.05 and protein ratio >2 were regarded as upregulation. p value <0.05 and protein ratio <0.5 were regarded as downregulation. In each experiment, three samples of SARS-CoV-2 were analyzed qualitatively and quantitatively to verify the biological reliability of the measurement.

Cell Culture and Virus Isolation

The SARS-CoV-2 virus isolate used in this study, GD01, was isolated by the Guangdong Provincial Center for Disease Control and Prevention. Vero E6 cells were grown in minimum essential medium (MEM) supplemented with 10% fetal bovine serum in a humidified atmosphere containing 5% CO2 at 37 °C. The viruses were isolated from the alveolar lavage fluid of a COVID-19 patient. A single layer of Vero E6 cells was inoculated with alveolar lavage fluid (10% of the total volume of the cell culture fluid) and incubated in a 37 °C incubator with 5% CO2 for 2 h. After the incubation, 1 mL of the MEM containing 2% fetal bovine serum was added. A negative control was set. Cells were observed for cytopathic effects (CPEs), and the viruses were collected when the CPE was visible in 75% of the cell layer. After repeated freezing and thawing three times, the supernatant was collected by centrifugation and stored at −80 °C.

Virus Inoculation

Monolayers of Vero E6 cells were grown to 80% confluence in 10 cm2 cell culture dishes. After removing the medium, the cells were gently washed with serum-free Dulbecco’s modified Eagle’s medium prior to inoculation. Then, the monolayers were inoculated with SARS-CoV-2 at a multiplicity of infection (MOI) of 5.0. Uninfected cells were used as the mock-infected group. The cells were collected from three independent experiments at 72 h post-infection (hpi). Viral propagation in Vero E6 cells was confirmed by the appearance of a virus-induced CPE and the one-step growth curve of SARS-CoV-2.

Protein Extraction and Trypsin Digestion

Uninfected cells and cells infected with SARS-CoV-2 were collected with cell scrapers. The cells were then lysed with lysis buffer (8 M urea and 1% protease inhibitor cocktail) and sonicated three times on ice using a high-intensity ultrasonic processor (Scientz). After removing the debris by centrifugation at 12,000g at 4 °C for 10 min, the protein concentration in the supernatant was determined using a BCA kit. Protein samples (2 mg) were reduced and digested as previously described.[17] Briefly, the samples were sequentially treated with 5 mM dithiothreitol and 11 mM iodoacetamide in the dark. For trypsin digestion, the protein sample was treated with trypsin overnight (at a trypsin/protein ratio of 1:50, wt/wt) and then for another 4 h (at 1:100, wt/wt).

Affinity Enrichment of Ubiquitinated Peptides

Anti-lysine ubiquitin (Kub) antibody beads (lot number 001, PTM Bio, China) were washed with cold PBS. Then, to enrich the Kub peptides, 2 mg of trypsin peptides was dissolved in NETN buffer (100 mM NaCl, 1 mM EDTA, 50 mM Tris-HCl, 0.5% NP-40, pH 8.0) and then gently shaken with pre-washed antibody beads at 4 °C overnight. Then, the antibody beads were washed four times with NETN buffer and twice with ddH2O in sequence. The peptides were eluted from the beads, and the eluted fractions were vacuum-dried and cleaned with a C18 ZipTip (Millipore, Billerica, MA).

Liquid Chromatography with Tandem Mass Spectrometry Analysis

Three parallel analyses for each fraction were conducted. The peptides (0.5 μg for global proteome profiling and all the peptide samples for ubiquitinome profiling) were dissolved in 0.1% formic acid (vol/vol) and loaded onto a reversed-phase analytical column (15 cm length and 75 μm i.d.) with an increased gradient of 6–23% solvent (0.1% formic acid in 98% acetonitrile) for 26 min, 23–35% for 8 min, increasing to 80% within 3 min, and then maintained at 80% solvent for 3 min on an EASY-nLC 1000 UPLC system at a constant flow rate of 400 nL/min.

The peptides were subjected to nanospray ionization and subsequently analyzed by tandem mass spectrometry (MS/MS) in a Q Exactive Plus (Thermo) coupled online to the ultra-performance liquid chromatography (UPLC). The electrospray voltage applied was 2.0 kV, and the scan range was set at 350–1800 m/z. Intact peptides were detected in the Orbitrap at a resolution of 70,000. The peptides were then selected for MS/MS using the NCE setting as 28 and the fragments were detected in the Orbitrap at a resolution of 17,500. Data were collected using a data-dependent procedure that alternated between a single MS scan followed by 20 MS/MS scans with 15.0 s dynamic exclusion.

Database Search

The resulting data were analyzed using the MaxQuant search engine (v.1.5.2.8) (PTMs including N-terminal acetylation, methionine oxidation, asparagine/glutamine deamidation, lysine ubiquitination, and carbamidomethyl). Tandem mass spectra were searched against the human UniProt database (Chlorocebus_sabaeus, 20191104 downloaded from Uniprot) concatenated with a reverse decoy database. For trypsin, a maximum of four missed cleavages were allowed. The mass deviation was set to 20 and 5 ppm for precursor ions in the first and main searches, respectively. The false discovery rate was adjusted to less than 1%, and the minimum score for modified peptides was set to more than 40.

Bioinformatics Methods

Motif Analysis

Soft MoMo (a motif-x algorithm) was used to analyze the model sequences in which amino acids constituted at specific positions with modify-21-mers (spanning 10 amino acids upstream and downstream) but phosphorylation with modify-13-mers. All the database protein sequences were chosen as the background database parameters. The minimum number of occurrences was set to 20. Emulate original motif-x was ticked, and all other parameters were the defaults.

Functional Enrichment

Enrichment of Gene Ontology Analysis

The retrieved protein information was classified by gene ontology (GO) annotation according to three categories: biological process (BP), cellular compartment, and molecular function (MF). A two-tailed Fisher’s exact test was used to identify differentially modified proteins, and a p value less than 0.05 was considered significant.

Enrichment of Pathway Analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to identify enriched pathways by a two-tailed Fisher’s exact test, and p values less than 0.05 were considered significant. Pathways were hierarchically organized into categories according to the KEGG website.

Enrichment of Protein Domain Analysis

For each protein category, the InterPro database was used to classify the differentially modified proteins, which were the proteins with a corrected p value less than 0.05, using a two-tailed Fisher’s exact test.

Enrichment-Based Clustering

Following the above functional classifications of differentially modified proteins, we performed a hierarchical clustering analysis. The filtered p value matrix was transformed by using the function x = −log10 (p value). Then, the x values for each functional category were z-transformed, and the z-scores were clustered by one-way hierarchical clustering (Euclidean distance, average linkage clustering) in Genesis. A heat map for cluster membership was generated with the heatmap.2 function of the “gplots” package in the R software package.

Protein–Protein Interaction Networks

Protein interactions between differentially expressed modified proteins were analyzed using STRING database version 10.1. Only interactions between the proteins belonging to the searched data set were selected, thereby excluding the external candidates. The confidence score was used to define the interaction confidence, and interactions with a confidence score ≥0.7, which were considered to have high confidence, were fetched. The interaction network form STRING was visualized in the R package “networkD3.”

Luciferase Reporter Assays

The dual luciferase reporter assay was performed using the Promega luciferase reporter system. HEK-293T cells were seeded at 1 × 105 per well in a 24-well plate and were co-transfected with a reporter plasmid encoding either NF-κB, IFN-β, IRF3, or ISRE-Luc, along with the phRL-TK vector and expression plasmids, to evaluate the ISRE and NF-κB pathway. Twenty-four hours later, the cells were treated with Sendi virus (SEV), and luciferase activity was detected after 16 h. All reporter assays were repeated at least three times. Data are presented as mean ± standard deviations.

The method of “Protein Extraction and Trypsin Digestion”, “Liquid Chromatography with Tandem Mass Spectrometry Analysis”, “Database Search”, and “Bioinformatics Methods” are used for both global ubiquitinome profiling and global proteome profiling experiments.

Results

Host Cell Infection and Proteome Analysis of SARS-CoV-2

An experimental protocol was designed to analyze protein ubiquitylation by liquid chromatography with MS/MS (LC–MS/MS). SARS-CoV-2-infected Vero E6 cells were used as the experimental group and mock-infected Vero E6 cells were used as the control group. After infection, total cellular proteins of cell samples were extracted and decomposed into small peptides by protease treatment, and then the peptides were enriched with ubiquitinated resin, separated from the resin, and analyzed by LC–MS/MS (Figure A). To identify the optimal time point after SARS-CoV-2 infection for the ubiquitome analysis, Vero E6 cells (most of the cells used in SARS-CoV-2 research are Vero E6 and Caco2)[18,19] were inoculated with SARS-CoV-2 and microscopically monitored for CPE at 24, 48, 72, and 96 hpi. As shown in Figure B, a minimal CPE was observed at 48 hpi and then readily became apparent as the infection progressed at 72 hpi. A large number of cells died after 72 h, which were not suitable for proteomics. Therefore, we selected 72 hpi as the time point for ubiquitome analysis (Figure B).

Figure 1

Experimental strategy and virus infection. (A) A workflow diagram for the identification of Kub sites in Vero E6 cells. (B) Vero E6 cells infected with SARS-CoV-2 at an MOI of 5.0 or mock-infected. 40× magnification.

Analysis of our proteomics data from virus- and mock-infected cells identified 5876 differential proteins, and 4948 of these proteins were quantified. Based on a fold-change >2 and a t-test p-value less than 0.05 as the threshold for significance, 46 proteins were upregulated and 141 were downregulated (Tables S4 and S5, Supporting Information Figures S1 and S2).

4D LFQ of Lysine Ubiquitination in SARS-CoV-2-Infected Vero E6 Cells and Quality Validation of the Ubiquitome Data

As shown in Figure A, the length of the trypsin peptides was predominantly distributed between 7 and 20 amino acids, which was in accordance with the guideline for trypsin digestion high-energy collisional dissociation (HCD) fragmentation. Peptides of less than 5 amino acids cannot yield an effective sequence identification because of the small number of fragment ions, and peptides longer than 20 amino acids are not suitable for the fragmentation of HCD due to their high mass and charge. Thus, the distribution of peptide lengths as examined by MS met the quality control requirements. From our ubiquitome data, a total of 424,771 secondary spectra were obtained by MS analysis. After searching the protein data library, the secondary spectrum of the mass spectrometer obtained an available effective spectrum number of 65,777 for a spectrum utilization rate of 15.5%. A total of 16,271 peptides and 8822 ubiquitinated modified peptides were identified by the spectrum analysis segment.

Figure 2

Overview of ubiquitome data in SARS-CoV-2-infected cells. (A) Distribution of tryptic peptide lengths. (B) Number of Kub sites per protein. (C,D) Numbers of proteins and sites that were identified and quantified by MS/MS. (E) Numbers of sites (left) and proteins (right) that met the threshold criterion (i.e., a fold change >2 or <0.5). (F) Crosstalk between proteomics data (LQ) and protein ubiquitination-modified omics data (LP), with Log2 protein virus-infected/mock-infected ratio as the vertical (y) axis and log2 Kub virus-infected/mock-infected ratio as the horizontal (x) axis.

A total of 8943 ubiquitination sites were mapped to 3086 proteins, and 3129 ubiquitination sites in 1371 proteins were quantified (Figure C, Table S1). The number of modification sites in each protein ranged from 1 to 30 (Figure B). Using a quantification ratio >2 as the threshold for upregulation and <0.5 as the threshold for downregulation, the results showed that 138 Kub sites on 104 proteins were upregulated and 828 Kub sites on 447 proteins were downregulated (Figure D).

We performed an analysis of 2720 lysine sites that had acceptable ratios and p values in both the proteomics data (LP) and protein ubiquitination-modified omics data (LQ), as shown in Figure E, with the log2 protein virus-infected/mock-infected group ratio as the vertical (y) axis and the log2 Kub virus-infected/mock-infected group ratio as the horizontal (x) axis. Based on the data in the figure, there were 224 K sites with upregulated protein expression and ubiquitination in the first quadrant, 1029 K sites with upregulated protein expression but downregulated ubiquitination in the second quadrant, 1083 K sites with downregulated protein expression and ubiquitination in the third quadrant, and 384 K sites with downregulated protein expression and upregulated ubiquitination in the fourth quadrant. In addition, there were two proteins with decreased protein expression and unchanged ubiquitination modification, which were located on the negative half of the y-axis (Figure E). We then screened out proteins with the same trend for all modification sites on a single protein, including 78 proteins in the first quadrant, such as CAPN5, ATP6V0D1, and RPS3; 319 proteins in the second quadrant, such as LAMTOR1, PMSD7, and TSG101; 362 proteins in the third quadrant, such as CCDC47, DDX5, and NAP1L1; and 124 proteins in the fourth quadrant, such as KARS, PRKCA, and POLR2B.

Characteristics of the Lysine Ubiquitination Motifs and Functional Categories of Host Proteins Ubiquitinated in Response to SARS-CoV-2 Infection

The characteristic sequence form is considered as a motif of the modified peptide if the number of peptides in this sequence is greater than 20 and the p value is less than 0.000001. In total, 17 motifs, named KubEL, KubXL, AKub, LXKub, IXKub, LKub, AXKub, VXKub, FXKub, LXXKub, KubL, VKub, KubXXXL, KubXV, VXXKub, IXXXKub, and KubI, were identified with different abundances (Kub = ubiquitinated lysine and X = any amino acid, Figure A and Table S3). The results showed both enrichment and depletion of some motifs, as shown in the heat map of the flanking amino acids from the −10 to +10 positions (Figure B). Valine (V) and leucine (L) were the most frequent amino acids occurring upstream of the Kub site, and leucine (L) was the most common amino acid occurring downstream of the Kub site. In infected samples, there were 12 ubiquitination site motifs in ligases and deubiquitinases (DUBs). There were 6 ubiquitination site motifs for ubiquitin ligases and 10 ubiquitination site motifs for DUBs (Table S7).

Figure 3

Sequence analysis of the ubiquitinated proteins in Vero E6 cells infected with SARS-CoV-2. (A) Enriched ubiquitination motif logos were produced by the Motif-X program and MoMo program. The size of the single letters refers to the frequency of the amino acid residue at that position. (B) Heat map of the amino acids surrounding Kub sites.

To gain a deeper understanding of the subcellular localization of proteins with differential ubiquitination during SARS-CoV-2 infection, COG/KOG and GO function classification analyses were performed. Using GO analysis, proteins with significant ubiquitination alterations were classified according to three categories: “BP”, “cellular component (CC)”, and “MF” (Supporting Information Figure 3). For the BP analysis, upregulated and downregulated Kub proteins were enriched in 28 processes, including protein localization or targeting to the endoplasmic reticulum, nuclear-transcribed mRNA catabolic process, mRNA catabolic process, regulation of cellular protein catabolic process, regulation of proteolysis or DNA metabolic process, and others. In the MF analysis, proteins with binding activity, catalytic activity, hydrolase activity, and others were enriched among proteins containing upregulated or downregulated Kub sites.

Based on the subcellular location analysis, among the 104 proteins containing upregulated Kub sites, 44.23% were localized in the cytoplasm, 22.12% in the plasma membrane, 9.62% in the nucleus, 7.69% in the extracellular space, 5.77% in the cytoplasm and nucleus, 2.88% in the mitochondrion, 2.88% in the cytoskeleton, 2.88% in the endoplasmic reticulum, and 1.92% in others (Figure A). Of the 447 proteins containing downregulated Kub sites, 46.09% were localized in the cytoplasm, 24.38% in the nucleus, 11.19% in the plasma membrane, 5.82% in the cytoplasm and nucleus, 4.92% in the mitochondrion, 4.47% in the extracellular space, and 3.13% in others (Figure B). In the COG/KOG analysis, we found that many upregulated proteins were associated with “translation, ribosomal structure and biogenesis” and “intracellular trafficking, secretion, and vesicular transport” (Figure C), while the downregulated proteins were associated with “post-translational modification, protein turnover, and chaperone”; “intracellular trafficking, secretion, and vesicular transport”; and “signal transduction mechanisms” (Figure D). These data reveal that SARS-CoV-2 may have a big impact on intracellular protein transport (Table S2).

Figure 4

Classification of the identified proteins containing upregulated and downregulated Kub sites. (A,B) Subcellular localization of proteins and (C,D) COG annotation of proteins.

Ubiquitinated Proteins Are Enriched in Functional Categories and Clusters Related to Host Responses

The identified proteins were subjected to functional enrichment (GO term, domain, and KEGG pathway) analyses. The GO terms were divided into three sub-ontologies (Figure A,B). In the “BP” categories, the most enriched among the upregulated and downregulated proteins were “protein localization to endoplasmic reticulum” and “regulation of cellular protein catabolic process”. The “CC” category analysis showed that upregulated proteins were most enriched in “endoplasmic reticulum sub-compartment”, while the downregulated proteins were most enriched in “myelin sheath”. For the “MF” categories, the terms “ubiquitin binding” and “nucleoside-triphosphatase activity” ranked first for the upregulated and downregulated proteins, respectively.

Figure 5

Functional enrichment analysis of identified proteins containing upregulated and downregulated Kub sites. BPs, CCs, and MFs.

KEGG pathway enrichment analysis showed that the upregulated proteins were associated with an “NOD-like receptor signaling pathway”, “ribosome”, “cholesterol metabolism”, and “drug metabolism-cytochrome P450” (Supporting Information Figure 4A), and the downregulated proteins were enriched in the pathways of “endocytosis”, “pancreatic secretion”, and “glycolysis/gluconeogenesis” (Supporting Information Figure 4B).

Protein domain enrichment analysis of the proteins revealed that upregulated proteins were enriched in “tubulin/FtsZ family and GTPase domain”, “fatty acid desaturase”, “cytochrome P450”, and “cytochrome b5-like heme/steroid binding domain” (Supporting Information Figure 4C), while the downregulated proteins were significantly enriched in “KH domain”, “Ras family”, and “ThiF family” (Supporting Information Figure 4D).

A heat map of the clustering analysis was generated to present the correlations between differentially ubiquitinated sites (Figure ). For the KEGG enrichment analysis (Supporting Information Figure 5A), upregulated proteins in category Q1, with a fold change <0.5, were enriched in “cardiac muscle contraction” and “endocrine and other factor-regulated calcium reabsorption”, whereas downregulated proteins were enriched in “ribosome”. In category Q2, which contained proteins with fold changes from 0.5 to 0.667, upregulated proteins were enriched in “pentose phosphate pathway”, and downregulated proteins were enriched in “phagosome” and “tight junction”. Upregulated proteins in Q3, with a fold change between 1.5 and 2, were enriched in a variety of important pathways, including “phototransduction”, “Chagas disease (American trypanosomiasis)”, and “toxoplasmosis”, while downregulated proteins were enriched in “protein processing in endoplasmic reticulum”. Finally, in category Q4, with a fold change >2, upregulated proteins were enriched in “Huntington disease”, “fluid shear stress and atherosclerosis”, and “NOD-like receptor signaling pathway”, while downregulated proteins were enriched in “antigen processing and presentation”. Notably, the Nod-like receptor family, containing pyridine domain 3 (NLRP3), regulates the secretion of the proinflammatory cytokines interleukin 1 beta (IL-1β) and IL-18.[20] The related virus SARS-CoV is known to initiate a cytokine/chemokine storm.[21] Therefore, we speculate that SARS-CoV-2 infection may induce an inflammatory cytokine storm by upregulating the NOD-like signaling pathway in Vero E6 cells.

Figure 6

Clustering analysis of identified proteins containing differentially expressed Kub sites. Modified sites were classified into four classes based on fold change. Heat maps showing the results of a cluster analysis of data from the BPs, CCs, and MFs. Red indicates upregulated and blue indicates downregulated.

SARS-CoV-2 Infection Causes Ubiquitination Changes in Important Protein–Protein Interaction Networks

The protein–protein interaction networks are composed of individual proteins that interact with each other. Systematic analysis of the protein interactions in biological systems is important for understanding the working principles, response mechanisms, energy metabolism, and functional links between proteins. In Figure , we screened the top 50 most closely interacting proteins and plotted the protein interaction network. Our KEGG analysis results showed that the differential proteins were involved in the process of intracellular transport, while GO analysis showed that these proteins were located in the endoplasmic reticulum. In the protein–protein interaction analysis, we also found an interaction with NAC-A/B domain-containing protein (NACAD). NACAD can prevent inappropriate targeting of non-secreted peptides to the endoplasmic reticulum. When they emerge from ribosomes, they may bind to nascent polypeptide chains and prevent them from interacting with signal recognition particles that normally target nascent secreted peptides to the endoplasmic reticulum. It may also reduce the intrinsic affinity of ribosomes to protein transport sites (M sites) in the endoplasmic reticulum membrane. Our later experiments will also study the effect of SARS-CoV-2 on this aspect of intracellular transport.

Figure 7

Protein–protein interaction networks. Different colors represent the differential modification of proteins (blue means downregulated protein and red means upregulated protein). The circle size represents the number of differentially modified proteins and their interacting proteins. The larger the circle, the more proteins it interacts with, indicating that the protein is more important in the network. To clearly show the protein interactions, the top 50 proteins were screened out, and the interaction network with the closest interactions was mapped.

Ubiquitination of SARS-CoV-2 Proteins and Functional Domains

Modification of ubiquitination can be important for antiviral defense and virus replication, and viral proteins can also be modified.[17,22,23] The extent of ubiquitination of SARS-CoV-2 proteins is unclear. Therefore, we analyzed SARS-CoV-2 proteins for the presence of ubiquitinated residues. A total of 259 peptides were identified by spectrum analysis, and 135 ubiquitinated residues on 11 viral proteins were detected (Table ). The 11 viral proteins had at least one ubiquitinated residue (envelope (E), nonstructural protein 8, and ORF14) and up to 66 ubiquitinated residues (rep). Of these 11 proteins, four have defined functional domains that play important roles in virus replication and host cell interaction. Interestingly, ubiquitinated sites were identified in important functional domains. The receptor binding domain (RBD) is involved in viral interaction with ACE2, and neutralizing antibodies usually target these sites. Seven lysines in the RBD (319–541), K386, K417, K424, K458, K462, K529, and K535 were ubiquitinated. Nucleoprotein is another protein with ubiquitination (K61, K65, K143, K169, K233, K266, K342, K347, K355, K361, K375, K388, and K405); four of these sites were in the N-terminal domain (NTD, 46–174) and five were in the C-terminal dimerization domain (CTD, 247–365) (Figure and Tables and S6).

Table 1

Ubiquitinated SARS-CoV-2 Proteins and Sites

identified SARS-CoV-2	identified sites
spike glycoprotein (S)	41, 97, 150, 187, 195, 202, 206, 304, 310, 386, 417, 424, 458, 462, 529, 535, 558, 776, 786, 790, 795, 811, 814, 835, 933, 986, 1028, 1038, 1086, 1157, 1181,1191, 1255, 1266, 1269
protein 3a (3a)	16, 21, 132
envelope small membrane protein (E)	63
membrane protein (N/A)	162, 166, 180, 205
non-structural protein 6 (6)	42, 48
protein 7a (7a)	32, 85
non-structural protein 8 (N/A)	53
nucleoprotein (N)	61, 65, 143, 169, 233, 266, 342, 347, 355, 361, 375, 388, 405
protein 9b (N/A)	4, 40, 59, 67, 80, 97
uncharacterized protein 14 (ORF14)	10
replicase polyprotein 1ab (rep)	72, 210, 225, 247, 258, 292, 458, 497, 500, 510, 518, 527, 564, 625, 634, 636, 669, 672, 684, 699, 701, 714, 734, 851, 1124, 1233, 1305, 1315, 1387, 1396, 1529, 1745, 1753, 1795, 1855, 1860, 1878, 2148, 2191, 2200, 2478, 2533, 2610, 2757, 2849, 3213, 3215, 3324, 3573, 3839, 3843, 3850, 3861, 3988, 4221, 4433, 4465, 4552, 4818, 4822, 4943, 4969, 5455, 5738, 5908, 6874

Figure 8

Ubiquitination sites in the functional domains of SARS-CoV-2 proteins. The red lines represent the ubiquitinated lysine sites.

SARS-CoV-2 Infection Can Regulate IFN-I Signaling Pathway by Modulating the Ubiquitination of Important Proteins

Type I IFN plays an important role in antiviral infection process, and USP5 can inhibit the production of type I IFN.[24] USP5 can increase Lys11-linked ubiquitination of RIG-I and interacts with the E3 ligases of RIG-I (such as STUB1) to promote the ubiquitination of Lys48 in RIG-I. In this scenario, USP5 recruits STUB1 to RIG-I, thereby inhibiting type I IFN signaling and enhancing viral replication, resulting in increased infection. USP5 also affects downstream IRF3-induced signaling pathways[25,26] In our study, ubiquitination levels of K184, K423, K720, and K770 in USP5 were downregulated following SARS-CoV-2 infection, which was accompanied by K558 upregulation and upregulated expression of USP5 (1.2 fold) (Table S4). Through verification experiment, we found that the ubiquitination of USP5 was significantly downregulated (2054.87333 vs 2962.94, p = 1.06 × 10–3) and the expression of USP5 was increased (121.491 vs 74.31267, p = 5.59 × 10–3) in SARS-CoV-2-infected cells (Figure A,B), which were consistent with our data. Then, we used the dual luciferase system to explore the effects of USP5 on IFN-I. SEV significantly enhanced IFN-β, IRF3, NF-κB, and ISRE promoter activity, while overexpression of USP5 significantly inhibited the promoter activity in a dose-dependent manner (Figure C–F), suggesting USP5 as a potent inhibitor of virus-induced type I IFN signaling pathway. Therefore, SARS-CoV-2 might modulate K48-linked polyubiquitination to increase the expression of USP5 and inhibit type I IFN signaling and as a result may enhance viral replication (Figure A).

Figure 9

SARS-CoV-2 infection regulates type I IFN immune response by modulating USP5 ubiquitination. (A) Ubiquitination of USP5 was significantly downregulated and the expression of USP5 was upregulated following SARS-CoV-2 infection. (B) Integrated band density analysis on ubiquitination and expression of USP5, and the p value was calculated by the two-sample two-tailed T-test method. ** indicates p ≤ 0.01. (C–F) USP5 suppresses the virus-induced type I IFN signaling pathway.

Figure 10

Schematic diagram of the possible regulatory mechanism. (A) Downregulated USP5 ubiquitination leads to increased USP5 expression, which promotes RIG-I K11 ubiquitination modification and inhibits the IFN-I pathways by recruiting STUB1. (B) K48 ubiquitin of IQGAP1 is downregulated, resulting in increased IQGAP1 expression and inhibition of IFN-I production. (C) The increase in Akt may affect the downstream mTOR protein through a multistep stimulatory modification process, which in turn acts on STAT3 and reduces the ubiquitination of STAT3. The reduction of JAK1 ubiquitination may increase the expression of Akt protein and reduce the ubiquitination of STAT3 protein.

Discussion

By functioning as a molecular regulator, ubiquitination is involved in a variety of physiological processes, including kinase activation, protein interaction, cell cycle regulation, and cell signal transduction.[27,28] More relevant to the current study, the UPS, ubiquitination, and deubiquitination are necessary for various stages of infection for numerous viruses, including coronavirus.[15,29] For instance, evidence has shown that coronavirus has evolved various strategies to antagonize the IFN-induced JAK/STAT pathway, leading to inhibition of the host innate immune response.[30−32] SARS-CoV protein 6 can mediate ubiquitin-dependent proteasomal degradation of N-Myc and STAT interactors, thus suppressing the IFN signaling in host cells to promote SARS-CoV survival.[33] In this study, we described for the first time the global profiling of ubiquitinated proteins in SARS-CoV-2-infected Vero E6 cells. We identified 966 ubiquitination sites in 551 proteins that were differentially expressed following SARS-CoV-2 infection, 138 Kub sites on 104 proteins, and 828 Kub sites on 447 proteins were upregulated and downregulated. Further analysis of the functions of these proteins and sites can help to better understand the pathogenesis of SARS-CoV-2.

Type I IFNs play important roles in antiviral innate immune responses and in the modulation of the adaptive immunity. The host response to and elimination of viral infections largely depend on the expression of type I IFNs.[34] SARS-CoV-2 effectively suppresses the activation of TRAF 3 and 6, thereby limiting the activation of the transcription factors NF-κB and IRF 3 and 7, resulting in suppression of early pro-inflammatory responses through type I IFN.[35] Therefore, restriction of the anti-viral response is a contributing factor for regulating the host immune response during SARS-CoV-2 infection. Our results showed that after SARS-CoV-2 infection, the ubiquitination level of USP5 is downregulated, resulting in increased expression of USP5, and the increase in USP5 expression inhibits IFN-I signaling. The analysis of virus–host interactome by Li et al.[36] showed that the nsp1 protein of SARS-CoV-2 could affect melanoma differentiation-associated protein 5 and RIG-I RNA sensor signaling. This is consistent with our results, which is important to further understand how SARS-CoV-2 manipulates the host cell machinery.

In addition, IQ motif containing GTPase activating protein 1 (IQGAP1) and mitochondrial protein phosphoglycerate mutant family member 5 (PGAM5) have a similar function to USP5 in regulating type I IFN. The scaffolding protein IQGAP1 is a negative regulator of the type I IFN signaling pathway.[37] In our study, SARS-CoV-2-infected cells showed remarkable reductions in the ubiquitination of K348, K1445, K1465, K1475, K1477, and K1541 in IQGAP1. Proteomics data showed that the expression of IQGAP1 was increased, indicating that SARS-CoV-2 infection might decrease the degradation of IQGAP1 by altering its K48-linked ubiquitination modification, thereby inhibiting the production of type I IFN (Figure B). PGAM5 is a positive regulator of IFN-β production in response to viral infection, which is mediated by the TBK1/IRF3 signaling pathway.[38] Increased ubiquitination of K84 and K137 on PGAM5 was observed, while protein expression was reduced by 0.72-fold, suggesting that SARS-CoV-2 might increase the K48-linked polyubiquitination on PGAM5 to target its proteasome degradation, thus inhibiting IFN-β production. Therefore, we speculate that SARS-CoV-2 may manipulate the ubiquitination of important regulators to suppress the host immune responses.

Cytokines play essential roles in the human immune response against pathogenic infection.[39] However, excessive or uncontrolled levels of cytokines may result in severely maladjusted immune responses.[40−42] Cytokine storms are thought to occur in a wide variety of infectious diseases caused by viruses such as SARS-CoV and H5N1.[43−45] Patients with SARS have higher levels of IFN-α, IFN-γ, and IFN-stimulated chemokines as well as higher expression levels of antiviral IFN-stimulated genes than healthy controls.[46] In patients with COVID-19, the initial plasma concentrations of IL1B, IL1RA, IL7/8/9/10, basic FGF, IFN-γ, IP10, MCP1, MIP1A/B, PDGF, TNFα, and VEGF were higher than those in healthy adults.[47]

When SARS-CoV-2 infects the upper or lower respiratory tract, it stimulates the body to release pro-inflammatory cytokines, including IL-1β and IL-6.[48] Abnormalities in the IL-6 signaling pathway are involved in the pathogenesis of autoimmune inflammatory diseases.[49] It has been reported that endogenous TRIM28 negatively regulates the production of inflammatory factor IL-6, which is induced by TNF-α.[50,51] IL-6 signaling reduces JAK1/STAT3 ubiquitination, affecting the IL-6/STAT3 signaling pathway cascade. In addition, we found that the Akt protein expression was significantly upregulated, while the ubiquitination level of LAMTOR1 (late endosomal/lysosomal adapter, MAPK, and MTOR activator 1) decreased. These changes might lead alterations in the downstream protein mTOR, which in turn acts on STAT3, resulting in decreased STAT3 ubiquitination. However, there were no changes in the expression levels of JAK1, STAT3, and LAMTOR1 proteins. These data suggest that SARS-CoV-2 may affect this signaling pathway by regulating K63-linked polyubiquitination. We speculate that the decrease in JAK1 ubiquitination might increase Akt protein expression and reduce STAT3 protein ubiquitination. Moreover, the increase in Akt expression and the decrease in LAMTOR1 ubiquitination might have an impact on mTOR, which in turn acts on STAT3, leading to a decrease in STAT3 ubiquitination. All these activities lead to abnormal IL-6 signal transduction and inflammatory responses (Figure C).

It has been reported that the mechanisms underlying the SARS-CoV-induced cytokine storm include activation of NF-κB by the ORF3a protein of SARS-CoV and activation of the NLRP3 inflammasome by promoting TRAF3-dependent ubiquitination of ASC.[14] Activation of the NRLP3 inflammasome requires two steps. The first step, known as the initiation step, is induced by PRR or TNF receptor (TNFR). During this step, NF-κB is activated to promote the expression of pro-IL-1β, IL-18, and NLRP3. The second step, known as the activation step, is triggered by a series of stimuli, including ATP,[52] that are induced during infection, tissue injury, or metabolic imbalance. The leucine-rich repeat domain in the C-terminus associates with Hsp90, SGT1, and PML and is believed to be involved in the regulation of NLRP3 inflammasome activity.[53,54]

In our study, the quantitative proteomics data showed that the expression levels of ATP increased significantly, while the expression levels of PRR and TNFR decreased. However, no changes were observed in the expression levels of NLRP3, pro-IL-1β, and IL-18, which indicated that NLRP3 was not affected by changes in PRR and TNFR expression at the initial stage. These data suggest that SARS-CoV-2 may stimulate the activation of NLRP3 inflammasome through itself or the increase of ATP expression in cells during the activation phase of NLRP3 inflammasome. In addition, KEGG pathway analysis revealed that SARS-CoV-2 infection caused upregulation of the NOD-like receptor signaling pathway. Further analysis of the data revealed that HSP90 ubiquitination was significantly downregulated at the K69, K107, K219, K275, K284, K354, K399, and K607 sites and was significantly upregulated at the K411 site. In addition, ubiquitination of SGTA at the K160 and K200 sites significantly decreased. The SGTA and SGT1 proteins belong to the SGT family and contain tetratricopeptide repeat domains. However, the expression levels of HSP90 and SGTA were not significantly affected. These data indicated that SARS-CoV-2 might activate the NOD-like receptor signaling pathway by regulating K63-linked polyubiquitination, thereby activating the NLRP3 inflammasome and leading to an inflammatory cytokine storm. Based on these results, we speculate that SARS-CoV-2 infection might trigger the K63-linked polyubiquitination of HSP90 and SGTA by increasing ATP to activate the NOD-like receptor signaling pathway and facilitate activation of the NLRP3 inflammasome, thus inducing an inflammatory cytokine storm.

The highly glycosylated trimeric S protein mediates virus entry into host cells. The S protein is composed of an N-terminal receptor-binding S1 domain and a C-terminal S2 domain. The S1 domain is subdivided into an NTD (amino acids 13–304), a RBD (amino acids 319–541), and two conserved subdomains (SD1 and SD2). When the RBD binds to the host cell receptor ACE2, it destabilizes the pre-fusion trimer, causing shedding of the S1 subunit and transition of the S2 subunit to a stable post-fusion conformation.[55−57] We found that the K386, K417, K424, K458, K462, K529, and K535 sites in the RBD domain were ubiquitinated (Figure and Tables and S6), suggesting that these ubiquitination sites might affect the structure of the S1 subunit and promote binding of the RBD to ACE2.

Coronavirus N proteins are typically composed of three distinct but highly conserved components, an N-terminal RNA-binding domain (NTD), a CTD, and a poorly structured central Ser/Arg (SR)-rich linker. Previous studies have revealed that NTD is responsible for RNA binding, CTD is responsible for oligomerization, and SR-rich linker is responsible for primary phosphorylation.[58−60] Our data indicated that K61, K65, K143, K169, K233, K266, K342, K347, K355, and K361 were ubiquitinated (Figure and Tables and S6), four of these sites are in the NTD (amino acids 46–174), and five were in the CTD (amino acids 247–365). We speculate that ubiquitination of K61, K65, K143, and K169 may promote binding of the NTD to RNA, assembling the viral RNA genome into a ribonucleoprotein complex, packaging the complex into virions, and facilitating replication of the virus in the body through interactions with M protein.

SARS-CoV-2 encodes three putative ion channels: E, 8a, and 3a. Studies of SARS-CoV-1 implicated 3a in viral release, inflammasome activation, and cell death, and in animal models, its deletion reduces the viral titer and morbidity. Therefore, 3a is considered a potential therapeutic target for SARS and COVID-19. SARS-CoV-2 protein 3a is predicted to have three transmembrane helices followed by a cytoplasmic domain with multiple β-strands in each protomer chain, which generally exists as a dimer or tetramer.[61] We found that the K16, K21, and K132 sites in the 3a protein were modified by ubiquitination (Figure and Tables and S6). Although the function of the N-terminus (amino acids 1–39) is unknown, it is impossible to speculate on the effect of ubiquitination of K16 and K21 on the structure and function. However, since K132 is located in the transmembrane region of the protein, these ubiquitination sites are thought to be associated with the formation of ion channels or the changes in ion channel pores.

Coronavirus employs multi-subunit machinery for replication and transcription. Coronavirus replication involves a set of nonstructural proteins encoded by open reading frame 1a (ORF1a), which are initially translated as polyproteins and then undergo proteolytic cleavage to become mature proteins. RNA-dependent RNA polymerase (RdRp, also known as nsp12) is a key component in the process that catalyzes the synthesis of viral RNA, albeit with extremely low efficiency, Nsp7 and nsp8 function as cofactors that significantly stimulate polymerase activity. Thus, the nsp12–nsp7–nsp8 complex may play an essential role in the replication and transcription cycle of SARS-CoV-2. In the polymerase complex, nsp12 binds to a nsp7–nsp8 heterodimer and an additional nsp8 subunit.[62,63] In this study, we found that the K53 site of the nsp8 protein underwent ubiquitination (Figure and Tables 1 and S6). We speculate that K53 ubiquitination of nsp8 may promote the binding of nsp8 to nsp12 and nsp7–nsp8 to enhance the enzymatic activity of nsp12 and facilitate rapid viral replication in host cells. Proteins that function in both virus replication and inhibition of the host innate immune response have also been identified.[64] Studies have shown that ubiquitination of viral proteins can regulate the processes of replication and other functions. Thus, we speculate that ubiquitination is very significant in controlling the life cycle of SARS-CoV-2 and the functions of viral proteins.

In conclusion, we described for the first time the global ubiquitination changes in the cellular proteins of SARS-CoV-2-infected Vero E6 cells. These data elucidate the role of ubiquitination in the pathogenesis of SARS-CoV-2 infection. This study provides new insights into the interaction between SARS-CoV-2 and the host and a theoretical basis for the prevention and treatment of COVID-19.