Main

The global coronavirus disease 2019 (COVID-19) pandemic caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has resulted in over 70 million recorded cases with a mortality rate of approximately 2.3%1. To date, the main strategy for ending the pandemic involves the development of effective vaccines, many of which are being evaluated2,3. Although most leading vaccine candidates can elicit virus-neutralizing antibodies, there is considerable uncertainty on how protective or durable these responses might be. Thus, a greater understanding of the innate and adaptive immune responses of individuals convalescing from COVID-19 is needed, particularly knowledge on whether immunological memory is established.

To date, studies on the immune responses to SARS-CoV-2 have mainly focused on the aspects of pathological inflammation3,4,5. According to single-cell transcriptome analyses, monocytes from the peripheral blood or the bronchoalveolar lavage fluid of COVID-19 patients are characterized by high levels of pro-inflammatory cytokines5,6,7,8 such as IL-6, IL-1β and TNFα. In addition, most of the patients with COVID-19 develop lymphocytopenia—especially reduced numbers of CD4+ and CD8+ T cells—within two weeks of disease onset, contributing to increased serum IL-6 levels. Lymphocytopenia is also predictive of the disease severity9,10,11, indicating that CD4+ and CD8+ T cells are critical for viral control and disease recovery. Moreover, single-cell transcriptomic analyses showed an increased ratio of plasma cells in patients infected with COVID-19 and early convalescent individuals7,12.

The epigenomic regulation of immune responses to primary SARS-CoV-2 infection remains unknown at present. Here we applied single-cell transposase-accessible chromatin with sequencing (scATAC-seq) to analyse chromatin remodelling in the peripheral immune cells of six matched and uninfected healthy donors (five donors collected in this study and an additional dataset from a public database13) and ten individuals convalescing from moderate or severe COVID-19 infection at 4–12 weeks following recovery. The approach of 10x-based scATAC-seq enables the production of high-quality single-cell accessible chromatin profiles at massive scale across all major immune-cell types in peripheral blood mononuclear cells (PBMCs). For the T-cell analysis, we developed a single-cell method that combines single-cell T cell-receptor (TCR) sequencing (scTCR-seq), fluorescence-activated cell sorting (FACS) with index sorting and ATAC-seq—termed TCR–FACS–index–ATAC sequencing (Ti-ATAC-seq)—thereby allowing us to obtain expression data for cell-surface markers, paired TCR sequences and the chromatin accessibility landscape for each cell analysed. This single-cell chromatin landscape approach facilitates a comprehensive understanding of how effector or memory cells are established in both the innate and adaptive immune responses of individuals convalescing from COVID-19.

Results

Overview of the peripheral immune-cell profiling in individuals convalescing from COVID-19

We mapped the transposase-accessible chromatin at the single-cell level using 10x chromium scATAC-seq (10x Genomics; Fig. 1a and Supplementary Table 1). We collected data from 72,318 PBMCs from eight individuals convalescing from COVID-19 as well as 24,997 PBMCs from five healthy donors and one additional dataset from a public database13 as controls (Fig. 1b and Extended Data Fig. 1a). Using uniform manifold approximation and projection (UMAP), we analysed the distribution of the immune-cell populations based on the gene scores of canonical lineage markers, which demonstrate the aggregate accessibility of several enhancers linked to the indicated genes—including CD14, CD16 (FCGR3A) and CEBPB for monocytes; PAX5 and MS4A1 for B cells; CD8A and CD3G for CD8+ T cells; CD4 and CD3G for CD4+ T cells; and KLRB1 for natural killer (NK) cells (Fig. 1c–e). The frequencies of monocyte-lineage and effector and memory CD8+ T cells increased in individuals convalescing from COVID-19, whereas those of the B-cell lineages decreased (Fig. 1f,g). We further analysed and clustered the cells of each lineage separately.

Fig. 1: Overview of the immune-cell epigenomic landscape of the blood of individuals convalescing from COVID-19.
figure 1

a, Outline of the two approaches used for scATAC-seq: 10x-based scATAC-seq (top) and Ti-ATAC-seq (bottom). b, Cell numbers (log10-transformed) in each sample for 10x-based scATAC-seq. c, UMAP plot showing 97,315 10x-based scATAC-seq profiles of immune cells in peripheral blood—including B, NK and T cells, monocytes and other clusters. The dots indicate individual cells and the cell-type identity is indicated by colour. d, UMAP plot showing the sample origin. e, Canonical markers overlaid on UMAP embedding, including CD14 and CEBPB for monocytes; PAX5 and MS4A1 for B cells; CD3G, CD8A and CD4 for T cells; and KLRB1 for NK cells. The UMAP plot is coloured based on the log-transformed normalized gene scores, which exhibited the accessibility of peaks linked to the indicated genes. The gene scores were calculated as log2(normalized count + 1). f, Cell-type frequencies in each sample, determined using 10x-based scATAC-seq. The colours indicate the cell type. g, Differences in the proportion of monocyte and B-cell lineages as well as effector and memory CD8+ T cells in the samples from individuals convalescing from COVID-19 (COV; n = 8 samples collected from eight individuals, one sample per individual) and healthy donors (HD; n = 4 samples, three samples collected from five HD in this study (including a pooled sample from three HD) and an additional sample from published data (GSE139369) as indicated in b and Extended Data Fig. 1a). A two-sided unpaired Student’s t-test was performed to determine the P values. The boxplots denote the median with the quartile range (25–75%), and the length of whiskers represents 1.5× the interquartile range (IQR). h, Number of T cells in each sample that had both scTCR-seq and ATAC-seq data.

