Persistent serum protein signatures define an inflammatory subcategory of long COVID

The study cohort consisted of 55 adults (21 men, 34 women; age 22–82 years) with persistent symptoms lasting ≥60 days after an acute, PCR-confirmed SARS-CoV-2 infection (termed "PASC"), 24 (9 men, 15 women; age 20–79 years) who symptomatically recovered after a PCR-confirmed SARS-CoV-2 infection (termed "Recovered"), and 22 (12 men, 10 women; age 29–77) who had a negative nasopharyngeal PCR test (termed "Uninfected"). Participants were enrolled during the ancestral strain infection of the COVID-19 pandemic. The majority of participants had mild COVID symptoms during their acute SARS-CoV-2 infection (World Health Organization (WHO) ordinal scale 2 or 3)¹¹. Only 3 participants were hospitalized and required oxygen therapy (WHO ordinal scale 5). Two of these received Remdesivir. One additional participant received steroids. No participants in this cohort required mechanical ventilation or underwent chest computed tomography (CT). All were unvaccinated when enrolled in the study. Summary demographics for the cohort are shown in supplemental Table S1. The uninfected individuals had blood drawn once at study entry while the PASC and recovered participants had one or more blood draws at timepoints ≥60 days and up to 379 days post-symptom onset (PSO) of acute COVID (see Supplemental Fig. S1). Symptoms and other clinical metadata for each PASC participant are provided in supplemental Data S1. Previous studies have divided PASC participants into subsets based on either type, number, or severity of clinical features^7,12–14. PASC participants in our cohort reported multiple symptoms at ≥60 days PSO ranging from fatigue, fever, chills to more clinically concerning symptoms like arrhythmia or brain fog. These individual PASC symptoms were combined into organ-related symptom groups like pulmonary, cardiovascular, neurologic, etc. (see Methods) (Supplemental Data S1). For our cohort, hierarchical clustering on PASC symptomatology alone at ≥60 days PSO did not clearly drive significant participant clustering (Supplemental Table S1 & Fig. S2A). We next attempted to cluster based on serum proteins that were significantly associated with reported symptoms, but no single symptom or combination of symptoms was able to clearly distinguish participant groups (Supplemental Fig. S2B–D) suggesting that symptoms alone are unable to differentiate subsets of PASC and that additional biologic measures are needed to clarify the groups.

We analyzed the serum proteome using the Olink Explore 1536 panel (see Methods, supplemental Data S2) to first interrogate specific proteins that distinguish PASC (at their first time point available ≥60-days post symptom onset), recovered (at their last time point available ≥60-days post symptom onset), and uninfected participants in our cohort. We identified 275, 25, and 14 proteins that were significantly differentially expressed (p < 0.05) between PASC, recovered, and uninfected groups respectively when each group was compared to the other groups (Supplemental Data S3). We noted that within the PASC group, there was variation in expression of the serum proteomic signature, suggesting that some have an inflammatory signature and others do not (Supplemental Fig. S2E), and highlighting the heterogeneous nature of the disorder.

It is possible that a poor immune response to SARS-CoV-2 may allow delayed clearance of viral particles that could increase the risk for developing PASC with persistent inflammation. To address this question, we evaluated SARS-CoV-2 receptor binding domain (RBD)-specific IgG titers in serum samples obtained 60 days PSO from PASC and recovered participants in each of the 5 clusters identified above.

There was no significant difference in RBD-specific IgG responses between previously infected COVID participants in any of the clusters (Fig. 1B). We also compared SARS-CoV-2-specific CD4+ and CD8+T cell frequencies^17,18 between PASC and recovered participants from all clusters and did not identify any significant differences (Supplemental Fig. S3C & S3D).

