Immunological and Antigenic Signatures Associated with Chronic Illnesses after COVID-19 Vaccination

Cohort Description

All the blood samples were collected between December 2022 and November 2023 from the Listen to Immune, Symptom and Treatment Experiences Now (LISTEN) study²³. The PVS cohort consisted of a total of 42 participants, including 29 females and 13 males with no preexisting comorbidities, whereas the control cohort consisted of 22 participants with 11 from each biological sex (Figures S1A and 1A). Upon recruitment of 44 PVS participants, two had to be excluded from the analyses due to evidence of pharmacological immunosuppression. Information on index vaccine types was reported by 39 out of 42 PVS participants included in the analyses, and they were Comirnaty (Pfizer) (n=14), Spikevax (Moderna) (n=21), and Jcovden (J&J) (n=4). The most frequent symptoms reported by participants were excessive fatigue (85%), tingling and numbness (80%), exercise intolerance (80%), brain fog (77.5%), difficulty concentrating or focusing (72.5%), trouble falling or staying asleep (70%), neuropathy (70%), muscle aches (70%), anxiety (65%), tinnitus (60%) and burning sensations (57.5%). Further, pairwise Euclidean distances were calculated in a sex-segregated manner based on the presence of symptoms at recruitment and two distinct clusters of symptoms were identified in both (Table S1, Figure S1B).

Between the case and control cohorts, a total of 15 (35.7%) and 10 (45.5%) reported having a history of one or more previous SARS-CoV-2 infections, respectively (Table S2). However, upwards of 40% of SARS-CoV-2 infections are asymptomatic²⁴. To further investigate prior history of SARS-CoV-2 infections, plasma specimens were analyzed using the EUA-cleared Elecsys® anti-SARS-CoV-2 immunoassays, which measure the presence of high affinity IgM, IgA, and IgG anti-N antibodies. A cut-off index ≥ 1 was defined as reactive based on previous literature²⁵. In the non-reactive group, the antibody indices varied between 0.09 and 0.17, whereas in the reactive group, they ranged between 1.37 and 94.4. Among participants with PVS, 26 (61.9%) were found to be reactive compared with 10 (45.5%) among controls. One participant from each cohort with self-reported history of infection had non-reactive test results. Based on both self-reports and serological analyses the two cohorts were further classified into four subgroups, PVS with no history of infection (PVS-I, n= 15), PVS with a history of infection (PVS+I, n= 27), controls with no history of infection (Control-I, n= 11) and controls with a history of infection (Control+I, n= 11) (Figure 1A, Table S2). Even though all PVS participants developed chronic symptoms following vaccination and not infection, it was important to consider the impact of a subsequent SARS-CoV-2 infection on immune phenotypes analyzed in this study.

• Download figure
• Open in new tab

Figure 1: Study overview and cohort information.

A. Summary of the study design, cohorts and assays included in this study B. Stacked bar plots showing the number of males and females in the four subgroups. C. Scatter plots showing the distribution of age across subgroups. D. Boxplots showing self-reported General Health Visual Analogue Scales (GHVAS) scores self-reported on the day of collection. E. Scatter plots showing the Patient-Reported Outcomes Measurement Information System 29-Item Profile Measure (PROMIS-29 v2.0) score differences across domains. Differences in biological sex among the four groups was assessed using Fisher’s exact test. For the rest, significance was assessed using either Mann-Whitney U and Kruskal-Wallis tests with Benjamini–Hochberg false-discovery rate (FDR) correction for multiple comparisons wherever necessary. Only significant differences are highlighted (p values ≤ 0.05). F. A pie chart showing the proportion of index vaccine doses within the PVS cohort. G. Lines plots showing the distribution of days post-vaccination that the participants develop any or severe symptoms associated with PVS.

Among the demographic variables, there were no significant differences in the number of males and females between cases and controls (Fisher’s Exact test, p= 0.58) or among the four groups (Kruskal-Wallis test, p = 0.25; Figure 1B). Similarly, no significant age differences were observed between cases and controls, median age (42.5 years, PVS; 38 years, controls, Mann-Whitney U test p= 0.27) and among the four groups (p = 0.17; Figure 1C).

