key: cord-0845222-x24g2q6v
authors: Li, Xu; Garg, Manik; Jia, Tingting; Liao, Qijun; Yuan, Lifang; Li, Mao; Wu, Zhengyu; Wu, Weihua; Bi, Yalan; George, Nancy; Papatheodorou, Irene; Brazma, Alvis; Luo, Huanle; Fang, Shisong; Miao, Zhichao; Shu, Yuelong
title: Single-Cell Analysis Reveals the Immune Characteristics of Myeloid Cells and Memory T Cells in Recovered COVID-19 Patients With Different Severities
date: 2022-01-03
journal: Front Immunol
DOI: 10.3389/fimmu.2021.781432
sha: 98b871d571bb196283ea39198fb0f76913d3c411
doc_id: 845222
cord_uid: x24g2q6v

Despite many studies on the immune characteristics of Coronavirus disease 2019 (COVID-19) patients in the progression stage, a detailed understanding of pertinent immune cells in recovered patients is lacking. We performed single-cell RNA sequencing on samples from recovered COVID-19 patients and healthy controls. We created a comprehensive immune landscape with more than 260,000 peripheral blood mononuclear cells (PBMCs) from 41 samples by integrating our dataset with previously reported datasets, which included samples collected between 27 and 47 days after symptom onset. According to our large-scale single-cell analysis, recovered patients, who had severe symptoms (severe/critical recovered), still exhibited peripheral immune disorders 1–2 months after symptom onset. Specifically, in these severe/critical recovered patients, human leukocyte antigen (HLA) class II and antigen processing pathways were downregulated in both CD14 monocytes and dendritic cells compared to healthy controls, while the proportion of CD14 monocytes increased. These may lead to the downregulation of T-cell differentiation pathways in memory T cells. However, in the mild/moderate recovered patients, the proportion of plasmacytoid dendritic cells increased compared to healthy controls, accompanied by the upregulation of HLA-DRA and HLA-DRB1 in both CD14 monocytes and dendritic cells. In addition, T-cell differentiation regulation and memory T cell–related genes FOS, JUN, CD69, CXCR4, and CD83 were upregulated in the mild/moderate recovered patients. Further, the immunoglobulin heavy chain V3-21 (IGHV3-21) gene segment was preferred in B-cell immune repertoires in severe/critical recovered patients. Collectively, we provide a large-scale single-cell atlas of the peripheral immune response in recovered COVID-19 patients.

Coronavirus disease 2019 , caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has resulted in more than 261 million confirmed cases and more than 5.2 million deaths according to the statistics of the World Health Organization (WHO) as of November 30, 2021 . The impact of this disease has led to extensive research work to quickly understand, control, and treat COVID-19. Most previous single-cell studies (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) focused on the COVID-19 progression stage and have provided important immune cellular and molecular characteristics. The large-scale integrated analysis of single-cell data by Ren et al. (9) included single-cell sequencing data of 140 different types of samples from 104 COVID-19 convalescent patients, showed the immune cell proportions of peripheral blood mononuclear cells (PBMCs) and T-cell receptor (TCR) clone diversity and other characteristics in convalescent patients. Zhang et al. (8) profiled adaptive immune cells of PBMCs from recovered COVID-19 patients with varying disease severity using single-cell RNA sequencing (scRNA-seq), single-cell TCR sequencing (scTCR-seq), and single-cell BCR sequencing (scBCR-seq). However, these studies cannot explain the phenomenon recently discovered by Amanat et al. (14) that convalescent patients with high serum anti-spike titers produce a higher proportion of non-neutralizing antibodies. Data from several studies suggested that patients with severe COVID-19 had higher serum anti-spike titers (15) (16) (17) (18) . In addition, according to a large-scale population survey in Denmark, the protection rate of individuals under 65 years of age against SARS-CoV-2 reinfection is higher than 80%, while patients 65 years of age and older have only 47.1% protection (19) . There is currently no clear explanation for this. Hence, there is an urgent need for a deeper and more comprehensive understanding of the immune characteristics of COVID-19 patients during the recovery stage.

Here we analyzed the scRNA-seq data of 41 individuals of more than 260,000 cells, including 16 mild/moderate recovered patients, 6 severe/critical recovered patients, and 19 healthy controls. CD14 monocytes (CD14 mono), CD4 T cells, and CD8 T cells in severe/critical recovered patients were still in a disordered state 27-47 days after symptom onset, accompanied by a high expression of cytokines and interferon-stimulated genes (ISGs). The percentages of CD4 T cells and CD8 T cells in mild/moderate recovered patients were comparable to healthy controls, but showed significant transcriptome changes. Our data and findings may have important implications for revealing the relationship between the immune response of patients recovering from COVID-19 and the immune protection against SARS-CoV-2 reinfection.

