Integrative systems immunology uncovers molecular networks of the cell cycle that stratify COVID-19 severity

Abstract

Several perturbations in the number of peripheral blood leukocytes, such as neutrophilia and lymphopenia associated with Coronavirus disease 2019 (COVID-19) severity, point to systemic molecular cell cycle alterations during severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) infection. However, the landscape of cell cycle alterations in COVID-19 remains primarily unexplored. Here, we performed an integrative systems immunology analysis of publicly available proteome and transcriptome data to characterize global changes in the cell cycle signature of COVID-19 patients. We found significantly enriched cell cycle-associated gene co-expression modules and an interconnected network of cell cycle-associated differentially expressed proteins (DEPs) and genes (DEGs) by integrating the molecular data of 1469 individuals (981 SARS-CoV-2 infected patients and 488 controls [either healthy controls or individuals with other respiratory illnesses]). Among these DEPs and DEGs are several cyclins, cell division cycles, cyclin-dependent kinases, and mini-chromosome maintenance proteins. COVID-19 patients partially shared the expression pattern of some cell cycle-associated genes with other respiratory illnesses but exhibited some specific differential features. Notably, the cell cycle signature predominated in the patients' blood leukocytes (B, T, and natural killer cells) and was associated with COVID-19 severity and disease trajectories. These results provide a unique global understanding of distinct alterations in cell cycle-associated molecules in COVID-19 patients, suggesting new putative pathways for therapeutic intervention.

Keywords: COVID-19 severity; SARS-CoV-2; cell cycle associated molecules; proteomics; systems immunology; transcriptomics.

Figure 1

Interactive network of multi‐omics data integration of cell cycle‐associated molecules. (A) Upset plot showing the intersection of transcriptomic (swab and blood) and proteomic data obtained from all studies included in this work. (B) Network of cell cycle‐associated molecules from proteome and transcriptome data sets from COVID‐19 patients. The color of the nodes represents gene ontology (GO) biological processes according to the figure legend. The genes labeled in red represent the differentially expressed genes (DEGs) from blood transcriptome data sets. DEGs from swab transcriptomes are denoted in the dark blue, and differentially expressed genes (DEPs) from the proteome data sets are labeled in purple. Molecules overlapping in the transcriptomics and proteomics data sets are exhibited in black. Triangles pointing up represent upregulated molecules, triangles pointing down indicate downregulated molecules, and the diamond shape shows DEGs and DEPs upregulated or downregulated across different data sets. The whole network comprises 429 proteins and 10 516 direct physical interactions.

Figure 2

Systems‐level view of proteome changes related to the cell cycle. (A) Schematic overview showing the number and classification of patient and control cohorts of each data set used for the proteomic data analysis. (B and C) Network of cell cycle‐associated proteins (gray nodes) obtained across the study cohorts and biological processes (BPs; blue nodes) enriched by (B) upregulated or (C) downregulated differentially expressed genes (DEPs). Gray edges reflect the association between DEPs and BPs. The interaction network was visualized using NAViGaTOR. More prominent nodes represent larger gene sets. The size of the squares increases according to the number of proteins enriching each BP.

Figure 3

Cell cycle‐associated gene co‐expression modules present in blood transcriptomes of COVID‐19 patients. (A) Schematic figure showing the number of samples of the peripheral blood leukocyte (PBL) transcriptome data set (GSE157103) and (D) the peripheral blood mononuclear cells (PBMCs) transcriptome data set (GSE152418). Bubble heatmap showing the results of gene set enrichment analysis, indicating the module (M) activity in (B) COVID‐19 patients versus non‐COVID‐19 (patients with other respiratory illnesses) and (E) COVID‐19 patients compared with healthy controls. Circle size and color reflect the normalized enrichment score (NES), as determined by CEMiTool. (C) Interaction plot for M3 (GSE157103) and (F) M1 (GSE152418), which contains genes enriching different cell cycle‐associated pathways, exhibited by the bar plot at the right side of figure (C) and (F). The most connected genes (hubs) are highlighted inside rectangles. The node size is proportional to its degree of interactivity. The bar plot indicates the top 10 enriched pathways from the over‐representation analysis of module M3 (C) and M1 (F).

Figure 4

Integrative meta‐analysis revealed the predominance of cell cycle enriched pathways in blood leukocytes of COVID‐19 patients. This analysis integrated the studies of peripheral blood leukocyte (PBL) (GSE157103) and peripheral blood mononuclear cells (PBMCs) (GSE152418). (A) Principal component analysis (PCA) and (B) density plots show the batch effect adjustment for GSE157103 + GSE152418 through empirical Bayes regression (using ComBat). (C) Meta‐analysis results displayed by the volcano plot, which was based on average Log2 fold change (FC) and combined p‐value obtained using the Fisher's method. Small blue and red circles denote up‐ and downregulated genes, respectively. Black circles show meta‐significant genes associated with the cell cycle process, indicating that they are predominantly upregulated. (D) Ridgeline chart denoting the fold‐change distribution of top 20 enriched pathways by the meta‐significant genes, suggesting the predominance of cell cycle enriched pathways across the studies. Gene (small gray circles) expression distribution is based on average Log2 FC across the enriched biological process. (E) Interactive heatmap exhibits the expression pattern of significant genes associated with the cell cycle process obtained by the meta‐analysis.