The self-reported General Health Visual Analogue Scale (GHVAS) scores on the day of biospecimen collection differed significantly among the four groups (Kruskal-Wallis test, p = <0.01). The controls in both subgroups had significantly higher median scores compared with the PVS subgroups (64, PVS-I; 60, PVS+I; 95, Control-I; 90, Control+I; Figure 1D). The PROMIS-29 physical function, fatigue, pain interference, depression, anxiety, sleep disturbance, and pain interference scores were compared independently among the four groups to gauge the physical and mental health status of the participants. The physical function scores were significantly higher among the controls than the cases irrespective of infection status, median scores (13, PVS-I; 14.5, PVS+I; 20, Control-I; 20, Control+I). The anxiety (9.5, PVS-I; 10, PVS+I; 5, Control-I; 4, Control+I), depression (8.5, PVS-I; 8, PVS+I; 4, Control-I; 4, Control+I), fatigue (16, PVS-I; 15, PVS+I; 6, Control-I; 7, Control+I) and pain scores (9.5, PVS-I; 12, PVS+I;4, Control-I; 4, Control+I) were significantly lower in controls compared with participants with PVS irrespective of infection status. Further, significantly higher sleep disturbance scores were observed only among infection-negative cases compared to control-I participants (13, PVS-I; 10, Control-I) (Figure 1E).

Most individuals in each cohort completed the primary series of vaccines based on WHO recommendations (83.3%, PVS; 100%, Controls; Fisher Exact test, p = 0.09). Participants with PVS received significantly fewer COVID-19 vaccine doses compared with controls, median vaccine numbers (2, PVS; 4, controls; Fisher Exact test, p = <0.01). On similar lines, the median number of days post the latest vaccination was significantly higher among cases with a median of 585 days (±190) compared with 199 days (±217) among controls (Mann-Whitney U test, p= <0.01). In 85% of the cases, participants identified the index vaccine dose as being part of the primary series [dose 1(45%) and dose 2(40%); Figure 1F]. The median number of days for the development of any symptom was 4 [Interquartile range (IQR): 23 days], while for severe symptoms, it was 10 (IQR: 44 days) post-vaccination. A high proportion of participants with PVS developed any symptoms (70%) or severe symptoms (52.2%) within 10 days of vaccination (Figure 1G).

Differences in circulating immune cell populations

