Thymosin alpha-1 Protected T Cells from Excessive Activation in Severe COVID-19

Two typical features of uncontrolled inammation, cytokine storm and lymphopenia, are associated with the severity of coronavirus disease 2019 (COVID-19), demonstrating that both innate and adaptive immune responses are involved in the development of this disease. Recent studies have explored the contribution of innate immune cells to the pathogenesis of the infection. However, the impact of adaptive immunity on this disease remains unknown. In order to clarify the role of adaptive immune response in COVID-19, we characterized the phenotypes of lymphocytes in PBMCs from patients at different disease stages using single-cell RNA sequencing (scRNA-seq) technology. Dynamics of the effector cell levels in lymphocytes revealed a distinct feature of adaptive immunity in severely affected patients, the coincidence of impaired cellular and enhanced humoral immune responses, suggesting that dysregulated adaptive immune responses advanced severe COVID-19. Excessive activation and exhaustion were observed in CD8 T effector cells, which might contribute to the lymphopenia. Interestingly, expression of Prothymosin alpha (PTMA), the proprotein of Tα1, was signicantly increased in a group of CD8 T memory stem cells, but not in excessively activated T cells. We further showed that Tα1 signicantly promoted the proliferation of activated T cells in vitro and relieved the lymphopenia in COVID-19 patients. Our data suggest that protection of T cells from excessive activation might be critical for the prevention of severe COVID-19.


Introduction
The current outbreak of the new coronavirus SARS-CoV-2 has resulted in a global pandemic. The World Health Organization (WHO) has reported 4,444,670 con rmed cases and 302,493 deaths in 212 countries and territories as of 15 May 2020. Symptoms of COVID-19 caused by an infection of this virus include fever, cough, fatigue, chills, diarrhea, muscle pain, headache, sore throat, shortness of breath, and pneumonia 1,2 . The illness can range from mild to severe, while some infected people are asymptomatic. About 15% of the con rmed cases advanced to severe diseases 3 . Patients developed acute respiratory distress syndrome (ARDS) and multiple organ failure have a high risk of mortality 4 . Due to the lack of a speci c drug against this virus, the current clinical management of this disease mainly depends on supportive care to reduce in ammatory responses and to keep the lung functioning 5 .
Understanding the underlying immunopathology of COVID-19 is therefore of paramount importance for improving the current treatment.
The outcome of a viral infection is mainly determined by the host immune system which is consisted of innate and adaptive immune cells. Appropriate immune responses eliminate invading viruses and establish immune memory for future protection 6 . However, an uncontrolled immune response may lead to severe immunopathology 7 . It is believed that the illness of coronavirus infection is resulted from abnormal host immune responses. The SARS-CoV-2-caused coronavirus disease is characterized by uncontrolled in ammation, which may be similar to diseases resulted from the infections of other pathogenic coronaviruses such as the severe acute respiratory syndrome coronavirus-1 (SARS-CoV-1) and the Middle East respiratory syndrome coronavirus (MERS-CoV) 2,8,9 . In ammatory cytokine storm and lymphopenia were frequently observed in COVID-19 1,3 , indicating that both innate and adaptive immune responses are involved in the development of this disease.
Most of the severe COVID-19 patients experienced cytokine storms which were characterized by the elevated plasma levels of pro-in ammatory cytokines, chemokines, and other in ammatory molecules such as C-reactive protein and D-dimer 2,3 . Increased numbers of In ammatory monocytes and neutrophils with upregulated expression of IFN-stimulated genes (ISGs) were also observed in these patients 2 . High levels of pro-in ammatory cytokines and chemokines might facilitate the in ltration of these in ammatory immune cells into infected organs, and thus mediating the immunopathology. These observations suggested that an aggressive innate immune response to the infection of SARS-CoV-2 was involved in the pathogenesis of this disease. Clinical trials of therapies reducing in ammatory cytokines were thereby quickly initiated 10-12 .
The adaptive immune response plays an essential role in viral clearance. SARS-CoV-2-speci c T cells and antibodies were detected in COVID-19 patients, indicating that the cellular and humoral immune responses were developed after the infection 13,14 . Some monoclonal antibodies cloned from B cells of infected patients could block the binding of the spike protein to the ACE2 receptor 15 . Nevertheless, recent studies found that the adaptive immune response was impaired in severely affected patients 16 . In addition to the decrease in lymphocyte number in blood, a common COVID-19 abnormality associated with the disease severity and poor clinical outcome, increased expression of exhaustion markers such as PD-1 and TIM3 on T cells was observed 17,18 . Dysfunction of adaptive immunity may not only impair the elimination of the infecting viruses but also lose the surveillance on other pathogens, which may lead to exacerbated infections and complicating diseases such as bacterial infections 19,20 . These ndings suggested that an insu cient adaptive immunity facilitated the progression of the severe disease.
Indeed, high rate of severe COVID-19 was found in immunocompromised patients 2,3,21 , supporting that an insu cient rather than an overactive immunity could be the basis for the development of this disease.
On the other side, a strong humoral immune response could lead to antibody-dependent enhancement (ADE) of viral entry into immune cells [22][23][24][25] . ADE-mediated hyperactivation of macrophages and monocytes may exacerbate the disease by secreting more pro-in ammatory molecules 24 . Therefore, it is unknown if the virus-speci c antibodies generated in COVID-19 patients are protective or pathogenic.
Recent studies have focused on the mechanisms of innate immune response to the infection of SARS-COV-2. However, there is limited understanding on the role of the adaptive immunity in the immunopathology of COVID-19. In order to comprehensively understand the pathogenesis of this disease, it is necessary to characterize the adaptive immune response during the course of the disease, which may provide critical information for improving the clinical management of infected patients and facilitating the development of vaccine.
In this study, we evaluated the adaptive immune response in COVID-19 by characterizing phenotypes of lymphocytes in patients' peripheral blood mononuclear cells (PBMCs) collected at different disease stages. In consistent with other reports, signi cant reduction in T cells and increased T cell exhaustion were found in patients with severe disease. We further found that CD8 T effector cells displayed phenotypes of excessive activation. Strikingly, high levels of plasma cells and activated B cells were observed at the mean time. T cell-mediated cellular immune response plays an essential role in the clearance of infected host cells, which produce most of the viruses after infection. Failed elimination of these infected cells could escalate the viral production. Without su cient T cell immunity, a strong antibody response may trigger ADE of viral infection. We therefore revealed an important feature of the adaptive immunity in severe disease, the coincidence of an impaired cellular immune response and an enhanced humoral immune response, which might facilitate the progression of the severe disease. We further found that Tα1 slightly reduced T cell activation but promoted proliferation of activated T cells. Most importantly, administration of Tα1 relieved the lymphopenia in COVID-19 patients. Taken together, our data revealed a dysregulated adaptive immune response in severe COVID-19, and suggested that early protection of T cells from excessive activation might be critical for the prevention of severe disease.