Source data

In addition, we explored Ti-ATAC-seq to further investigate special clonally expanded T-cell clones (Fig. 1a and Extended Data Fig. 2). We collected 2,880 αβ T cells from ten individuals convalescing from COVID-19, among which 1,574 quality-control-positive T cells had both paired TCRαβ sequences and ATAC-seq (Fig. 1h). Ti-ATAC-seq had similar performance to previously published microfluidic platform-based T-ATAC-seq14 and other joint ATAC-RNA approaches15,16 (Extended Data Fig. 2c).

Sequential trained and activated monocytes in COVID-19

We identified and annotated 13 subclusters of monocytes (Fig. 2a). Clusters of CD14+ monocytes (clusters 11, 12 and 13) showed a high CD14 gene score. Clusters of CD16+ monocytes (clusters 3 and 4) were characterized by a high CD16 gene score. CD1c+ conventional dendritic cells (clusters 1 and 8) were defined by their high CD1c gene score. CLEC4C+ plasmacytoid dendritic cells showed high gene accessibility in the CLEC4C locus. Clusters 9 and 10, with their predominantly high KLF4 and CD68 gene scores, probably represented monocyte-derived dendritic cells (Fig. 2b and Extended Data Fig. 3a,b).

Fig. 2: Epigenomic signatures of trained and activated monocytes in individuals convalescing from COVID-19.
figure 2

a, Subclustering UMAP of all monocytes (see Fig. 1c). The five clusters of CD14+ and CD16+ monocytes indicated in the legend were annotated according to sequential differentiation states; cDC, conventional dendritic cells; pDC, plasmacytoid dendritic cells and moDC, monocyte-derived dendritic cells. b, Surface markers overlaid on UMAP embedding coloured according to the log-transformed normalized gene scores. c, Differences in the proportions of clusters 3, 4, 11 and 12 in the samples of individuals convalescing from COVID-19 (COV; n = 8 samples collected from eight individuals, one sample per individual) and healthy donors (HD; n = 4 samples, three samples collected from five HD in this study (including a pooled sample from three HD) and an additional sample from published data (GSE139369)). A two-sided unpaired Student’s t-test was performed to determine the P values. The boxplots denote the median with the quartile range (25–75%), and the length of whiskers represents 1.5× the IQR. dg, Volcano plots showing the differential TF motif accessibility using the mean TF motif accessibility in the chromVAR TF bias-corrected deviation in clusters 3 and 4 (d), 12 and 11 (e), 13 and 11 (f), and 13 and 12 (g). The P values were calculated using two-sided pairwise Wilcoxon test and the false discovery rate (FDR) was corrected using the Benjamini–Hochberg procedure.

Source data

We focused on CD14+ and CD16+ monocytes because the individuals convalescing from COVID-19 had a higher abundance of two clusters of CD14+ monocytes (clusters 11 and 12) and one cluster of CD16+ monocytes (cluster 4) than healthy individuals (Fig. 2c). We found and annotated sequential differentiation states from homeostasis to the mature immune inflammatory effector state in the COVID-19 group through differential cis- and trans-element analysis within the CD14+ and CD16+ clusters (Fig. 2a). For example, cluster 4 (CD16-activation) was associated with leukocyte activation and the regulation of cytokine production, and showed high levels of activity of transcription factors (TFs) involved in IFNγ-induced rapid signatures and myeloid differentiation, including TBET, JUN and FOSB17,18,19,20, whereas cluster 3 (CD16-maturation) was associated with haemopoiesis and the regulation of the immune response (Fig. 2d and Extended Data Fig. 3c). Interestingly, cluster 4 had increased chromatin accessibility in the CCL5, PRF1 and GZMB loci (Extended Data Fig. 3a), demonstrating that IFNγ drives the differentiation and activation of CD16+ monocytes in the hyperinflammatory microenvironment of patients with COVID-19. For the CD14+ monocytes, the clusters 12 and 13 (CD14-activation and CD14-maturation, respectively) showed increased chromatin accessibility at inflammatory cytokine genes—including IL1β, IL6, IL8, CCL2, CCL3 and CCL7 (Extended Data Fig. 3a)—and enriched TFs that represent a mature state—including FOS, JUN and MAF (Fig. 2e–g)19,21. Meanwhile, cluster 11 (CD14-trained) was characterized by the accessibility of genes involved in leukocyte activation (Extended Data Fig. 3c) and showed high activity of TFs involved in haematopoietic commitment and survival of monocytes, including HOXA9 and NR4A1 (Fig. 2e,f)22,23. In addition, NR4A1 TF activity, which was enriched in cluster 11 (Fig. 2e,f), was also enriched in the activation stage of the CD16+ monocyte (Fig. 2d).

Innate immune cells have recently been shown to display immunological memory after certain infections or vaccines, also termed trained immunity24,25,26. To understand the developmental dynamics of these monocytes, we constructed a lineage trajectory of CD14+ based on sequential differentiation states, which progressed from trained to activation and finally to maturation, and generated ordered single cells (termed as ‘pseudotime’) based on their epigenetic similarity (Fig. 3a). The dynamic cis-element and TF motif accessibility across the trajectory were consistent with the sequential differentiation states (Fig. 3b). For instance, the trained, activated and mature monocytes showed sequential chromatin accessibility in the IL1β and CCL5 loci (Fig. 3b,c). The trained stage of the CD14+ lineage trajectory was accompanied by increased accessibilities of IRF1, IRF3 and IRF8 TF motifs (Fig. 3b,d,e), thus enhancing susceptibility to regulate rapid IFN-β induction and host inflammatory defences in human blood monocytes20,27. The activation and maturation stages of these trajectories shared accessibility at TFs involved in monocyte activation and maturation for the activator protein 1 (AP-1) factors FOS and JUN (Fig. 3b,d,e)19,28,29.