To determine immune signatures of PVS, peripheral blood mononuclear cells (PBMC) were analyzed using flow cytometry. Among the cell populations of myeloid lineage, proportions of non-classical monocytes (CD14^lowCD16^high; Mann-Whitney U test, p= 0.03) were significantly higher in the PVS cohort compared to the controls without significant differences in the percentage of total monocytes despite greater median values in PVS (Figure 2A). The median percentage of conventional type 2 dendritic cells (cDC2; CD304^-/HLA-DR⁺/CD1c⁺) was significantly lower among the participants with PVS compared to the controls (p= 0.02) while no differences were observed in the proportions of conventional type 1 dendritic cells (cDC1; (CD304^-/HLA-DR⁺/CD141⁺) (Figure 2B). Pairwise comparisons were also executed to understand the differences among the PVS subgroups with or without a history of infection. Among the low-density granulocytes, no differences were observed between cases and controls in the proportions of eosinophils (CD66b⁺CD56^-CD16^-) or between the cases and controls in the infected or uninfected subgroups but the proportion of neutrophils (CD66b⁺/CD56^-/CD16⁺) was significantly higher (p= 0.02) in infection positive PVS subgroup (PVS+I) compared to the convalescent controls (control+I) (Figure S2A). Significantly lower and higher proportions of classical and non-classical monocytes, respectively, were observed among the PVS+I compared to the control+I subgroup (p(cMonocytes)= <0.01; p(ncMonocytes)= 0.03) with no differences between the cases and controls without prior history of SARS-CoV-2 infection (Figure S2B). Next, significantly higher proportions of both cDC1 and cDC2 cells were observed in the control+I subgroup compared to the PVS+I subgroup (p= 0.03 and p= <0.01), respectively (Figure S2C).

• Download figure
• Open in new tab

Figure 2: Immune cell feature of myeloid and lymphoid cells in PVS patients.

A-B. Violin plots of myeloid peripheral blood mononuclear populations (PBMCs) plotted by groups as percentages of respective parent populations (live cells). C. Violin plots of B lymphocyte subsets from PBMCs plotted as percentages of respective parent populations (total B cells). D. Violin plots of various CD4 T cell subsets. E. Violin plots of various cytokine-producing CD4 T cell subsets. F. Violin plots of various CD8 T cell subsets and cytokine-producing CD8 T cell subsets. Significance differences were assessed using Mann-Whitney U tests with Benjamini–Hochberg false-discovery rate (FDR) correction for multiple comparisons. G. Linear regressions of TNFα producing CD8 T cells and IFNγ producing CD8 T cells. Spearman’s correlation was calculated with corresponding p-values. Dotted lines depict linear regressions, with the area inside representing 95% CI.

Among the B cell populations, relative proportions of unswitched memory B cells (US memory B cells; CD19⁺/CD27⁺/IgD⁺) were significantly higher (p= 0.02) while the proportion of double negative B cells (DN B; IgD^-/CD27^-/CD24^-/CD38^-) was observed to be lower (p= 0.01) in the PVS cohort compared with controls (Figure 2C).

Significant differences in subsets of circulating immune cell populations were observed across T cell lineages. Upon assessment of the T cell populations, notably higher proportions of effector memory CD4 T cell subsets (CD4⁺Tem; CD45RA^-/CD127⁺/CCR7^-; p= 0.01) and resting natural CD4⁺ Treg; CD45RA⁺/CD25⁺/CD127^-/HLA-DR^-; p= 0.05) were observed among the controls (Figure 2D). However, the PVS cohort had significantly higher proportions of exhausted CD8 T cell (CD8⁺ Tex; PD-1⁺/TIM3⁺; p= 0.02) (Figure 2F) with no observed differences in the CD4⁺ central memory (CD4⁺ cm; CD45RA^-/CD127⁺/CCR7^-) and exhausted (CD4⁺Tex; PD-1⁺/TIM3⁺) CD4 T cell populations (Figure 2D). Upon in-vitro stimulation, the expression of CXCR3 on the cell surface (Mann-Whitney U test, p= <0.01), intracellular IL-4 levels (p= 0.04) and IL-4, IL-6 in combination were found to be significantly lower in the CD4 T of PVS cohort (p= <0.01), with no differences were observed in IFNγ & TNF⍺ levels (Figure 2E). Significant increases in intracellular TNF⍺ levels (p= <0.01) with non-significant increases in IFNγ in the stimulated CD8 T cells were observed in PVS cohort (Figure 2F). Only a total of 32.23% of the variability in intracellular TNF⍺ could be explained by IFNγ levels in the CD8⁺ T cell populations (R²= 0.32; Figure 2G).

In the subgroup analyses, no differences in proportions of DN B cell subpopulations were observed (Figure S2D). Proportions of CD4⁺ CD45RA⁺ effector memory T cells (CD4⁺ T_EMRA; CD45RA⁺/CD127^-/CCR7^-; p= 0.02) and rnTregs (p= 0.03) were both observed to be significantly lower in PVS+I compared to the control+I subgroup (Figure S2E). Proportions of CXCR3 expressing stimulated CD4 T cells was much lower in PVS+I cases (p= <0.01) and the proportions of both IL-4⁺ (p= 0.04), and IL-4⁺/IL-6⁺ (p=< 0.01) double-positive cells were also lower compared to the controls (Figure S2E). Higher proportion of CD8⁺ Tcm cells was retained in the PVS+I subgroup compared to the control+I group (Figure S2F; p= 0.02). No differences were observed in immune cell populations between the infection-negative cases (PVS-I) and controls (Controls-I).

Lower levels of spike-specific antibody responses in PVS

Given the differences in the number of vaccine doses received between participants in the PVS cohort and the control group, we compared spike-specific immunoglobulin G (IgG) levels in relation to the number of vaccine doses administered. Correlation analyses revealed a significant positive correlation between the number of vaccine doses and plasma anti-S IgG levels (Spearman’s Rank Correlation Coefficient, ⍴ = 0.85, p= <0.01) as well as anti-RBD IgG levels (⍴ = 0.83, p= <0.01) in the PVS-I subgroup. In the PVS+I subgroup, only anti-S IgG levels showed a significant correlation (⍴ = 0.55, p = 0.01) with the number of doses (Figure 3A). Next, correlation analyses were performed to assess the relationships between plasma anti-S, anti-RBD, and anti-N IgG levels with the number of days post last vaccination among the four groups. No significant changes in anti-S and anti-RBD antibody levels were observed with increasing days since vaccination in the control group, regardless of infection history, and in the PVS+I subgroup (Fig 3B). In contrast, significant negative correlations were found in the PVS-I subgroup between the number of days post-vaccination and both anti-S (⍴ = -0.87, p= <0.01) and anti-RBD (⍴ = -0.83, p= <0.01) IgG levels, indicating a decline in these antibodies over time (Fig. 3B). Additionally, as expected, no correlations were observed between anti-N IgG levels and days post vaccination across the infection-positive subgroups.