The dataset generated in this study was termed as the "Li dataset". To generate this dataset, human blood samples were collected by Shenzhen Center for Disease Control and Prevention, Shenzhen, China. We collected PBMC samples at 27-47 days after onset of symptoms or tested with SARS-CoV-2 nucleic acid positive. All patients were in recovery stage or had no clinical symptoms at sample collection (Table S1A) . PBMCs were isolated immediately, using lymphocyte separation fluid under the enhanced biosafety level 2 facility. Then, we used a freshly prepared freezing solution [fetal bovine serum (FBS) containing 10% dimethyl sulfoxide (DMSO)] to freeze the PBMCs. In addition, we collected two fresh PBMC samples from two healthy individuals for sequencing as healthy control (Table S1A) .

For the Li dataset, after sample collection, the PBMCs were stored in liquid nitrogen. Cell suspensions were barcoded through the 10x Chromium Single Cell platform using the Chromium Single Cell 5′ Library, Gel Bead and Multiplex, and Chip kits (10x Genomics). Twenty thousand PBMCs were loaded on each 10X Chromium A Chip. Single-cell lysis and RNA firststrand synthesis were performed using the 10X Chromium Single Cell 5′ Library and Gel Bead Kit according to the manufacturer's protocol. RNA and VDJ library preparation were performed according to the manufacturer's instructions, using the Chromium Single Cell 5′ v3 Reagent (10x Genomics) and Chromium Single Cell V(D)J Reagent kits (10x Genomics). Each sequencing library was generated with a unique sample index. All libraries were quantified by Qubit 3.0 (Thermo Fisher), Agilent 2100, and Qsep 100. Sequencing was performed on a Hiseq4000 platform with a paired-end 150 sequencing strategy.

For the Li dataset, we used the Cell Ranger single-cell software suite (version 3.0.0, 10x Genomics) to compare and quantify the single-cell sequencing data against the GRCh38 human reference genome. Firstly, the cells in each sample were screened, and cells expressing at least 200 genes were kept. Next, cells were filtered according to three criteria: (1) cells must have a proportion of mitochondrial gene counts (UMIs from mitochondrial genes/ total UMIs) of less than 15%; (2) cells must have a total number of unique molecular identifiers (UMI) counts per cell (library size) of more than 500; (3) genes must be expressed in more than two cells. Doublets were identified using Scrublet (20) and were removed from the analysis. After quality control filtering, a total number of 86,650 cells were retained for downstream analysis (details are shown in Table S2A, B) .

