Long COVID involves activation of proinflammatory and immune exhaustion pathways

Assessment of SARS-CoV-2 neutralizing antibody (nAb) responses using luciferase (LUC)-based pseudovirus neutralization assays and T cell responses using pooled peptide interferon-γ (IFNγ) enzyme-linked immunospot (ELISpot) assays against SARS-CoV-2 WA1/2020, Delta and Omicron BA.1 in the CC and the LC groups did not detect differences in SARS-CoV-2 nAb titers and identified a trend toward higher Spike-specific IFNγ ELISpot responses in the LC compared with the CC group (Fig. 1b and Extended Data Fig. 1a–c). We did not detect plasma SARS-CoV-2 viral loads in any CCs or individuals with LC using PCR with reverse transcription genomic or subgenomic viral load assays (Fig. 1c).

Bulk RNA-seq in PBMCs from individuals with LC (n = 26) and CCs (n = 21) at day 90–180 and in uninfected (n = 35) and acute (n = 54) individuals identified reads that mapped to the human genome and multiple common viruses (SARS-CoV-2, varicella zoster virus, Western equine encephalitis virus, Epstein–Barr virus, human cytomegalovirus, herpes simplex virus 1 and 2, influenza virus, rotavirus). We did not detect significant differences in viral read counts for these common viruses in the LC compared with the CC group (Extended Data Fig. 1d). Unsupervised clustering of bulk RNA-seq transcriptomic data revealed a distinct stratification between the LC group and the other cohorts, whereas the CC group clustered with the UCs (Fig. 1d). Differential bulk RNA-seq gene expression analysis in PBMCs across these groups showed upregulation of multiple proinflammatory markers in the LC compared with the CC group and uninfected individuals, including chemokines and cytokines (CXCL2, CXCL3, CCL3, IL10, IFNG, IL6, TNF, IL1B, IL1A, NFKBIZ), the NLRP3 inflammasome and the complement and coagulation genes C5, F3 and THBS1 (Fig. 1e,f). In contrast, downregulation of activating (for example, KLRC2) and inhibitory (for example, KLRC1, KIR3DL2) natural killer (NK) cell receptors and T cell activation markers in the LC compared to the CC group was observed (Fig. 1f).

To evaluate whether activation of proinflammatory pathways during acute COVID-19 infection correlated with the subsequent development of LC, we performed an exploratory analysis of transcriptomic and proteomic changes in the subset of participants in this cohort (LC: n = 8; CC: n = 5) who had PBMC and plasma samples both during acute COVID-19 (< 30 days) and at day 90–180 after infection. Transcriptomic profiling of PBMCs indicated that acutely infected participants who subsequently developed LC had higher levels of proinflammatory pathways such as IFNβ and IFNγ, JAK-STAT and IL-6 signaling, as well as innate immune cell signatures of monocytes neutrophils and complement and coagulation cascades (CCL3, CCL20, CD160, F13A1, F3, IL6, NR4A1, NLRP3, THBS1) during acute infection compared with acutely infected patients who fully recovered (CCs) (Extended Data Fig. 4a,b). Plasma proteomics profiling validated the significant increase of proinflammatory pathways (IL-6 signaling, complement cascade, leptin signaling pathway) during acute infection in participants who subsequently developed LC compared to those who recovered (CC) (Extended Data Fig. 4c,d).

To explore further the potential association between early activation of proinflammatory pathways during acute infection and the subsequent development of LC, we used a supervised random forest (RF) algorithm to identify key blood and plasma features during acute infection that predict the development of LC. This model revealed that gene expression and protein levels of complement activation, proinflammatory response, JAK-STAT, IL-6, IL-6-JAK-STAT3, IFNβ and IFNγ signaling pathways during acute infection were among the top predictors for the development of LC (Extended Data Fig. 4e,f). A feature importance analysis revealed upregulation of the IL-6 and JAK-STAT signaling, corticotropin-releasing hormone, IL-10 and TNF signaling pathways during the acute phase were among the top predictors for the development of LC. Together, these observations suggested that early activation of proinflammatory pathways strongly predicted LC development based on gene and protein profiling.