Results
Elderly patients have a high risk of lymphopenia.
A high rate of severe COVID-19 was reported in immunocompromised patients, suggesting that an insu cient rather than an overactive antiviral immunity could be the basis for the development of this disease 2,3,21 . Meanwhile, lymphopenia, a reduction in the number of lymphocytes in the blood, was associated with the disease severity of COVID-19 2,16 . We analyzed the incidence of lymphopenia in 284 patients infected with SARS-CoV-2 collected from designated hospitals in Fujian province of China (basic information of the patients was listed in Table S1), and found that a reduction of lymphocytes was more frequently observed in aged patients (Fig. 1a). These ndings denote the pivotal role of the adaptive immunity for the viral clearance and disease control.
scRNA-seq analysis identi ed main clusters in lymphocyte.
We hypothesized that single-cell transcriptomic analysis of peripheral lymphocytes during the course of the disease might clarify the role of the adaptive immune system in this disease. Therefore 13 samples of PBMCs were collected from 10 patients at different stages of the disease, namely pre-severe disease (PR, 1 sample), severe disease (SD, 3 samples), post severe disease (PS, 3 samples), post mild disease (PM, 3 samples) and convalescence of mild disease (CM, 3 samples), which also included seven samples from our previous report focusing on innate immune cells (samples without enough number of lymphocytes were not included) 11 . Among all enrolled patients, four of them experienced severe disease under intensive care and the other six showed fever, cough, chill, and fatigue etc mild symptoms. Results of COVID-19 examination and times of sample collection were shown in Fig. 1b. Patient information was listed in Table S1. Samples of PM were collected from patients during hospitalization whose viral testing turned negative after treatment. Convalescent samples were collected within the rst 5 days after hospital discharge. Three normal PBMCs from healthy donors (HD) were used as controls (named HD-1, HD-2, and HD-3).
We performed single-cell mRNA sequencing (scRNA-seq) of these samples on the 10x genomics platform. 243 million RNA transcripts in 91,649 cells were obtained after ltering cells with low quality, in which 61,386 cells were from COVID-19 patient samples. Canonical correlation analysis (CCA) was performed in all immune cells to reduce batch effect (Fig. S1a). We then used t-Distributed Stochastic Neighbor Embedding (t-SNE) to visualize the clusters of all the cells identi ed by a shared nearest neighbor (SNN) modularity optimization based clustering algorithm implemented in Seurat ( Fig. 1c) 26 .
Three major groups of immune cells in PBMCs, namely myeloid, B and T cells were clearly identi ed by speci c gene expression signatures (Fig. 1d). Classic lineage markers further con rmed the identities of these clusters (Fig. 1e). Consistent with reported studies, unique cluster with features of NK cells was not identi ed 27 . We then focused on adaptive immune cells (CD19+, CD79A+, CD79B+, CD3D+, CD3E+, and CD3G+ cells), and grouped the lymphocytes into B cells, Plasma Cells, CD4 and CD8 T cells ( Fig. 1f and S1b). Analysis of the lymphocyte composition at different disease stages revealed a signi cant reduction in T cells along with an increase in B cells in SD and PS samples (Fig. 1h), indicating the abnormal dynamics of adaptive immune cells in patients with severe disease.
Humoral and cellular immune responses in patients with severe disease.
Effector cells carry out the major functions of the adaptive immune system, and decide the e cacy of adaptive immunity against viral infection. We checked the distribution of effectors in the clusters of lymphocytes. Despite the small numbers of cells identi ed, the plasma cells were prominently grouped ( Fig. 1f). The T effectors were identi ed by the expression of Granzyme B (GZMB), an effective molecule speci cally expressed in activated T cells. As shown in Fig. 1g, most of the effector T cells were located in the CD8 cluster, while CD4 effectors were barely detected.
In order to examine changes of effectors in different types of adaptive immune cells, we performed separate cluster analysis on the B, CD4 and CD8 T cells. B cells were grouped into 4 major clusters, namely mature, memory, activated B cells and plasma cells (Fig. 2a, S1c and S1f). CD4 T cells consisted of naïve, memory, effector and regulatory T cell clusters ( Fig. 2c and S1d), while clusters of naïve, memory, effector and mucosal-associated invariant T cells (MAIT) were found in CD8 T cells ( Fig. 2e and S1e). We found that the levels of effectors in these cells displayed different dynamics during the disease. The frequency of effectors in B cells was signi cantly increased at the SD stage but substantially reduced at the PS stage (Fig. 2b). Among CD4 T cells, low levels of effectors were observed in all of the samples except for one from PR (Fig. 2d). Consistent with the above observation (Fig. 1g), a large portion of effectors was found in CD8 T cells (Fig. 2e), but a signi cant low level of CD8 effectors was found at the PS stage (Fig. 2f).
Next, we examined changes of effectors in the humoral and cellular immune responses by measuring the frequencies of B (plasma cells plus activated B cells) and T effectors (CD4 plus CD8 effectors) in total lymphocytes respectively. In consistent with the observation in B cells, high levels of B effectors were observed in lymphocytes at the SD stage (Fig. 2g). For T effectors, however, low frequencies were found in lymphocytes at both the SD and PS stages (Fig. 2h), which was different from the observation in CD8 T cells among which low levels of effectors were only observed at the PS stage (Fig. 2f). This could be explained by the reduced portion of total T cells in lymphocytes at these two stages (Fig. 1h).
Taken together, the dynamics of effector levels in lymphocytes revealed a coincidence of high humoral and low cellular immune responses in patients with severe disease. The drastically damaged T cell immunity might impair the elimination of infecting viruses, while a strong antibody response with insu cient T cell support could mediate the ADE of viral infection. The disordered adaptive immune response likely contributes to the rapid progression of severe COVID-19.
Excessive activation led to exhaustion of CD8 T effector cells in sever COVID-19.
In order to understand the reduction of effector cells, we increased the modularity to group more distinct subsets, which allow us to explore the changes inside the major clusters. In CD8 T cells, two naïve, three memory and two effector clusters were identi ed after stringent grouping (Fig. 3a). No sub-clusters were detected in the MAIT cells. Gene expression signatures and pseudotime analysis demonstrated unbiased cellular dynamic processes of these sub-clusters ( Fig. S2a and S2b). We then compared the two groups of effector T cells. Differentially expressed gene (DEG) and pathway enrichment analysis showed that one group was excessively activated, which was then named as Tea (excessively activate T cells). Compared to the normal effectors (Te), Tea cells had higher expression of activation genes (FCER1G, KLRB1, KLRF1, NKG7 and IGFBP7), co-inhibitory receptors (LAG3, CD300A, CD244, CD160 and HVCR2), IFN downstream signaling molecules (JAK1, IFITM2, IFITM3, TYK and IRF8), and chemokines (CCL3 and CCL5) ( Fig. 3b and 3c). Signi cant reductions in TCR signaling molecules (CD3D, CD3G, CD8A, CD8B and ZNF683), and cytokine/receptor genes (IL7R, IL2RB, IL32 and IFNG) were found in Tea cells. Increased expression of NKT cell genes (CD160, TRROBP, FCER1G and XCL1) was also observed in Tea cells, suggesting that excessive activation might happen in NKT cells as well ( Fig. S2a and S2c). For effector molecules, Tea cells expressed less GZMB but more PRF1 than Te cells. Although it is known that activation induces cell death in T cells 28-31 , we did not observe any signi cant difference in expression levels of apoptosis genes or pathways between these two clusters. This might have been due to the sample process of single cell analysis which excluded apoptotic cells to ensure the accuracy of clustering. Indeed, although the Tea cells expressed slightly more FAS and FASLG, both Tea and Te clusters expressed extremely low levels of these two genes (Fig. S2d). IFN signaling might be involved in the differentiation of Tea cells. Although CCL3 and CCL5 secreted by Tea cells may promote innate immunity, the marginal increase of these two chemokines suggested that they had a limited impact on innate immune cells. As shown in Fig. 3d, CD8 effectors were dominated by Tea cells in SD samples, indicating that the adaptive immune system had been overwhelmed by these dysfunctional cells. Together, these data support that the Tea cells are a group of excessively activated T cells with exhausted phenotype and diminished function of antigen recognition. Excessive activation-induced continuous expansion of Tea cells may waste the majority of T cells, and lead to lymphopenia which paralyzes the adaptive immune system.
We then checked the CD8 memory T cells. The three sub-clusters were named as Tm-1, Tm-2 and Tm-3. Signi cant accumulation of Tm-3 was observed in the SD and PS stage (Fig. 3e), suggesting that this group of cells was speci cally generated during severe disease. The percentage of Tm-3 in CD8 memory T cells was strongly correlated with the proportion of Tea in the CD8 effector cells (R 2 = 0.7813) (Fig. 3f), suggesting that Tm-3 might be a group of memory cells derived from the Tea cells during excessive activation. Differential gene expression analysis and pathway enrichment analysis showed that Tm-3 was a cluster of highly proliferating cells with features of memory stem cells ( Fig. S2e and S2f). Compared to Tm-1 and Tm-2, Tm-3 had higher expression levels of T memory stem cell markers (SELL, CXCR3, CCR7, FAS, CD27 and CD28), and the proliferation gene Ki67 (Fig. 3g) 32 . Furthermore, a slightly higher expression of GZMA and lower level of IL7R in the Tm-3, as compared to Tm-1 and Tm-2, also showed that this was a group of recently developed memory T cells. The low level of Tm-3 at the PR stage further con rmed its differentiation during severe disease stages (Fig. 3e).
Although both Tea and Tm-3 cells experienced excessive activation, exhaustion was observed in the Tea but not Tm-3 cells. We then performed DEG and pathway-enrichment analysis to understand the difference between these two groups of cells. The Tea cells had comparatively higher expression level of T cell activation and effector genes (GNLY, NKG7 and GZMB) but almost no expression of Ki67 ( Fig. 3h  and 3i), suggesting that they were terminal-differentiated cells and more prone to death. In contrast, the Tm-3 cells had high expression levels of genes involved in epigenetic modi cation (DNMT1 and EZH2), oxidative phosphorylation (NDUB6 and NDUFV2), regulation of telomerase, cell cycle and proliferation etc ( Fig. 3h and 3i). These distinct signaling transduction and epigenetic changes might have shaped the stem-like memory phenotype of Tm-3. Interestingly, we found the expression of PTMA was also signi cantly increased in Tm-3 ( Fig. 3h and 3j). Thymosin alpha-1 (Tα1), the rst 28 amino acids of PTMA, was found to be involved in T cell development, and has been used for the treatment of certain infection diseases including COVID-19 33-35 . In addition, Tea expressed less PTMA than Te, and Tm-3 expressed the highest level of PTMA in all the sub-clusters of CD8 T cells (Fig. 3j). We thus suspected that Tα1 might protect T cells from excessive activation.
We further analyzed the subsets in naïve T cells but did not observe any signi cant change during the different disease stages (Data not shown). Since no sub-clusters were identi ed in the effector clusters of B cells and CD4 T cells, we did not analyze the subsets in these two types of adaptive immune cells (Data not shown).
Thymosin alpha-1 protected T cells from excessive activation.
Next, we tested the effect of Tα1 on T cell activation. PBMCs from healthy donors were activated by anti-CD3/CD28 antibodies in vitro and treated with 200ng/ml Tα1 for 3 days, followed by cultures with IL-2 (200U/ml) and Tα1 (200 ng/ml) for 6 more days. Compared to the control group, Tα1 had signi cantly increased T cell numbers at day 6 and 9 although a slight decrease in T cell numbers was observed at day 3 (Fig. 4a), indicating that Tα1 promoted the proliferation of T cells after activation. After 3 days of activation, we found that cell size measured by FSC and SSC was slightly reduced in the group with Tα1 (Fig 4b), suggesting that these T cells were less activated. Tα1 treatment did not change the proportion of CD4 and CD8 T cells (Fig 4c), but reduced the production of IFNγ and TNFα, although no signi cant statistical difference was observed (Fig 4d). Signi cant reduction of GZMB was observed in both CD4 and CD8 T cells, while no change in the expression of PD-1 was observed. These data showed that Tα1 could reduce T cell activation and promote the proliferation of T cells after activation, indicating it may protect the T cells from excessive activation.
The use of Tα1 in some COVID-19 patients may allow the evaluation of its impact on lymphocytes. For this purpose, the data from 25 severe and critical COVID-19 cases treated at the Huoshenshan hospital (Wuhan, China) were collected (basic information of the patients was listed in Table S1). Of them, 11 patients received daily Tα1 treatment for at least one week, while the other 14 patients were not treated with Tα1 during the hospitalization period. Compared to the non-treated patients, the lymphocyte counts of the treated patients were signi cantly increased after one week of Tα1 treatment (Fig. 4e). The fold change of lymphocytes counts in each patient was shown in Fig. 4f. These data showed that Tα1 enhanced the number of lymphocytes in patients with severe and critical disease. Due to the limited number of patients in this retrospective analysis, we were not able to clearly evaluate the clinical bene ts of Tα1 treatment. Nevertheless, our data suggest that the administration of Tα1 could be a potential approach to protect T cells from excessive activation in COVID-19.

