Participants were enrolled in LIINC (www.liincstudy.org; NCT04362150)7, a prospective observational study enrolling individuals with prior nucleic acid-confirmed SARS-CoV-2 infection, regardless of the presence or absence of postacute symptoms. At each study visit, participants underwent an interviewer-administered assessment of 32 physical symptoms that were newly developed or had worsened since the COVID-19 diagnosis. Detailed data regarding medical history, COVID-19 history, SARS-CoV-2 vaccination and SARS-CoV-2 reinfection were collected. Two participants had biospecimens collected via the COVID-19 Host Immune Response Pathogenesis (CHIRP) study5. For the present study, we selected participants who consistently met a case definition for LC based on the presence or absence of at least one symptom attributable to COVID-19 for the 8 months following SARS-CoV-2 infection (Fig. 1a). The LC group (n = 27) had a median age of 46 years, and was comprised of 63% females and 26% of whom were previously hospitalized for COVID-19. The R group (n = 16) had a median age of 45.5 years, and was comprised of 44% females and 12.5% of whom were previously hospitalized for COVID-19 (Supplementary Table 1). Participants were deliberately not matched by age and sex, but we ensured that there was overlap in the groups. Blood samples were collected between September 16, 2020 and April 6, 2021. All participants provided a post-COVID blood sample before a SARS-CoV-2 vaccination to exclude the potential effects of SARS-CoV-2 vaccination on our study. Specimens were collected 8 months postinfection from individuals. All assays were performed from the same parent set of n = 27 LC and n = 16 specimens. All participants provided written informed consent.
Whole blood was collected in EDTA tubes followed by isolation of PBMCs and plasma as described in ref. 15. Serum was obtained concomitantly from serum-separator tubes.
Antibody responses against SARS-CoV-2 spike RBD were measured on sera using the Pylon COVID-19 total antibody assay (ET Health) and reported as relative fluorescence units (RFUs).
Peptides used for T cell stimulation comprised a mix of overlapping 15-mers spanning the entire SARS-CoV-2 spike protein (PM-WCPV-S-1, purchased from JPT), and peptides corresponding to CD8+ T cell epitopes identified by T-scan16 synthesized in-house (Supplementary Table 12). Final peptide concentrations were 300 nM for the 15-mers and 450 nM for the T-scan peptides.
Sample preparation was performed similar to methods described2,3,4,5. Upon revival of cryopreserved PBMCs, cells were rested overnight to allow for antigen recovery17 and then divided equally into two aliquots. To the first aliquot, we added 3 µg ml−1 brefeldin A (BFA; to enable intracellular cytokine detection), the costimulation agonists anti-CD28 (2 µg ml−1; BD Biosciences) and anti-CD49d (1 µg ml−1; BD Biosciences), and the SARS-CoV-2 peptide pool prepared as described above. To the second aliquot, we added 1% DMSO (Sigma-Aldrich) and 3 µg ml−1 BFA. Cells from both treatments were incubated at 37°C for 6 h. Cells were treated with cisplatin (Sigma-Aldrich) as a live/dead distinguisher and fixed in paraformaldehyde (Electron Microscopy Sciences) as described2,3,4,5. CyTOF antibody conjugation was performed using the Maxpar X8 Antibody Labeling Kit (Standard BioTools) according to the manufacturer’s instructions. CyTOF staining was performed as described2,3,4,5, but using the CyTOF panel created for this study (Supplementary Table 4). Stained samples were washed with CAS buffer (Standard BioTools), spiked with 10% (vol/vol) EQ Four Element Calibration Beads (Standard BioTools) and run on a Helios CyTOF instrument (UCSF Parnassus Flow Core).
CyTOF data analyses
EQ bead-normalized CyTOF datasets were concatenated, de-barcoded and normalized using Standard BioTools Software version 6.7. Following arcsinh transformation of the data18, cells were analyzed by FlowJo (version 10.8.1, BD Biosciences). Intact (Ir191+Ir193+), live (Pt195−), singlet events were identified, followed by gating on CD3+ T cells, and sub-gating on CD4+ T cells and CD8+ T cells (Extended Data Fig. 1g,h).
CyTOF antibody validation
CyTOF antibodies in our panel (Supplementary Table 4) were validated using methods previously described, including the use of human lymphoid aggregate cultures generated from tonsils2,3,4,5,18,19. The observed expression patterns among tonsillar T and B cells (Extended Data Fig. 10a) were similar to those previously observed18. To validate the detection of cytokines and other effectors, we stimulated PBMCs with 16 nM phorbol 12-myristate 13-acetate (PMA) (Sigma-Aldrich) and 1 μM ionomycin (Sigma-Aldrich), or 1 μg ml−1 lipopolysaccharides (LPS; eBioscience), for 4 h in the presence of 3 μg ml−1 BFA solution (eBioscience), combined the cells and prepared them for CyTOF as described above. We observed the expected induction of cytokines or cytolytic markers (Extended Data Fig. 10b)2,3,4,5 and preferential expression of Treg lineage marker Foxp3 among CD3+CD4+CD45RO+CD45RA−CD127−CD25+ Treg cells (Extended Data Fig. 10c). We also observed preferential expression of CD30 and Ki67 in CD4+ TM as compared to CD4+ TN cells (Extended Data Fig. 10d). Examples of pTFH and TFH gates are depicted in Extended Data Fig. 10e.
Identification of SARS-CoV-2-specific T cells
For identification of SARS-CoV-2-specific T cells, we compared unstimulated specimens to their peptide-stimulated counterparts. Effector cytokines (IFN-γ, TNF, IL-2, IL-4, IL-6, IL-17 and CCL4) and cytolytic effectors (granzyme B and perforin) were assessed for the ability to identify antigen-specific T cells at the single-cell level. The following criteria were established to identify effector molecules appropriate for identifying SARS-CoV-2-specific T cells: (1) counts of positive cells in unstimulated sample (not receiving peptide) was less than 5 events, or the frequency of positive cells was lower than 0.1%; (2) counts of positive cells in the peptide-stimulated sample was not less than 5, or the frequency was higher than 0.1%; (3) differences in frequencies of positive cells between unstimulated and peptide-stimulated samples cells was not less than 0.01%; (4) fold change in frequencies of positive cells between unstimulated and peptide-stimulated samples cells was greater than 10 and (5) the aforementioned four criteria could identify SARS-CoV-2-specific T cells among >50% of participants. Effectors that fulfilled all five criteria were IFN-γ, TNF and IL-2 for CD4+ T cells and IFN-γ, TNF and CCL4 for CD8+ T cells. For a sub-analysis to identify responding cells that may only exist in a small subset of individuals, we removed criterion 5 and reduced the positive cell counts to number 3 within criteria 1 and 2. This approach allowed us to determine that SARS-CoV-2-specific CD4+ T cells producing IL-6 were exclusively detected from LC (Extended Data Fig. 2f). SARS-CoV-2-specific T cells were detected at a median of 163 cells (134 for CD4+ T cells and 29 for CD8+ T cells) and a mean of 221.7 cells (185.2 for CD4+ T cells and 36.4 for CD8+ T cells), per participant. SARS-CoV-2-specific T cells, once identified, were analyzed by Boolean gating20 and exported for further analyses.
SPICE analyses were performed using version 6.1 software21. CD4+ and CD8+ T cells were subjected to manual gating based on the expression of cytokines used to define SARS-CoV-2-specific T cells (IFN-γ, TNF, IL-2 and CCL4, see above) using operations of Boolean logic. The parameters for running the dataset were as follows: iterations for permutation test = 10,000 and highlight values = 0.05. The parameters for the query structure were set as follows: values = frequency of single cytokine positive cells in total CD4+/CD8+ T cells; category = IFN-γ, TNF, IL-2 and CCL4; overlay = patient type (LC versus non-LC); group = all other variables in the data matrix.
T cell subsetting
Manual gating was performed using R (version 4.1.3). Arcsinh-transformed data corresponding to total or SARS-CoV-2-specific CD4+ or CD8+ T cells were plotted as 2D plots using the CytoExploreR package. Visualization of datasets by t-SNE was performed using methods similar to those described2,3,4,5. CytoExploreR and tidyr packages were used to load the data, and t-SNE was performed using Rtsne and RColorBrewer packages on arcsinh-transformed markers. Total CD4+/CD8+ T cells were downsampled to n = 8,000 (maximal cell number for individual samples) before t-SNE analysis. The parameters for t-SNE were set as iteration = 1,000, perplexity = 30 and θ = 0.5.
T cell clustering analysis
Flow cytometry standard (FCS) files corresponding to total and SARS-CoV-2-specific CD4+ and CD8+ T cells were imported in R for data transformation. Packages of flowcore, expss, class and openxlsx were loaded in R. Arcsinh-transformed data were then exported as CSV files for clustering analyses. Biological (LC status, biological sex and hospitalization status) and technical (batch/run of processing) variables were visualized using the DimPlot function of Seurat22. Batch correction was performed by RunHarmony23. Optimal clustering resolution parameters were determined using Random Forests24 and a silhouette score-based assessment of clustering validity and subject-wise cross-validation, as detailed in ref. 25. A generalized linear mixed model (GLMM, implemented in the lme4 (ref. 26) package in R with family argument set to the binomial probability distribution) was used to estimate the association between cluster membership and LC status and the sex of the participant, with the participant modeled as a random effect. For each individual, cluster membership of cells was encoded as a pair of numbers representing the number of cells in the cluster and the number of cells not in the cluster. Clusters having fewer than three cells were discarded. The sex-specific log odds ratio of cluster membership association with LC status was estimated using the emmeans27 R package using the GLMM model fit. The estimated log odds ratio represented the change (due to LC status) in the average over all participants of a given sex in the log odds of cluster membership. The two-sided P values corresponding to the null hypothesis of an odds ratio value of 1 were computed based on a z statistic in the GLMM model fit. These P values were adjusted for multiple testing using the Benjamini–Hochberg method.
Flow cytometry was performed on PBMCs from 25 LC and 15 R individuals from our cohort, obtained from aliquots of specimens analyzed by CyTOF. Cells were stained with the panel shown in Supplementary Table 13, using Zombie UV or Zombie NIR (BioLegend) as viability indicators. All cells were analyzed on a Fortessa X-20 (BD Biosciences). FCS files were exported into FlowJo (BD, version 10.9.0) for further analysis. Flow cytometric data were arcsinh-scaled before analyses. In flow cytometric experiments, SARS-CoV-2-specific CD8+ T cells were defined as those specifically inducing IFN-γ and/or TNF in response to SARS-CoV-2 peptide stimulation, as the CCL4 antibody exhibited background staining in flow cytometry and could not be used to define SARS-CoV-2-specific T cells.
RNA-seq was performed on PBMCs from 23 LC and 13 R individuals from our cohort, obtained from aliquots of specimens analyzed by CyTOF. Samples were prepared using the AllPrep kit (Qiagen) per the manufacturer’s instructions. RNA libraries, next-generation Illumina sequencing, quality control analysis, trimming and alignment were performed by Genewiz (Azenta). Briefly, following oligo dT enrichment, fragmentation and random priming, cDNA syntheses were completed. End repair, 5′ phosphorylation and dA-tailing were performed, followed by adaptor ligation, PCR enrichment and sequencing on an Illumina HiSeq platform using PE150 (paired-end sequencing, 150 bp for reads 1 and 2). Raw reads (480 Gb in total) were trimmed using Trimmomatic (version 0.36) to remove adapter sequences and poor-quality reads. Trimmed reads were mapped to Homo sapiens GRCh37 using star aligner (version 2.5.2b)28. log2 fold changes were calculated between LC versus R individuals. Two-sided P values corresponding to a null hypothesis of fold change of 1 were calculated using DESeq2’s (ref.29) Wald test and were adjusted for multiple testing using false discovery rates. Genes with an adjusted P value < 0.05 and absolute log2(fold change) > 1 were considered significant DEGs. Clustered heatmaps of DEGs were constructed with groups of genes (rows) defined using the k-means algorithm to cluster genes into k clusters based on their similarity. K = 4 was determined using the Hierarchical Ordered Partitioning and Collapsing Hybrid (HOPACH) algorithm30, which recursively partitions a hierarchical tree while ordering and collapsing clusters at each level to identify the level of the tree with maximally homogeneous clusters.
scRNA-seq was performed on PBMCs from 8 LC and 4 R individuals from our cohort, obtained from aliquots of specimens analyzed by CyTOF. Library preparation was performed using the Chromium Next GEM Single-Cell 5′ Reagent Kits v2 (10x Genomics) and sequenced on the Illumina NovaSeq 6000 S4 300 platform. Samples were sequenced at a mean of >50k reads per cell (minimum 51k, maximum 120k and median 83k). A median of 7,888 cells was analyzed per donor (minimum 4,189 and maximum 9,511). Demultiplexed fastq files were aligned to human reference genome GRCh38 using the 10x Genomics Cell Ranger v7.1.0 count pipeline31. The include-introns flag for the count pipeline was set to true to count reads mapping to intronic regions. The filtered count matrices generated by the Cell Ranger count pipeline were processed using Seurat22. Each sample was preprocessed as a Seurat object, and the top 1% of cells per sample with the highest numbers of unique genes, cells with ≤200 unique genes and cells ≥10% mitochondrial genes were filtered out for each sample. The samples were then merged into a single Seurat object, and normalization and variance stabilization were performed using sctransform86 with the ‘glmGamPoi’ method32 for initial parameter estimation.
Graph-based clustering was performed using the Seurat22 functions FindNeighbors and FindClusters. First, the cells were embedded in a k-nearest neighbor graph (with k = 20) based on the Euclidean distance in the principal component analysis (PCA) space. The edge weights between the two cells were further modified using Jaccard similarity. Next, clustering was performed using the Louvain algorithm33 implementation in the FindClusters Seurat function. Clustering with 15 principal components (PCs, determined based on the location of the elbow in the plot of variance explained by each of the top 25 PCs) and 0.1 resolution (determined using the resolution optimization method described above for CyTOF data clustering) resulted in 11 distinct biologically relevant clusters (clusters 0–11), which were used for further analyses. Marker genes for each cluster were identified using the FindAllMarkers Seurat function. Marker genes were filtered to keep only expressed genes detected in at least 25% of the cells, with at least 0.5 log2 fold change. Cluster annotation was performed according to subset definitions previously established34,35,36. Classification markers included CD19, MS4A1 and CD79A for B cells; CD3D, CD3E, CD5 and IL7R for CD4+ T cells; CD3D, CD3E, CD8A, CD8B and GZMK (CTL subset) for CD8+ T cells; CD14, CD68, CYBB, S100A8, S100A9, S100A12 and LYZ for monocytes; CSF2RA, LYZ, CXCL8 and CD63 for granulocytes and PF4, CAVIN2, PPBP, GNG11 and CLU for platelets.
The counts-per-million reads for ALAS2 and OR7D2 were assessed using edgeR37, and associations with group status were made using the two-sample Welch t test, followed by multiple correction testing using the Holm38 procedure. For establishing associations between clusters and group status, GLMM implemented in the lme4 R package was used. The model was performed with the family argument set to the binomial probability distribution and with the ‘nAGQ’ parameter set to 10 corresponding to the number of points per axis for evaluating the adaptive Gauss–Hermite approximation for the log-likelihood estimation. Cluster membership was modeled as a response variable by a two-dimensional vector representing the number of cells from a given sample belonging or not to the cluster under consideration. The corresponding sample from which the cell was derived was the random effect variable, and the group (R, LC, OR7D2high LC, or ALAS2high LC) was considered the fixed variable. The log odds ratio for all pairwise comparisons was estimated using the model fits provided to the emmeans function in the emmeans R package27. The resulting P values for the estimated log odds ratio and clusters were adjusted for multiple testing using the Benjamini–Hochberg method39. For associations of gene expression with group status, raw gene counts per cell were loaded as a SingleCellExperiment object. Cells from clusters 9 and 10 were not included in this analysis as the median number of cells across samples was less than 20 per cluster. The aggregateData function in the muscat bioconductor package40 was used to pseudo-bulk the gene read counts across cells for each cluster group. Genes with raw counts less than ten in more than eight samples were removed from the analyses. The pbDS function implementing the statistical methods in the edgeR package37 was used to assess associations of gene expression with group identity. Results from the cluster-specific pseudo-bulked gene expression association analyses were visualized as volcano plots using EnhancedVolcano41,42. Select genes of interest or genes that passed a multiple testing-adjusted P value threshold of 0.05 or 0.1 as indicated were indicated in the volcano plots. For gene set enrichment analyses, the raw P values for each gene derived from hypothesis tests for associations of interest were combined with a list of genes annotated with each of the gene sets in the biological processes domain of GO43 and analyzed via the simultaneous enrichment analysis method44 using the rSEA R package45. The family-wise error rate-adjusted P values for cluster-specific associations of interest with each of the annotated gene sets were used to identify significant associations.
The Olink EXPLORE 384 inflammation protein extension assay was performed per manufacturer’s protocol as published in ref. 46.
HOPACH30 was used to find the best cluster number. Gene expression values were log-transformed and centered using the average expression value. Clustering was performed by running the k-means algorithm using the best cluster number k found, and the results were plotted using the pheatmap package47. For gene network analyses, the STRING interaction database was used to reconstruct gene networks using stringApp48 for Cytoscape49. For the network, the top 50 genes or 25 proteins with the lowest P values were selected from the RNA-seq data and Olink data, respectively. They were then subjected to stringApp with an interaction score cutoff = 0.5 and the number of maximum additional indirect interactors cutoff = 10.
Unless otherwise indicated, permutation tests, two-tailed unpaired Student’s t tests and Welch’s t test were used for statistical analyses. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001 and NS. Error bars corresponded to s.d. Graphs were plotted by GraphPad Prism (version 9.4.1). All measurements were taken from distinct samples, no samples were measured repeatedly to generate data. Where appropriate, P values were corrected for multiple testing (across three pairwise comparisons) using the Holm procedure38. Tests involving cluster membership differences assumed a binomial probability distribution, and those involving RNA expression differences assumed a negative binomial probability distribution, but these were not formally tested. All other tests were based on the normality assumption but this was not formally tested.
Statistics and reproducibility
No statistical method was used to predetermine the sample size. Samples were chosen based on the availability of specimens meeting our LC criteria. No samples were excluded from the analyses. Randomization was not implemented as the study compared LC to R individuals. Data collection and analysis were not performed blind to the conditions of the experiments.
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.