Previous studies have suggested that severity of the acute infection may lead to a higher incidence of PASC^19–21. We hypothesized that an inflammatory serum protein signature may also correlate with being more symptomatic during acute infection. However, because our cohort primarily experienced mild acute COVID-19 symptoms (WHO ordinal scale 2 or 3), commonly used COVID severity indices did not capture the range of heterogeneity in symptomatology that we observed clinically. We therefore developed a clinical activity score for the acute phase of mild COVID. The clinical activity score for each participant reflects the impact on Activities of Daily Living (ADLs) for each day of illness. The score evaluates symptom activity for each day of illness ranging from no symptoms to life threatening illness with profound effect on ADLs and accounts for the time spent at each level (see Methods). Inflammatory PASC participants in clusters 4 & 5 had a significantly higher pre-PASC clinical activity score or higher impact on ADLs (two-sided Wilcoxon test p = 0.002) compared to non-inflammatory PASC participants in clusters 2 & 3 (Fig. 1C).

We assessed if other covariates or clinical parameters such as age, gender, BMI, or underlying comorbidities were associated with having an inflammatory serum proteomic signature. While there were no overall significant differences in age or BMI at enrollment between PASC, recovered and uninfected participants (Supplemental Figs. S4A & S4B), BMI was significantly higher in the inflammatory PASC groups (clusters 4 and 5) compared to participants in the other non-inflammatory clusters (Fig. 1D). We noted that well-described BMI associated proteins like (Leptin (LEP) and Fatty Acid Binding Protein (FABP4)) and the Leptin signaling module were significantly increased in participants with high BMI (clusters 4 and 5) in our cohort (Fig. 1E, F)^22–24. We also noted that interleukin-6 (IL-6), identified as part of the inflammatory signature associated with PASC, was positively correlated with BMI (Fig. 1E). These data suggest that BMI may be a risk factor for exhibiting an inflammatory PASC phenotype but does not account fully for the increased IL-6 observed in inflammatory PASC because a small subset of participants whose BMI was within the "healthy" range also had increased IL-6 as part of their inflammatory proteomic signature.

We next evaluated the impact of age on the inflammatory serum protein signature because we noted that PASC participants in cluster 4 were significantly older (Fig. 1G). Specifically, CXCL9, CXCL10, and IL18BP in addition to the type II interferon signaling module that have been associated with the process described as "inflammaging"²⁵, were positively correlated with higher age in cluster 4 in our cohort (Fig. 1H, I). Inflammatory proteins present in the published proteome data associated with obesity or age^22–24,26 overlapped with some of the signature proteins we found associated with inflammatory PASC but did not account for the breadth of the signature we identified. We have summarized the numerical overlap of protein markers reported in obesity, aging, and the PASC inflammatory signature in a Venn diagram (Supplemental Fig. S4C, Data S5) to demonstrate the similarities and differences between these states. Together, these data suggest that increased BMI and older age may be associated with a risk for persistent inflammation in participants with PASC but do not account for the full inflammatory signature observed in inflammatory PASC since several identified inflammatory markers were specific to PASC. We also noted a non-significant (p = 0.088) trend toward higher blood pressure in the inflammatory clusters (4 & 5) of PASC participants (Supplemental Fig. S4D) but the known association between increased blood pressure, high BMI, and older age confounds interpretation of this finding. Very few participants had other comorbidities and the presence or absence of other comorbidities was not significantly associated with the inflammatory serum proteomic signature. All p-values tested for the presence or absence of a comorbidity between the inflammatory and non-inflammatory groups with a Fisher Exact test were >0.05 (Fig. S4D). While these results are interesting, additional larger studies will be needed to confirm these observations.

To determine whether IFN-γ and IFN-γ driven cytokines and chemokines remained persistently elevated over time in inflammatory PASC, we compared expression levels in all inflammatory PASC individuals (clusters 4 + 5) with those observed in non-inflammatory PASC or in individuals that were infected with SARS-CoV-2 but recovered. Proteins were evaluated longitudinally in all available samples beginning from early acute infection to 275 days post-symptom onset (PSO). IFN-γ, IL-12 p40, and IFN-γ-driven chemokines were consistently elevated in inflammatory PASC relative to the other 2 groups (Fig. 3C and supplemental Fig. S9). To strengthen this observation, we extended our analysis to include all proteins that are part of IFN-γ related signaling modules and show that over time, the differences between inflammatory PASC vs. non-inflammatory PASC and the recovered COVID groups are even more distinct (Fig. 3D and Supplemental Fig. S10A).