Discussion
Overall, our results revealed an unexpected elevation of effector B cells and the early absence of CD4 effector T cells in severely affected patients, while excessive activation-mediated exhaustion of CD8 T cells might be an important mechanism of lymphopenia found in patients with severe disease. These ndings suggested that a disordered adaptive immune response might facilitate the infection of SARS-CoV-2 and thus cause the outbreak of severe disease, highlighting that the early protection of T cell from excessive activation may be a key step to prevent severe disease.
The dynamics of B effector levels in lymphocytes observed in this study is similar to the reported antibody responses in patients infected with SARS-CoV-1, the virus that caused the severe acute respiratory syndrome (SARS). High activities of neutralizing antibody were quickly developed but diminished rapidly later in deceased SARS patients, while recovered patients slowly developed these activities and maintained them for a long time 20,22,24 . In addition, recent reports also observed increased IgG and total antibody titers in severe COVID-19 patients, and found that high antibody titer was associated with disease severity and worse outcome 23,25 . These ndings suggest a possible role of ADE in the progression of severe COVID-19, which is a potential challenge for vaccine development as well, and needs to be further addressed. Since the ADE of viral infection depends on the binding of virusspeci c antibodies to Fc receptors, administration of the Intravenous Immunoglobulin (IVIg), a blood product containing polyclonal immunoglobulin G from healthy donors, may block the Fc receptors and thus reducing the entry of virus into immune cells. Interestingly, a recent report showed that administration of high-dose IVIg reversed the deteriorating course of disease in 3 severe COVID-19 patients 36 . Although the immunomodulatory mechanisms of IVIg are not fully understood, it has been used in treating in ammatory diseases, prophylaxis and severe infections, and exhibited good tolerance 37,38 . It is worth testing if this therapy is able to block the ADE of viral infection in the future.
CD4 T cells contribute to various activities in immune responses, for instance, helping antibody responses by B cells, the expansion of CD8 effector T cells, and recruitment of immune cells to lymphoid tissue or sites of pathogen infection. The early absence of CD4 effector T cells in patients with severe disease may impede both the adaptive and innate immune responses against the SARS-CoV-2.
Despite the promising preliminary results observed in Tα1 treated clinical patients, more novel therapies that are able to synergize with the patients' adaptive immune responses are needed in future. Enhanced activation-induced cell death has been reported in patients with advanced age 29 . Future evaluation of excessive activation in elderly patients may clarify if this is the cause for the high rate of severe disease observed in this group of patients. Since the Tea cells identi ed in this study could re ect the damage of the adaptive immunity, they might serve as a good prognostic marker to identify patients at high risk of developing severe disease.