Fig. 3: Epigenomic differentiation trajectory of CD14+ and CD16+ monocytes.
figure 3

a, UMAP showing the lineage trajectory of CD14+ monocytes ordered based on trained, activation and maturation states. Pseudotime values were overlaid on the UMAP embedding; the smoothed line and arrow represent the visualization of the trajectory path from the spline fit. b, Heatmaps of the ordered cis-element accessibility (left) and TF motif accessibility (right) across pseudotime in the CD14+ monocytes (see Fig. 3a). The cis-element and TF motif accessibilities are indicated by the gene score and chromVAR TF-motif bias-corrected deviation, respectively. c, Aggregated single-cell genome tracks for the indicated clusters at the IL1β (left) and CCL5 (right) loci with peak co-accessibility (Co-access). The Co-access is indicated by the inferred peak-to-gene links for distal regulatory elements. Green shading indicates differential peaks within clusters. d, Violin plots showing the ChromVAR TF-motif bias-corrected deviation scores of the indicated TF regulators across clusters 3 (n = 1,003 cells), 4 (n = 1,472 cells), 11 (n = 988 cells), 12 (n = 5,719 cells) and 13 (n = 4,537 cells). The boxplots denote the medians and the quartile range (25–75%), and the length of whiskers represents 1.5× the IQR. e, TF footprints of the TBET motif in clusters 3 and 4 (left) as well as the FOS (middle) and IRF1 (right) motifs in clusters 11, 12 and 13. The Tn5 insertion bias track is shown. f, Levels of IL-1β and IL-6 secreted by virus-exposed PBMCs (COV-S, n = 3 samples from three individuals that recovered from severe COVID-19; COV-M, n = 7 samples from seven individuals that recovered from moderate COVID-19; HD, n = 10 samples from ten healthy donors). The error bars indicate the mean ± s.e.m. A two-sided unpaired Student’s t-test was performed to determine the P values.

Source data

To validate this, we re-challenged PBMCs with spike-nCoV pseudovirus or hepatitis B virus for 24 h in vitro, as monocytes predominantly showed high gene scores for IL1β and IL6 (Fig. 1e). The PBMCs from individuals convalescing from COVID-19 secreted higher levels of IL-1β and IL-6 than those from the healthy donors (Fig. 3f), which is consistent with trained monocytes having enhanced sensitivity in response to different pathogens after an initial challenge25. There was no significant difference in the levels of cytokine secretion between the patients that were tested early (4 weeks) or late (10–12 weeks) following discharge (Fig. 3f).

Collectively, although systemic inflammation had returned to a normal and healthy level (Supplementary Table 1), the CD14+ and CD16+ monocytes maintained chromatin reprogramming that promotes trained immunity and thus enables a rapid inflammatory response against subsequent infection.

A facilitated B-cell developmental programme in individuals convalescing from COVID-19

We identified ten B-cell clusters: immature B cells (clusters 6 and 7) with a high SDC1 score; naive B cells (clusters 4 and 5) with high TCL1A, CD19 and CD20 scores; memory B cells (clusters 1, 2 and 3) with high CD27 and CD38 scores; and plasma cells (clusters 8, 9 and 10) with a high XBP1 score (Fig. 4a,b and Extended Data Fig. 4a).

Fig. 4: Accelerators facilitated the B-cell developmental programme in individuals convalescing from COVID-19.
figure 4

a, Subclustering UMAP plot of all B cells (see Fig. 1c). b, Gene markers overlaid on UMAP coloured according to the log-transformed normalized gene scores for SDC1, TCL1A, CD27 and XBP1. c, Relative frequencies of the B-cell subclusters in the different samples. d, Integrative lineage trajectory of B-cell states in healthy donors (HD) and individuals convalescing from COVID-19 (COV). The smoothed line and arrow represent the visualization of the trajectory path across different states (immature, naive, memory and plasma), and the sample origins are denoted by colour. e, Heatmaps showing the ordered gene-score trajectory across pseudotime for B-cell differentiation. f, Heatmaps showing the positive TF regulators obtained from the integration of ordered TF gene scores (right) with ordered TF motif accessibility (left) across pseudotime for B-cell differentiation. Positive TF regulators are TF motifs that show high bias-corrected chromVAR TF-motif deviations that also exhibit similarly dynamic gene scores across differentiation states. g, Volcano plots demonstrating the differential TF motif accessibility using the mean TF motif accessibility in the chromVAR TF bias-corrected deviation between COV and HD individuals in the indicated B-cell states. The P values were calculated using a two-sided pairwise Wilcoxon test and the false discovery rate (FDR) was corrected using the Benjamini–Hochberg procedure. h, Schematic of differential positive regulatory TFs driving B-cell differentiation of immature B cells to plasma cells in HD and COV individuals.

Source data

Although the individuals convalescing from COVID-19 and healthy donors had comparable frequencies of B-cell subsets (Fig. 4c), the UMAP projection revealed preferential cell distribution in each B-cell subset between healthy donors and individuals convalescing from COVID-19 (Fig. 4d and Extended Data Fig. 4b). The transcriptional regulation of the B-cell developmental programme in healthy donors has been extensively studied30,31. However, the B-cell lineage trajectories based on accessible chromatin differed between the healthy donors and individuals recovering from COVID-19 (Extended Data Fig. 4b–d), suggesting that SARS-CoV-2 infection induced different B-cell developmental programmes.