We used enzyme-linked immunosorbent assay (ELISA) and Meso Scale Discovery (MSD) assays to evaluate the plasma levels of selected proinflammatory markers at day 90–180 (LC, n = 19; CC, n = 13), more than 180 days (LC, n = 19; CC, n = 13) and UCs (n = 13) in the 2023–2024 cohort (Fig. 5c,d and Extended Data Fig. 6). Plasma levels of IL-6R were significantly elevated in the LC group compared to the CC and UC groups at day 90–180 and more than 180 days using ELISA (Fig. 5c) and MSD (Fig. 5d), suggesting a link between chronic inflammation and LC^8,12,16,33.

Finally, in an exploratory sex-stratified analysis at day 90–180 after infection, using transcriptomic profiling in peripheral blood in the combined cohort, females with LC (n = 29) showed stronger enrichment of inflammatory pathways than females who recovered (CC, n = 14), including monocyte signatures, proinflammatory cytokine/chemokine signaling, complement and coagulation cascades, IL-6, JAK1, JAK-STAT signaling, and T cell exhaustion (Extended Data Fig. 8a). In the LC group, these pathways were generally higher in females than males (Extended Data Fig. 8b), whereas no significant sex differences were detected within the CC group (Extended Data Fig. 8c). These data demonstrate persistent activation of proinflammatory pathways in a validation LC cohort, with possibly greater inflammation in females compared with males, although these observations require confirmation in larger studies.

The 2020–2021 cohort samples were collected between April 2020 and October 2021 as part of the MassCPR COVID-19 biorepository. This discovery cohort of 142 individuals included UCs (n = 35, 51% female and 49% male), acutely infected individuals (n = 54, 50% female and 50% male), CCs (n = 24, 60% female and 40% male) and patients with LC (n = 28, 86% female and 14% male). Samples were collected during the acute phase (< 30 days after infection) and during the chronic phase between 30 and 300 days after infection.

The 2023–2024 validation cohort samples were collected between October 2022 and December 2024 as part of the BIDMC biorepository jointly with the BIDMC clinical site of the NIH RECOVER study with the complete clinical meta-data and structured symptom questionnaires. This validation cohort of 38 individuals included CCs (n = 20, 55% female and 45% male) and patients with LC (n = 18, 78% female and 22% male). Samples were collected during the acute phase (< 30 days after infection) and during the chronic phase between 30 and 700 days after infection.

Symptoms in the 2020–2021 cohort were self-reported, whereas symptoms in the 2023–2024 cohort were collected with the RECOVER standardized symptom questionnaire. All samples were collected as part of the MassCPR, RECOVER and BIDMC biorepository studies with institutional review board approval and participant informed written consent. De-identified and cryopreserved samples were used in this study.

The nAb titers against the SARS-CoV-2 variants were determined using pseudotyped viruses expressing a LUC reporter gene. In brief, a LUC reporter plasmid pLenti-CMV Puro-Luc (Addgene), the packaging construct psPAX2 (AIDS Resource and Reagent Program) and Spike protein expressing pcDNA3.1-SARS-CoV-2 SΔCT were cotransfected into human HEK 293T cells (CRL-3216, ATCC) with lipofectamine 2000 (Thermo Fisher Scientific). Pseudotyped viruses of the SARS-CoV-2 variants were generated using the Spike protein from WA1/2020 (Wuhan/WIV04/2019, GISAID accession ID: EPI_ISL_402124), B.1.617.2 (Delta, GISAID accession ID: EPI_ISL_2020950) and Omicron BA.1 (GISAID accession ID: EPI_ISL_7358094.2). Forty-eight hours after transfection, the supernatants containing the pseudotyped viruses were collected and purified using filtration with 0.45-μm filter. To determine nAb titers in human sera, HEK-293T-hACE2 cells were seeded in 96-well tissue culture plates at a density of 2 × 10⁴ cells per well overnight. Three-fold serial dilutions of heat-inactivated serum samples were prepared and mixed with 60 μl of pseudovirus and incubated at 37 °C for 1 h before adding to the HEK-293T-hACE2 cells. Forty-eight hours later, cells were lysed in Steady-Glo Luciferase Assay (Promega Corporation) according to the manufacturer's instructions. SARS-CoV-2 neutralization titers were defined as the sample dilution at which a 50% reduction (NT₅₀) in relative light units was observed relative to the average of the virus control wells. Titers greater than 1:20 were considered positive.