In addition to the IFN-γ related signature noted above, we also observed that TNF, TNF-driven cytokines and chemokines (including IL-6 and CCL7 (MCP3)), and several TNF receptor superfamily members were also increased in both inflammatory PASC clusters but most extensively in cluster 4 (Fig. 3A, E). TNF, IL-6, and CCL7 remained persistently elevated in all inflammatory PASC (clusters 4 + 5) over time when compared to non-inflammatory PASC or individuals who were infected but recovered (Fig. 3F). In addition to elevated cytokine levels, there is evidence for persistent inflammatory cytokine signaling based on the enrichment of proteins involved in TNF signaling, the IL-18 pathway, and the NF-κB canonical signaling pathway (Fig. 3G and Supplemental Fig. S10B).

In inflammatory PASC cluster 5, we also noted a protein signature suggestive of persistent type I interferon signaling including elevation of proteins that are induced by type I IFNs during acute SARS-CoV-2 infection (SAMD9L, DDX58, MNDA, and LAMP3)^28,29. Interestingly, these type I IFN associated proteins were increased at the earliest sampling timepoint available for analysis in inflammatory PASC and remained elevated for approximately 180 days post infection (Fig. 3H and Supplemental Figs. S8 & S9). Similarly, ssGSEA analysis showed enrichment for the pathway associated with regulation of IFN-α signaling in inflammatory PASC that followed a similar kinetic (Fig. 3I and supplemental Fig. S10C). We do not have a specific readout for systemic type I IFN levels since these are notoriously difficult to accurately quantify in circulation and the Olink assay only measures IFN-γ and IFN-λ1, however the accumulated evidence points to the persistent activity of type-I interferons in PASC individuals that exhibit signs of persistent inflammation. This is notable considering recent studies reporting detection of SARS-CoV-2 RNA and spike protein in gastrointestinal and hepatic tissue of convalescent participants up to 180 days after acute infection, diverse extrapulmonary tissues including brain up to 230 days after acute symptom onset, and SARS-CoV-2 spike protein in serum up to a year post-acute infection^30–32. Whether residual viral RNA and/or protein may serve as a driver of the phenotype in inflammatory PASC remains to be investigated more thoroughly.

To summarize, broad-based serum proteomic screening identifies a subset of PASC individuals that have evidence of persistent inflammation that is sustained over the time course of this study. The signatures group inflammatory PASC into two major clusters (identified here as 4 & 5). One group (cluster 4) has prominent IFN-γ and evidence of both IFN-γ and NF-kB driven inflammation with upregulation of several inflammation-driven cytokines, chemokines, receptors, and immune checkpoint proteins. The second group (cluster 5) has evidence of ongoing IFN-γ and NF-κB driven inflammation as well but with a more prominent signature of neutrophil activation, degranulation, and/or NETosis paired with evidence for persistent type I IFN signaling. It is unclear whether these represent distinct inflammatory states or are part of a continuum of inflammation and resolution. To address this, we evaluated whether participants who were initially in one cluster at the earliest post-60 day PSO timepoint, transitioned to a different cluster at later timepoints, suggesting a transition to a different inflammatory or non-inflammatory state. We found that most individuals remain in the same cluster over the time course of observation in this study. There were a small number that transitioned within the inflammatory clusters or from inflammatory to non-inflammatory but there was no clear pattern (Fig. S). Longer follow-up and a larger cohort may be needed to understand clearly how the different immune states are related in PASC. 3B

To determine whether these observations could be extended to an independent cohort of PASC participants collected across a broader range of acute COVID severities, we used the INCOV cohort, the only other recent published dataset using Olink technology to evaluate over 442 proteins in blood samples from PASC^36,37.

The INCOV cohort includes data from 204 SARS-CoV-2-infected participants and 289 healthy controls. Of the 204 INCOV participants, 75 met the criteria used for our cohort (Olink data available from sample obtained ≥60 days after acute infection + clinical data available). Forty-three percent (43%) of these had 1 or more PASC symptoms like the PASC participants in our cohort, and the remainder had no recorded PASC symptoms, similar to the "recovered" group in our cohort. The Olink panel employed in the INCOV study measured only 442 of the 1463 proteins measured in our study but 163 proteins overlapped with the inflammatory signatures that significantly defined the two inflammatory clusters (4 & 5) in our cohort (Supplemental Fig. S). There are two other differences worth noting: First, in our study, the Olink assay was performed on serum while the INCOV cohort used plasma. Second, while both studies utilized the Olink proximity extension assay (PEA)® technology to detect the serum proteins, this was quantified using a next-gen sequencing approach in our study whereas it was quantified by multiplex quantitative polymerase chain reaction (qPCR) in the INCOV study. 11