To investigate this, we reconstructed the lineage trajectory of B cells from the healthy donors and individuals convalescing from COVID-19 (Fig. 4d). We observed cis-elements near known regulators of every stage of B-cell development—such as TCL1A, CXCR4, MEF2C, BHLHE41, BAFF, PLCG2, PAX5, EBF1 and RORA (Fig. 4e and Extended Data Fig. 5a)30,31,32,33,34—suggesting that the integrated B-cell lineage trajectory was a well-defined developmental programme that could be used to compare regulatory mechanisms between the healthy donors and individuals convalescing from COVID-19. Next, through integrated chromVAR TF deviations with similarly dynamic gene scores across differentiation states, we identified positive TF regulators with sequential activities of BCL11A, IRF8, PAX5, REL, BATF, IRF4, EBF1, POU2F2, TBET and LEF1 that promote B-cell commitment, differentiation, maintenance and class-switch recombination30,35,36,37,38,39, thus resolving the integrated timing of TF activity for comparison (Fig. 4f and Extended Data Fig. 5b–d).

We then measured the TF deviation scores and variations to identify the differential positive TF regulators of four B-cell subsets of the healthy donors and individuals convalescing from COVID-19 (Fig. 4g). Interestingly, NF-kB subunits, including REL, RELA and RELB, which are involved in germinal centre B-cell maintenance and homeostasis40,41, were enriched in four subsets of B cells of the healthy donors (Fig. 4g and Extended Data Fig. 5d). Meanwhile, AP-1 factors, including FOS and JUN, which are involved in the B-cell receptor signalling pathway and indicate B-cell differentiation and activation42, were enriched in naive, memory and plasma B cells of the individuals convalescing from COVID-19. Transcription factors—including SPI1, EBF1, IRF4 and POU2F2—that are vital for B-cell survival, differentiation and receptor signalling responses were enriched in the naive B cells of the individuals convalescing from COVID-1936,38,43,44,45,46. In comparison to the memory and plasma B cells of healthy donors that showed a high REL deviation score, the plasma cells and a subset of memory B cells of the COVID-19 group showed high activity of TFs that are involved in class-switch recombination and promote specialized immune function in class-specific IgG+ memory B cells, including TBET and BATF39,47 (Fig. 4g and Extended Data Fig. 5d). We also detected specific IgG antibodies of the spike protein of SARS-CoV-2 in the blood samples of the individuals convalescing from COVID-19 (Extended Data Fig. 5e). Collectively, according to the results of integrated timing analysis of TF activity, differential positive TF regulators promoted B-cell maintenance and homeostasis in healthy donors, whereas they facilitated B-cell activation, differentiation and IgG class-switch recombination in the individuals convalescing from COVID-19 (Fig. 4h).

CD8+ T cell-fate decisions in individuals convalescing from COVID-19

We first sub-grouped NK and T cells into 14 subclusters (Fig. 5a and Extended Data Fig. 6a–c). The CD8+ T-cell states consisted of naive, intermediate, memory and effector T cells as well as mucosal-associated invariant T cells (MAIT), whereas the CD4+ T-cell states were composed of naive, central memory (TCM) and effector memory (TEM) T cells. In contrast to healthy donors enriched in naive CD8+ T cells, the individuals convalescing from COVID-19 showed increased numbers of effector and memory CD8+ T cells (Figs. 1h and 5b). Thus, clusters 6 and 8 may represent SARS-CoV-2-induced effector and memory CD8+ T cells, respectively.

Fig. 5: Single-cell epigenomic profiles of T cells in individuals convalescing from COVID-19.
figure 5

a, Subclustering UMAP of NK and T cells (see Fig. 1c). For the CD8+ T cells, naive cells had high TCF7 and CCR7 gene scores; effector cells had high CD8A, TBX21 (TBET) and IFNG gene scores; memory cells had high CD8A, TBX21 and KLF2/13 gene scores but low effector gene scores; and MAIT cells had high SLC4A10 gene scores. For the CD4+ T cells, naive cells had high TCF7 and CCR7 gene scores; TCM cells had high CD4, CCR7, AQP3 and SELL gene scores; and TEM cells had high CD4 and effector gene scores but low CCR7 gene scores. The NK cells had high NCAM, KLRC1, KLRD1 and FCGR3A gene scores; and NKT cells had high NCAM and CD3E gene scores. b, Cell-type frequencies in the different samples with the cluster identities indicated. c, Lineage trajectory of CD8+ T-cell states, which included naive, intermediate, effector and memory CD8+ T cells. Pseudotime values overlapped on the UMAP. d,e, Pseudotime heatmap showing the ordered cis-element accessibility (left) and TF motif accessibility (right) in effector (d) and memory (e) CD8+ T-cell lineage trajectories (see c). f, Aggregated single-cell genome tracks for the indicated clusters and states of CD8+ T cells at the GZMB (top) and IFNG (bottom) loci with peak co-accessibility (Co-access). The regions shaded in blue indicate differential peaks within clusters. g, Ridge plots showing TF deviation scores across different clusters and states of CD8+ T cells. h,i, Gene scores of IFNG (left) and GZMB (right; h), as well as the gene score (left) and TF deviation score of TBET (right; i) in the indicated cell types. COV, individuals convalescing from COVID-19; HD, healthy donors; COV-CD8-effector, n = 12,353 cells; COV-CD8-memory, n = 1,313 cells; COV-CD8-intermediate, n = 2,737 cells; COV-CD8-naive, n = 1,858 cells; HD-CD8-effector, n = 2,052 cells; HD-CD8-memory, n = 279 cells; HD-CD8-intermediate, n = 1,408 cells; HD-CD8-naive, n = 1,272 cells. The boxplots denote the medians and the quartile range (25–75%), and the length of whiskers represents 1.5× the IQR.

