Immune cell communication networks and memory CD8+ T cell signatures sustaining chronic inflammation in COVID-19 and Long COVID

Single-cell RNA sequencing (scRNA-seq) data (27) of PBMCs from COVID-19 patients with varying disease severities were retrieved from Data of 197 patients admitted to Yale New Haven Hospital with COVID-19 between 18 March and 5 May 2020 which were previously described (28). This dataset includes samples from four groups: Healthy controls, Exposed individuals (close contacts not yet diagnosed), Outpatients (Infected), and Hospitalized patients (Severe), with two samples per group. The validation Bulk RNA-seq transcriptomic cohort was previously described and generated on the GPL24676 Illumina NovaSeq 6000 platform, encompassing 100 COVID-19 samples and 26 non-COVID-19 controls (29).

The raw unique molecular identifier (UMI) count matrix was processed using the R package Seurat (version 5.1.0) and converted into a Seurat object (30). Cells were filtered based on the following criteria: fewer than 1,000 detected RNA molecules, fewer than 200 or more than 10,000 detected gene features, mitochondrial gene proportions (percent.mt) exceeding 20%, and hemoglobin gene proportions (percent.hb) exceeding 90%. After filtering, 73,110 cells were retained for downstream analysis. The dataset was normalized using the ScaleData() function, and the top 2,000 highly variable genes were identified with FindVariableFeatures(). Principal component analysis (PCA) was performed using RunPCA(), selecting the top 10 principal components for dimensionality reduction. t-SNE and UMAP analyses were then conducted using RunTSNE() and RunUMAP(), respectively (31). Clustering was performed with FindNeighbours() and FindClusters(), with the resolution set to 0.3. Cluster-specific marker genes were identified using FindAllMarkers(), applying thresholds of |log2FC| > 0.25 and p-value < 0.05 (32). Finally, cluster marker genes were annotated using the scMayoMap package, integrating scMayoMap Database and lung tissue-specific data (tissue = 'lung') (33, 34).

To evaluate the distribution preferences of cell types across different tissue states (Healthy, Exposed, Infected, and Hospitalized), the calTissueDist() function from the sscVis package (version 0.1.0) was used to calculate the R_o/e ratio, which quantifies differences in cell type distribution among tissues. The statistical significance of these differences was assessed using the chi-squared test (method = "chisq"), and a cell type ratio matrix associated with tissue states was extracted. The results were visualized as a heatmap created with the ComplexHeatmap package (version 2.18.0). Color gradients were used to intuitively represent changes in R_o/e values, and symbolic annotations were added to heatmap cells to indicate the degree of preference: "+++" for R_o/e > 1, "++" for 0.8 < R_o/e ≤ 1, "+" for 0.2 ≤ R_o/e ≤ 0.8, "+/-" for 0 < R_o/e < 0.2, and "-" for R_o/e = 0. This approach provided a clear visualization of the significant differences in cell type distribution across tissue states (35).

To assess the activity of the Cytokine and Inflammatory gene sets at the single-cell level, the AddModuleScore() function was applied to calculate module scores for each cell, which were then mapped onto UMAP plots. Using the ggplot2 package (version 3.5.1), pie charts were generated to illustrate the activity proportions of these two gene sets across different cell types. To compare Cytokine Score and Inflammatory Score among different tissue states (Exposed, Infected, Hospitalized, and Healthy), boxplots were created, and statistical significance was evaluated using the Wilcoxon test. Additionally, T cell subpopulations were analyzed for their scores in Cytokine, Exhaustion, Inflammatory, and Regulatory Effector gene sets, with boxplots used to visualize score comparisons among different T cell subtypes (30).

T cell subpopulations were extracted using the subset function and reclustered into eight clusters. The sctype package (version 1.0) was used to annotate these clusters based on gene set scoring of specific T cell markers, identifying Effector CD8+ T cells, Memory CD4+ T cells, Memory CD8+ T cells, Naive CD4+ T cells, and Naive CD8+ T cells (36, 37).

To explore the functional characteristics of different T cell subpopulations, differential expression analysis and functional enrichment analysis were performed based on single-cell transcriptomic data. Differentially expressed genes (DEGs) for each T cell subset were calculated using the FindAllMarkers function in the Seurat package, with the top 1,000 genes showing significant upregulation selected based on average fold-change values. The selected DEGs were then subjected to Gene Ontology (GO; biological processes and cellular components) and KEGG pathway enrichment analysis using the compareCluster function in the ClusterProfiler package (38). Gene names were converted to ENTREZ IDs using the org.Hs.eg.db database, and statistically significant enriched terms (p-value < 0.05) were identified. Finally, the enrichment results were visualized as dot plots using the enrichplot package.