Informed consent was obtained from all participants in the Seattle COVID-19 Cohort Study to Evaluate Immune Responses in Persons at Risk and with SARSCoV-2 Infection and to publish individual's indirect identifiers such as exact age, sex, and BMI. The Fred Hutch Cancer Center Institutional Review Board approved this study (IR10440).

Serum was collected from participants enrolled in the longitudinal study, "Seattle COVID-19 Cohort Study to Evaluate Immune Responses in Persons at Risk and with SARS-CoV-2 Infection"⁴³. Eligibility criteria included adults in the greater Seattle area at risk for SARS-CoV2 infection or those diagnosed with SARS-CoV-2 by a commercially available PCR assay. Study data were collected and managed using REDCap electronic data capture tools hosted at Fred Hutch Cancer Center, including detailed information on symptoms during acute infection and longitudinal follow-up ranging from 33–379 days post symptom onset. All but 2 persons in the "uninfected" group had at least 1 symptom of SARS-CoV-2 infection within the 14 days prior to study screening but had a negative SARS-CoV-2 nasopharyngeal PCR test. Sex of participants was determined based on self-reporting. Participants were not compensated for being in this study.

At the time of study design, the definition of PASC was evolving including some literature that utilized a persistence of symptoms lasting as few as 30 days from symptom onset to define Long COVID. However, we noticed that for many participants, the loss of smell and taste were the sole persistent symptoms in individuals that had otherwise recovered symptomatically after their acute SARS-CoV-2 infection. In most cases, this is resolved by 60 days post symptom onset. To prevent any confounding from acute recovery, we, therefore, included study participants whose symptoms had lasted for at least 60 days from symptom onset. To assure that our cohort was also representative of current definitions of PASC, we determined whether or not they also met the recent WHO Delphi Consensus criteria (at least 60 days of symptoms persisting for a minimum of 90 days post-symptom onset) (https://www.who.int/publications/i/item/WHO-2019-nCoV-Post_COVID-19_condition-Clinical_case_definition-2021.1↗). We confirmed that all but 1 PASC participant in our cohort continued to experience symptoms for more than 90 days post-symptom onset. The individual without 90-day follow-up left the study and could not be contacted.

We collected symptom information from each donor over multiple visits. Participants reported symptoms ranging from fatigue and fever to more clinically concerning symptoms like arrhythmia and brain fog. Each individual PASC symptom that was reported at ≥60 days PSO were merged into six major categories including fatigue/malaise, pulmonary, cardiovascular, gastrointestinal, musculoskeletal, and neurologic. Other mild symptoms were combined into a single category as "any mild symptoms" (Supplemental Data S). We acknowledge that these symptoms are reported and collected in snapshots of time and may not fully capture the day-to-day variability. 1

The more clinically concerning symptoms were combined into organ-related symptom groups as follows:

Fever/Chills = Fever 101 F or greater and/or Chills

Fatigue/Malaise = Fatigue or malaise

Pulmonary = Dry cough and/or SOB or DIB and/or Chest congestion and/or Chest tightness and/or Oxygen supplementation

Cardiovascular = Arrhythmias and/or Palpitations but not Tachycardia

Gastrointestinal = Diarrhea and/or Nausea or vomiting

Musculoskeletal = Muscle aches or pains and/or arthralgias and/or arthritis

Neurologic = Brain fog and/or Memory loss

Tinnitus = Tinnitus

All other non-specific symptoms were combined into the "any mild symptom" category, which include the following:

Sore throat

Scratchy throat

Low fever under 101 F

Runny nose

Loss of smell

Loss or change in taste

Headache

Itchy eyes

Dry lips in corners (angular cheilitis)

Tachycardia

Dizziness

Lightheadedness

Hair loss/Alopecia

Hands or feet tingling

Diaphoresis or sweating

Symptom information was converted to binary format where yes = 1 and no = 0. Missing symptom information is denoted by NA. For each symptom category we identified symptom-specific differential serum proteins using a linear mixed model utilizing all longitudinal timepoints available for each participant post ≥60 days PSO. We used the lme4 package (v1.1) to carry out linear mixed model analysis where age and sex were fixed variables and donor information was a random variable⁴⁴.1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{\rm{NPX}}}}}} \sim {{{{{\rm{Symptom}}}}}}\,{{{{{\rm{status}}}}}}+{{{{{\rm{Age}}}}}}+{{{{{\rm{Sex}}}}}}\,+(1|{{{{{\rm{Donor}}}}}})$$\end{document}NPX~Symptomstatus+Age+Sex+(1∣Donor)