For the preprocessing of the Ren dataset (9), we downloaded the scRNA-seq expression profile from NCBI Gene Expression Omnibus (GEO) database: GSE158055 (https://www.ncbi.nlm. nih.gov/geo/query/acc.cgi?acc=GSE158055). Only healthy control samples and samples of 27 and 47 days after recovery (hospital discharge) were included in our analysis. As the downloaded data had already been analyzed and annotated, all the cells in the dataset were used in our analysis without quality control. The two datasets were merged together based on the raw counts using the concatenate function in Scanpy (21) . An overlap of 22,094 genes was found in the merged data. To integrate the Li dataset and the Ren dataset, we selected the top 2,000 highly variable genes for the integrated dataset using the "seurat_v3" flavor in the scanpy.pp.highly_variable_genes() function in Scanpy (21) .

Scanpy (21) was used to analyze the data, including normalization, transformation, highly variable gene selection, and dimension reduction. The expression profile was first normalized to counts per ten thousand (CPTT) by scanpy.pp.normalize_total() function and then log-transformed by scanpy.pp.log1p(). Highly variable genes were selected with scanpy.pp.highly_variable_genes() according to the mean expression and dispersion of the genes. Principal Component Analysis (PCA) was performed on the expression profile of the highly variable genes. Harmony (version 0.0.5) (22) was used to integrate data on the latent space from different samples on the PCA space. Dimension reduction was performed with Uniform Manifold Approximation and Projection (UMAP) (23) . To present cell clustering using Louvain method (24) on the Harmony corrected latent space. The cluster stabilities were assessed by self-projection (25) . A machine learning-based method (25) was used to infer cell types according to two annotated reference datasets (3, 6) .

We calculated immune cell proportions for each major cell type and cell subtype. For each sample, cell type proportion was calculated by the number of cells in a certain cell type divided by total number of cells. To identify changes in cell proportions between samples in different disease severity states and sex, we performed T tests and non-parametric tests on the proportions of each major cell type and cell subtype across different groups. We performed Spearman's correlation analysis to assess the association between cell type proportion and patient age. A value of p less than 0.05 is regarded as significant.

To further evaluate the influence of different sample technical factors, patient phenotypes, and their potential interactions to cell type proportions, we performed One-way ANOVA on cell type proportions based on different phenotypes (Figures S2A), including disease severity, sex, and sample type (fresh or frozen). In Figure S2A , we included the sample data from the Ren dataset with 21-90 days sampling time (days after symptom onset) for One-way ANOVA and linear regression analysis. Cell type proportions were used as the outcome in a regression analysis with age and sampling time (days after symptom onset) as predictors, respectively. Following a multiple testing correction, phenotypes were regarded as significantly associated with cell type proportions when q value is less than 0.05. Statistical analyses were performed using R software (v 4.1.1; The R Foundation).

To model the differential expression of genes between cells together with technical effects, we used the hurdle model in MAST (26) , while accounting for the covariates of sex and age. In the model, the effects of covariates are regressed out such that the differential expression represents the effect of the disease condition. Differentially expressed genes with a false detection rate (FDR) lower than 0.01 were used for volcano plots and pathway/gene ontology analysis. Upregulated genes were defined as the ones with a positive log fold change value. GProfiler (27) was used to analyze and visualize the regulated pathways based on the Gene Set Enrichment Analysis (GSEA) database of hallmark genes (28) , while the Kyoto Encyclopedia of Genes and Genomes (KEGG) hallmark gene set was used in the analysis.

The gene module score was calculated as the average expression of a set of genes in a given gene module subtracted with the average expression of a reference set of genes. The latter were randomly sampled from the gene_pool for each binned expression value. For a given set of genes belonging to each module (such as HLA class II, cytokine module), the scores were generated using scanpy.tl.score_genes() function of Scanpy (v1.8.1) (21) with the parameter ctrl_size=100.

The Ren TCR/BCR data (9) were first downloaded from NCBI Gene Expression Omnibus (GEO) database: GSE158055 (https:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE158055). These data have been preprocessed by Ren et al. (9) , and only those T cells with at least one TCR a-chain (TRA) and one TCR b-chain (TRB) were provided in this dataset. Similar preprocessing was done for BCR data where cells with at least one heavy chain (IGH) and one light chain (IGK/IGL) were provided. The clonotype frequencies for each sample were also provided. We preprocessed our TCR/BCR datasets in a similar manner before integrating them with Ren's dataset and calculated the clonotype frequencies per sample. Please note that only those B-cell/T-cell clones with corresponding cell-type annotations obtained from scRNA-seq analysis (see Methods section Single-Cell RNA-Seq Cell Type Annotation) were considered for further analysis. The final numbers are given in Table S2A . The UMAPs were plotted using the scanpy.pl.umap() function in (21) , and the rest of the analysis was done using Rpackage Immunarch (v0.6.6) (29, 30).

The Scanpy AnnData (21) containing all cells was subsetted to T cells, Natural killer (NK) cells, monocytes, dendritic cells (DCs), and plasmacytoid dendritic cells (pDCs) belonging to the 41 samples with a sampling date 27-47 days. This AnnData was further divided into three subsets corresponding to the three severities. The R-package CellChat (v1.1.3) (31) was used to infer cell-cell interaction networks where these annData objects used by Scanpy were converted to separate cellchat objects and then merged for comparison. To identify the upregulated and downregulated signaling pathways using differential expression analysis, two cellchat objects were analyzed at a time. Those differentially expressed signaling ligands with a Bonferroni corrected p-value lower than 0.05 and a log fold change higher than 0.01 were considered upregulated in the second dataset, while those ligands and receptors with a Bonferroni corrected pvalue lower than 0.05 and a log fold change higher than 0.01 in the first dataset were considered downregulated in the second dataset. Only those ligand/receptor genes expressed in at least 10% of cells in the respective datasets were considered for visualization and analysis.

To profile the transcriptional immune landscape of COVID-19 patients in the recovery stage, we collected 10 PBMC samples from recovered COVID-19 patients 27-47 days after symptom onset, and 2 PBMC samples from healthy control. Using singlecell sequencing technologies, we performed single-cell RNA sequencing as well as single-cell immune profiling (both singlecell B-cell receptor sequencing, scBCR-seq, and single-cell T-cell receptor sequencing, scTCR-seq) on these 12 samples (we term the dataset we generated as the "Li dataset"; see Methods). To improve the reliability and reproducibility of the data analysis, we added 12 PBMC samples from COVID-19 patients in the recovery stage from the previously reported dataset (9) along with 17 PBMC samples from healthy control (this added dataset is termed as the "Ren dataset"). These samples were selected to match the sampling time (according to the day after symptom onset) of the Li dataset, i.e., both datasets had a sampling day of 27-47 days after symptom onset (Tables S1A, B). The overall integrated data included 19 healthy control samples (HC), 16 mild/moderate recovered (MR) samples, and 6 severe/critical recovered (SR) samples, which were classified according to WHO criteria (https://www.who.int/publications/i/item/covid-19-therapeutic-trial-synopsis) ( Figure 1A and Tables S1A, B).

In Table S1C , we compare the detailed characteristics of patients and controls (including median sampling day and sample type, patient median age, gender, comorbidities, and outcome). Consistent with other reports (9, 16), we found that the median age of severe/critical patients is greater than that of mild/moderate patients. The median sampling day of the mild/moderate and severe/critical recovery groups are 33.5 and 34.5 days (days after symptom onset), respectively (Table S1C) .

This basically eliminates the effect of detection time on the result of scRNA-seq of patients, and these two groups are comparable in our data.

After quality control, we obtained transcriptomes of 86,650 cells, 3,693 productive BCR clones, and 15,717 productive TCR clones from the Li dataset (Table S2A) . From the Ren dataset we retrieved 178,239 cells, 13,899 productive BCR clones, and 41,676 productive TCR clones (Table S2A) . A Unified Manifold Approximation and Projection (UMAP) based on the Harmony-corrected latent space was generated ( Figure 1B ; see Methods). We identified 28 distinct cell populations using a machine learning-based approach (25) in comparison with two annotated reference datasets (3, 6) ( Figure 1B) . These cell populations were further confirmed with known marker gene expression ( Figures S1A, B) .

We first analyzed the compositional changes of the broad categories of immune cells in PBMC ( Figure 1C) . Notably, in severe/critical COVID-19 patients in the recovery stage, the proportion of CD14 monocytes (based on the expression of the marker genes CST3, LYZ, CD14) in PBMCs were elevated, but CD4 T cells (IL7R, LTB) were decreased ( Figure 1C ), consistent with a previous report (9) . CD8 T cells (which express CD8A, CD8B) decrease in severe/critical recovered patients in comparison with healthy controls.

Multiple subtypes of myeloid cells significantly changed in cell proportions and genes transcription during the progression of COVID-19 dependent on symptom severity (6, 9) . To understand the immune characteristics of these myeloid cell subtypes during the recovery stage of COVID-19, we analyzed the relationship between patient age, sex, symptom severity, and PBMC compositions ( Figure 2 and Figures S2, S3) . The percentages of CD14 Mono and CD14 Mono (NFKBIA) cells were significantly higher in severe/critical recovered patients compared with healthy controls. We analyzed the correlation between all COVID-19 convalescent patient samples and age, and found CD14 Mono increases with age in the convalescent COVID-19 patients (r=0.7294, p=0.0019), but CD14 Mono (NFKBIA) has no significant correlation with age (r=0.3153, p=0.1530) ( Figure 2C ). Sex has no significant effect on the proportions of these two cell subtypes (p>0.05) ( Figure 2B) . Surprisingly, pDCs were significantly elevated in mild/moderate recovered patients but were comparable with healthy controls in severe/critical recovered patients (Figure 2A ). Sex has no significant effect on the percentage of pDCs, but age is negatively correlated with the percentage of pDCs (r=−0.4825, p=0.023) ( Figures 2B, C) .

Next, we defined an inflammatory score and HLA class II score for each cell based on the expression of the reported inflammatory response genes (32) and HLA class II genes (6) ( Table S3) , respectively. We used these two scores to evaluate the inflammation status and antigen presentation ability for each cell (Figures 2D, E) . The expression of HLA class II genes was significantly reduced in severe/critical recovered patients, and this reduction was more significant in monocyte and dendritic cell populations (Figures 2D, E) . Similarly, we calculated the cytokine score and ISG score of each cell according to the expression of the collected cytokine genes and ISG genes ( Table S3 ; see Methods). The expression of genes in cytokine module was significantly increased in severe/critical recovered patients but was comparable to healthy controls in mild/ moderate recovered patients or showed a downward trend ( Figures S2B, C) . This phenomenon was observed in most cell subtypes. The expression of genes in the ISG module was essentially restored to healthy levels in recovered COVID-19 patients, but remained high in a subset of cell subtypes in severe/ critical recovered patients, for example Prolif T, CD8m T (GZMH), and NK(GZMH) (Figures S2B, C) .

Through the KEGG (33) pathway analysis, we further detected the differences in cell function of CD14 Mono, CD14 Mono (NFKBIA), and pDCs in recovered patients with different clinical severity ( Figure 2F and Figure S2D ). In CD14 monocytes, the nucleotide-binding oligomerization domain (NOD)-like receptor signaling pathway was enriched in severe/critical recovered patients, but the antigen processing and presentation and Intestinal immune network for IgA production were downregulated ( Figure 2F ). In pDCs, Lysosome pathway was downregulated in severe/critical recovered patients. Together with the significant decrease in HLA class II score, it suggested that the severe/critical recovered patients' antigen processing and presentation ability was reduced ( Figures 2E, F) . Single-Cell Transcriptional Landscape of T Cells

To further clarify the T-cell immune characteristics of COVID-19 patients during the recovery stage, we analyzed each T-cell subtype. According to the cell proportion analysis, CD4m T, CD4m T (GZMK), CD8m T (IL7R), gd T, and Vg9Vd2 T (34, 35) cells significantly decreased in severe/critical recovered patients, but were comparable with healthy controls in mild/moderate recovered patients ( Figure 3A and Figure S3A ). Unlike these Tcell subtypes, proliferating T (prolif T) cells were increased in severe/critical recovered patients compared to healthy controls ( Figure 3A) , which was consistent with the results reported by Ren et al. (9) . In the correlation analysis with age, prolif T cells were positively correlated with age (r=0.5206, p=0.0005), while gd T cells and Vg9Vd2 T cells were negatively correlated with age (r=−0.4498, p=0.0032; r=0.4344, p =0.0045, respectively) as shown in Figures  S3B , C. Similar to myeloid cells, sex had no significant effect on all T-cell subtypes (data not shown). Interestingly, in NK and T cell types, the transcriptional activator genes and functional immune genes, such as FOS, JUN, CD69, CXCR4, NFKBIA and TNFAIP3, were generally elevated in mild/moderate recovered patients ( Figure 3B) , while the expression of GIMAP4, SELPLG, S1PR1, and STAT1 genes decreased ( Figure 3B ). In the KEGG pathway analysis, CD4m T cells exhibited enriched Th17 cell differentiation in mild/moderate recovered patients, while severe/critical recovered patients lacked IL-17 and TNF signaling pathway ( Figure 3C ). IL−17 signaling pathway was also downregulated in severe/ critical recovered patients for CD8m T (GZMK) ( Figure 3C ). In CD4m T (GZMK), the cytotoxic activity-related PRF1 and SLC9A3R1 gene expression level was significantly reduced, KLF6 and CD83 were significantly increased in mild/moderate recovered patients, and the expression level of these genes in severe/critical recovered patients was comparable to healthy controls ( Figure 3D) . These four functional genes had similar expression pattern in the different severity for three CD8 memory T cell types ( Figure S3D ). IFNG was significantly increased in the mild/moderate recovered patients for the three CD8 memory T cell types ( Figure S3D ). These results suggested that the percentages of CD4 T cells and CD8 T cells was significantly reduced in severe/critical recovered patients, but the expression level of transcriptional activator genes and functional immune genes was comparable to healthy controls or had a slightly downward trend. Although the percentage of CD4m T, CD4m T (GZMK), and three CD8 memory T cell types in mild/moderate recovered was equal to or slightly higher than that in healthy controls. The memory T cells differentiationrelated genes were significantly upregulated. These different immune characteristics may produce different SARS-CoV-2specific memory T cells in different recovered patients.

Our sequencing data also included scTCR-seq and scBCR-seq dataset to investigate the characteristics of TCR/BCR immune repertoires in recovered COVID-19 patients. UMAP results showed the distribution of TCR/BCR clone size in T/B cell subpopulations. The expanded clonotypes of the TCR clonal size >=5 were mainly from CD8m T (GZMH), CD8m T (GZMK), and Vg9Vd2 T cells (Figures 4A, B) . Also, the expanded BCR clonotype ratio was very small. The clonotypes of BCR clonal size >=5 mainly came from plasma and memory B cells, and the clonotypes of BCR clonal size 2-4 mainly came from naïve B(TCL1A) and memory B cells (Figures 4A, B) . Surprisingly, we found comparable T-/B-cell clonal expansion in healthy controls, mild/moderate, and severe/critical recovered patients ( Figure 4C and Figure S4A ). This was not consistent with the results reported by Zhang et al. (8) . Whether this inconsistency was caused by the difference of sample types (the sample types used by Zhang et al. were peripheral blood CD3+ T cells and CD3−CD19+CD20+CD27+ antigen-experienced B cells) or the number of productive TCR/BCR clones obtained by sequencing or number of samples, it needs to be further verified. Further, in these cell types, clinical severity did not affect the TCR/BCR diversity ( Figure 4E ). The TCR/BCR diversity of CD8m T (GZMH), CD8m T (GZMK), and memory B was negatively correlated with age (R<0), while that of prolif T and plasma had a positive trend with age with no significant (p>0.05; Figure 4D ). In addition, sex does not affect the TCR/BCR diversity of these cells ( Figure S4B) .

Next, to reveal the unique gene patterns and preferences of BCR or TCR in recovered COVID-19 patients, we compared the usage rate of immunoglobulin variable (V) genes. The severe/critical recovered patients use IGHV3-21 more frequently, compared to healthy controls ( Figure S4C) . Patients in the recovery stage had a certain bias towards the usage of other V genes as well, but there were large variations among individuals, and there was no statistically difference ( Figure S4C ). Finally, we analyzed the distribution of the heavy chain CDR3 (HCDR3) amino acid sequence length in the BCR repertoires of memory B cells. There was no significant difference in the distribution of HCDR3 lengths in memory B cells between the recovery stage of COVID-19 and the healthy control ( Figure S4D ).

In myeloid cell subtypes, the HLA class II score significantly decreased in CD14 monocytes and DCs of severe/critical recovered patients. We also found that the antigen processing or presentation pathways of CD14 monocytes and pDCs were downregulated in severe/critical recovered patients. We further investigated whether these changes contributed to the differences in CD4m T and CD8m T cell proportions, and cell-cell interaction analysis on the main subgroups of monocytes, DCs, and CD4m T cells and CD8m T cells was performed ( Figure 5 ; see Methods).

For each signaling pathway considered for the cell-cell interaction analysis (see Methods), we compared the aggregated incoming and outgoing signaling for each cell population. While mild/moderate recovered and severe/critical recovered patients presented similar overall signaling patterns with different patterns to healthy controls ( Figure 5A ). Based on these signaling patterns, we focused on the interaction of MHC-II, TNF, IL1, IL16, INF-II, and CD48 and other ligand-receptor pairs in monocytes, DCs, and CD4 and CD8 memory T cells ( Figures 5B, C) . In the Circos plot, we found in an unbiased manner that the MHC II signaling was downregulated in CD14 monocytes, DCs, pDCs in severe/critical recovered patients compared to healthy controls. The TNF_TNFRSF1B ligandreceptor pair related to inflammation expression was upregulated in mild/moderate recovered patients compared to healthy controls ( Figure 5B) . These results were consistent with , mild/moderate recovered (n = 16), severe/critical recovered (n = 6) patients. T tests (and non-parametric tests), *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. our previous observations with HLA class II score and inflammatory score trends ( Figures 2D, E) . Surprisingly, HLA-DRA, HLA-DRB1, and other HLA class II genes were downregulated in CD14 monocytes, DCs, pDCs in severe/ critical recovered patients, but they were upregulated in mild/ moderate recovered patients ( Figures 5B, C) . 

The immune status of COVID-19 patients during the recovery stage is the key to whether they can obtain immune protection against SARS-CoV-2 reinfection. Basic immune memory and high serum antibody levels can be obtained about 28 days after foreign antigens invade the human body (36) . Therefore, collecting peripheral blood of recovered patients 27-47 days after the onset of SARS-CoV-2 infection is one of the best choices to understand the immune characteristics of COVID-19 patients during the recovery stage, and find molecular markers related to disease severity or protection rate against reinfection. In this study, we generated single-cell sequencing data from the blood of SARS-CoV2 recovered patients and performed an integrated analysis with published data. We constructed a single-cell transcriptome landscape map of peripheral immune cells of recovered COVID-19 patients.

The results of the proportions and cytokine scores of cell subtypes in recovered patients in our study are consistent with those reported by Ren et al. (9) . However, due to differences in cell population annotations and the number of samples included, the results of some cell subtypes cannot be directly compared. The samples we analyzed only included samples with a sampling time of 27-47 days (days after symptom onset). To determine if this impacted on results, we extracted sample data from the Ren dataset with 21-90 days sampling time (days after symptom onset) and assessed how the sampling time, severity, sample type, gender, and age influence cell subtype proportions ( Figure S2A ). We found gender has no significant effect on all cell subtypes, and age has significant effects on CD4m T, CD4m T (GZMK), CD8m (IL7R), prolif T, gd T, and Vg9Vd2 T ( Figure S2A and Figure S3C ), all consistent with our subpopulation results. This shows that these findings were verifiable and stable in recovered COVID-19 patients for at least 2-3 months after symptom onset. Compared with Zhang et al.'s study on isolated T cells and B cells (8) , our study obtained fewer TCR and BCR clones in each sample. Therefore, we have obtained very few characteristics of the T-and B-cell immune repertoire of recovered COVID-19 patients.

Our findings on HLA class II are different from those of Wilk et al. The Wilk et al. study focused on the progression stage of COVID-19 patients, while our study discusses the recovery stage. In the study of Wilk et al. (6) , the HLA class II score was downregulated in patients with COVID-19 in the progression stage and the most decreased in severe/critical patients. In samples with sampling date 27-47 days after symptom onset, severe/critical recovered patients maintained lower HLA class II scores and higher ISG scores, compared to healthy control, while mild/moderate recovered patients were comparable to healthy controls. Furthermore, through cell-cell interaction analysis, we found several HLA class II genes, which are downregulated in severe/critical recovered patients but upregulated in mild/ moderate recovered patients. These results suggested that the immune response of mild/moderate patients could gradually return to normal levels during the recovery stage, while severe/ critical patients could remain in an immune disorder state. This might cause lower protection rate of severe/critical recovered patients against SARS-CoV-2 reinfection.

We also found that there was a higher proportion of CD14 monocytes in severe/critical recovered patients. These CD14 monocytes exhibited phenotypes such as upregulation of cytokine expression, downregulation of HLA class II genes, and antigen processing and presentation signaling pathway. Although the proportion of dendritic cells did not change, they also showed upregulation of cytokines and ISGs expression, and downregulation of HLA class II genes in severe/critical recovered patients. We observed a decreased proportion of CD4m T and CD8m T cells and showed a phenotype of downregulation of IL-17 and TNF signaling pathways in severe/critical recovered patients. This could suggest that the low antigen processing of dendritic cells and monocytes might negatively affect the memory T cell differentiation necessary to provide protection against SARS-CoV2 reinfection. The downregulation of HLA class II genes may be one of the reasons why convalescent patients with high serum anti-spike titers produced a higher proportion of non-neutralizing antibodies (15, 16, 18) .

The proportions of most peripheral immune cell types in mild/moderate recovered patients were equivalent to that of healthy controls, but CD14 monocytes exhibited upregulated expression of inflammatory genes such as TNF. The expression of HLA-DRA, HLA-DRB1, and other HLA class II genes were also upregulated in CD14 monocytes and DCs. Similarly, T-cell differentiation regulation and memory T cell-related genes FOS, JUN, CD69, CXCR4, and CD83 were upregulated. Therefore, we believe that the recovery of CD14 monocytes includes the return of cytokine expression and cell proportions to healthy levels. The upregulation of HLA-DRA, HLA-DRB1, and other HLA class II genes in CD14 monocytes and DCs may promote the change of CD4m T cell and CD8m T cell transcriptomes, helping the formation of T cell immune memory, thereby providing effective cellular immunity against SARS-CoV-2 reinfection. Further experiments are required to validate this hypothesis. Recent studies from Public Health England, London, UK, Andrews et al. found that the effectiveness of the Vaxzevria and Comirnaty vaccines against symptomatic diseases has greatly waned in people over 65 years of age (Unpublished, https:// twitter.com/jburnmurdoch/status/1438100712441974786?s=19). This outcome seems to be similar to that found in SARS-CoV-2infected people. This suggests that the differential immune characteristics we found in mild/moderate recovered and severe/critical recovered patients may be the keys to the development of an effective SARS-CoV-2 vaccine for the elderly.

SARS-CoV-2-specific memory T cell responses are long-lasting in recovered COVID-19 patients (37-39) even though SARS-CoV-2-specific antibody response may decrease (15, 37, (40) (41) (42) . Previous studies have shown that T cell responses to SARS-CoV-1 and Middle East Respiratory Syndrome Coronavirus (MERS-CoV) are long-lasting, up to >17 years (43) (44) (45) . Recent studies have also demonstrated that SARS-CoV-2-specific memory T cell response lasts for more than 10 months in recovered COVID-19 patients, and these cells are stem-like memory T cells with multifunctionality and proliferation ability (38) . In addition, Katherine et al. found memory T cells contribute to protection against SARS-CoV-2 rechallenge in a rhesus monkey model (29, 30) . These results support that the generation and persistence of memory T cells in recovered COVID-19 patients are essential for humans to prevent SARS-CoV-2 reinfection. Currently, there is no report on myeloid cell immunity in recovered patients 1 year or more after infection, and their role in preventing SARS-CoV-2 reinfection remains unknown. In addition, it is still unclear whether the T cell responses in severe COVID-19 patients who have recovered a year or longer after infection have returned to normal.

In summary, our analysis of a large-scale scRNA-seq dataset covering diverse disease severity has revealed multiple immune characteristics during the recovery stage of COVID-19 that have not been adequately studied previously. Such results provided a critical resource and important insights in dissecting the human body's immune protection mechanism against SARS-CoV-2 reinfection and may help to develop effective SARS-CoV-2 vaccines for the elderly.

Experimental protocols and the data analysis pipeline used in our work follow the 10X Genomics, Scanpy, Immunarch, and CellChat official documentation. The analysis steps, functions, and parameters used are described in detail in the Methods section. Custom scripts for analyzing data are available upon reasonable request.

The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive of the Beijing Institute of Genomics (BIG) Data Center, BIG, Chinese Academy of Sciences, under accession code HRA000864 and are publicly accessible at http://bigd.big.ac.cn/gsa-human.

The studies involving human participants were reviewed and approved by the Research Ethics Committee of the School of Public Health (Shenzhen) (2020-007), Sun Yat-sen University, China. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin. Written informed consent was obtained from the individual(s), and minor(s)' legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

YS designed the study. XL performed experiments, analyzed the data, and wrote the manuscript. ZM and MG analyzed the data and wrote the manuscript. QL, LY, ML, and ZW contributed to the data collection, data analysis, and data interpretation. HL, TJ, WW, and SF contributed to the samples collection and experiments. YS, ZM, and SF supervised the study and revised the manuscript. YB, NG, IP, and AB participated in revising and improving the English language of the manuscript. All authors contributed to the article and approved the submitted version. 

Potent Neutralizing Antibodies Against SARS-CoV-2 Identified by High-Throughput Single-Cell Sequencing of Convalescent Patients' B Cells

COVID-19 Severity Correlates With Airway Epithelium-Immune Cell Interactions Identified by Single-Cell Analysis

Single-Cell Landscape of Bronchoalveolar Immune Cells in Patients With COVID-19

Multi-Omics Resolves a Sharp Disease-State Shift Between Mild and Moderate COVID-19

Immune Cell Profiling of COVID-19 Patients in the Recovery Stage by Single-Cell Sequencing

A Single-Cell Atlas of the Peripheral Immune Response in Patients With Severe COVID-19

Single-Cell Landscape of Immunological Responses in Patients With COVID-19

Adaptive Immune Responses to SARS-Cov-2 Infection in Severe Versus Mild Individuals

COVID-19 Immune Features Revealed by a Large-Scale Single-Cell Transcriptome Atlas

Single-Cell Transcriptome Analysis Highlights a Role for Neutrophils and Inflammatory Macrophages in the Pathogenesis of Severe COVID-19

A Human Circulating Immune Cell Landscape in Aging and COVID-19

T-Cell Hyperactivation and Paralysis in Severe COVID-19 Infection Revealed by Single-Cell Analysis

Imbalance of Regulatory and Cytotoxic SARS-CoV-2-Reactive CD4(+) T Cells in COVID-19

SARS-CoV-2 Mrna Vaccination Induces Functionally Diverse Antibodies to NTD, RBD, and S2

Antibody Responses to SARS-Cov-2 in Patients With COVID-19

The Characterization of Disease Severity Associated Igg Subclasses Response in COVID-19 Patients

Seroprevalence and Humoral Immune Durability of Anti-SARS-CoV-2 Antibodies in Wuhan, China: A Longitudinal, Population-Level, Cross-Sectional Study

Antibody Response to SARS-CoV-2 Infection in Humans: A Systematic Review

Assessment of Protection Against Reinfection With SARS-Cov-2 Among 4 Million PCR-Tested Individuals in Denmark in 2020: A Population-Level Observational Study

Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data

SCANPY: Large-Scale Single-Cell Gene Expression Data Analysis

Sensitive and Accurate Integration of Single-Cell Data With Harmony

Dimensionality Reduction for Visualizing Single-Cell Data Using UMAP

Local Leaders in Random Networks

Putative Cell Type Discovery From Single-Cell Gene Expression Data

MAST: A Flexible Statistical Framework for Assessing Transcriptional Changes and Characterizing Heterogeneity in Single-Cell RNA Sequencing Data

Profiler: A Web Server for Functional Enrichment Analysis and Conversions of Gene Lists

Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles

Correlates of Protection Against SARS-CoV-2 in Rhesus Macaques

Immunomind/ Immunarch: 0.6.5: Basic Single-Cell Support (0.6.5). Zenodo

Inference and Analysis of Cell-Cell Communication Using Cellchat

The Molecular Signatures Database (Msigdb) Hallmark Gene Set Collection

New Perspectives on Genomes, Pathways, Diseases and Drugs

Butyrophilin 2A1 is Essential for Phosphoantigen Reactivity by Gammadelta T Cells

Different B Cell Populations Mediate Early and Late Memory During an Endogenous Immune Response

Persistence of SARS-Cov-2-Specific B and T Cell Responses in Convalescent COVID-19 Patients 6-8 Months After the Infection

SARS-CoV-2-Specific T Cell Memory is Sustained in COVID-19 Convalescent Patients for 10 Months With Successful Development of Stem Cell-Like Memory T Cells

T Cell and Antibody Kinetics Delineate SARS-Cov-2 Peptides Mediating Long-Term Immune Responses in COVID-19 Convalescent Individuals

The Dynamic Changes of Serum Igm and Igg Against SARS-CoV-2 in Patients With COVID-19

SARS-CoV-2 Immunity and Functional Recovery of COVID-19 Patients 1-Year After Infection

Clinical and Immunological Assessment of Asymptomatic SARS-CoV-2 Infections

SARS-CoV-2-Specific T Cell Immunity in Cases of COVID-19 and SARS, and Uninfected Controls

Recovery From the Middle East Respiratory Syndrome is Associated With Antibody and T-Cell Responses

Lack of Peripheral Memory B Cell Responses in Recovered Patients With Severe Acute Respiratory Syndrome: A Six-Year Follow-Up Study

We would like to thank Professor Xiangjun Du for providing part of the computer resources needed for data analysis for this research, and Dr. Guoqin Mai for her guidance in BCR data analysis. Thanks to Jiaying Yang, Feng Zhang, Xiaohao Xu, Yani Peng, Liangliang Wang, and Sidian Xie from professor YS's team for their help in the preliminary analysis of BCR data. We would also like to thank Dr. Alvis Brazma and Dr. Irene Papatheodorou for their guidance and support. We would also like to thank Jonathan R. Manning for helping to improve the English of the manuscript.

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.781432/ full#supplementary-material Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.Publisher's Note: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.Copyright © 2022 Li, Garg, Jia, Liao, Yuan, Li, Wu, Wu, Bi, George, Papatheodorou, Brazma, Luo, Fang, Miao and Shu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.