Antibodies and reagents
The antibodies for ow cytometry were purchased from Biolegend or BD. Tα1 (ZADAXIN) was purchased from SciClone Pharmaceuticals. Chromium TM Single cell 5' Library & Gel Bead Kit (1000006), Chromium TM Single cell 3'/5' Library Construction Kit (1000020) was purchased from 10x Genomics. Cytokines were purchased from Peprotech. All cell culture reagents were purchased from Gibco unless otherwise indicated.

Patients and study approval
This study enrolled a total of 10 COVID-19 patients. The 4 patients with severe disease were reported in our previous study 11 . The 6 patients with mild disease (3 at post mild disease stage and 3 at convalescence stage) were recruited from the Fifth A liated Hospital (Zhuhai, China) and the Third A liated Hospital (Guangzhou, China) of Sun Yat-sen University. Lymphocyte count data from 25 SARS-CoV-2 patients were collected in Huoshenshan hospital (Wuhan, China) by managing doctors.
Lymphocyte count data from 256 SARS-CoV-2 patients were collected from designated hospitals in Fujian province during January 3 2020 to march 13 by managing doctors. Informed consent has been obtained from all participated patients. This study received approval from the Research Ethics Committee of the Sun Yat-sen University Cancer Center, the Fifth A liated Hospital, the Third A liated Hospital of Sun Yat-sen University, the Huoshenshan hospital and the Fujian Provincial Hospital, China (K176-1, HSSLL030, K2020-03-01).