The p-value is obtained from chi-square statistics. The specific symptom category associated with differential plasma proteins selected using p < 0.05. The identified differential proteins from six symptom specific categories were merged together and their expression visualized in a heatmap using package ComplexHeatmap (v2.4).

Symptom activity in mild to moderate acute COVID was classified by participant report of impact on Activities of Daily Living (ADLs) for each day of illness. Days hospitalized were recorded as were any treatment or therapies received. Participants were scored according to their maximum symptom activity for each day: 0, no symptoms; 1, mild impact on ADLs reported; 2, moderate impact on ADLs reported; 3, severe illness without hospitalization; 4, severe illness with hospitalization; 5, hospitalized with ICU care, or 6, life threatening illness. Duration was assigned for days spent at each level of symptom activity. A clinical activity score was calculated for each participant by multiplying the symptom activity score by the number of days spent at each level, then summing all values.

Blood was drawn into a serum separator tube and serum samples were processed, aliquoted and frozen within 4 h of blood draw.

Serum samples were inactivated with 1% Triton X-100 for 2 h at room temperature according to the Olink COVID-19 inactivation protocol. Inactivated samples were then run on the Olink Explore 1536 platform, which uses paired antibody proximity extension assays (PEA) and a next generation sequencing (NGS) readout to measure the relative expression of 1472 protein analytes per sample. Analytes from the inflammation, oncology, cardiometabolic, and neurology panels were measured. For plate setup, samples were randomized across plates to achieve a balanced distribution of age and gender. Longitudinal samples from the same participant were run on the same plate. To ensure consistent results between batches, 42 previously run serum samples from the initial batch were included on the second Olink batch as bridging controls. Per protocol, samples were run alongside a negative control (buffer), plate control (pooled serum), and a sample control (pooled serum).

Olink's standard data normalization was performed on this dataset. Protein expression values were first normalized across wells using an internal extension control (IgG antibodies conjugated with a matching oligo pair). Plates were then standardized by normalizing to the inter-plate pooled serum controls run in triplicate on each plate. Data were then intensity normalized across all cohort samples. Final normalized relative protein quantities were reported as log2 normalized protein expression (NPX) values.

Olink results and QC flags were reviewed for overall quality. Internal controls added at the incubation (non-human antigen+antibody pair) and detection (double-stranded DNA) steps were used to ensure consistency in the assay procedure among samples. There were no sample or assay failures flagged in the dataset. Warned samples were observed on some marker subsets (internal controls differ +/- 0.3 NPX from plate sample median) but were retained in the analysis. Results for TNF, IL-6 and CXCL8, which were measured on all 4 Olink panels, were reviewed prior to averaging to a single NPX value for analysis. Two samples had discrepant cross-panel measurements on these proteins. The results that trended most consistently with the participant's longitudinal measurements were kept and averaged. Serum samples were analyzed in two batches. Following the method recommended by Olink, results of the later batch were bridged to those of the earlier batch using a set of 42 cohort samples that were tested in both batches. A batch offset for each analyte was calculated as the median batch-to-batch difference on the 42 serum samples, excluding samples with QC warning flags. The analyte-specific offsets were then added to the raw NPX values of the later batch.