Source data

We next orchestrated the effector and memory CD8+ T-cell trajectories to explore CD8+ T cell-fate decisions. The effector and memory trajectories almost overlapped at the early state and diverged into two distinct branches, that were predominantly comprised of cells from the individuals convalescing from COVID-19 (Fig. 5c). The analysis of cis-elements near effector genes revealed distinct regulatory patterns of accessibility across the pseudotime of the effector and memory CD8+ T cells (Fig. 5d,e). For example, the cis-element accessibility in the promoter and distal enhancers in the GZMB and IFNG loci gradually increased from the naive to the intermediate and then to the memory and effector states (Fig. 5f), demonstrating that the expression of effector genes was accurately regulated by a distinct cis-element network in the different states.

The results of TF-activity analysis demonstrated the shared and unique TF programmes across the pseudotime of effector and memory CD8+ trajectories in the individuals convalescing from COVID-19 (Fig. 5d,e). For instance, in the first stage of trajectories from the naive to intermediate state, effector and memory CD8+ T cells shared accessibility at the AP-1 factors FOS and JUNB, consistent with their roles in T-cell differentiation and activation14,31; however, effector CD8+ T cells also showed accessibility at NFKB1/2 and RUNX3, whereas memory CD8+ T cells were characterized by accessibility of BATF. IRF4 activity gradually increased across the effector CD8+ T-cell trajectory but was depleted in that of memory CD8+ T cells (Fig. 5g and Extended Data Fig. 6d), which is consistent with their role in limiting the development of memory T cells48. In contrast, BATF motifs showed accessibility earlier in the trajectory of memory CD8+ T cells than that of effector cells (Fig. 5d,e). Similarly, the second stage of trajectories from the intermediate state to either the effector or memory state identified the pivotal roles of EOMES and TBET motifs in each pathway (Fig. 5d,e,g). Effector CD8+ T cells were also characterized by accessibility of SREBF1/2 (Fig. 5d,g), which is involved in the metabolic reprogramming of effector T cells during extensive clonal expansion49. Memory CD8+ T-cell commitment was accompanied by the accessibility of motifs involved in promoting memory T-cell trafficking and maintenance of Krüppel-like factor 2/13 (KLF2/13; Fig. 5e,g)50,51,52,53. Notably, compared with healthy donors, we also observed gradual increased gene accessibility of GZMB, IFNG and TBET in the intermediate, memory and effector CD8+ T cells of individuals convalescing from COVID-19 (Fig. 5h,i).

Furthermore, we performed TCR clonality analysis using the paired TCRαβ sequences generated from Ti-ATAC-seq. The clonal expansion rates of CD8+ T cells were significantly higher than those of CD4+ T cells of individuals convalescing from COVID-19 (Fig. 6a,b). Interestingly, CD8+ T cells were strongly dominated by one or two T-cell clones in each patient and the clonal expansion rate of the largest clone ranged from 2.6% to 41.1% of the total CD8+ T cells (Fig. 6a). These results suggest that these putative SARS-CoV-2-specific CD8+ T-cell clones, rather than CD4+ T cells, play a critical role in viral control and long-term immune protection. We next integrated TCR clonality with the single-cell epigenomic profiles generated from Ti-ATAC-seq. Notably, on comparison with unexpanded CD8+ T-cell clones, the two largest clonally expanded CD8+ T-cell clones were particularly enriched in TBET and EOMES TFs (Fig. 6c). This is consistent with our 10x-based scATAC-seq data showing that these two TFs are important for the development of effector and memory CD8+ T cells in the individuals convalescing from COVID-19 (Fig. 6d).

Fig. 6: Single-cell epigenomic and TCR profiling of CD8+ T cells in individuals convalescing from COVID-19.
figure 6

a, Distribution of CD4+ (top) and CD8+ (bottom) T-cell TCR clonotypes, according to size, from ten individuals convalescing from COVID-19 (COV). The full-length paired TCR sequences were obtained from the Ti-ATAC-seq platform-based scTCR-seq. For the expanded TCRs with clonotype size > 1, each slide indicates a unique clonotype. Single viable DAPICD3+TCRαβ+ T cells were gated and sorted, and the T-cell types were identified by FACS-indexed sorting. The CDR3α-CDR3β sequences and frequencies of the largest CD8+ T-cell clonotypes are shown. b, Clonal expansion rate of CD4+ and CD8+ T cells in ten individuals convalescing from COVID-19. Each data point indicates a single individual and the matched data of the CD4+ and CD8+ T cells are linked by distinct lines. A two-sided unpaired Student’s t-test was performed to determine the P value. c, Ranked TF deviation enrichment values in the aggregated largest (left; TCR clone: CLVGGGDNTDKLIF-CASSQDRFYEQYF, n = 50 cells) and second-largest clonally expanded CD8+ T cells (right; TCR clone: CAVGDTDKLIF-CASSLGSLGGGELFF, n = 24 cells) versus aggregated unexpanded CD8+ T cells (n = 298 cells). The T-cell clonality and single-cell epigenomic profiles were obtained from the Ti-ATAC-seq platform. The TF enrichment values were calculated as the difference in the mean TF deviation between the two single-cell populations. d, Anti-SARS-CoV-2 CD8+ T-cell response. Schematic of the differential regulatory TFs driving CD8+ T-cell differentiation from the naive state to the effector-cell or memory state in individuals convalescing from COVID-19.

