Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the pathogen responsible for the global pandemic of coronavirus disease-2019 (COVID-19), has emerged as a significant threat worldwide.^1,2 It was first reported towards the end of 2019,³ and since then, it has rapidly spread across more than 150 countries on every continent except Antarctica. Although COVID-19 is relatively less contagious than its predecessor, SARS-CoV-1, it is associated with a higher mortality rate compared to the seasonal flu, estimated to be around 3.4%. The average incubation period of the disease is approximately 5.2 days, and in severe cases, individuals succumb to the illness within a median duration of 14 days after the onset of symptoms.⁴ The leading causes of death in patients with COVID-19 are acute respiratory distress syndrome, hypoxic respiratory failure, and a cytokine storm triggered by a prolonged infection.⁵

While COVID-19 is generally characterized by asymptomatic or mild to moderate symptoms such as fever and dry cough, it can progress into acute pneumonia, particularly in susceptible populations such as the elderly with metabolic, cardiovascular, and pulmonary conditions.^6,7 In terms of diagnosis, RT-PCR testing is the gold standard for early confirmation of COVID-19 infection, with rapid antigen tests developed for large-scale screening. Additionally, imaging techniques such as chest CT scans are used to assess the severity of lung involvement. Treatment primarily involves supportive therapies, including oxygen therapy and mechanical ventilation.⁸ Current treatment guidelines recommend the administration of anticoagulants, anti-inflammatory drugs, antiviral medications, broad-spectrum antibiotics, and antifungal agents for patients suspected or confirmed to have COVID-19.⁹ However, despite implementing similar clinical interventions, patients may exhibit significant differences in recovery periods, clinical outcomes, or mortality rates.¹⁰ The molecular mechanisms underlying these different outcomes remain largely unknown.¹¹

Upon entering host cells, SARS-CoV-2 exploits the angiotensin-converting enzyme 2 (ACE2) receptor, binding to it through the receptor-binding domain of the spike (S) protein.¹² This interaction triggers fusion between the viral envelope and the host cell membrane via a Toll-like receptor 7 (TLR-7)-dependent pathway, allowing the release of the viral genome into the cytoplasm.¹³ The virus uses its own RNA-dependent RNA polymerase to replicate its genome.¹⁴ The replication-transcription complexes are formed within double-membrane vesicles, mediated by two large polyproteins (such as pp1a and pp1b), which encode non-structural proteins.¹⁵ This process generates multiple subgenomic RNAs that encode structural proteins and accessory proteins.¹⁶ Viral assembly and budding occur within smooth-walled vesicles of the endoplasmic reticulum-golgi intermediate compartment. Finally, the vesicles containing viral particles fuse with the cellular membrane, leading to the release of the virus through exocytosis.¹⁷

Numerous studies have established a strong correlation between innate immune response, adaptive immune development, and the severity of COVID-19.¹⁸ Excessive inflammatory responses have been observed in various samples obtained from patients with COVID-19, including blood, nasopharynx, and bronchoalveolar lavage fluid. Elevated levels of cytokines such as IL6, tumor necrosis factor-alpha (TNF-α), and monocyte chemoattractant protein-1 have been detected, which can contribute to the development of SARS, coagulation disorders, and multi-organ failure.^4,19 IL6 is closely associated with respiratory failure and acute kidney injury, with its levels positively correlated with disease severity.²⁰ IL8 is linked to 30-day mortality.²¹ IL-10, as an anti-inflammatory cytokine, plays a pivotal role in modulating immune responses. While IL-10 typically protects tissues by dampening excessive immune reactions, its abnormal elevation in COVID-19 patients may contribute to disease complexity and progression. Specifically, elevated IL-10 levels might impair the body’s ability to clear the virus by suppressing the function of critical immune effector cells, while also potentially sustaining a chronic inflammatory state.

Although several studies have established prognostic models for patients with COVID-19 based on transcriptomic levels,^22–24 the overall effectiveness of these models has been found to be unsatisfactory. Although transcriptome sequencing data analysis can capture subtle biological changes that may be difficult to detect through clinical features, this type of data presents challenges due to its high dimensionality and complexity. It is particularly susceptible to various factors, including batch effects, which complicate data standardization. Additionally, obtaining gene expression data is generally more expensive and requires more time for experimentation and analysis compared to collecting clinical features. In contrast, clinical rapid testing is typically easier to perform and often exhibits high accuracy. As a result, models based on clinical features are usually more straightforward to understand and interpret, making them more accessible for clinicians in practical applications.

