Variant- and vaccination-specific alternative splicing profiles in SARS-CoV-2 infections

Abstract

The COVID-19 pandemic, driven by the SARS-CoV-2 virus and its variants, highlights the important role of understanding host-viral molecular interactions influencing infection outcomes. Alternative splicing post-infection can impact both host responses and viral replication. We analyzed RNA splicing patterns in immune cells across various SARS-CoV-2 variants, considering immunization status. Using a dataset of 190 RNA-seq samples from our prior studies, we observed a substantial deactivation of alternative splicing and RNA splicing-related genes in COVID-19 patients. The alterations varied significantly depending on the infecting variant and immunization history. Notably, Alpha or Beta-infected patients differed from controls, while Omicron-infected patients displayed a splicing profile closer to controls. Particularly, vaccinated Omicron-infected individuals showed a distinct dynamic in alternative splicing patterns not widely shared among other groups. Our findings underscore the intricate interplay between SARS-CoV-2 variants, vaccination-induced immunity, and alternative splicing, emphasizing the need for further investigations to deepen understanding and guide therapeutic development.

Keywords: Immunology; Molecular biology; Virology.

Graphical abstract

Figure 1

SARS-CoV-2 infection leads to aberrant global alternative splicing (A) Principal component analysis (PCA) plot of exon inclusion level from 190 RNA-seq dataset. (B) Number of differential alternative splicing events (DASEs) of five alternative splicing types. The DASEs were defined as following criteria: absolute PSI differences value >0.1 and corrected p value <0.05. (C) Results of the top 5 terms of GO biological process and KEGG pathways enriched with the 1,928 genes from 3,381 DASEs (HC vs. All patients). (D) Venn diagrams displaying the gene overlap between DASEs in the SE and MXE categories, respectively, and the set of SARS-CoV-2 susceptible genes. (E) Heatmap of differentially alternative spliced SARS-CoV-2 susceptible genes between HC and COVID-19 patients. Z score indicates relative exon inclusion levels. Hierarchical clustering of DASEs was performed with a Euclidean distance matrix of relative exon inclusion levels. (F) Significant differential spliced events of TLR4 (chr9:117,704,403-117,724,735) and JAK3 (chr19:17,824,782-17,847,982) showing two skipped exon (SE) events and six mutually exclusive exon (MXE) events, respectively. Exon inclusion levels represent the usage of spliced exons in the case of SE events, while in the context of MXE events, they indicate the ratio of the second mutually exclusive exon within each event. The expression level refers to the overall gene expression level. (G) Correlation analysis of expression differences between the inclusion level differences of total genes.

Figure 2

Dysregulated alternative splicing-related genes in COVID-19 patients (A) PCA plot of gene expression levels from 190 samples. (B) Results of the top 10 pathways enriched with the 7,529 DEGs. (C) Heatmap of 92 DEGs involved in RNA splicing (GO:0008380). Z score indicates relative gene expression levels. Hierarchical clustering of DEGs was performed with a Euclidean distance matrix of relative gene expression levels. (D) Protein-protein network of significant DEGs related to RNA splicing. The color of each gene indicates statistically significant and fold change to HC. (E) Gene expression levels of six hub genes which are depressed in COVID-19 patients. The error bar indicates standard deviation of the mean.

Figure 3

Differential regulation of alternative splicing in T and B cell receptor signaling genes among vaccinated and unvaccinated Omicron-infected patients (A) PCA plot of exon inclusion levels from 190 samples, as presented in Figure 1A. The samples are marked differentially based on the vaccination status of Omicron-infected patients. Group 1 and 2 represent clusters established through k-means clustering. (B) Vaccination status ratio of Omicron-infected patients within group 1 and group 2. A chi-squared test was conducted to confirm a significant difference in vaccination status proportions between the two groups. (C) Volcano plot showing the DASEs between Omicron groups. The x- and y axis indicate inclusion level difference (ΔPSI) and negative log₁₀ transformed p value. The q-value indicated corrected p value. (D) Results of the top 5 terms of GO biological process and KEGG pathways enriched with the 1,789 genes from 2964 DASEs (Omi.group 1 vs. Omi.group 2). (E) A subset of shared pathways based on 16 common DASEs identified in the T cell receptor signaling pathway (hsa04660) and the B cell receptor signaling pathway (hsa04662). The coloration of each gene box signifies the expression difference between the two groups, with solid lines indicating a direct effect and dashed lines representing an indirect effect. (F) Significant differential spliced events of KRAS (chr12:25,205,246-25,250,929) and NFATC2 (chr20:51,386,963-51,542,719) showing two MXE events, respectively. Inclusion levels represent the usage of spliced exons. The error bar indicates standard deviation of the mean.

Figure 4

Patients infected with different variants exhibit differential regulation of splicing machinery genes (A) PCA plot of 304 genes related to RNA splicing across 190 samples. Each sample in HC and COVID-19 patients (left panel) and Omicron patients with vaccination status (right panel) is marked as dots. (B) Biplot results to identify representative genes with high contributions in each group. Blue arrows represent genes representing the Omicron patient group, red indicates genes representing the Alpha and Beta infection patient groups, and black represents genes representing the HC group. Each gene marked with an arrow on the Biplot is displayed as a heatmap, with values indicating log2 fold change of gene expression levels against HC. Hierarchical clustering was performed with a Euclidean distance matrix of fold changes.