SARS-CoV-2 variants evolve convergent strategies to remodel the host response

Abstract

SARS-CoV-2 variants of concern (VOCs) emerged during the COVID-19 pandemic. Here, we used unbiased systems approaches to study the host-selective forces driving VOC evolution. We discovered that VOCs evolved convergent strategies to remodel the host by modulating viral RNA and protein levels, altering viral and host protein phosphorylation, and rewiring virus-host protein-protein interactions. Integrative computational analyses revealed that although Alpha, Beta, Gamma, and Delta ultimately converged to suppress interferon-stimulated genes (ISGs), Omicron BA.1 did not. ISG suppression correlated with the expression of viral innate immune antagonist proteins, including Orf6, N, and Orf9b, which we mapped to specific mutations. Later Omicron subvariants BA.4 and BA.5 more potently suppressed innate immunity than early subvariant BA.1, which correlated with Orf6 levels, although muted in BA.4 by a mutation that disrupts the Orf6-nuclear pore interaction. Our findings suggest that SARS-CoV-2 convergent evolution overcame human adaptive and innate immune barriers, laying the groundwork to tackle future pandemics.

Keywords: SARS-CoV-2; innate immunity; protein-protein interactions; proteomics; systems biology; transcriptomics; variants; virus-host interactions.

Figure 1.. VOCs impact RNA and protein landscape during infection.

(A) Global COVID-19 case numbers (log₂ transformed) over time, annotated for each SARS-CoV-2 variant of concern (VOC) based on sequences from GSAID. Thin lines reflect raw counts, and thick lines represent a LOESS fit. Inset: y-axis is displayed in linear scale. (B) Number of protein coding and non-coding mutations for each VOC. Data were from covariants.org on Jan 5, 2022, and correspond to Alpha 20I V1, Beta 20H V2, Delta 21J, Gamma 20J V3, and Omicron 21K. (C) Number of protein coding mutations per protein in the VOCs. (D) Experimental workflow. Infected Calu-3 cells were harvested at 10 and 24 hours post-infection (hpi) and processed for bulk mRNA sequencing and global mass spectrometry abundance proteomics and phosphoproteomics. We additionally performed affinity purification mass spectrometry (AP-MS) on individually overexpressed VOC and W1 viral proteins in HEK293T cells to quantify changes in virus-host protein interactions. (E) Quantitative reverse-transcription PCR (qRT-PCR) for viral E gene copies to quantify viral replication over time for each experiment. Calu-3 lung epithelial cells were infected with 2000 E copies/cell of each SARS-CoV-2 VOCs, VIC, IC19 or mock. Stars indicate significant difference relative to time-matched VIC (adjusted p-value < 0.05). (F) Viral replication over time in experiments 1 and 2 based on Orf1a leader sequence-containing counts from bulk mRNA sequencing. (G) Quantification of the sum of non-structural protein intensities from abundance proteomics for each virus in experiment 1 or 2 at 10 and 24 hpi in Calu-3 cells. (H) Flow cytometry assessing the percentage of cells positive for nucleocapsid (N) staining for each virus at 10 and 24 hpi in Calu-3 cells. Error bars represent standard error (SE). (I) Fraction of mRNA, protein, or phosphorylation sites that change (black; abs[log2FC]>1 & q<0.05) in response to infection with at least one virus, for at least one time, in at least one experiment, relative to mock. Fraction of protein-protein interactions from AP-MS data that change between VOC and W1 viral proteins (abs[log2FC]>0.5 & p<0.05; right bar). (I) Number of mRNA transcripts, proteins, or phosphorylated peptides that increase or decrease for each condition and time, compared to mock. Numbers in parentheses indicate experiment number.

Figure 2.. Convergent molecular strategies of VOCs.