Cellular immune responses specific to SARS-CoV-2 were assessed using IFNγ ELISpot assays using pools of overlapping 15-amino-acid peptides for the Wisconsin delta and omicron variants (21st Century Biochemicals). Ninety-six-well multiscreen plates (Merck Millipore) were coated with 1 µg per well of mouse antihuman IFNγ (MabTech) overnight in endotoxin-free Dulbecco's PBS (DPBS) overnight at 4 °C. Plates were washed with DPBS three times and blocked using Roswell Park Memorial Institute 1640 media containing 10% FCS for 2–4 h at 37 °C. Peptides pools were prepared at a concentration of 2 µg per well, and 200,000 cells per well were added. Peptides and cells were incubated for 15–20 h at 37 °C. The plates were washed with DPBS-Tween seven times and then incubated with 1 µg ml⁻¹ per well of biotinylated antihuman IFNγ (MabTech) for 2–4 h at room temperature, followed by four washes with DPBS-Tween and 1.33 µg ml⁻¹ per well of alkaline phosphatase-conjugated anti-biotin (Rockland) for 2–3 h at room temperature. Plates were developed with nitroblue tetrazolium-5-bromo-4-chloro-3-indolyl-phosphate chromogen (Pierce), stopped by washing with tap water, and read using an ELISpot reader (KS ELISPOT Reader, Carl Zeiss). The number of spot-forming cells per 10⁶ cells were calculated, subtracted over background (PBMCs incubated with medium and dimethylsulfoxide without peptide).

PBMCs were lysed in 700 μl of TRIzol and then extracted using the miRNeasy Mini Kit (QIAGEN) with on-column DNase digestion. RNA quality was assessed using a TapeStation 4200 (Agilent Technologies) and then 10 ng of total RNA was used as input for complementary DNA (cDNA) synthesis using the Clontech SMART-Seq v4 Ultra Low Input RNA Kit (Takara Bio) according to the manufacturer's instructions. Amplified cDNA was fragmented and appended with dual-indexed barcodes using the Nextera XT DNA Library Preparation Kit (Illumina). Libraries were validated using capillary electrophoresis on a TapeStation 4200, pooled at equimolar concentrations and sequenced with PE100 reads on an Illumina NovaSeq 6000, yielding ~30 million reads per sample on average. Alignment was performed using STAR v.2.7.3a⁴⁷; transcripts were annotated using a composite genome reference that included the GRCh38 Ensembl release 100 and SARS-CoV-2 (GCF_009858895.2, ASM985889v3, MN985325.1). Transcript abundance estimates were calculated internal to the STAR aligner using the htseq-count algorithm. Transcript abundance estimates were calculated internal to the STAR aligner using the htseq-count algorithm. DESeq2 (ref. ⁴⁸) was used for normalization, producing a normalized read count.

To assess for viral transcription, sample reads were aligned to a reference of 59 complete viral genome sequences using the Burrows–Wheeler Aligner (v.0.7.17)⁴⁹ to produce a sorted alignment (.bam) file. SAMtools (v.1.3.1) was used to summarize the number of reads aligned to each genome. The complete FASTA file is provided along with the bulk RNA-seq raw data.

A total of 55 ml serum or plasma from all participants, five pooled plasma controls and one buffer control were analyzed using the SomaScan Assay Kit for human plasma V4.1 (cat. no. 900-00020), measuring the expression of 6,596 unique human protein targets using 7,596 slow off-rate modified aptamer reagents (SOMAmer), single-stranded DNA aptamers, according to the manufacturer's standard protocol (SomaLogic). The modified aptamer binding reagents, SomaScan assay, and its performance characteristics and specificity to human targets, have been described previously. The assay used standard controls, including 12 hybridization normalization control sequences used to control for variability in the Agilent microarray readout process, and five human calibrator control pooled plasma replicates and three quality control (QC) pooled replicates used to mitigate batch effects and verify the quality of the assay run using standard acceptance criteria. The readout was performed using the Agilent microarray hybridization, scan and feature extraction technology.