This study employed the CellChat R package, a specialized tool for inferring, analyzing, and visualizing cell-cell communication from single-cell RNA sequencing data (39). The analysis focused on T cell subsets, including Effector CD8+ T cells, Memory CD8+ T cells, Memory CD4+ T cells, and Naive CD8+ T cells. Overexpressed genes and ligand-receptor pairs were identified to calculate intercellular communication probabilities and infer signaling pathway networks. Visualizations, such as circular and bubble plots, were used to depict communication frequencies and strengths among subsets. Key pathways, particularly MHC-I, were analyzed in-depth, revealing the central role of Memory CD8+ T cells in immune modulation. These findings provide critical insights into the mechanisms of T cell-mediated immune regulation and their implications in COVID-19 progression (40).

Differential expression analysis was conducted on Memory CD8+ T cells to investigate gene expression differences across various pathological states (Healthy, Exposed, Infected, and Hospitalized). The analysis employed the Seurat package and the MAST method, comparing each of the three experimental groups to the Healthy control group. The process involved extracting the Memory CD8+ T cell subset, identifying differentially expressed genes (DEGs) using the FindAllMarkers function in Seurat, and filtering genes with adjusted p-values (p_val_adj) less than 0.05 and average log fold-change (avg_log2FC) greater than 0. Results were visualized using the ggplot2 package, generating volcano plots, bubble plots, and lollipop plots to depict the distribution, expression proportions, and commonly upregulated genes among groups (41).

Nine commonly used machine learning methods were applied to model COVID-19-related data, including Linear Discriminant Analysis (LDA), Flexible Discriminant Analysis (FDA), Logistic Regression, Naive Bayes, Support Vector Machine (SVM), Random Forest, Gradient Boosting Machine (GBM), Mixture Discriminant Analysis (MDA), and XGBoost (42). Through multiple randomized experiments assessing AUC in both training and testing datasets, XGBoost demonstrated the best stability and performance. Consequently, XGBoost was selected as the final diagnostic model (43). The xgboost package was employed to train the model, using 70% of the data for training and 30% for testing. Additionally, the rmda package was utilized for Decision Curve Analysis (DCA) to evaluate the model's net benefits at varying risk thresholds, highlighting its clinical application potential.

The shapviz package was employed to analyze SHAP values, providing interpretability for the XGBoost model. Key diagnostic genes, such as RPS26, RPS29, and RPL41, were identified through feature importance bar plots and beeswarm plots, quantifying their contributions to the prediction results. Furthermore, a force plot was generated for the highest-scoring positive sample, visualizing the positive and negative contributions of individual genes to the prediction outcome.

We obtained peripheral blood mononuclear cells (PBMCs) from eight individuals with confirmed COVID-19, covering a range of disease severities. After quality control procedures, including filtering based on mitochondrial gene content, transcript counts, and the number of detected genes, we retained 73,110 high-quality cells for downstream single-cell RNA sequencing (scRNA-seq) analysis (Figures 1A, B). Uniform Manifold Approximation and Projection (UMAP) embedding revealed 15 distinct cell clusters (Figure 1D). These clusters were assigned to 10 major immune and epithelial cell types based on canonical marker gene expression (Figures 1C, E), including T cells (CD3D, CD3G, CD3E), B cells (MS4A1, CD79A, CD79B), myeloid cells (CD14, LYZ, ITGAM), neutrophils (S100A8, S100A9, CSF3R), macrophages (CD68, CD163, MRC1), secretory cells (SCGB1A1, SCGB3A2), alveolar epithelial type I cells (AGER, HOPX, PDPN), alveolar epithelial type II cells (SFTPC, SFTPB, SFTPA1), plasma cells (MZB1, JCHAIN, TNFRSF17), and basal cells (KRT5, KRT14, TP63).

UMAP analysis compared T cell distributions across four disease states (Healthy, Exposed, Infected, Hospitalized), revealing a progressive decline in T cell proportions with increasing disease severity (Figure 2A). Analysis of cell-type enrichment ratios (R_o/e) showed significant enrichment of macrophages, neutrophils, and secretory cells in severe disease states, while T cells were notably depleted in the Infected and Hospitalized groups (Figure 2B). T cell enrichment values exhibited a marked decline from Healthy to Hospitalized states (Figure 2C), underscoring their potential role in disease progression. Cytokine expression analysis identified IL32, LTB, and MIF among the top 10 cytokines highly expressed in pathological conditions, with IL32 showing the highest expression across all groups, suggesting its pivotal role in inflammatory responses (Figure 2D). These findings collectively demonstrate dynamic changes in cell-type composition and highlight their functional implications in disease progression.