(A) Comparative systems omics revealed VOCs converge on key molecular strategies to remodel the host environment by altering viral gene expression, viral protein phosphorylation, and virus-host protein complexes relative to W1 viruses. (B-C) Normalized read counts from bulk mRNA sequencing containing the leader sequence and mapping to a portion of the SARS-CoV-2 genome (B) or normalized protein intensities from abundance proteomics (C) at 24 hpi in Calu-3 cells. RNA quantities are normalized to Orf1a genomic (log2(counts/genomic/W1)) per virus and viral protein intensities to the summed intensity of Nsps (log2(intensity/sum of nsps/W1)) to control for differences in viral replication and defined relative to W1 virus VIC. Experiments 1 and 2 were integrated after normalization. Stars indicate abs(log2FC)>1 & adjusted p-value < 0.05. (D) Absolute value log₂ VOC- and VIC-normalized viral protein quantities from (C) for each VOC. Colored dots indicate proteins regulated more than twofold in expression, compared to VIC. Dashed line denotes twofold expression change. All VOCs encode at least one viral protein expression change relative to VIC and IC19. (E) SARS-CoV-2 Orf9b viral RNA and protein for Alpha and Delta isolated from (B) and (C). (F) Genomic sequence of the region surrounding the start of the N gene. “−3” indicates three nucleic acid positions upstream of the N translation start (“N start”). Orf9b translation start (“Orf9b start”) is indicated within the N coding region. Mutations in Alpha (−3 deletion and N protein D3L) and Delta (−3 deletion only) are indicated in red. The Orf9b transcriptional regulator sequence (TRS) is also indicated and is thought to be enhanced by the GAU→CUA. (G) Depiction of how mutations colored red in (F) may affect transcription and translation of Orf9b. (H) Expression of viral sgRNA (from bulk mRNAseq) and protein (bulk proteomics) from mutant viruses (derived via reverse genetics approaches) at 48 hpi in Calu-3 (MOI=0.01): one with a N D3L mutation (GAU→CUA), one with a −3 deletion upstream of N start, and one with both. Quantifications are normalized to genomic (RNA, log2(counts/genomic/W1)) or summed intensity of all Nsps (protein, log2(intensity/sum of nsps/W1)) per virus to control for differences in viral replication and defined relative to wave 1 virus from Washington, USA (USA-WA1/2020) (I) Absolute value t-statistic from t-tests comparing phosphorylated peptide intensities between all possible pairs of viruses (see Table S2 and Methods). Each dot represents one phosphopeptide compared between two viruses. Intensities are normalized by corresponding total protein abundance. Comparisons were restricted to viral peptides with identical sequences between virus pairs, given that peptide intensities of peptides with different sequences are not directly comparable using mass spectrometry. If p≤0.001, dots are colored black, otherwise they are grey. (J) The t-statistic from t-tests in (I) restricted to phosphorylation sites on N and relative to VIC. If p≤0.001, dots are colored black, otherwise they are grey. All VOCs, except Gamma, show evidence of remodeling N protein phosphorylation relative to VIC. (K) Significantly changed phosphorylation sites (p≤0.001) on N protein between pairs of viruses. Length of each lollipop depicts the abs(t-statistic) between the pairs of viruses. (L) In vitro ADP-GLO kinase activity assay of 122 recombinant kinases (predicted to phosphorylate N protein sequence based on GPS 5.0) incubated with full-length recombinant W1 SARS-CoV-2 N protein.Y-axis depicts the log2 fold-change between the kinase incubated with N or alone. Red indicates our positive control, SPRK1 and GSK3β co-incubated with full-length N . Labeled kinases indicate those with greater activity against N than our positive control and at least half the magnitude of the canonical positive control substrate for that kinase. (M) Quantified changes in protein-protein interactions for 127 protein-coding mutations in 16 mutated viral proteins across all 5 VOCs. Left, number of high confidence virus-human interactions for all mutant VOC and wave 1 (W1) viral proteins using APMS. Of 1746 interactions, 1473 were unchanged (gray), 150 increased in binding affinity with at least one mutant (red), and 123 decreased in binding affinity with at least one mutant (blue). Significant increase in mutant binding affinity (red) is defined as log₂ fold-change>0.5 & p<0.05. Significant decrease in mutant binding affinity (blue) is defined as log₂ fold-change<−0.5 & p<0.05. Right, same but broken down by VOC. (N) Virus (“bait”, diamonds) host (“prey”, circles) protein-protein interaction (PPI) map for N*, Orf9b, and Orf6 depicting significantly changing interactions (absolute value log₂ fold-change>0.5 & p<0.05) comparing VOC to W1 forms (see Fig 3 and Table S3).