Figure 5

Alterations of cell cycle signatures in the blood are not present in transcriptomes from nasopharyngeal swabs of COVID‐19 patients. (A) Meta‐analysis results displayed by the volcano plot, which is based on average Log2 fold change (FC) and p‐value from the meta‐analysis performed to obtain combine p‐values from swab transcriptomes (GSE152075 and GSE156063 data sets) using the Fisher's method. Small blue and red circles denote up‐ and downregulated genes, respectively. Black circles show significant genes that are associated with the cell cycle, indicating that these genes are primarily downregulated in swab transcriptomes. (B) Ridgeline chart denoting the fold‐change distribution of the top 20 enriched pathways by the meta‐significant genes, indicating the predominance of the inflammatory response, immune and defense activation, while no cell cycle enriched pathways across the swab studies. Gene (small gray circles) expression distribution is based on average Log2 FC across the enriched pathways. (C) Circular heatmaps of a set of genes involved in each phase of the cell cycle in swab and blood (peripheral blood leukocytes [PBLs] and peripheral blood mononuclear cells [PBMCs]) transcriptomes. Color scale refers to up‐ (red) and downregulated (blue) genes. Gray fields indicate genes not differentially expressed in the data sets (list of all genes described in Supporting Information: Table S9). (D) The number of upregulated cell‐cycle associated DEGs across different leukocyte subpopulations from single‐cell RNAseq data sets (Supporting Information: Table S11). Each data set's GSE number and cell subpopulations are shown below the bar plot.

Figure 6

Expression of cell cycle‐associated genes upregulates during severe COVID‐19 disease phases and normalizes with disease recovery. (A) Heatmap showing the average expression of 126 cell cycle‐associated genes (list of genes [see Supporting Information: Table S13] obtained by meta‐analysis of peripheral blood leukocyte [PBL]: GSE157103 and peripheral blood mononuclear cells [PBMCs]: GSE152418 data sets) in different disease stages of COVID‐19 patients (pseudotimes 0–7 according to figure legend) using the previously published gene expression data from Bernardes et al. Row‐wise z‐scores of the scaled mean expression per pseudotime are plotted and hierarchically clustered in the heatmap. Differentially expressed genes (DEGs) in the longitudinal analysis are denoted with an asterisk (*), and DEGs resulting from the comparison of COVID‐19 versus healthy controls are denoted by the + sign according to the adjusted p‐value (*/+p ≤ 0.05, **/++p ≤ 0.01, ***/+++p ≤ 0.001). (B) Box plots illustrate normalized counts distribution across pseudotimes (0–7) for seven significant genes (significance level shown in (A)). Each box plot shows the median with the first and third interquartile range (IQR), whiskers representing minimum and maximum values within IQR, and individual data points.

Figure 7

Cell cycle‐associated genes stratify COVID‐19 from other respiratory illnesses. (A) Principal component analysis (PCA) of meta‐significant cell cycle‐associated genes (list of genes [see Supporting Information: Table S14] obtained by meta‐analysis of peripheral blood leukocyte [PBL]: GSE157103 and peripheral blood mononuclear cells [PBMCs]: GSE152418 data sets) showing the stratification of COVID‐19 from non‐COVID‐19 patients with other respiratory illness (patient cohorts from data set GSE157103 are described in Supporting Information: Figure 4); both groups containing patients admitted to an intensive care unit (ICU) or not (non_ICU). Confidence ellipses are shown for each group. Density plots associated with the PCA indicate the sample distribution across the PCA axes. (B and C) Variable importance score plot obtained by random forest classification analysis. The importance score plot is based on the Gini decrease and number (no.) of nodes for each variable (cell cycle‐associated meta‐DEGs). It indicates the top 10 variables predicting COVID‐19 severity when comparing (B) COVID‐19_ICU versus COVID‐19_non_ICU and (C) COVID‐19_ICU versus non‐COVID‐19_ICU. Small blue circles show the top 10 gene rankers of COVID‐19 severity while the black circles represent those below (with less importance) the top 10. (D and E) Stable curve indicating the number of trees and out‐of‐bag (OOB) error rate of the random forest analysis. (F and G) Receiver operating characteristic (ROC) curves with an area under the curve (AUC) exhibiting the relationship between true and false positive classification rates.

Figure 8

Cell cycle‐associated genes show a strong correlation in severe COVID‐19. (A) Correlation analysis of top cell cycle‐associated rankers (obtained from random forest analysis: Figure 7B,C) according to COVID‐19 severity (COVID‐19 and non‐COVID‐19 [other respiratory illness] cohorts in an intensive care unit [ICU] or not [non_ICU]: from data set GSE157103). The Spearman's rank correlation coefficient is shown according to the color scale bar ranging from −1 to 1. (B) Box plots illustrating the distribution of absolute values (considering all correlation coefficients despite the correlation direction [+ or − sign]) of Spearman correlation coefficients across study cohorts. Each box plot shows the median with the first and third interquartile range (IQR), whiskers representing minimum and maximum values within IQR, and individual data points. Significance was determined using two‐sided Wilcoxon rank‐sum test and is indicated by asterisks (*p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001, and ****p ≤ 0.0001). (C) Heatmap of correlation index. The bars aside from the heatmap represent the sum of correlation indexes. The color scale bar represents the correlation index for each gene.