• Download figure
• Open in new tab

Figure 3: Plasma reactivity to SARS-CoV-2 antigens.

A. Correlation comparisons of virus-specific ancestral anti-S and anti-RBD IgG levels by number of COVID vaccine doses. B. Correlation and linear regression comparisons of virus-specific ancestral anti-S, anti-RBD and Anti-N IgG levels by days post last vaccination. C. Correlation and linear regression comparisons of virus-specific ancestral anti-S, anti-RBD and Anti-N IgG levels by days post last exposure. Regression lines are shown colored by groups Control-I, Control+I, PVS-I, and PVS+I as indicated in the figure legend. Spearman’s ρ coefficients and linear regression significance are colored; accordingly, shading represents 95% confidence interval. D. Plasma reactivity to ancestral S, RBD, and N proteins measured by ELISA are shown by groups Control-I, Control+I, PVS-I, and PVS+I. Significance of difference in group median values was assessed using Kruskal– Wallis with Benjamini–Hochberg false-discovery rate (FDR) correction for multiple comparisons. The central lines indicate the group median values, and the whiskers show the 95% CI estimates. E. Generalized linear model analysis for virus-specific ancestral anti-S, anti-RBD and Anti-N IgG levels. Model predictors are indicated on the x axis and include days from vaccination (DFV) among others. Predictors with p ≤ 0.05 are highlighted in pink to indicate significance, while non-significant predictors are displayed in black. Detailed model results are shown in table S3.

The next step was to evaluate if the most recent exposure to SARS-CoV-2 or vaccination correlated with the observed differences in waning patterns. No significant changes in anti-S, anti-RBD and anti-N antibody levels were observed with an increase in the number of days from self-reported viral infection dates among the Control+I and PVS+I subgroups (Figure 3C). In addition, the plasma titers of anti-S IgG were significantly lower among the PVS-I cases compared to the Control-I subgroup (p= <0.01) (Figure 3D). However, no differences were observed in the anti-RBD IgG levels across the four subgroups (Figure 3D). As expected, the uninfected PVS-I and the Control-I subgroups had much lower anti-N IgG levels as detected by in-house ELISAs (Figure 3D). To further account for variations in vaccine doses and infection, we developed linear models. Those models indicated that both prior SARS-CoV-2 infection and the number of vaccine doses were significantly associated with higher levels of anti-RBD and anti-S IgG (Figure 3E, Table S3).

Serological evidence of recent EBV reactivation in PVS

Many human pathogens are ubiquitous, opportunistic, and capable of establishing lifelong infections with alternate latency and reactivation cycles²⁶. These cycles can be triggered by physiological perturbations and can contribute to systemic inflammation²⁷. Therefore, we used serum epitope repertoire analysis (SERA) to evaluate seropositivity against a range of pathogens, including five bacterial, seven parasitic,14 viral and one fungal species. On performing two group analyses, no significant differences were observed for all pathogens, indicating similar levels of prior exposure. (Fig. 4A).

• Download figure
• Open in new tab

Figure 4: Elevated responses to Epstein Barr Virus in PVS patients.