Figure 3.. VOC evolution rewires virus-host protein-protein interactions.

Virus (“bait”, diamonds) host (“prey”, circles) protein-protein interaction (PPI) map for significantly changing interactions (absolute value log₂ fold change>0.5 & p<0.05) comparing VOC to W1 forms (see Table S3). Color of edge represents log₂ fold change in the abundance of each human prey protein in the affinity purification as determined by mass spectrometry, comparing VOC and W1 forms. Multiple edges are displayed when the same affinity purification was performed multiple times and both results were significantly differentially interacting. Black edges indicate human-human protein complexes annotated by CORUM , also highlighted and annotated using yellow shading. Biological processes are indicated using cyan shading.

Figure 4.. Integrative computational analysis reveals conservation and divergence of host response to variants.

(A) Number of RNA transcripts, proteins, or phosphorylation sites that significantly changed during VOC infection, compared to mock at 24 hpi. For transcriptomics and phosphoproteomics, we required absolute value log₂ fold-change log2FC >1 and p<0.001. For abundance proteomics, we required abs(log2FC)>log₂(1.5) & q<0.05. For each dataset, a molecule had to pass the threshold twice at either times, viruses, or experiments. (B) Flowchart of computational pipeline. Host genes regulated during infections from (A) were extracted from the STRING network and clustered into 85 pathway modules based on a diffusion measure of network node proximity (see Methods). (C) Average absolute value log₂ fold-change versus mock for each module (gray lines) using RNA or phosphorylation data. Virus conditions are ordered by their timeline of emergence around the world. The purple line defines the average across the module intensities to define the AHR. (D) Pearson’s R correlation between viral genomic (leader + Orf1a) RNA counts (R=0.59) or sum of non-structural protein (Nsps) intensities (R=0.77) versus VIC plotted against the AHR defined in (C). Each dot represents one virus condition at 24 hpi. (E) Pearson’s R correlation between the average log₂ fold-change of each module, per dataset, and the RNA-derived AHR. Red dots indicate 10 modules most correlated with the AHR, based on a geometric mean across RNA, abundance proteomics, and phosphoproteomics datasets (“composite R”). Blue dots indicate 10 modules least correlated with AHR. (F) 10 most (red; highest correlation coefficients) and 10 least (blue; lowest correlation coefficients) correlated modules with AHR, annotated by the most prevalent GO Biological Process term, module number, and number of genes within the module that connect to the top GO term. The x-axis depicts the composite R value (defined in E). Colored numbers indicate ranking of modules based on composite R. Terms in bold highlight prevalent biological categories: translation related terms in top 10 and innate immune/inflammation related terms in the bottom 10. Red modules represent pathways similarly regulated by all variants (“less variant specific”). Blue modules represent pathways differently regulated across the variants (“more variant specific”). (G) t-Distributed Stochastic Neighbor Embedding (t-SNE) plot representing the STRING network proximity between genes, colored according to the module annotation. Top (red) and bottom (blue) 10 modules are bolded, and their locations are annotated using contours. (H) Innate immune- and inflammation-related modules within the 10 least correlated with the AHR. (I) RNA and protein expression of interferon stimulated genes (ISGs) and RNA expression of proinflammatory genes at 24 hpi in Calu-3 cells. Expression is defined as the average log2FC of ISGs or proinflammatory genes (see Table S4 for list of genes) for each virus compared to VIC and averaged across batches. Error bars depict SE. Proinflammatory genes were sparsely detected at the protein level and excluded. (J) Western blot of MX1, IFIT1, IFITM3, and SARS-CoV-2 nucleocapsid (N) protein expression in Calu-3 cells infected with 2000 E copies/cell at 24 hpi. Protein quantification over β-actin is shown normalized to IC19 levels. (K) Correlation of viral genomic (leader + Orf1a) counts and ISG expression (versus VIC) over time. Viral counts are represented as transcripts per million. Times 10 and 24 hpi are represented by circles and triangles, respectively. (L) Viral replication of Alpha, Delta, or Omicron BA.1 in Calu-3 cells at 48 hpi with and without ruxolitinib, a JAK/STAT inhibitor. Fold-change between conditions is noted. (M) Relationship between ISG and proinflammatory gene expression for each virus at 24 hpi in Calu-3 cells relative to VIC for RNA. (N) Pearson’s R correlation between average log2FC for each module and levels of each viral protein across the viruses. Innate immune and inflammatory modules in the 10 least correlated category, including one additional innate immune term ranked 11th (m43, regulation of immune response), are colored. Viral proteins are ranked from left-to-right according to the average R values across the five inflammation-related terms.