In this study, we conducted a comprehensive analysis of data from 189 patients with COVID-19 hospitalized in real-world settings, examining 149 clinical features to identify those most closely associated with clinical prognosis. Using these features, we developed a robust prognostic model with strong predictive capabilities. Our analysis identified IL-10 as a key factor influencing the prognosis of COVID-19 patients. Importantly, we validated these findings across multiple bulk transcriptomic and single-cell sequencing datasets, underscoring the robustness and reliability of our results. The identification of IL-10 as a significant prognostic biomarker has important clinical implications. Elevated IL-10 levels could serve as an early indicator of high-risk patients, allowing for timely intervention and more personalized treatment strategies. This biomarker may also guide clinicians in adjusting therapeutic approaches, such as intensifying immunomodulatory treatments for patients identified as high-risk. By providing a novel insight into the immune pathogenesis of COVID-19, our findings could inform clinical decision-making and contribute to improved patient outcomes.

Materials and Methods

Collection and Cleaning of Clinical Features

We retrieved data from the electronic medical record database of the Second Hospital of Hebei Medical University, consisting of 189 hospitalized patients admitted between November 1, 2022, and March 31, 2023. All patients were diagnosed with novel coronavirus pneumonia according to the “Diagnosis and Treatment Plan for Novel Coronavirus Pneumonia (Trial Tenth Edition).” The collected patient data included gender, age, smoking history, body mass index, symptoms at admission, vital signs, underlying diseases, blood gas analysis, laboratory examinations, presence of pathogenic bacteria, imaging examinations, electrocardiograms, presence of lower limb deep venous thrombosis, oxygen administration method, treatment, length of hospital stay, and final outcome. Regarding cytokine protein expression levels, we used the flow fluorescence technique to detect them. For the imaging examinations, we followed established protocols and divided each lung into three regions: the upper region (above the clavicle), the middle region (from the clavicle to the lower pulmonary vein), and the lower region (below the lower pulmonary vein). Each lung region was assigned a quantitative score based on the extent of involvement: 0 indicated no involvement; 1 indicated involvement of <25% of the area; 2 indicated involvement of 25–50% of the area; 3 indicated involvement of 50–75% of the area; and 4 indicated involvement of >75% of the area. Each patient had six lung regions, resulting in a maximum score of 24 (4 x 6 regions).²⁵ We performed data cleaning procedures to ensure data quality. First, we removed clinical features with >30% missing values, resulting in a final set of 149 valid indicators. For handling missing values, we used the mice package in R to perform imputation. Multivariate linear regression-based imputation was applied for numerical indicators, while multiple imputations by chained equations were used for categorical indicators.

Prognostic Model Building Based on Machine Learning Integration Framework

In this study, we aimed to develop a prognostic model using a machine learning integration framework. Specifically, we incorporated seven different machine learning algorithms, namely random survival forest (RSF), elastic network (Enet), least absolute shrinkage and selection operator (Lasso), Ridge, stepwise Cox, CoxBoost, and partial least squares regression for Cox (plsRcox). To evaluate their performance, we employed a ten-fold cross-validation strategy. To ensure model consistency and avoid over-reliance on a single feature, we required at least two features for each model, as models with only one feature may perform well on the training cohort but have limited generalization ability to new, unseen data. Through our integration process, we generated a total of 77 different combinations of machine learning algorithms. The performance of each model was assessed using Harrell’s concordance index (C-index), a widely used measure of predictive accuracy for survival models. The model with the highest average C-index value was identified as having the best signature and selected for further analysis.

Gathering and Processing of Bulk Transcriptomic Data

For this study, we utilized bulk transcriptomic data obtained from the Gene Expression Omnibus (GEO) database, accessible at https://www.ncbi.nlm.nih.gov/geo/. Specifically, we analyzed two login identification numbers: GSE217948 and GSE179850. The GSE217948 dataset consisted of transcriptomic data from peripheral blood samples collected from 71 healthy individuals and 203 patients with COVID-19. Similarly, the GSE179850 dataset included transcriptomic data from 16 healthy individuals and 31 patients with COVID-19. Each dataset was analyzed individually to extract meaningful insights. To ensure the accuracy and reliability of our analysis, we employed specific methods in an R environment. Firstly, we used the “avereps” function available in the “limma” package to calculate average gene expression levels for replicate genes. Furthermore, we employed the “normalizeBetweenArrays” function to correct for any potential variability that may exist between the datasets.

Enrichment Analysis

The enrichment analysis of differential genes was conducted using the “ClusterProfiler” and “org.Hs.eg.db” packages. The EnrichGO and EnrichKEGG functions were used for enrichment analysis. Pathways with a P <0.05 were considered to be significantly enriched. Visualization of the results was accomplished using the “ggplot2” package.

Collection and Processing of Single-Cell Transcriptome Data