A. Proportion of each group (PVS: n = 42, control: n = 22, pre-pandemic healthy control: n = 3448) seropositive for each of 31 common pathogen panels as determined by SERA, grouped by pathogen-type. Statistical significance determined by Fisher’s exact test corrected with FDR (Benjamini Hochberg). Star indicates panels for which pre-pandemic healthy controls were not analyzed. B. Heatmap showing supervised clustering of SERA-determined seropositivity to EBV, CMV, HSV-1, and HSV-2 across samples. Clusters were named for their herpesvirus dominance and are labeled accordingly. C. Herpes seropositivity composition for each cohort. Significance of relative enrichment for each cluster was calculated using Chi-square test of observed composition vs. expected composition. D, E Plasma reactivity to EBV gp42 protein measured by ELISA shown by cohort, PVS and Control (D) and by groups Control-I, Control+I, PVS-I, and PVS+I (E). F. SERA-derived z scores for the gp42 motif [VI]XLPHW among EBV-seropositive individuals only, plotted by cohort, n = 20 (Control), n = 38 (PVS) (F) and group, n = 11 (Control-I), n = 9 (Control+I), n = 12 (PVS-I), n = 26 (PVS+I) (G). The dashed line represents the z-score threshold for epitope positivity defined by SERA. H. Three-dimensional mapping of the PVS-enriched linear peptide sequence [VI]XLPHW (purple) onto EBV gp42 (blue) in a complex with gH (light grey) and gL (dark grey) (PDB: 5T1D). I. Relationship between EBV gp42 [VI]XLPHW SERA z score and plasma concentration of anti-gp42 IgG. J. SERA-derived z scores for the gp350 motif KXRX[RQ]WXF among EBV-seropositive individuals only, plotted by cohort. The dashed line represents the z-score threshold for epitope positivity defined by SERA. K. The relationship between plasma concentration of IgG against EBV gp42 and the percentage of TNFα CD8+ T cells (of total CD8+ T cells). For all box plots, the central lines indicate the group median values, the top and bottom lines indicate the 75th and 25th percentiles, respectively, the whiskers represent 1.5× the interquartile range. Each dot represents one individual. Statistical significance of the difference in median values was determined using Kruskal–Wallis tests with Post hoc Dunn’s test and Bonferroni–Holm’s method to adjust for multiple comparisons. Correlation was assessed using Spearman’s correlation. The black line shows linear regression, and shading shows the 95% CIs.

Moreover, the seropositivity for each pathogen did not significantly differ from seropositivity in 3448 healthy controls collected before the COVID-19 pandemic (Figure 4A). Given the high seropositivity rates for herpesviruses, we further analyzed the seropositivity patterns in combination, for cytomegalovirus (CMV), Epstein-Barr Virus (EBV), Herpes Simplex Virus Type 1 (HSV-1) and Herpes Simplex Virus Type 2 (Figure 4B). Significant differences were observed between cases and controls (Mann-Whitney U test, p = 0.01; Figure 4C), where the participants with PVS had higher prevalence of EBV and HSV coinfection, and lower prevalence of EBV and CMV coinfection. There are reports of similarities in symptom phenotypes between PVS and long COVID, as well as evidence of EBV reactivation in long COVID cases, including elevated antibodies against EBV surface protein gp42^{22, 28}. Therefore, we further investigated the prevalence of antibodies against EBV gp42 and identified significantly elevated antibodies in the plasma of PVS participants compared with controls (Kruskal-Wallis test, p = <0.01, Figure 4D, E). As an orthogonal validation, we tested the distribution of linear peptide reactivities across the EBV proteome. Greater reactivities to two peptides corresponding to two envelope glycoproteins necessary for B cell infection, gp42 and gp350 were observed. For the gp42 protein, the antibody reactivity to peptide ([VI]XLPHW) was significantly higher among the PVS participants irrespective of the SARS-CoV-2 infection status (Mann Whitney U test, p= <0.01; Figure 4F) and across the four subgroups (Kruskal-Wallis test, p= 0.02; Figure 4G). Greater reactivities were also observed for the gp350 peptide (KXRX[RQ]WXF) among the PVS participants compared to controls (Kruskal-Wallis test, p= 0.04; fig 4J) and across the four subgroups (Kruskal-Wallis test, p= 0.03; Figure S3C). Further, anti-gp42 ([VI]XLPHW) reactivity by SERA significantly correlated with anti-gp42 ELISA measurements thus validating the finding (R = 0.37, p= <0.01; Figure 4I). We also mapped this motif onto available structures of gp42 complexed with EBV gH/gL (PDB: 5T1D), demonstrating its location close to the transmembrane (TM) domain of gp42 and surface-exposed (Figure 4H). Study participants with greater antibody reactivity to gp42 as assessed by ELISA also exhibited higher percentages of TNFα-producing CD8+ T-cells (R = 0.47, p= <0.01, Figure 4K). This correlation was not observed for IL-4, IL-6 double-positive CD4 T cells (Figure S3D) as was previously reported for long COVID²².

Participants with PVS have a distinct set of autoantibodies