Twelve hybridization control SOMAmers were added alongside SOMAmers to be measured from the serum samples and controls of each well during the SOMAmer elution step to control for readout variability. The control samples were run repeatedly during assay qualification and robust point estimates were generated and stored as references for each SOMAmer result for the Calibrator and QC samples. The results are used as references for the SomaScan v.4.1 Assay. Plate calibration was performed by calculating the ratio of the calibrator reference relative fluorescence unit (RFU) value to the plate-specific calibrator replicate median RFU value for each SOMAmer. The resulting ratio distribution was decomposed into a plate scale factor defined by the median of the distribution and a vector of SOMAmer-specific calibration scale factors. Normalization of QC replicates and samples was performed using adaptive normalization by maximum likelihood with point and variance estimates from a normal US population. Post-calibration accuracy was estimated using the ratio of the QC reference RFU value to the plate-specific QC replicate median RFU value for each SOMAmer. The resulting QC ratio distribution provides a robust estimate of accuracy for each SOMAmer on every plate. SomaScan RFU values and clinical information were obfuscated to protect personally identifiable information while preserving biologically relevant biomarkers. These reference datasets were provided by SomaLogic. We used the limma R package to identify differentially expressed proteins in the LC and CC groups. The method involves fitting a linear model to the data and then performing a t-test to identify proteins that are differentially expressed between two or more groups. P values were corrected for multiple testing using the Benjamini–Hochberg method. The R packages ggplot2 and ComplexHeatmap were used to generated the figures. Pathway enrichment analysis was performed using GSEA (https://www.gsea-msigdb.org/gsea↗). Genes were preranked according to the fold change from the highest to the lowest; GSEA was used to assess the enrichment of selected gene sets. Cytokine signaling, immune cell signatures and molecular pathways were compiled from the MSigDB Hallmark, C2, C7 and C3 gene sets (www.gsea-msigdb.org/gsea/msigdb/collections.jsp↗) and the blood transcriptional modules. The GSEA Java desktop program was downloaded from the Broad Institute (www.broadinstitute.org/gsea/index.jsp↗) and used with GSEA preranked module parameters (number of permutations: 1,000; enrichment statistic: weighted; 10 ≤ gene set size ≤ 5,000)⁵⁰. Sample-level enrichment analysis⁵⁰ was used to investigate the enrichment of pathways in each animal. Briefly, the expression of all genes in a specific pathway was averaged across samples and compared to the average expression of 1,000 randomly generated gene sets of the same size. The resulting z-score was then used to reflect the overall perturbation of each pathway in each sample.

We estimated pathway activity per sample using ssGSEA⁵¹. Leading-edge gene lists were parsed from the prior GSEA output and used as the gene sets for each pathway. For every sample, genes were ranked according to normalized expression and an enrichment score was computed as the difference between the empirical cumulative distributions of ranks for leading-edge genes versus all other genes; higher scores indicate greater pathway activation. Scores were normalized (per sample and then z-scored across samples) to enable comparisons, yielding a samples × pathways ssGSEA matrix for the downstream analyses.

To explore associations between symptom presence and gene expression signatures, we merged the one-hot symptom mapping with the original cohort data containing pathway-level gene expression profiles. For each of the 18 symptom categories, we then identified all patients who had that symptom and averaged their gene expression values across all measured pathways. In other words, for each symptom like fatigue or brain fog, we looked at all the patients who had that symptom and calculated the average expression level of each gene pathway in that subgroup. This produced an 18 × n matrix (where n is the number of gene pathways), with each cell representing the average expression in patients with that symptom. The resulting correlation matrix was visualized using a heatmap with color mapping reflecting average expression levels. This can be formally represented as:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{j,k}=\frac{1}{|{S}_{j}|}\mathop{\sum }\limits_{i\in {Sj}}{G}_{i,k}$$\end{document}Cj,k=1|Sj|∑i∈SjGi,k

Where C_j,k is the average expression of gene pathway (k) for symptom (j), S_j is the set of patients with symptom (j) and G_i,k is the expression value of gene pathway (k) for patient (i).

Cytokines were assessed using ELISA. An antihuman cytokine coating antibody was adsorbed onto 96-microwell plates. Microwells were washed with wash buffer; human serum samples and human cytokine standards were prediluted in assay buffer and added to each plate. Plates were then incubated for 1 h before adding antihuman cytokine horseradish peroxidase (HRP). After an additional 1-h incubation, plates were again washed with wash buffer. SeraCare KPL TMB SureBlue Start solution was added to each well; plate development was halted by adding SeraCare KPL TMB Stop solution to each well. Absorbance at 450 nm was recorded with a VersaMax Microplate Reader (Molecular Devices). The standard curve was prepared from human cytokine standard dilutions and the human cytokine concentration was determined. For each sample, the cytokine concentration was calculated using a four-parameter logistic curve fit; cytokine concentrations of the unknown samples were interpolated from the linear portion of the standard curve generated from the human cytokine standards of known concentration.

Serum levels of human IL-6R were tested using the R-PLEX Human IL-6R Kits from Meso Scale Discovery (cat. no. K1510GR-2) by the Metabolism and Mitochondrial Research Core (BIDM) according to the manufacturer's instructions. Briefly, the plate was coated using the provided biotinylated capture antibody (1:100 dilution). The highest Calibrator standard is 2,000 pg ml⁻¹. Then a fourfold serial dilution was done to generate seven calibrator curve using a four-parameter logistic model. The detection limit is 0.4 pg ml⁻¹. Samples were thawed on ice and diluted at a 1:200 ratio using the Diluent 7 provided by the kit. The assay plate was read by a MESO QuickPlex SQ 120 instrument and data were analyzed using the Discovery workbench 4.0 software.

We applied ssGSEA, which transforms gene expression profiles into enrichment scores for predefined gene sets (for example, pathways), thereby producing a matrix in which each row corresponds to a sample and each column represents a pathway activity score. To assess the contribution of individual pathways, we trained an RF classifier using the ssGSEA scores as input features. The dataset was randomly split into a training set (70%) and a testing set (30%) using stratified sampling to preserve the class distribution. The RF model was trained with 500 trees using the randomForest R package, where the outcome variable was CLASS (LC or CC) and the predictors were the ssGSEA-derived pathway scores. Model performance was evaluated using the held-out test set, with classification accuracy and a confusion matrix computed to assess prediction quality. To determine the relative importance of each pathway in predicting the LC outcome, we extracted the mean decrease in Gini impurity (MDG) from the trained RF model. The MDG quantifies how much each feature (pathway/gene) contributes to reducing classification uncertainty. A scatter plot was created to display the relationship between prediction accuracy and MDG for each pathway using transcriptomic and proteomic pathways. We added bar plots to visualize the top pathways/genes ranked according to importance generated by the RF model. Each bar represents a pathway; its length reflects its importance in the RF model. This plot highlights the most influential pathways in distinguishing LC from CC based on ssGSEA scores. All figures were generated using R v.4.4.2, using standard visualization packages for statistical and transcriptomic data analysis.

The data from the RNA-seq and proteomics were analyzed with RStudio (v.4.4.0) in PRISM 9.0 (GraphPad Software). Pathways enrichment analysis was performed using the GSEA software. Sample-level ssGSEA was performed using the GSVA (v.1.40) R package. The ssGSEA analysis delineates a comparative exploration of the groups' median values, visually represented through the data distribution and median lines. Statistical significance, denoted by the calculated P value, is highlighted above the comparison line, where a P < 0.05 is considered statistically significant. The two-sided Mann–Whitney U-test was used to identify significant differences between groups. For more than one group comparison, we used a Kruskal–Wallis test to compare all three groups; when significant, we performed post hoc pairwise comparisons with the Dunn's test with Bonferroni correction for multiple comparisons. nAb titers and IFNγ ELISpot counts were analyzed after log₁₀ transformation; to avoid undefined values, measurements <1 were set to 1 before transformation. For each viral variant (WH/2020, Delta, BA.1), CC and LC groups were compared using a two-sided Wilcoxon rank-sum (Mann–Whitney U-test) test on the log₁₀ values. Multiple testing across variants within each assay (nAb and ELISpot analyzed separately) was controlled using the Benjamini–Hochberg procedure; adjusted qvalues are reported in the subpanel titles; nominal P values are shown above the parentheses. A nonparametric test was selected a priori because immune readouts are typically right-skewed with potential outliers and unequal variances, and the groups are independent. Individual data and group medians are displayed in the plots. Analyses were performed in R using the base wilcox.test with Benjamini–Hochberg adjustment. Because of the limited sample availability, no statistical methods were used to predetermine sample sizes and we used all the samples that were made available to us. Where applicable, data distribution was assumed to be normal but this was not formally tested.

Further information on research design is available in thelinked to this article. Nature Portfolio Reporting Summary