Source data

In summary, our in-depth analysis of the epigenomic landscape and single-cell TCR clonality in individuals convalescing from COVID-19 revealed pivotal roles for effector CD8+ T cells in the initial viral control and the formation of memory CD8+ T cells via global chromatin accessibility remodelling and accurate regulatory programmes.

Discussion

SARS-CoV-2 causes severe pulmonary disease and complications with significant morbidity and mortality1. Our current understanding of the epigenomic regulatory mechanisms in the host immune response to SARS-CoV-2 infection and long-term immune protection is limited; thus, it is challenging to develop urgent therapeutics or assess the effect of different vaccine candidates.

Here we utilized high-throughput 10x-based scATAC-seq technology to measure the chromatin accessibility of all PBMCs, including monocytes and B, NK and T cells. We also used Ti-ATAC-seq to further investigate special clonally expanded T-cell clones with paired TCRαβ sequencing and responded chromatin accessibility from each T cell. Chromatin remodelling was significantly altered in almost all immune-cell compartments in the individuals convalescing from COVID-19. Our epigenomic profiles revealed a sequential differentiation state from homeostasis to the mature immune inflammatory effector response in CD14+ and CD16+ monocytes of individuals convalescing from COVID-19, which featured trained, activation and maturation states. The trained and activation states of CD14+ and CD16+ monocytes were dominantly enriched in the individuals convalescing from COVID-19. Although immune memory is a well-known feature of the acquired immune system, the activation of the innate immune system can also heighten responsiveness to subsequent triggers, which may provide protection during the early stage of reinfection. In B-lineage cells, we found substantial differences in the TF regulators of each state between the individuals convalescing from COVID-19 and healthy donors. Furthermore, B-cell lineage trajectories revealed an accelerated developmental programme from immune B cells to antibody-producing plasma cells in COVID-19. Integrated analysis of single-cell TCR clonality with the single-cell chromatin accessibility landscape showed the dramatic clonal expansion of CD8+ T cells in individuals convalescing from COVID-19 as well as the bifurcation of cell-fate decisions for the putative SARS-CoV-2-specific clonally proliferating CD8+ T cells. This is consistent with effector CD8+ T cells in initial viral control and memory CD8+ T cells in long-term immune protection.

Consistent with our data, the reported scRNA-seq profiles of PBMCs identified the abundant IL-1+CD14++ monocytes and activated monocytes in individuals convalescing from COVID-19 (ref. 7). However, our scATAC-seq profiles revealed further epigenetic changes in individuals convalescing from COVID-19, which can be attributable to the fact that epigenetic changes are more long-lived compared with transcriptional differences. Trained immunity caused by epigenetic and metabolic reprogramming is a ‘double-edged sword’. Although it offers broad benefits for the host immune defence against pathogens, it identifies potentially detrimental outcomes in immune-mediated and chronic inflammatory diseases24,25. In the context of COVID-19, persistent infection with SARS-CoV-2 results in long-term chromatin accessibility reprogramming, which permits monocytes to remain in a ‘trained’ functional state. Well-controlled trained immunity is likely to protect against subsequent infection25. The effect of cytokine storms on trained immunity should be studied in the future.

Previous studies of immune responses to SARS-CoV-1 and Middle East respiratory syndrome coronavirus (MERS-CoV) can partially offer insights into SARS-CoV-2 memory immunity. Similar to SARS-CoV-2, SARS-CoV-1 and MERS-CoV induce strong inflammatory responses and associated lymphopenia9. Furthermore, long-term immune protection has been observed in patients that have recovered from SARS-CoV-1 and MERS-CoV and is attributed to virus-specific antibody responses and long-lived virus-specific memory T-cell responses2,54,55, which is consistent with our observations on the chromatin reprogramming of B and CD8+ T cells. Recent studies on memory immune protection against SARS-CoV-2 reinfection, including a rhesus macaque reinfection model and COVID-19 vaccines under phase I, II and III clinical trials, suggest that primary infection with SARS-CoV-2 and vaccines against SARS-CoV-2 can provide partial immune protection56,57,58,59.

The mechanisms of memory immunity—including that in trained and activated monocytes, accelerated B cells and putative SARS-CoV-2-specific clonally expanded effector and memory CD8+ T cells—to SARS-CoV-2 may be induced by cytokine storms, such as the high production of IFNγ in the serum. As we observed, IFNγ-induced TBET was significantly enriched in the activation state of trained monocytes, B-cell acceleration and CD8+ T cell-fate decisions, which is consistent with its role in monocyte activation and maturation, IgG+ memory B-cell development and effector and memory T-cell differentiation17,20,24,31,39,60,61,62. Therefore, the epigenomic regulation of innate and adaptive immune memory responses we demonstrated may not be specific for SARS-CoV-2 and could be elicited following other infections, such as SARS-CoV-1 and MERS. Future studies comparing epigenetic changes among convalescing individuals infected with SARS-CoV-2 and other viruses would be valuable.

Overall, our broad analysis of the epigenomic landscape and TCR profiling demonstrated that individuals convalescing from COVID-19 established immune memory formation and trained immunity via the global remodelling of the chromatin accessibility landscape. The stability of these changes over a longer period requires further study.

Methods

Subjects and specimen collection