To evaluate differences in immunoglobulin isotypes and IgG subtypes in the plasma, Luminex assays were performed. No significant differences were observed between the PVS cohort and the controls (Figure S4A). Next, to determine the presence of autoantibodies in PVS, we screened for reactivities across a range of 120 known autoantigens using microarrays for three different immunoglobulin isotypes, IgM, IgG, and IgA. We observed significant increases in IgM reactivities against 65 antigens, IgG reactivity against 1 antigen and IgA reactivities against 39 antigens in PVS compared to controls after multiple testing corrections (Table S4). Among these antigens, two showed log₂fold change of greater than 2: anti-nucleosome IgM (Mann-Whitney U test, p= <0.01) and anti-AQP4 IgA (p= <0.01) (Figure S4B). Conversely, control participants exhibited higher reactivities against a total of 21 antigens, 18 of which were of the IgG isotype and five were of IgA isotype with two common antigens between the two isotypes (Table S3). Among these autoantibodies, anti-histone H1 IgG differed by greater than log₂fold change (p= <0.01, Figure S4B). Infection-positive subgroups had a higher number of reactivity differences between cases and controls (Figure S4C). Among the PVS-I participants, anti-calprotectin/S100 IgM, anti-genomic DNA IgA and anti-ssDNA IgA reactivities were significantly higher while anti-histone H3 IgG, anti-MBP IgA, and anti-PR3 IgA reactivities were higher among the controls-I (Figure S4C, Table S5).

Circulating hormones and immune modulators in PVS

Two group analyses of circulating hormones and immune modulators revealed significantly lower levels of fetuin A26 and neurotensin (Mann-Whitney U tests, p= 0.01 and p= 0.03; Figure S5A) in participants with PVS with fold differences of 1.3 and 1.9 respectively. Additionally, four group analysis was performed to evaluate the impact of infection on PVS. Given the smaller number of samples in the four group analyses, each panel of analytes was independently evaluated. Significantly lower levels of circulating fetuin A36, and neurotensin were also observed among participants with PVS with a history of SARS-CoV-2 infection compared to convalescent participants (p= 0.01 for both analytes; Figures S5B-C). No differences were observed for other factors across the subgroups except for β endorphin which was significantly lower in PVS+I compared to the control+I group (p= 0.01; Table S7) without any significant differences in the two group analyses. No differences were observed in the uninfected subgroups.

Increase in circulating SARS-CoV-2 Spike protein in participants with PVS

It has been reported that the BNT162b2 or mRNA-1273 derived S proteins circulate in the plasma of those vaccinated as early as one day after the vaccine and interactions of the circulating protein¹⁶. Hence, we next sought to investigate whether the S1 subunit of the SARS-CoV-2 S protein could be detected in the plasma. For this, we used an anti-S1 Successive Proximity Extension Amplification Reaction (SPEAR) immunoassay. This method can detect S1 levels as low as 5.64 fM. We conducted a one-sided Kolmogorov–Smirnov test with 1000 permutations to see if the participants with PVS had higher circulating S1. The results indicated that participants with PVS had significantly higher circulating S1 levels compared with the control group (p = 0.01). However, circulating S1 was found in only a subset of participants with PVS at varying concentrations while the control group mostly exhibited a bimodal distribution of zero and non-zero values (Fig. 5A, Table S2). Detectable S1 was found in participants’ plasma ranging from 26 to 709 days from the most recent known exposure (Figure 5B). To fully account for the width of this dataset, we included all non-detectable values in the analysis and applied a generalized regression model accounting for zero-inflation. We found that both PVS-I and PVS+I groups displayed significantly elevated S1 levels than the Control-I group (p= <0.01 and p= 0.02, respectively) (Figure 5C).

• Download figure
• Open in new tab

Figure 5: Circulating SARS-CoV-2 Spike protein.