The single-cell transcriptome data used in this study were obtained from the GEO dataset, specifically the GSE217948 and GSE158055 datasets. For the GSE217948 dataset, quality control measures were performed within the R environment (version 4.1.2) using standard single-cell processing procedures. To process the count matrix, we used the Read10X function from the Seurat package (version 4.0.4) and converted it to the dgCMatrix format. To integrate all individual objects into an aggregate object, we employed the merge function. Additionally, the RenameCells function ensured unique cell labels. To filter out low-quality cells, we applied specific criteria. Genes expressed in <3 cells were removed, and cells expressing <200 genes were deleted. To equalize the total gene expression in each cell, we performed global-scaling normalization using the “LogNormalize” method with a scale factor of 10,000. For downstream analysis, we focused on the top 2000 variably expressed genes, which were identified using the FindVariableFeatures function. Unwanted sources of variation, such as unique molecular identifiers and percent mitochondrial content, were regressed using the ScaleData function with the “vars.to.regress” option. Principal component (PC) analysis was utilized to reduce the dimensionality of the dataset by incorporating highly variable features. The first 30 PCs were selected for further analysis. The Harmony method²⁶ was used to address the batch effects between samples. To visualize the cells in a lower-dimensional space, we applied the uniform manifold approximation and projection method, which effectively reduces the dimensionality of the data while preserving its local structure. Clustering analysis was performed based on edge weights between cells. We constructed a shared nearest-neighbor graph using the Louvain algorithm, implemented in the FindNeighbors and FindClusters functions. The resolution parameter in the FindClusters function was systematically adjusted between 0.1 and 1 to find the optimal clustering resolutions. We used the clustree function to evaluate the clustering trees at different resolutions and ultimately selected a resolution of 0.5, which provided clear and meaningful clustering results. To identify differentially expressed markers associated with the resulting clusters, we used the FindAllMarkers function. This function utilizes a default nonparametric Wilcoxon rank sum test with Bonferroni correction. To annotate all cells, we used cell surface markers and known genes from relevant literature, as well as the CellMarker database²⁷ (available at http://xteam.xbio.top/CellMarker/). For the GSE158055 dataset, we performed a rapid analysis using the online web version available at http://covid19.cancer-pku.cn/#/differential-expression.

Cellular Communication Network

For a comprehensive investigation of cell-cell interactions, we employed the “CellChat” (v1.0.0) R package.²⁸ This package provides us with a curated database of ligands, receptors, cofactors, and their respective interactions, which is publicly available at http://www.cellchat.org/. To analyze cell-cell communication, we projected the expressed ligand and receptor genes of each cell onto a reference communication network. By utilizing gene expression data, we inferred the probability of communication in each pathway. The CellChat package proved to be a versatile and user-friendly tool, enabling us to easily infer, analyze, and visualize cell-cell communication within our single-cell RNA sequencing data. For visualization purposes, we utilized the default parameters of the netVisual_bubble function, which allowed us to visually represent the communication pathways with bubble plots.

Statistical Analysis

To perform our statistical analyses, we used R (version 4.1.2). We considered a P-value of <0.05 to be statistically significant.

Results

Prognostic Model Building Based on Machine Learning Integration Framework

A total of 189 COVID-19 patients were included in the study. Patient clinical and demographic characteristics are shown in Table 1. Our workflow, shown in Figure 1, began with conducting univariable Cox regression analysis on all 149 features to identify prognostic indicators. Among these features, 43 were found to be prognostically relevant (P <0.05), comprising 13 protective features (hazard ratio [HR] <1) and 30 hazard features (HR >1) (Table S1). To construct an optimal prognostic model, we integrated these 43 features into our machine-learning process. For validation purposes, we randomly divided the 189 case samples into two cohorts, namely the training and validation cohorts, at a ratio of 1:1. The comparative 43 features of patients in training cohorts, validation cohorts and all patients are described in Table S2. The model training was exclusively conducted in the training cohort, and we subsequently evaluated the accuracy and validity in both the validation and entire cohorts. Using ten-fold cross-validation, we used a combination of seven machine learning algorithms, namely RSF, Enet, Lasso, Ridge, stepwise Cox, CoxBoost, and plsRcox. One algorithm was used for feature filtering, while another was utilized for constructing a prognostic model. Through this approach, we successfully generated 77 different model combinations within the training cohort. For each model, we calculated the C-index across all validation datasets, selecting the model with the highest average C-index as the optimal choice. Among all the model combinations, the combination of Lasso and RSF demonstrated superior prognostic prediction performance (Figure 2A). Consequently, we named this model the COVID-19-related prognostic model (CRM). The C-index values of the CRM were 0.918, 0.911, and 0.914 in the training cohort, validation cohort, and entire cohort, respectively, indicating the robustness and reliability of the CRM. Using the predict function within the CRM, we calculated a risk score for each patient. The Lasso regression identified four features for inclusion in the CRM (Figure 2B and C). Subsequently, RSF assessed the importance of these four features, ranking them in the following order of importance: IL-10, alpha-hydroxybutyric dehydrogenase, pro-brain natriuretic peptide, and unconsciousness (Figure 2D and E).

Table 1 Clinical Characteristics of Study Participants

Figure 1 Overview of all analytical processes conducted in this study.

Figure 2 Establishment of the machine learning-integrated prognostic model. (A) C-index values of 77 machine learning combinations represented by a color spectrum from green to yellow, indicating increasing C-index values. Each row represents a signature, and each column represents a cohort. (B and C) Lasso regression used for signature selection analysis. (D) RSF analysis to select minimum number of trees. (E) Assessment of signature importance using the RSF algorithm. Abbreviations: RSF, random survival forest; Lasso, least absolute shrinkage and selection operator; C-index, concordance index.

Accuracy of CRM in Predicting the Prognosis of Patients with COVID-19

In order to stratify the patients into high-risk and low-risk groups, we used the “srvminer” package to determine the optimal threshold value based on individual risk scores in the training cohort. Our analysis revealed that the high-risk patients had significantly reduced overall survival (OS) compared to the low-risk patients (Figure 3A). This finding was consistent and replicated in both the validation cohort (Figure 3B) and the entire cohort (Figure 3C), providing strong evidence for the association between high-risk scores and poorer prognostic outcomes. To further assess the predictive performance of the risk scores, We performed a receiver operating characteristic (ROC) analysis, and the results demonstrated the excellent discriminative ability of the risk scores at different time points. In the training cohort, the area under the curve (AUC) values for risk scores at 10, 25, and 40 days were 0.932, 0.969, and 0.949, respectively (Figure 3D). Similarly, the AUC values for risk scores at the same time points were 0.953, 0.960, and 0.975 in the validation cohort (Figure 3E) and 0.932, 0.960, and 0.960 in the entire cohort (Figure 3F), respectively. These results demonstrate the robust effectiveness and accuracy of the CRM in predicting patient outcomes.

Figure 3 Evaluation of the effectiveness and accuracy of the CRM. (A) Survival analysis comparing high- and low-risk patients in the training cohort. (B) Survival analysis comparing high- and low-risk patients in the validation cohort. (C) Survival analysis comparing high- and low-risk patients in the entire cohort. (D) ROC curve analysis of risk scores for the training cohort. (E) ROC curve analysis of risk scores for the validation cohort. (F) ROC curve analysis of risk scores for the entire cohort. Abbreviations: CRM, coronavirus disease 2019-related prognostic model; ROC, receiver operating characteristic.

Stable Diagnostic and Prognostic Value of IL-10

The CRM identified the expression of IL-10 as an important prognostic factor for patients. IL-10, a cytokine regulator, plays a role in modulating immune responses by exerting anti-inflammatory and immunosuppressive effects. In the case of patients with COVID-19, high levels of IL-10 can hinder the host immune response, preventing viral clearance and worsening the condition of the patient. Based on these findings, we conducted a comprehensive analysis of IL-10. First, our analysis revealed that the expression of IL-10 was consistently higher in the high-risk group compared to the low-risk group within training (Figure 4A) and validation cohorts (Figure 4B), as well as across the entire cohort (Figure 4C). Moreover, patients with elevated expression of IL-10 in these cohorts exhibited significantly reduced OS compared to patients those with lower expression of IL-10 (Figures 4D–F). To further assess the diagnostic potential of IL-10 in distinguishing healthy individuals from patients with COVID-19, we examined the GSE217948 and GSE179850 transcriptomic datasets. GSE217948 consisted of transcriptomic data from 71 healthy samples and 203 patients with COVID-19, while GSE179850 included transcriptomic data from 16 healthy individuals and 31 patients with COVID-19. Our analysis of both datasets revealed a significant elevation in the expression of IL-10 in patients with COVID-19 compared to healthy individuals (Figure 4G). Furthermore, the ROC curve showed strong diagnostic performance for IL-10, with an AUC of 0.906 in the GSE217948 dataset and 0.915 in the GSE179850 dataset (Figure 4H). In summary, our findings highlight the critical role of elevated levels of IL-10 in patients with COVID-19, indicating its association with a poor prognosis.

Figure 4 Analysis of the diagnostic and prognostic value of IL-10. (A–C) Expression of IL-10 in high- and low-risk patients in the training, validation, and entire cohorts. (D–F) Survival analysis comparing groups expressing high and low levels of IL-10 in the training, validation, and entire cohorts. (G) Box plots depicting differential expression of IL-10 between healthy individuals and patients with COVID-19 in the GSE217948 and GSE179850 datasets. (H) ROC curves illustrating the differentiation of IL-10 between healthy individuals and patients with COVID-19 in the GSE217948 and GSE179850 datasets. Abbreviations: IL, interleukin; COVID-19, coronavirus disease-2019; ROC, receiver operating characteristic.

IL-10-Related Functional Analysis

To further explore the functional pathways potentially influenced by IL-10, we conducted a correlation analysis between IL-10 and the expression levels of all other genes in patients with COVID-19 using the GSE217948 dataset. Our analysis revealed a positive correlation between the expression of IL-10 and several genes associated with the inflammatory response and immune regulation (Figure 5A). For instance, MS4A4A, which is involved in signal transduction and cellular vitality, could potentially influence immune cell activity. IL-19 and SOCS3, both anti-inflammatory cytokines, may work in conjunction with IL-10 to modulate immune responses. TLR5, responsible for pathogen recognition, exerts immunomodulatory effects. COX6B2, a protein related to oxidative phosphorylation, is involved in cellular energy metabolism. Additionally, CD163, a receptor expressed on macrophages and monocytes, suggests that IL-10 may impact these cells, potentially promoting the formation of M2 macrophages. To gain deeper insights into the functional impact of IL-10, we performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the top 100 co-expressed genes with IL-10. GO enrichment analysis (Figure 5B) revealed enrichment in various metabolic and developmental processes. Metabolic pathways such as the ribonucleotide metabolic process, the ribose phosphate metabolic process, and the long-chain fatty-acyl-CoA metabolic process suggest that IL-10 may regulate immune responses and antiviral capabilities through the regulation of metabolic pathways. Moreover, pathways related to the development of the placenta and the morphogenesis of the embryonic placenta shed light on the potential role of IL-10 in governing cell proliferation and differentiation processes. In terms of cellular components, the significant enrichment of pathways like the mitochondrial outer membrane, peroxisome, trans-golgi network, and nucleolar exosome further supports the notion that IL-10 may be involved in cellular defense against viruses, as well as cell proliferation and differentiation processes, by regulating and repairing cellular organelles. In the molecular function category, enriched activities such as carbohydrate kinase activity, fatty acid ligase activity, acid-thiol ligase activity, cysteine-type endopeptidase inhibitor activity, and calmodulin-dependent protein kinase activity suggest that IL-10 could participate in metabolic pathways and cell proliferation processes by modulating enzyme activity. KEGG enrichment analysis (Figure 5C) highlights pathways such as thermogenesis, peroxisome proliferator-activated receptor (PPAR) signaling, fatty acid metabolism, and inflammatory bowel disease, all of which are relevant to metabolic disorders. This finding aligns with the results of the GO analysis. Additionally, enriched pathways such as Salmonella infection and ferroptosis provide insights into the potential role of IL-10 in modulating immune responses. Collectively, our findings suggest that IL-10 may influence the prognosis of patients with COVID-19 by affecting pathways related to metabolism, proliferation and differentiation of cells, as well as the regulation and repair of cellular organelles and immune responses. These results establish a foundation for future studies exploring the immunopathological mechanisms of COVID-19.

Figure 5 Functional analysis related to IL-10. (A) Pearson correlation analysis of genes co-expressed with IL-10. (B) GO enrichment analysis. (C) KEGG enrichment analysis. Abbreviations: IL, interleukin; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Single-Cell Sequencing Analysis: Monocytes as the Primary Source of the Elevated IL-10 in Patients with COVID-19

Although our clinical data and transcriptomic analysis have demonstrated the robust diagnostic and prognostic value of IL-10, the specific source of IL-10 in patients with COVID-19 has remained elusive. Bulk transcriptomic analysis tends to mask crucial differences among cellular subpopulations. However, the emergence of single-cell sequencing technology presents an opportunity to uncover the intrinsic features of patients with COVID-19 at the single-cell level. In this study, we obtained single-cell transcriptome data from peripheral blood samples of 11 patients with COVID-19 and four healthy individuals, sourced from the GSE149689 dataset. Following quality control measures, we included a total of 70,559 cells for subsequent analyses. Using cell surface markers, we annotated all cells (Figure 6A) and identified nine distinct cell types (Figure 6B). These cell types encompassed B cells (CD79A and MS4A1), dendritic cells (DC) (CLEC10A, CST3, HLA-DPA1, and HLA-DQA1), plasmacytoid DC (pDC) (ALOX5AP, APP, CCDC50, IRF7, IRF8, and LILRA4), erythroid cells (HBA1, HBA2, HBB, and AHSP), innate lymphoid cells (ILC) (GNLY, GZMA, GZMM, and NKG7), megakaryocytes (PPBP, GNG11, MYL9, and PF4), monocytes (S100A9, S100A11, AIF1, and S100A8), plasma cells (IGHA1, IGHG1, IGHG3, IGHG4, and IGKC), and T cells (CD3D, CD3E, and IL7R). Consistently, our analysis revealed significantly elevated expression levels of IL-10 in patients with COVID-19 compared to healthy samples (Figure 6C), aligning with the findings from bulk transcriptomic data analysis. Notably, we observed that monocytes exhibited the highest expression of IL-10 compared to other cell types (Figure 6D), suggesting that monocytes serve as the primary source of elevation of IL-10 in patients with COVID-10. To validate our findings externally, we further analyzed the GSE158055 dataset, which included 284 single-cell transcriptomic profiles collected from 196 individuals, including 22 patients with mild or moderate symptoms, 54 hospitalized patients with severe symptoms, and 95 recovered convalescent individuals, as well as 25 healthy individuals. The dataset encompassed samples from various tissue types, such as peripheral blood mononuclear cells, bronchoalveolar lavage fluid, and pleural effusion sputum. Similarly, our analysis demonstrated that the expression of IL-10 was highest at the monocyte level (Figure 6E). Additionally, expression levels of IL-10 were elevated in nucleic acid test-positive specimens compared to negative ones (Figure 6F). Furthermore, the healthy individuals exhibited the lowest expression level of IL-10, while the progression group showed the highest expression, with levels of IL-10 partially recovering in the convalescent group (Figure 6G). The levels of IL-10 gradually increased from the normal level to mild, moderate, severe, or critical levels (Figure 6H). We also discovered a gradual increase in the expression levels of IL-10 with age (Figure 6I), which may contribute to the relatively severe condition of COVID-19 in elderly individuals. In conclusion, our single-cell analysis findings indicate that monocytes are the primary source of increased IL-10 in patients with COVID-19.

Figure 6 Differential expression analysis of IL-10 at the single-cell level. (A) A heatmap showing marker gene expression used for cell annotation in each cell in the GSE149689 dataset. (B) Dimensional landscape showing cell annotation in the GSE149689 dataset. (C) Differential expression of IL-10 in patients with COVID-19 and healthy individuals within the GSE149689 dataset. (D) Differential expression of IL-10 in different cell types in the GSE149689 dataset. (E) Differential expression of IL-10 in different cell types in the GSE158055 dataset. (F–I) Differential expression of IL-10 in different clinical feature subgroups within the GSE158055 dataset. Abbreviation: COVID-19, coronavirus disease-2019.

Monocytes Expressing High Levels of IL-10 are Involved in Immune Cell Communication in Patients with COVID-19

Cellular communication plays a crucial role in understanding the complexity of the immune microenvironment in patients with COVID-19, as changes in this microenvironment can influence disease progression. However, bulk transcriptomic data fails to capture this aspect adequately. To address this, we divided all monocytes into two groups based on the expression of IL-10: high- and low-expressing monocytes (Figures 7A–C). Subsequently, we examined changes in cellular communication between these two types of cells and other cells (Figures 7D and E). Our analysis revealed several specific changes in signaling between monocytes expressing high levels of IL-10 and other cells compared to those monocytes expressing low levels of IL-10. Notably, the activation of ANXA1-FPR2 signaling between monocytes expressing high levels of IL-10, DCs, ILCs, pDCs, T cells, and other immune cells suggests a critical role in modulating the immune response, as this pathway is associated with the resolution of inflammation. The activation of MIF-(CD74+CXCR4) signaling between B and T cells as ligand cells and monocytes expressing high levels of IL-10 as receptor cells indicates a potential involvement in cellular communication. Similarly, the activation of NAMPT-INSR signaling between DCs as ligand cells and monocytes expressing high levels of IL-10 as receptor cells highlights another interaction. Additionally, the activation of transforming growth factor (TGF) B1-(TGFBR1 + TGFBR2) and TGFB1-(ACVR1B + TGFBR2) signaling between other cells and monocytes expressing high levels of IL-10 is consistent with previous studies suggesting a role for TGF-β signaling in promoting immunosuppression. The activation of tumor necrosis factor superfamily (TNFSF) 12-TNFRSF12A signaling between DCs, ILCs, and monocytes expressing high levels of IL-10 may contribute to the activation and differentiation of T cells. Finally, we observed enhanced autocrine signaling of CD99-CD99 in monocytes expressing high levels of IL-10 compared to those expressing low levels of IL-10, which may impair cellular immune function and hinder humoral immune responses, potentially increasing susceptibility to infection. CD99-CD99 signaling is involved in vital cellular processes such as cell adhesion, cell-mediated toxicity, and intercellular interactions. Consequently, enhanced autocrine signaling of CD99-CD99 in monocytes expressing high levels of IL-10 could be related to the disease deterioration and poor prognosis in patients with COVID-19. In summary, these findings shed new light on the intricate cellular communication among immune cells in patients with COVID-19 and altered immune microenvironments. Our study provides a comprehensive understanding of the specific changes in cellular communication associated with monocytes expressing high levels of IL-10 in patients with COVID-19. These results have the potential to guide the development of novel therapeutic approaches and aid in the identification of biomarkers for predicting disease progression.

Figure 7 Analysis of cell communication. (A) The re-dimensionalization landscape of all monocytes after re-extraction. (B and C) Splitting of single monocytic cells into groups with high or low expression of IL-10 based on the expression of IL-10; (D and E) Communication analysis between the two groups expressing IL-10 using CellChat software, with “high” representing single monocytic cells with high expression of IL-10, and “low” representing single monocytic cells with low expression of IL-10. Abbreviation: IL, interleukin.

Discussion

In this study, our aim was to gain a deeper understanding of the disease status of hospitalized patients with COVID-19 and identify novel biomarkers that could assist in predicting patient prognosis. Initially, we analyzed real-world datasets to examine various clinical features of patients with COVID-19, including imaging, drug treatment, and peripheral blood characteristics. Subsequently, we developed a robust clinical feature model using consensus frameworks in machine learning to predict patient prognosis effectively. Notably, our model identified IL-10 as a crucial clinical indicator that significantly influences patient prognosis, a finding consistently validated in both bulk transcriptome and single-cell sequencing datasets. Moreover, our analysis of single-cell transcriptome data revealed that monocytes in the peripheral blood of the patients were the primary source of IL-10. These monocytes expressing high levels of IL-10 were involved in cell-to-cell communication within the immune microenvironment of the patient. Collectively, our findings highlight the importance of IL-10 as an important biomarker, emphasizing the need to test for levels of IL-10 early in inpatient admission and throughout the entire hospitalization period.

Previous studies have highlighted the significant changes in levels of various cytokines experienced by patients with COVID-19.²⁹ The infection triggers the immune system, resulting in variations in immune cell numbers and proportions, leading to the secretion of various cytokines that promote and regulate immune responses.³⁰ However, the specific cytokines that play a crucial role remain unclear. It has been observed that COVID-19 elevates the levels of cytokines in the body, including IL-10, IL-6, TNF-α, interferon-gamma (IFN-γ), and IL-1β, among others.^31–34 These cytokines have the capacity to regulate the inflammatory response of the body, potentially leading to pulmonary and cardiovascular complications.³⁵ While some studies suggest that IL-10 is not a risk factor for COVID-19, they emphasize the association of higher levels of IL-6 with the severity of the disease and mortality rates.³⁶ Conversely, some studies indicate that the combined action of TNF-α and IFN-γ can trigger extensive apoptosis, pyroptosis, and necroptosis of immune cells in the bodies of patients with COVID-19.³⁷ Furthermore, studies demonstrate that compared to previous outbreaks such as SARS and Middle East respiratory syndrome (MERS), COVID-19 exhibits higher levels of IL-10, indicating that the immune response to COVID-19 exhibits higher levels of IL-10, suggesting potentially distinctive immunomodulatory characteristics.⁸ Overall, COVID-19 infection is characterized by significant changes in the levels of cytokines, which have implications for the severity and prognosis of the disease.¹ Understanding these changes is crucial for the development of more effective treatment strategies and the prediction of the risk of infection following COVID-19.

In our analysis, we investigated the genes associated with the expression of IL-10 and identified the top six co-expressed genes: MS4A4A, IL-19, SOCS3, TLR5, COX6B2, and CD163. Notably, MS4A4A and IL-19 exhibited high co-expression coefficients with IL-10, reaching approximately 0.7. Previous studies have also confirmed the relevance of these genes to COVID-19. MS4A4A is a transmembrane protein widely expressed in immune cells, particularly monocytes and B cells.³⁸ It plays a role in immune regulation during viral infections. A large-scale gene association study targeting COVID-19 revealed that MS4A4A, located on chromosome 11q12.2, is associated with severe COVID-19 infection.³⁹ Expression of MS4A4A is affected by COVID-19 in monocytes and B-cells, with further downregulation observed in patients with severe infection.⁴⁰ Additionally, MS4A4A has been found to regulate cellular apoptosis and activation of T cells, indicating its significance in the immune response to COVID-19.⁴¹ Decreased expression of MS4A4A in patients with COVID-19 may lead to a compromised antiviral immune response, potentially exacerbating the disease.⁴² Moreover, other studies have suggested that MS4A4A can regulate inflammatory responses, which are major contributors to pathological damage caused by lung inflammation in COVID-19, leading to death.⁴² Therefore, MS4A4A may play a crucial role in suppressing pathological inflammatory responses, and its decreased expression may be associated with the severity of the disease.

Similarly, IL-19 is also regulated during severe COVID-19 infections and plays a vital role in regulating immune and inflammatory responses, proliferation, and differentiation.⁴³ During COVID-19 infection, immune cells produce excessive inflammatory cytokines, leading to an overactive immune response. IL-19 may act as a negative regulator, mitigating this excessive immune response.⁴⁴ For example, IL-19 can inhibit the recognition of pathogens by alveolar macrophages and reduce macrophage activity.⁴⁵ While this may increase the risk of bacterial infections in certain cases, it could help alleviate lung inflammation and atelectasis in COVID-19. Furthermore, studies suggest that IL-19 is involved in immune response and antiviral immunity.⁴⁶ IL-19 can promote viral clearance, facilitate pathogen-mediated cross-interference, and regulate the polarization state of immune cells, leading to reduced inflammatory responses.⁴⁷ Patients with severe COVID-19 infection often experience a cytokine storm, which is characterized by heightened inflammation and tissue damage due to an overactive immune response.³⁷ IL-19 may contribute to the occurrence and development of cytokine storms,⁴⁸ making it an important regulator of immune responses in COVID-19. By inhibiting excessive immune reactions and alleviating inflammatory responses, IL-19 plays a crucial role in mitigating the severity of COVID-19 infections. In summary, IL-19 holds significant potential in the context of severe COVID-19 infections. While our understanding of its precise mechanism is still evolving, studying the role of IL-19 in COVID-19 provides valuable insights for better prevention and treatment of COVID-19.

We investigated the impact of COVID-19 on cellular functions and signaling pathways related to IL-10. Our enrichment analysis revealed that IL-10 may influence various metabolic pathways and immune responses. These findings are consistent with previous studies that have explored the role of IL-10 function in similar contexts.⁴⁹ Specifically, in the case of pulmonary viral infections, IL-10 can inhibit the development of IL-17-producing cells that can cause tissue damage, suppress innate inflammatory responses, and potentially reduce the production of cytokines.⁵⁰ The increased levels of IL-10 observed in patients with COVID-19 may directly result from the progression of lung inflammation and serve as a part of feedback inhibition circuits that counteract inflammatory responses.⁵⁰ Furthermore, elevated levels of IL-10 can serve as blood-based biomarkers for patients with COVID-19 who experience worsening lung damage.⁵¹ However, individuals with higher levels of IL-10 may exhibit adaptive immune deficiencies, which have been associated with long-term, persistent symptoms of viral infection.⁵² This overproduction of IL-10 could compromise the immune response to recurrent bacterial infections, leading to secondary pneumonia and playing a significant role in the progression of COVID-19 in patients.⁵³

In subsequent studies, we categorized monocytes from the single-cell data into two groups based on their expression levels of IL-10. The analysis revealed notable differences in cell communication between cells with high expression of IL-10 and those with low expression, particularly in the TNFSF and TGFβ signaling pathways. This finding is consistent with a previous study demonstrating the involvement of various TNFSF proteins, such as CD40L, FasL, gitr, and OX40, in COVID-19. CD40L, for example, has been associated with severe COVID-19 infections and is involved in B cell and multiple antibody responses.⁵⁴ FasL, expressed by CD8 T cells, may contribute to COVID-19 by inducing cell apoptosis.⁵⁵ Additionally, GITR and OX40, expressed in T cells and plasma cells, may also have implications for COVID-19.⁵⁶ Consequently, further extensive and in-depth research is needed to comprehensively understand the role of TNFSF in COVID-19, offering valuable insights for the treatment and prevention of COVID-19. TGFB1, another significant factor, participates in the regulation of the immune response in patients with COVID-19. Studies have shown that in COVID-19, virus-induced immune cells release excessive inflammatory cytokines, leading to an extreme activation of the immune system, characterized by the activation and infiltration of inflammatory cells, which results in an excessive inflammatory response and lung damage.⁵⁷ TGFB1 can regulate the activation of inflammatory cells and suppress the inflammatory response, thus playing a protective role.⁵⁸ A study observed a decrease in the levels of TGFB1 in the serum of patients with severe COVID-19, possibly related to the excessive activation of immune cells caused by the viral infection.⁵⁹ Furthermore, TGFB1 also plays an important role in the onset and progression of COVID-19-related lung injury.

Studies have shown that the COVID-19 virus impacts the upregulation of ACE2 and TMPRSS2, receptors associated with virus infection. This leads to the release of inflammatory mediators, which further activate the immune response, ultimately resulting in damage to lung parenchymal tissue.⁵⁷ Within this process, the regulation of the expression of TGFB1 influences the functionality of immune cells, including inflammatory cells, T cells, B cells, macrophages, and fibroblasts, thereby affecting their involvement in the onset and progression of lung injury.⁶⁰ Overall, TGFB1 plays a crucial regulatory role in COVID-19, particularly in severe cases, by regulating various aspects such as immune and inflammatory cells, lung injury, and prognosis. However, there is currently limited research on the specific mechanisms and therapeutic effects of TGFB1 in COVID-19. Therefore, further investigation is necessary to enhance our understanding of the association between this protein and severe COVID-19 infection, enabling better guidance for the treatment and prevention of COVID-19 in patients.

In conclusion, we established a prognostic model for patients with COVID-19 using a machine-learning consensus framework, incorporating real-world clinical features. Notably, our study identified IL-10 as a critical prognostic factor. These findings offer valuable insights into the intricate pathogenesis and progression of COVID-19. They have the potential to aid in the identification of high-risk patients with COVID-19 at fatal risk and may enable the implementation of personalized strategies for risk assessment and therapeutic intervention, ultimately improving patient outcomes. However, our study still has certain limitations. Our prognostic model needs to be validated in larger patient cohorts.