Half-well area plates (Greiner) were coated with purified RBD protein at 16.25 ng/well (kind gift from Dr.Leo Stamatatos⁴⁵) in PBS (Gibco) for 14–24 hr at room temperature. After 4 x 150 ul washes with 1X PBS, 0.02% Tween-20 (Sigma) using the BioTek ELx405 plate washer, the IgA and IgG plates were blocked at 37 °C for 1-2 h with 1X PBS, 10% non-fat milk (Lab Scientific), 0.02% Tween-20 (Sigma); IgM plates were blocked with 1X PBS, 10% non-fat milk, 0.05% Tween-20. Serum samples were heat inactivated by incubating at 56 °C for 30 minutes, then centrifuged at 10,000 x g for 5 min, and stored at 4 °C before use in the assay. For IgG ELISAs, serum was diluted into a blocking buffer in 7–12 1:4 serial dilutions starting at 1:50. For IgM and IgA ELISAs, serum was diluted into 7 1:4 serial dilutions starting at 1:12.5 to account for their lower concentration. A qualified pre-pandemic sample (negative control) and a standardized mix of seropositive sera (positive control) were run in each plate and used to define passing criteria for each plate. All controls and test sera at multiple dilutions were plated in duplicate and incubated at 37 °C for 1 h, followed by 4 washes in the automated washer. 8 wells in each plate did not receive any serum and served as blocking controls. Plates then were plated with secondary antibodies (all from Jackson ImmunoResearch) diluted in blocking buffer for 1 h at 37 °C. IgG plates used donkey anti-human IgG HRP diluted at 1:7500 (cat # 109-035-098); IgM plates used goat anti-human IgM HRP diluted at 1:10,000 (cat # 109-035-043); IgA plates used goat anti-human IgA HRP at 1:5000 (cat # 109-035-011). After 4 washes, plates were developed with 25 ul of SureBlock Reserve TMB Microwell Peroxide Substrate (Seracare) for 4 min, and the reaction stopped by the addition of 50 ml 1 N sulfuric acid (Fisher) to all wells. Plates were read at OD_450nm on a SpectraMax i3X ELISA plate reader within 20 min of adding the stop solution.

OD_450nm measurements for each dilution of each sample were used to extrapolate RBD endpoint titers when CVs were less than 20%. Using Excel, endpoint titers were determined by calculating the point in the curve at which the dilution of the sample surpassed that of 5 times the average OD_450nm of blocking controls + 1 standard deviation of blocking controls.

RBD titers at day 60 PSO were estimated by a linear mixed effects model of titers over time from day 30 PSO with random effects for the intercept and slope, using lme from the nlme R package.

Flow cytometry was used to examine SARS-CoV-2-specific CD4+ and CD8 + T-cell responses using a validated ICS assay. The assay was similar to a published report^43,46,47. Peptide pools covering the structural proteins of SARS-CoV-2 were used for the six-hour stimulation. Peptides matching the SARS-CoV-2 spike sequence (316 peptides, plus 4 peptides covering the G614 variant) were synthesized as 15 amino acids long with 11 amino acid overlaps and pooled in 2 pools (S1 and S2) for testing (BioSynthesis). All other peptides were 13 amino acids long, overlapping by 11 amino acids and were synthesized by GenScript. The peptides covering the envelope (E), membrane (M) and nucleocapsid (N) were initially combined into one peptide pool, but most of the assays were performed using a separate pool for N and one that combined only E and M. Several of the open reading frame (ORF) peptides were combined into two pools, ORF 3a and 6, and ORF 7a, 7b and 8. All peptide pools were used at a final concentration of 1 microgram/ml for each peptide. As a negative control, cells were not stimulated, only the peptide diluent (DMSO) was included. As a positive control, cells were stimulated with a polyclonal stimulant, staphylococcal enterotoxin B (SEB). Cells expressing IFNγ and/or IL-2 and/or CD154 were the primary immunogenicity endpoint for CD4 + T cells and cells expressing IFNγ were the primary immunogenicity endpoint for CD8 + T cells. The overall response to SARS-CoV-2 was defined as the sum of the background-subtracted responses to each of the individual pools. A sample was considered positive for CD4+ or CD8 + T cell responses to SARS-CoV-2 if any of the CD4+ or CD8 + T cell responses to the individual peptide pool stimulations was positive. Positive responses to a given peptide pool stimulation were determined using the MIMOSA (Mixture Models for Single-Cell Assays) method ⁴⁸. The MIMOSA method uses Bayesian hierarchical mixture models that incorporate information on cell count and cell proportion to define a positive response by comparing peptide-stimulated cells and unstimulated negative controls. MIMOSA estimates the probabilities that peptide-stimulated responses are responders and applies a false-discovery rate multiplicity adjustment procedure⁴⁹. Responses with false-discovery rate q-values < 0.05 were considered positive. The total number of CD4 + T cells must have exceeded 10,000 and the total number of CD8 + T cells must have exceeded 5000 for the assay data to be included in the analysis (supplemental Data S11 and Fig. S14).

Partial area under the receiver operating characteristic curve (pAUC)^50,51 was used to evaluate the rule-in performance¹⁵ of individual pathways in identifying PASC participants with respect to recovered and uninfected participants. The pAUC bounded by a specificity between 90–100% and the corresponding 99% confidence interval (two-sided) of each pathway were calculated using the "ci.auc" function in the R package pROC with the following parameters: partial.auc = c(0.9, 1), conf.level = 0.99, boot.n = 1000. A pathway was identified as significant with p < 0.01 if its pAUC lower confidence bound was above the corresponding pAUC of a random, non-performing classifier, i.e. 0.005.

We collected the canonical pathway "c2.cp.v7.2.symbols" genesets collection and associated gene information from MsigDB (v7.2). The 2871 canonical pathways were used to perform single sample GSEA (ssGSEA)⁵² using GSVA (v1.40) R package. Among 2871 pathways, 1960 pathways with overlapping serum proteins were used as input for GSVA with min.size 2 and max.size 2000 genes as parameters. The ssGSEA resulted in a normalized enrichment score (NES) for each pathway. One sample for each PASC donor was selected as follows: The first time point ≥60 days PSO for infected PASC donors (n = 55) and the last time point ≥60 days PSO for individuals that were infected but recovered (n = 24). Including the uninfected individuals (n = 22), a single sample from each of 101 donors was included in the biomarker analysis.

The rule-in approach implemented to identify pathways significantly associated with PASC donors was performed as follows: Parameters including confidence interval (CI), pAUC and bootstrap (boot.n) were used. Bootstrap analysis was performed using random seed over multiple processors using the R function mcapply (v3.4.1). Range of CI 0.8–0.99 and pAUC 0.8–0.95 was used to identify pathways associated with the PASC group. These pathways were used to differentiate the uninfected and PASC donors into separate clusters incorporating >50% of cluster size. The clustering was performed by the k-means approach implemented in ComplexHeatmap (v2.4) and visualized. The bootstrap analysis resulted in CI of 0.99 and pAUC of 0.95 which can differentiate uninfected and PASC donors in clusters. These parameters were used to identify pathways associated with PASC with a bootstrap of 1000 as noted above. The analysis resulted in 85 pathways. These 85 pathways were then collapsed into 54 modules.

A module is defined if pairwise genests had an overlap of at least 25% (jaccard index 0.25) of the genes/proteins between them using the Enrichment Map approach¹⁶. The 54 modules were then used to perform module enrichment at a single sample level using GSVA. The normalized enrichment score for each module was scaled and clustered using k-means clustering implemented in ComplexHeatmap (v2.4) with parameter row_km and column_km. The identified clusters are then visualized in a heatmap.

The INCOV proteomics data was generated using plasma samples and our proteomics data was generated using serum samples. Both used the PEA® technology but were quantified using two different assays, (NGS for our cohort and multiplex qPCR for the INCOV cohort). Olink advises against combining such datasets for analysis, so they were not combined. However, since both platforms measure relative protein abundance, the analysis was performed independently within each cohort to determine whether the same patterns and trends of protein abundance and a similar proteomic signature can be observed across groups. The Olink proteomic data from Su Y et al. consisted of 204 SARS-CoV-2 (INCOV) participants and 289 healthy controls. The INCOV participants were studied at clinical diagnosis (T1), acute disease (acute, T2), and 2–3 months post onset of initial symptoms (convalescent, T3). Olink plasma proteomic data was available for a total of 443 proteins. Among these, 163 proteins overlapped with the differentially expressed proteins found in inflammatory signatures that significantly defined clusters 4 & 5 in our cohort. K-means unsupervised clustering of the INCOV Olink proteomic data was performed on the 163 protein overlap. To remain consistent with our cohort, we used only those samples available for timepoints ≥60 days PSO per INCOV participant (74 INCOV participants met these criteria). If more than one sample was available ≥60 days PSO, we used the sample closest to the ≥60 day timepoint. The k-means function of the stats R package was used with k = 5, allowing 100 iterations.

A list of fifteen proteins (Supplemental Data S10) were selected as candidates for a protein panel to distinguish inflammatory versus non-inflammatory PASC. Several criteria had to be met for inclusion of proteins in this group: All proteins had to be measured in both our study and INCOV, were significant in distinguishing the two PASC groups (adjusted p < 0.008) in our study cohort and had strong support in the literature as biomarkers for COVID-19 and/or PASC. Logistic regression (LogReg) models⁵³ were used to combine protein measurements into LogReg probability scores as follows:2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W={c}_{0}+\mathop{\sum }\limits_{i=1}^{n}{ci}*{{NPX}}_{i},$$\end{document}W=c0+∑i=1nci*NPXi,3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${score}=\frac{1}{1+{\exp }^{-W}}$$\end{document}score=11+exp−Wwhere \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${c}_{0}$$\end{document}c0 is the intercept, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document}n is the protein number, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${c}_{i}$$\end{document}ci is the coefficient of protein \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i$$\end{document}i, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{NPX}}_{i}$$\end{document}NPXi is the NPX value of protein \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i$$\end{document}i. The models were assessed on our PASC data and evaluated in a stepwise backward elimination procedure: In each step, the approach of 10-fold cross validation was repeated 10,000 times to train LogReg models and evaluate the corresponding p values of individual proteins having a non-zero coefficient. The protein having the highest p-value in the step was eliminated. This procedure was repeated until all proteins had p < 0.05. Coefficients of the final LogReg model were then set to the corresponding median coefficients of the remaining proteins.