A. Density plots describing the distribution of circulating S1 levels across controls (n= 22) and PVS (n= 42) measured by SPEAR assays. B. Levels of circulating S1 in plasma days post last known self-reported exposure. C. Circulating S1 levels measured by SPEAR assay are shown across groups Control-I, Control+I, PVS-I, and PVS+I. A parametric test incorporating a zero-inflated Poisson model was used to account for the excess zeros in the data. D. Circulating S1 antigen levels above LLoD across cohort groups, MY-LC-HC/CVC (n= 41), MY-LC-LC (n=45), Control (n= 7), PVS (n= 15). E. Correlation between circulating Spike protein assays using antibodies for the S1 subdomains and S1& S2 subdomains (full length Spike). F. Correlation between circulating Spike protein assays among samples with values above LLoD and LLoQ. Correlations were assessed using Spearman’s rank correlation. G. A participant-specific graphic representation illustrating key events including vaccination, infection, sample collection and the presence or absence of anti-N antibodies and circulation spike protein. Each participant is represented by a single horizontal line and each vaccination event is marked by triangles, and the index doses are highlighted in blue. The number of days between latest exposure and biospecimen collection is also indicated. The abbreviations for the vaccine types are: J for Jcovden (Johnson & Johnson), M for Spikevax (Moderna), and P for Comirnaty (Pfizer-BioNTech). H. A classification tree of PVS participants based on infection status with or without detectable S1 in circulation and a heatmap of distinct demographic and immunological variables that differentiate PVS within the infected and uninfected subgroups based on Mann-Whitney U tests.

Given the similarities between PVS and long COVID symptoms, one hypothesis in the literature is that shared exposure to the S protein may play a role and several groups have independently reported the presence of circulating S1 & full-length S in long COVID using various detection methods^16,29. To further investigate this circulating S1 positivity percentages and levels in the LISTEN cohort subgroups were compared with an external cohort of healthy, convalescent controls and LC participants (MY-LC cohort). This external cohort, collected from Mount Sinai clinics, included 134 healthy/convalescent controls and 134 long COVID participants, which were all assayed together with the LISTEN cohort biospecimens.

Among the MY-LC healthy [HC(n= 62); no reported SARS-CoV-2 infection], convalescent controls [CVC(n=72); with reported SARS-CoV-2 infection], and long COVID participants, detectable S1 was observed in 30.6% of control participants (41/134; HC= 12.9%; CVC= 22.2%) and 33.6% (45/134) of individuals with long COVID, with mean S1-ln(x_fM +1) values of 3.72 and 3.85 respectively among those above the LLoD. These figures were comparable to the percentages observed among the LISTEN controls (31.8%) and PVS (35.7%) groups. S1 levels were moderately elevated in the MY-LC control group compared to the PVS control group (p=0.06), potentially reflecting differences in exposure timing or SARS-CoV-2 variant of concern (VOC). Despite this, the PVS group demonstrated significantly higher S1 levels compared to both control cohorts (LISTEN-control: p<0.01; MY-LC control: p=0.03; mean S1-ln(x_fM +1) = 6.24) (Fig. 5D, Table S8).

To further validate the findings and to investigate whether the presence of S1 reflects the presence of full-length S protein among the LISTEN participants, we next conducted a full-length S SPEAR assay. The calculated values for the LLoD and the LLoQ were 1.81 and 8.24 fM, respectively. The full-length S SPEAR assays showed a significant correlation between S and S1 across all samples (Fig. 5E), as well as for values above the LLoD and LLoQ (Fig. 5F, Table S2). Thus, based on SPEAR assays, the individuals with PVS exhibited elevated levels of circulating full-length S compared to healthy controls.

Immune signatures in PVS subgroups based on the presence of circulating S1 protein

To gain a clearer understanding of the variability of circulating S1 protein levels, we first compiled a structured timeline that summarizes the self-reported infection dates, vaccine numbers (including types and administration dates), and the number of days between the latest known exposure and the collection of biospecimens. This timeline was organized for both the PVS-I and PVS+I groups (Figure 5G). Notably, we observed that the highest levels of detectable S1 in the PVS-I group were the furthest away from the last known exposure and ranging between greater than 600-700 days (NI-1 & NI-5; Figure 5G). This suggested that prolonged antigen persistence might be associated with PVS in a subgroup of patients. Further, most of the PVS+I group participants experienced breakthrough SARS-CoV-2 infections with the exception of two cases, indicating that PVS symptoms started prior to infection (Figure. 5G).

Given the possible heterogeneity in immunological trajectories leading up to PVS and the lack of adequate sample numbers in each PVS subgroup, we next took a more descriptive approach to look for peripheral immune signatures stratified based on their infection status and detectible S1 above SPEAR assay’s lower limit of quantitation. In order to begin with a valid method of selection despite the small sample sizes, non-parametric Mann-Whitney tests were implemented without multiple testing corrections to look for differences in distributions of 547 independent variables including GHVAS scores, circulating modulator levels, anti-SARS-CoV-2 antibody titers and autoantibody scores within both the PVS-I and PVS+I subgroups. Variables showing significant differences were further filtered based on greater than 1.5 fold changes to identify the distinct determinants associated with each of the four trajectories.