Figure 5.. Omicron subvariants evolved innate immune antagonism by modulating Orf6.

(A) Experimental workflow. Calu-3 cells were infected with the indicated variants, and W1 control IC19. Cells were harvested at 10, 24, and 48 hpi for bulk mRNA sequencing. Cells were harvested at 24 and 48 hpi for abundance mass spectrometry proteomics analysis. (B) Counts of viral genomic (leader + Orf1a) RNA over time for each virus. (C) ISG (RNA and protein) and proinflammatory gene (RNA) expression at 48 hpi relative to BA.1. Proinflammatory genes were sparsely detected at the protein level and excluded. Error bars represent SE. (D) Heatmap of ISG protein expression for Omicron subvariants at 24hpi. Color indicates the log₂ fold-change in expression, relative to Omicron BA.1. (E) Same as D, but at 48hpi. (F) Relationship between ISG and proinflammatory gene expression for each virus at 48 hpi in Calu-3 cells, relative to BA.1 for RNA as in (C). (G) Expression of viral RNA for Omicron subvariants, normalized to viral genomic (leader + orf1a) counts for each virus and set relative to BA.1. Error bars represent SE. (H) Expression of viral protein for Omicron subvariants, normalized to the sum non-structural protein (Nsps) intensities for each virus and set relative to BA.1. Error bars represent SE. (I) ISG RNA (left) and protein (right) expression versus VIC for Alpha, Omicron BA.1, an Orf6 knock-out version of Alpha created using reverse genetics, and poly I:C. (J) AP-MS of Orf6 D61L (occurring in BA.2 and BA.4 but not BA.1 or BA.5) compared to wave 1 (W1) Orf6 in HEK293T cells. All detected proteins are plotted with high-confidence interactions that are also significantly differentially interacting (abs(log2FC)>0.5 & p<0.05) highlighted in blue. (K) X-ray crystallography structure (PDB 7VPG) of RAE1, NUP98, and Orf6. Orf6 D61 residue is pink and forms a hydrogen bond (gray sticks) with RAE1. Other Orf6 residues (E59, E56, D53, and E55) that participate in interactions with RAE1 are shown. M58 inserts into a RAE1 hydrophobic pocket. (L) Model of the effects of Orf6 levels and Orf6 D61L mutant status on the innate immune response. The nuclear pore (RAE1 and NUP98) physically interacts with Orf6, which suppresses the nuclear translocation of ISG-inducing transcription factors and export of ISG mRNAs during infection. This interaction is weakened, but not abolished, when the D61L mutation is present, resulting in diminished innate immune antagonism by Orf6. BA.1 and BA.2 downregulate Orf6 relative to early-lineage SARS-CoV-2 (IC19), which results in an increased innate immune response during infection, exacerbated by the presence of the D61L mutation in BA.2. Conversely, BA.4 and BA.5 upregulate Orf6 protein to similar levels. However, BA.5 more strongly antagonizes the innate immune response, which we speculate is due to the absence of the D61L mutation.

Figure 6.. VOCs balance adaptive and innate immune escape during the course of their evolution.

(A) Data on Spike sequence similarity relative to W1 virus (green) and ISG expression (this paper, red) are shown for each variant. Impact of viral protein expression and point mutations are indicated (red and green text). Importantly, the line does not imply sequential evolution. (B) Model for SARS-CoV-2 VOC convergent strategies to modulate the immune responses. Spike mutations increase evasion from the adaptive immune system (top). In this work, we discovered VOC mutations that enhance viral protein expression and rewire virus-host protein interactions that modulate the innate immune response. A coordinated balance between adaptive and innate immune evasion is required by successful variants. Combination therapeutic strategies that enhance adaptive and innate immunity may improve disease outcomes for patients.