PBMCs Isolation
All procedures were carried out within P2+ laboratory certi ed for studies of infectious materials. PBMCs were isolated from peripheral venous whole blood samples density gradient cell separation (Ficoll, TBD Science, China), and were preceded for Single cell transcriptome sequencing immediately.

Cell culture
Fresh PBMCs from healthy donors were cultured in IL-2 (200U/ml) containing media in 96-well ubottomed plates coated with anti-CD3/CD28 (5ug/ml respectively, Peprotech). To test the effect on T cell activation, Tα1 (200ng/ml) was added into T cells culture media while seeding. Cell numbers were recorded at day 3, 6 and 9. Cells harvested after 3 days activation were analyzed by ow cytometry.

Flow cytometry
The antibodies used for cell surface labeling were BUV737 anti-human CD4 (564305, BD), BUV395 antihuman CD8 (563795, BD), FITC anti-human PD-1 (329904, Biolegend). The antibodies used for intracellular staining were PE anti-human IFNG (506507,Biolegend), APC anti-human TNFA (502912, Biolegend), BV421 anti-human Granzyme B (563389, BD). All antibodies for ow cytometry staining were used at 1:200 dilution. Cells were harvested and washed once in 2ml PBS and labeled on ice with indicated antibodies for analysis of surface markers. To determine intracellular cytokine expression, the cells were pre-incubated with Brefeldin A (Biolegend) for 4h before staining. After incubation, cells were stained following surface staining, xation and permeabilization using the BD Transcription Factor Buffer Set (562574, BD) according to the manufacturer's instructions. All samples were analyzed on a LSRFortessa™ X-20 cell analyzer (BD Biosciences). Analysis of acquired data was performed with the FlowJo software (FlowJo LLC).
Single cell transcriptome sequencing and data preprocessing The single cell RNA libraries were prepared using the Chromium TM Single cell 5' Reagent Kit of Chromium platform (10x Genomics, USA) following the manufacturer's instruction. Generated single cell RNA libraries were sequenced on the Illumina HiSeq X Ten platform. The CellRanger software (version 3.1.0) was used for preprocessing of the PE150 Illumina sequencing reads. Brie y, raw reads in bcl format were converted to FASTQ format using "cellranger mkfastq", and then the reads in FASTQ format were aligned to human genome reference (hg38, GRCh38) using STAR, and then "cellranger count" was used to derive gene expression matrix for each sample.
Determination of cell types from single cell transcriptome sequencing data Seurat (v3.1.3) R toolkit was used to analyze the single cell transcriptome sequencing data. Firstly, cells with low quality were ltered out. Brie y, the dead or dying cells with more than 20% mitochondrial RNA content were removed, and the cells with too low number (less than 200) or too high number (more than 2500) of genes were also removed, and the cells expressed more than one markers among the three markers (CD2, CD79A, CD68) were de ned as doublets and removed. Then, the ltered gene expression matrix for each sample was normalized using "NormalizeData" function in Seurat, and only highly variable genes were remained using "FindVariableFeatures" function in Seurat. Next, "FindIntegrationAnchors" and "Integratedata" functions in Seurat were used to integrate the gene expression matrices of all samples, where batch effects between different samples have been adjusted. Next, "RunPCA" function was used to perform the principle component analysis (PCA) and "FindNeighbors" function was used to construct a K-nearest-neighbor graph. Next, the most representative principle components (PCs) selected based on PCA were used for clustering analysis with "FindCluster" function to determine different cell types. Lastly, tSNE was used to visualize the different cell types.
We annotated the cell types using the following rules: Based on the most 10 differentially expressed genes that were derived using "FindAllMarkers" function in Seurat, genes such as CD2, CD3D, CD3E, CD3G and CD247 were used as T cell markers, and genes such as CD19, CD79A, CD79B, BLNK, FCRL5, MS4A1 were used as B cell markers, and genes such as CD14, CD163, CD68, CSF1R, FCGR2A, and CD33 were used as Myeloid cell mark. The percentage of CD4 gene expression and CD8A was counted to de ne CD4+ T cells or CD8+ T cells.
The CD19+, CD79A+, CD79B+, CD3D+, CD3E+, and CD3G+ lymphocytes were further clustered using the single cell analysis pipeline as described above. To get higher resolution clusters, the "resolution" parameter used in FindCluster was set from 0.3 to 0.5.