We collected blood from five healthy donors, seven convalescing individuals who had recovered from moderate COVID-19 and three convalescing individuals who had recovered from severe COVID-19. They were enrolled at the Fifth Medical Center of PLA General Hospital in May 2020. We isolated PBMCs from the study participants using Ficoll solution according to the manufacturer’s instructions. We defined individuals convalescing from COVID-19 as patients who had been diagnosed as infected with COVID-19 at admission by real-time RT–PCR, and had recovered and been discharged from the hospital at least one month before participating in the study. The individuals convalescing from COVID-19 were classified into mild and severe groups based on the Fifth Revised Trial Version of the Novel Coronavirus Pneumonia Diagnosis and Treatment Guidance. The healthy donors—three women and two men, 18–55 years old—were defined by negative test results in real-time RT–PCR analysis of the SARS-CoV-2 gene and IgG antibody detection. This study and study protocol was approved by the Ethics Committees of the Fifth Medical Center of Chinese PLA General Hospital, Beijing, China (2020005D) on healthy volunteers or individuals convalescing from COVID-19 in Medical Research and written informed consent was obtained from all of the participants. All of the participants received participant compensation. The clinical features of these individuals convalescing from COVID-19 and healthy donors are listed in Supplementary Table 1.

Single T-cell index sorting

Human PBMCs were thawed after storage in liquid nitrogen and allowed to recover overnight in R10 culture medium (RPMI1640 medium supplemented with 10% fetal calf serum, 25 mM HEPES, 1×non-essential amino acids, 50 μM β-mercaptoethanol and 1×penicillin and streptomycin). Cell surface staining was performed with 4,6-diamidino-2-phenylindole (DAPI) and the following antibodies: Alexa Fluor 700 anti-human CD3 (clone SK7, 1:100; BioLegend, cat. no. 344822), FITC anti-human TCRαβ (clone IP26, 1:50; BioLegend, cat. no. 306706), PerCP/Cyanine5.5 anti-human CD8 (clone SK1, 1:200; BioLegend, cat. no. 344710), PE/Dazzle 594 anti-human CD4 (clone A161A1, 1:200; BioLegend, cat. no. 357411), PE anti-human CD25 (clone BC96, 1:200; BioLegend, cat. no. 302606), Brilliant Violet 650 anti-human CD127 (clone A019D5, 1:200; BioLegend, cat. no. 351326), APC anti-human CD137 (4-1BB) (clone 4B4-1, 1:100; BioLegend, cat. no. 309810) and Brilliant Violet 605 anti-human CD279 (PD-1) (clone EH12.2H7, 1:100; BioLegend, cat. no. 329923). Single viable DAPICD3+TCRαβ+ T cells were immediately index sorted into 96-well plates for Ti-ATAC-seq analysis using a BD FACS Aria III system.

Preparation of 10x-based scATAC-seq libraries

After storage in liquid nitrogen, PBMCs were thawed, allowed to recover overnight in R10 culture medium and washed three times with PBS buffer to remove debris and aggregated cells. The cell viability of each sample exceeded 90%. For the healthy donor controls, we performed 10x-based scATAC-seq on PBMCs from five healthy donors, including a sample of mixed PBMCs from three healthy donors. The cells were lysed, the nuclei were isolated, the nuclei suspensions were washed and counted, and the transposition reaction and nuclei barcoding were performed according to the manufacturer’s instructions. Approximately 5,000–10,000 nuclei were collected in each sample. The protocol for library preparation and the settings for the instrument and sequencing were provided by the manufacturer (Annoroad and Berry Genomics) and are available at https://support.10xgenomics.com/single-cell-atac. We also collected scATAC-seq data from a healthy donor control from a published dataset (GSE139369)13.

Preparation of single-cell Ti-ATAC-seq libraries

Single antibody-indexed cells were sorted into 96-well plates with 6 μl ATAC-RSB buffer (10 mM Tris–HCl at pH 7.5, 10 mM NaCl, 3 mM MgCl2 and 0.8 U μl−1 RNase inhibitor) per well. Next, 1.5 μl transposition reaction mix with 0.1 μl TTE Mix V5 (Vazyme Biotech), 50 mM TAPS–NaOH (pH 8.5), 25 mM MgCl2, 50% dimethylformamide, 0.8 U μl−1 RNase inhibitor and 0.2% NP40 was carefully added to each well. The transposition reaction was performed by incubating the mixture at 37 °C for 30 min and then stopped by the addition of 1.5 μl release buffer (125 mM EDTA and 10 mM Tris–HCl at pH 8.0), following which the reaction mixture was incubated at 50 °C for 30 min. Quenching buffer (1 μl; 187.5 mM MgCl2 and 10 mM Tris–HCl at pH 8.0) was then added for MgCl2 quenching. The TRA and TRB chains and transposed DNA fragments were reverse transcribed with 11 μl of the RT–PCR reaction mix (8.8 μl of 2×1 Step buffer, 1.5 μl of TRA and TRB RT primer mix, 0.2 μl of non-indexed custom Nextera ATAC-seq PCR primer mix (forward, 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-3′ and reverse, 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-3′) and 0.5 μl PrimeScript 1 Step enzyme mix (Takara Bio)) under the following RT–PCR conditions: 50 °C for 36 min; 95 °C for 15 min; and eight cycles at 94 °C for 30 s, 62 °C for 1 min and 72 °C for 1 min. The RT–PCR product was then immediately equally divided into two 96-well plates for amplification and the barcoding of transposed DNA fragments and TRA and TRB chains, which were performed as previously described.63 Briefly, 1.25 μM of dual-index Nextera ATAC-seq PCR primers with 1×TruePrep Amplify reaction mix (Vazyme Biotech) was prepared, and ATAC-seq amplification and barcoding were performed under the following PCR conditions: 72 °C for 5 min, 98 °C for 30 s, and 15 cycles of 98 °C for 15 s and 63 °C for 30 s. Finally, scTCR-seq and scATAC-seq libraries were respectively pooled and purified on a single Zymo DNA clean and concentrator 5 column. Fragment-length selection was performed (0.5×/1.5× for scATAC-seq and 300–400 bp for scTCR-seq). The scTCR-seq and scATAC-seq libraries were 150-bp paired-end sequenced on a HiSeq X Ten platform (Illumina) with the sequencing read length and dual indexing according to the manufacturer’s instructions (Annoroad).