To investigate the role of various cell types in immune regulation and inflammatory responses, we conducted an analysis of cytokine and inflammatory scores. T cells exhibited the highest cytokine scores, constituting 34.7% of the spatial distribution (Figure 3A), underscoring their pivotal role in immune responses. Myeloid cells (27.9%) and lymphocytes (18.3%) also displayed significant cytokine activity, indicating their critical involvement in the pathological states. Inflammatory score analysis (Figure 3B) revealed that myeloid cells had the highest contribution (35.4%), followed by T cells (31.2%). Comparative analysis across different disease states (Figures 3C, D) showed that T cells' cytokine and inflammatory scores were significantly elevated compared to the Healthy group (p < 0.05), with the highest increase observed in the Hospitalized group. Additionally, neutrophils, macrophages, and secretory cells showed a progressive rise in scores from Exposed to Infected and Hospitalized states, suggesting their amplified activation in exacerbating inflammation. Conversely, alveolar epithelial and basal cells demonstrated minimal variations.

Through UMAP analysis, T cells were subdivided into eight clusters (Figure 4A) and annotated into classical subtypes based on specific markers (Figure 4C), including Effector CD8+ T cells, Memory CD8+ T cells, Naive CD8+ T cells, Memory CD4+ T cells, and Naive CD4+ T cells (Figure 4B). Enrichment analysis across disease groups (R_o/e values) revealed a higher proportion of naive T cells in Healthy individuals, while effector and memory T cells were significantly enriched in diseased groups, with Memory CD8+ T cells showing the strongest enrichment in the Hospitalized group (Figure 4D).

Functional scoring of gene sets further elucidated subtype-specific characteristics. Effector CD8+ T cells displayed the highest cytotoxic scores, reflecting their critical killing functions during disease progression (Figures 5A, E). Memory T cells showed elevated exhaustion and inflammatory scores, suggesting functional impairment during prolonged immune responses (Figures 5B, F). Moreover, regulatory effector scores highlighted the role of Memory CD4+ T cells in immune modulation (Figures 5C, D, G, H). Collectively, these findings emphasize significant functional heterogeneity among T cell subtypes, particularly the roles of effector and memory T cells in immune responses and pathophysiology.

Interaction strength analysis revealed that Memory CD8+ T cells exhibited the highest interaction intensity with other subtypes, such as Effector CD8+ T cells and Memory CD4+ T cells (Figures 6A, B). Further exploration of MHC-I signaling pathways underscored the critical involvement of Memory CD8+ T cells (Figure 6C). Quantitative analysis of interaction strength (Figures 6D, E) highlighted their significant contribution to all pathways, particularly MHC-I signaling. At the transcriptional level, Memory CD8+ T cells demonstrated prominent expression of MHC-I-related genes (e.g., HLA-A, HLA-B, HLA-C) (Figure 6F), reinforcing their pivotal role in antigen presentation and immune response regulation.

KEGG and GOBP analyses revealed the crucial functional roles of Memory CD8+ T cells under pathological conditions. In KEGG pathways, these cells were significantly enriched in immune-related processes such as natural killer cell-mediated cytotoxicity, antigen processing and presentation, chemokine signaling, and COVID-19-related pathways (Figure 7A). Additionally, metabolic pathways like oxidative phosphorylation and ribosome-related pathways were prominently enriched, highlighting their potential role in energy metabolism under active immune states. GOBP analysis further emphasized their involvement in immune regulation, including leukocyte-mediated cytotoxicity, T cell activation, and immune response regulation (Figure 7B). These enriched pathways illustrate the essential functions of Memory CD8+ T cells in antigen recognition, immune responses, and inflammation control.

By comparing Exposed, Infected, and Hospitalized groups against the Healthy group, volcano plots depicted significant differential gene expression (Figures 8A–C). The Hospitalized group exhibited the highest number of upregulated genes, indicating that severe pathological states are associated with substantial transcriptional changes. A cross-group comparison identified seven key genes (RPS26, RPS29, RPL36, RPL39, CD3E, RPS28, RPS21) (Figure 8E). Their expression patterns across groups (Figure 8D) revealed substantial upregulation in the Hospitalized group, suggesting their involvement in immune responses and protein synthesis during disease progression.

Among nine machine learning models, the XGBoost algorithm demonstrated superior performance in both training and testing datasets, achieving the highest average AUC and stable results (Figures 9A, B). Decision Curve Analysis (DCA) validated the model's clinical utility, with standardized net benefits significantly surpassing baseline models ("All" and "None") across various risk thresholds (Figure 9C). SHAP-based analysis identified key predictive genes (e.g., RPS26, RPS29, RPL36), highlighting their significant contributions to disease diagnostics (Figure 9D). Additionally, force plots illustrated how these genes cumulatively influenced the prediction outcomes of specific positive samples (Figure 9E).