Differential gene expression analysis between different cell types
Seurat v3 was used to perform differential gene expression analysis between different cell types. For each cell type, the differentially expressed genes (DEGs) were obtained relative to all of the other cell types using "FindCluster" function in Seurat. DEGs between Tea cells and Te cells, and DEGs between Tm-3 and Tea were obtained using R package edgeR with log2 Fold change > 0.58 and P value < 0.5. Gene functional annotation For DEGs were performed using the R package clusterPro lter with default settings. The P values were corrected for multiple testing.

Single cell trajectories analysis
The R package Monocle 2 was used to analyze single cell trajectories in order to discover the cell-state transitions. We sampled 200 cells from each cell type and used genes with average expression > 0.1 to sort cells in pseudo-time order. "DDRTree" function in Monocle 2 was applied to reduce dimensions and "plot_cell_trajectory" function was used to visualize the minimum spanning tree on cells.
Statistical analysis.
All sample sizes were large enough to ensure proper statistical analysis. Statistical analyses were performed using GraphPad Prism (GraphPad Software, Inc.). P values < 0.05 were considered as statistically signi cant. All t-test analyses were one-tailed t-tests (paired or unpaired depending on the experiments). The number of replicates (n), number of independent experiments performed, and p values for each experiment are reported in the corresponding gure legends.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.