Among other factors, the infection-naïve PVS participants with quantifiable S1 had lower GHVAS scores indicative of poorer general health (GHVAS) and lower anti-S IgG titers, whereas higher circulating IL-7 and IL-21 levels were detected compared to other groups (Figure. 5H). Elevated growth hormone levels alongside low TSH levels were also observed among the PVS-I participants with S1 protein in circulation. By contrast, among the infection-positive subgroups, participants with circulating S1 were observed to have higher anti-N antibody titers based on the clinical COBAS assays and in-house ELISA indicative of the contribution of infection history. Moreover, anti-nucleosome IgA levels were higher among those in the PVS+I subgroup without detectable S1 (Figure. 5H).

Machine Learning-based identification of peripheral immune signatures of PVS

To establish a combined global immune signature for persistent symptoms following COVID-19 vaccination, we built machine learning models to predict PVS outcomes. The goal was to identify prominent features that could effectively distinguish PVS from the controls in a parsimonious manner. Given that autoantibodies are also common in the general population at low levels, we chose to exclude them from this analysis in the absence of further validation³⁰. Additionally, we excluded any variables with greater than 20% missing values for either the PVS or control groups and SERA variables because most of the dataset lacked a significant number of values above LLoD. A total of 193 variables were included.

Weighted Gene Co-Expression Network Analysis (WCGNA)³¹ was applied to this set to find groups of highly correlated variables (Figure S6B). The final feature set was then created by taking all variables that did not form a tight cluster (141, Module 1) and the eigengene, or first principal component, of every other set of variables (Modules 2-6).

This gave us a total of 146 features. Next, we performed classification utilizing Least Absolute Shrinkage and Selection Operator (LASSO) with nested cross-validation. We achieved an overall model accuracy of 78.1% on the validation folds (per-fold range 62.5% - 100%; Figure S6A), and an AUC of 0.80 (per-fold range 0.67-1.00, 95% CI = 0.67-0.92; Figures 6A and S6A). A permutation test further concluded that this performance was significantly above random (p = 0.02; Figure S6A). Segregating test data per-fold by infection status yielded accuracy on the infected population of 86.5% and the uninfected population of 73.3%. (Figure S6A) Per-class accuracy showed some divergence, with PVS accuracy of 85.7% vs. control accuracy of 63.6%. This was due to a subset of the controls clustering primarily with PVS samples, making control classification more difficult (Figure 6B).

• Download figure
• Open in new tab

Figure 6: Machine Learning Results and Prominent Feature Identification.

A. Confusion matrix inspired barplot describing actual and predicted labels for PVS and control with a classification threshold of 0.65. The p-value from a permutation test is also displayed to motivate model legitimacy. Boxplot displaying accuracy and AUC for each outer fold of nested-cross validation with both inner and outer loops set to 10. Mean and median accuracy and AUC stats are also shown. B. t-SNE plot showing a two-dimensional similarity-based representation of the LASSO variable space. C. Scatter plot containing all variables identified through either LASSO or marginal logistic regression models. Features only identified by marginal test, and not by LASSO are displayed in the center panel with coefficient of 0. Points with a shared color display Pearson correlation > 0.4. D. Correlation heatmaps of LASSO-selected and WGCNA module six features, with features clustered through hierarchical clustering.

The LASSO model selected 21 features using all data, consisting of CD4 T cell populations, immune modulators, neuropeptides, and antibodies (Figure 6C). Among the features selected, there were several negatively associated with PVS. These included circulating factors sIL-1R1, fetuin A36, granzyme A and B, FLT−3L and HMGB1, and subsets of circulating CD4 T cell populations (CXCR3⁺ CD4 T cells CD4⁺ T_EMRA cells, and IL-4⁺/IL-6⁺ CD4 T cells). Multiple hormones and neuropeptides synthesized by the hypothalamus, pituitary glands, and the peripheral nerves and involved in nociception and stress responses such as oxytocin, neurotensin, ꞵ endorphin, melanocyte-stimulating hormones (MSH), and substance P were also negatively associated with PVS and formed a single module (Module 6) (Figure 6D). The features that were positively associated with PVS were anti-EBV gp42 IgG titers, MMP1 levels, and TNFɑ⁺ CD8 T cells. We observed that no single variable or small subset of variables had a particularly strong differentiating power.