Challenges associated with sample collection, supply-chain difficulties, and a shifting research environment during the global COVID-19 pandemic required that the Olink analysis be performed in two batches that impacted batch design. The first batch was collected early in the pandemic (April 2020) with the aim of studying the difference between SARS-CoV-2 infected (n = 15) versus uninfected donors (n = 22). We later found out two of the fifteen infected donors had PASC and had persistent inflammatory signatures over time, which motivated us to expand our study by including a larger cohort of PASC. By the time the second batch was collected for analysis, PASC had become a major concern in COVID-19 patients, so we added an additional 53 PASC donors and 11 recovered COVID-19 donors in the second batch based on sample availability. Based on this, the two batches have a bias related to disease status (PASC, recovered, and uninfected; p < 0.05) but no biases on biological sex or age (p > 0.05). Longitudinal samples were collected from participants based on their availability, but no replicate samples were collected at any single timepoint. Hence, no statistical method was used to predetermine sample size. We were not blinded to allocation and outcome assessment. To minimize possible batch effects, we used 42 cohort samples as bridging controls for batch correction between the two batches as described in the 'Olink data preprocessing' section of the methods. No data were excluded from the analyses.

All statistical analyses were performed using the corresponding functions in RStudio (version 4.1). Comparisons of single protein Olink NPX or module ssGSEA scores between groups were tested using a two-sided Wilcoxon rank sum test and when appropriate, the Benjamini-Hochberg method was applied to adjust p-values in multiple-testing correction. Unless specified, an adjusted p-value of 0.05 was considered significant and p-values less than 0.0001 were reported as p-value < 0.0001. All figure panels were assembled using Adobe Illustrator (version 27.2; https://www.adobe.com/products/illustrator.html↗).

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

Supplementary Information Description of Additional Supplementary Files Supplementary Data 1 Supplementary Data 2 Supplementary Data 3 Supplementary Data 4 Supplementary Data 5 Supplementary Data 6 Supplementary Data 7 Supplementary Data 8 Supplementary Data 9 Supplementary Data 10 Supplementary Data 11 Reporting Summary