Data processing of scTCR-seq libraries

The scTCR-seq data were analysed as previously described63,64. Briefly, raw reads were unpacked and joined by shared regions and demultiplexed using a custom pipeline. The TCR V, D and J segments were analysed and assigned using MiXCR65.

Primary data pre-processing of scATAC-seq libraries

We used a previously described workflow with minor modifications14. Briefly, a sample name was added to the header of each read, adaptor sequences were trimmed using Cutadapt, reads were mapped to the hg38 human genome with Bowtie2 and those aligned to the mitochondria were removed. The resulting cleaned single-cell files (./bam files) of each individual convalescing from COVID-19 were merged and sorted by read name using SAMtools. Finally, the resulting merged file of each individual convalescing from COVID-19 was converted to a fragment file, the Tn5 insertion site was adjusted and unique nuclear fragments were retained using BEDtools. We then calculated the TF deviation of cells using chromVAR66.

For the pre-processing of the 10x-based scATAC-seq data, we used the ‘cellranger-atac count’ function (cellranger-atac, v1.2.0) to generate single-cell accessibility counts for each library (https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/using/count). Reads were aligned to the hg38 human genome. Similarly, we created an Arrow file for each resulting fragment file and included those single-cell libraries with at least 1,000 unique fragments and a transcription-start-site enrichment of eight using the ArchR package67 in the R statistical environment (v3.5.1). Finally, we created an ArchRProject for downstream analysis by combining all of the Arrow files of the 10x-based scATAC-seq.

Dimensionality reduction and clustering analysis

We used the ‘addIterativeLSI’ function of ArchR to perform iterative latent semantic indexing31,67. We then used the harmony algorithm to correct for batch-effect differences68 and added clusters using the ‘addClusters’ function. We ran UMAP using the ‘addUMAP’ function and plotted the results using the ‘plotEmbedding’ function in ArchR67.

Identification of marker features and differential analysis for clusters

We created pseudo-bulk replicates using the ‘addGroupCoverages’ function for each cluster and called peaks using the ‘findMacs2’ function. We used the ‘addDeviationsMatrix’ function to compute per-cell deviations across all motif annotations to create the deviation matrix ‘MotifMatrix’. Next, we used the ‘getMarkerFeatures’ function with ‘bias’ parameter to account for transcription-start-site enrichment and the number of unique fragments per cell to identify the markers of clusters, including peaks, genes (based on gene scores) and TF motifs (based on chromVAR deviations). We applied the ‘getMarkers’ function to obtain the marker list of each cluster, the ‘addImputeWeights’ function to impute the weights of marker features and the ‘plotEmbedding’ function to visualize the marker features. We performed pairwise comparisons of peaks, genes and TFs for the indicated two clusters using the getMarkerFeatures function and plotted an MA or volcano plot using the ‘markerPlot’ function in ArchR67.

Plotting browser tracks of cis-element co-accessibility

Co-accessibility represents the correlation between accessibility peaks across many single cells and is useful for identifying cell type-specific peaks. We used the ‘addCoAccessibility’ function to examine co-accessibility in ArchR67, which returned a loop track that represented the co-accessibility information using the ‘getCoAccessibility’ function. Finally, we plotted the genome browser tracks of peaks and co-accessibility using the ‘plotBrowserTrack’ function.

TF footprinting

We performed TF footprinting to precisely predict the binding location of a TF at a particular locus as previously described31,67 using ArchR. The Tn5 bias signal was normalized using the ‘Divide’ strategy.

Trajectory construction

We performed cellular trajectory analyses as previously described using ArchR67. First, we defined the trajectory backbone of cell groups or clusters based on the cell differentiation states, such as the order of naive, intermediate and effector CD8+ cells in the effector CD8+ T-cell trajectory. We then created a trajectory using the ‘addTrajectory’ function and plotted the pseudotime values on UMAP embedding using the ‘plotTrajectory’ function. We next plotted pseudotime heatmaps of TFs, gene scores and peak accessibility using the ‘plotTrajectoryHeatmap’ function. To identify positive TF regulators, we performed an integrative analysis of gene scores and motif accessibility across pseudotime using the ‘correlateTrajectories’ and plotTrajectoryHeatmap functions.

Statistics and reproducibility

Statistical analyses were performed in GraphPad Prism (v6) or R (version 3.5.1). For the frequencies of the identified clusters and the virus rechallenge assay, a two-sided unpaired Student’s t-test was performed and the P values are reported. For the differential TF motif accessibility, differential cis-element accessibility and Gene Ontology analyses, multiple-test corrections were performed, the P values were calculated using a two-sided pairwise Wilcoxon test and the false discovery rate was corrected using the Benjamini–Hochberg procedure. No data were excluded from the analyses.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.