Association of circulating biomarkers with illness severity measures differentiates myalgic encephalomyelitis/chronic fatigue syndrome and post-COVID-19 condition: a prospective pilot cohort study

A unicenter, prospective, cross-sectional, pilot cohort study was conducted involving 31 ME/CFS patients (mean age: 49.3 ± 3.1 years; 65% female), 23 individuals with long COVID (mean age: 48.7 ± 2.4 years; 65% female), and 31 matched healthy sedentary controls (mean age: 41.7 ± 1.8 years; 71% female) recruited consecutively between September 2020 and December 2022 from the largest outpatient tertiary referral center in Spain (ME/CFS Clinical Unit, Vall d'Hebron University Hospital, Barcelona, Spain). After receiving verbal and written information on the study protocol, they all signed informed consent to participate prior to enrollment. The study protocol was approved by the Institutional Review Board of Vall d'Hebron University Hospital (reference number PR/AG 201/2016), and all procedures were conducted in accordance with the ethical standards of the board and the 1964 Declaration of Helsinki including its later amendments.

ME/CFS patients were eligible if they were aged ≥ 18 years and had a confirmed diagnosis by a specialist physician based on international consensus criteria (2011 ICC), the recommended case criteria for ME/CFS research purposes being those set out in an updated EUROMENE report [19].

COVID-19 "long haulers" (a.k.a. long COVID) were eligible for enrollment if they were aged ≥ 18 years and, following a confirmed diagnosis of acute COVID-19 infection based on a positive SARS-CoV-2 RT-PCR test on nasal swab during the COVID-19 pandemic, were still suffering from persistent, unexplained symptoms/signs ≥ 3 months after the acute COVID-19 infection, as established in the 2021 WHO clinical case definition for post-COVID-19 condition [6].

Healthy sedentary volunteers were eligible if they had neither experienced any self-reported autonomic symptoms nor had been in recent contact with anyone who was infected with SARS-CoV-2 (COVID-19) within 90 days prior to the study, no significant neurological, cardiac, endocrine or neuroimmune disorders, no alcohol or drug dependence, and did not use daily prescribed medications. All healthy sedentary subjects were recruited through word-of-mouth from the local community and did not meet the case criteria for either ME/CFS or long COVID at the time of enrollment.

All participants with cardiovascular dysautonomia were diagnosed by a physician based on consensus criteria from the POTS Working Group for the U.S. NIH [20]. All participants were of Caucasian descent, from the same geographical area, and had a sedentary lifestyle at the time of the study. They were subject to stringent exclusion criteria, as previously described by our group [21, 22]. The major exclusion criteria were a relevant previous or current diagnosis of an autoimmune disorder, multiple sclerosis, psychosis, major depression disorder, heart disease, hematological disorders, infectious diseases, sleep apnea or metabolic disorders; pregnancy or breast-feeding; smoking habit; strong hormone-related drugs; and preexisting fatigue-associated symptoms or evidence of multi-organ failure that did not meet the case criteria for ME/CFS and long COVID used in this study. Baseline demographic and clinical characteristics of the study population are displayed in Table 1.

All participants attended the aforementioned hospital between 8:00 a.m. and 11:00 a.m. for a face-to-face clinical assessment by a specialist physician. Baseline demographic and clinical characteristics of each participant were recorded based on self-reported outcome measures. A blood sample was taken from all participants for routine biochemical analysis and circulating biomarkers to assess endothelial dysfunction and systemic inflammation was assayed. All participants underwent the same cardiovascular orthostatic challenge protocol (10-min NASA lean test) to evaluate orthostatic intolerance.

Orthostatic intolerance was evaluated using a standardized passive standing test (10-min NASA lean test, NLT), a simple and well-established non-invasive orthostatic challenge protocol used to assess cardiovascular compensatory autonomic responses to standing. The NLT enables detection of orthostatic intolerance (OI) phenotypes such as orthostatic hypotension (OH) and postural orthostatic tachycardia syndrome (POTS) by measuring cardiovascular hemodynamic parameters (SBP/DBP and HR), and it is suitable for both clinical and research purposes [23].

Briefly, the NLT was conducted in a consistent manner by the same examiner in the morning between 8:30 a.m. and 11 a.m., in a quiet room with an average relative temperature of 22.3 ± 1.5 °C and humidity of 55.7 ± 4.8%. Participants were first asked to lie down on an exam table for 5 min and then to stand and lean against a wall, with their heels 15–20 cm away from the wall. An automated BP cuff with a monitor (Beurer BM-26, Beurer GmbH & Co., Ulm, Germany) was placed on the left arm, recording the systolic BP (SBP), diastolic (DBP), and heart rate (HR) at 1-min intervals. Baseline SBP/DBP HR were consecutively recorded twice during the supine position and every minute during the full 10 min after attaining the upright position. Throughout the recording, participants stood with only their shoulder blades touching the wall and their heels were positioned 15–20 cm from the wall. Participants were asked to remain still, and any talking or movement was discouraged, except for reporting any symptoms of concern. The NLT was stopped early at the request of the subject, or in the event of severe pre-syncope. After 10 min upright, each participant was asked about the frequency/severity and impact of orthostatic symptoms on a 5-item OGS score [23].

Criteria for OI were based on the 2021 expert consensus statement and guidelines for the definition of OH and POTS, as follows: (1) OH was defined as a decrease in SBP of ≥ 20 mmHg, or a decrease in DBP of ≥ 10 mmHg in the first 3 min standing, compared with resting supine values; and (2) POTS was defined as an increase in HR ≥ 30 bpm and/or a current HR ≥ 120 bpm without BP changes, based on the average of the final 3 min standing. Orthostatic intolerance (OH and POTS) during the 10-min NLT was quantified as the difference between supine and standing hemodynamic changes (ΔBP and ΔHR), as previously described [20, 24, 25]. For the purposes of this study, SBP/DBP and HR were used as maximum raw values recorded during the orthostatic test. For correlation analysis and group comparison, we calculated mean values during the final 2 min in the supine position (sp), mean values during the first 3 min standing (3 F), and mean values of final 3 min (3 L). Changes in these variables compared with values in the supine position (baseline) were also calculated (Δ3F or Δ3L).

Participants were asked to fill out a set of validated self-administered screening questionnaires regarding their current health status one week after the first clinical assessment. Changes in perceived fatigue (Fatigue Impact Scale, FIS-40), sleep quality (Pittsburgh Sleep Quality Index, PSQI), anxiety and depression symptoms (Hospital Anxiety and Depression Scale, HADS), autonomic dysfunction symptoms (Composite Autonomic Symptom Score-31, COMPASS-31), frequency and severity of orthostatic symptoms (Orthostatic Grading Scale, OGS), and health-related quality of life (Short-Form Health Survey, SF-36) were rated by each study participant under the supervision of two trained investigators (J.A.-M. and J.C.-M.), who oversaw compliance as described in our previous study [22].

Fatigue severity was assessed using the 40-item FIS-40. This questionnaire assesses three domains reflecting the perceived feeling of fatigue: physical (10 items), cognitive (10 items), and psychosocial functions (20 items). Each item is rated from 0 (no fatigue) to 4 (severe fatigue), and a total score is calculated by summing the item scores (range 0 to 160). Higher scores indicate greater functional limitations in daily life due to fatigue [26].

Sleep quality disturbances were evaluated with the 19-item PSQI, which assesses seven domains: subjective sleep quality, sleep latency, sleep duration, habitual sleep efficiency, sleep perturbations, use of sleeping medication, and daytime dysfunction. The overall PSQI score ranges from 0 to 21 points, with scores ≥ 5 indicating poorer sleep quality [27].

The HADS is a self-report scale used to screen for the presence of anxiety/depression symptoms in people with chronic illnesses. It consists of a 14-item inventory scored on a four-point Likert-type scale (from 0 to 3) and divided into two subscales, anxiety (7 items) and depression (7 items), which are scored independently. Each subscale score ranges from 0 to 21, and the higher the score, the greater the level of anxiety or depression [28].

To assess dysautonomia, all participants were screened using the COMPASS-31, a questionnaire designed to evaluate the frequency and severity of autonomic function symptoms in six core domains: orthostatic intolerance (four items), vasomotor (three items), secretomotor (four items), gastrointestinal (twelve items), bladder (three items), and pupillomotor symptoms (five items). The six domain scores are summed to provide a global COMPASS-31 score ranging from 0 (no symptoms) to 100 (worst symptoms). Higher scores indicate more severe autonomic complaints [29].

The orthostatic grading scale (OGS) is a validated 5-item self-report questionnaire designed to assess symptoms of orthostatic intolerance due to cardiovascular orthostatic dysfunction. The five questions address the frequency/severity of and interference due to orthostatic symptoms in daily life activities. Respondents rate each item on a scale of 0 to 4. Item scores are summed to give an overall OGS score ranging from 0 (no orthostatic symptoms) to 20 (worst orthostatic symptoms). Higher scores indicate greater severity of orthostatic intolerance [30].

The SF-36 is a generic scale that provides a health status profile, and it was used here to assess quality of life in study participants. Its 36 items explore eight dimensions of health status (physical function, role limitations due to physical health, bodily pain, general health, vitality, social functioning, emotional role, and mental health), and it yields summary scores for the physical and mental components [31].

After a 12-hour overnight fasting, 20 ml of peripheral whole blood were collected from each participant by venipuncture (from an antecubital vein with a 19-gauge needle) into SST-tubes and K₂EDTA-containing tubes (BD Vacutainer, Becton Dickenson, Sarstedt, Barcelona, Spain). This was performed between 8:00 a.m. and 10:00 a.m. by a trained phlebotomy nurse (USIC Outpatient Clinical Unit, Vall d'Hebron University Hospital, Barcelona, Spain).

One tube was transported and delivered to the local core laboratory within 2 h of collection for routine blood tests/analyses, following standard and recommended procedures. All other blood samples were immediately centrifuged at 2,500 rpm for 15 min at 4 °C followed by a second centrifugation step under the same conditions (Thermo Scientific, Waltham, MA, USA), after which serum and plasma specimens were collected and stored in aliquots at − 80 °C until further assays. No sample was thawed more than twice. Repeated samples from each participant were measured in the same analytical batch.

Blood levels of vascular endothelium proteins, namely soluble endothelin-1 (ET-1), thrombospondin-1 (TSP-1), vascular cell adhesion molecule-1 (VCAM-1), intracellular adhesion molecule-1 (ICAM-1), plasminogen activator inhibitor-1 (serpin E1/PAI-1), endothelial selectin, and adiponectin, as well as nitric oxide (NO_X assessed as NO₂^– + NO₃^–), were measured to assess peripheral vascular endothelium function. Plasma concentrations of each protein were assayed in each participant with commercially available ELISA kits, according to the manufacturer's instructions (Quantikine R&D Systems, Minneapolis, USA), using a Synergy™ H1M, hybrid multi-mode microplate reader (BioTek Instruments, Inc., Winooski, USA) at O.D. 450 nm. Results were analyzed by comparison with standard calibration curves in each well, and are presented as averages of two duplicated plasma samples.

Twelve cytokines/chemokines, namely IL-1β, IL-4, IL-6, IL-8, IL-10, IL-13, TNF-α, monocyte chemoattractant protein-1 (MCP-1/CCL2), interferon-gamma (IFN-γ), interferon gamma-induced protein-10 (IP-10/CXCL10), macrophage inflammatory protein-1-alpha (MIP-1α), and leptin were measured simultaneously in duplicate plasma samples using a fluorescently labeled microsphere-based multiplex bead-array immunoassay, according to the manufacturer's instructions (Cat# HCCBP1MAG-58 K and Cat# HCYTA-60 K, Milliplex Map Human Cytokines/Chemokines Magnetic Bead Panel I/II, Linco Research Millipore, Billerica, MA, USA). Measurements were taken on a Luminex-100 ISv2 reader (Linco Research Millipore, Billerica, MA, USA). For each cytokine/chemokine the standard curve ran from 3.2 to 10,000 pg/mL. Results were recorded in pg/mL based on standard curve values. The intra- and inter-assay coefficients of variation for each cytokine/chemokine were, respectively: IL-1β: 7% and 12%; IL-4: 3% and 11%; IL-6: 9% and 5%; IL-8: 7% and 5%; IL-10: 2% and 11%; IL-13: 2% and 11%; TNF-α: 8.7% and 4.9%; MCP-1: 2% and 11%; IFN-γ: 8% and 13%; IP-10: 7% and 11%; MIP-1α: 7% and 13%; and leptin: 5% and 7%.

Data for continuous variables are reported as mean ± standard error of the mean (SEM), while for categorical variables we report numbers and percentages. Due to a limited sample size and skewed distribution, statistical comparisons were performed using non-parametric methods. Specifically, Fisher's exact test was used for categorical variables, while for continuous variables we used either the Kruskal-Wallis test with Dunn´s post-hoc multiple comparison test or the Mann-Whitney U (rank sum) test for two groups. The Kolmogorov-Smirnov and Shapiro-Wilk tests were used to assess normality in the data distributions of the studied parameters. Box plots show mean ± SEM, including 25th and 75th percentiles of triplicate sampling. Statistical analyses were performed using GraphPad Prism software (GraphPad Prism 10.1.0.316 for Windows, serial number: GPS-1,359,963-L###-###; machine ID: E32AAD88B97; Boston, MA, USA). A two-tailed p-value ≤ 0.05 was considered statistically significant.

For the principal component analysis (PCA) we used the prcomp and autoplot function from the ggfortify library in R [32]. Autoplotting was used for both cluster analysis and representation of the different ellipses with coverage probability for a 95% confidence interval. The overlap between ellipses was calculated using the shipunov package in R, again with a 95% confidence interval. For the correlation heatmap we used the R package for multi-environment trial analysis (metan), selecting the Spearman's method and a cut-off p-value < 0.05 [33].

Finally, in order to differentiate between the groups, we conducted a MANOVA followed by discriminant analysis. The discriminant function analysis evaluates canonical discriminant functions based on combinations of the selected markers which contribute maximally to group separation, and it assesses how well these functions discriminate the diagnosis. Examining Wilk's lambda values for each of the predictors reveals how important the independent variable is to the discriminant function, with smaller values representing greater importance. Data were analyzed using IBM SPSS Statistics for Windows 27.0 (IBM Corp., Armonk, NY, USA).

Table 1 summarizes the baseline demographic and clinical characteristics of study participants recruited in this study. Patients (two groups) and controls did not differ in terms of age, gender or BMI. The mean duration of illness was 7 years for ME/CFS and 2 years for long COVID patients. As shown in Table 1, all patients with ME/CFS and long COVID presented with a wide range of persistent symptoms, with no hospitalizations during the acute COVID-19 infection. Thirty of the 31 ME/CFS individuals and 22 of the 23 long COVID subjects (96.8% and 95.7%, respectively), as well as all matched sedentary healthy controls, were vaccinated with at least one dose against COVID-19 using Comirnaty (Pfizer-BioNTech) ≤ 8 months prior to study inclusion. The ME/CFS and long COVID groups differed from healthy subjects in their scores on self-reported outcome measures, with patients reporting more fatigue, more anxiety/depression symptoms, poorer sleep quality, greater dysautonomia, and lower health-related quality of life (all p < 0.001).

Table 2 summarizes data for the cardiovascular autonomic response variables assessed through the 10-min NLT. It can be seen that SBP and DBP in the supine position were significantly higher in the ME/CFS and long COVID groups than in matched sedentary healthy controls (all p < 0.05). Regarding HR, this was higher in ME/CFS patients than in matched sedentary healthy controls (p < 0.005), but no differences were observed between long COVID and healthy individuals on this variable. No differences by group were found for changes in cardiovascular hemodynamic variables. None of the participants had abnormal pulse pressure (less than 25% SBP), and this variable did not differ between groups. Assessment of orthostatic intolerance using the 10-min NLT indicated that 4 ME/CFS patients (13%), 1 with long COVID (4%), and 1 healthy control (3%) presented POTS. No differences between the groups in terms of distribution were found (Fig. 1A-C).

The principal component analysis (PCA) depicted in Fig. 5 is based on datasets from clinical and biochemical variables. The analysis of outcome measures revealed that the first two principal components explained 88.7% of the total variance within the dataset (PC1: 82.7% and PC2: 6.1%, as shown in Fig. 5A). This analysis yielded two significant clusters, effectively distinguishing the patient groups (ME/CFS and long COVID) from the healthy control group. At the group level, patients (with ME/CFS or long COVID) displayed a notably more scattered distribution in the PCA plot compared with matched healthy controls, confirming the greater data variance. However, there seemed to be no distinctive difference between the patient groups (Fig. 5A). Indeed, the PCA plot showed substantial overlap (around 80%) between the two patient groups, rendering them indistinguishable in terms of clinical questionnaires.

As an additional and independent approach, we explored whether a combination of endothelial and/or inflammatory biomarkers could aid in distinguishing between patients with ME/CFS and those with long COVID. To assess the impact of phenotypes, we conducted PCA to maximize the separation between the ME/CFS and long COVID groups, using the results of the clinical questionnaires as a starting point and introducing the studied biomarkers (Fig. 5B-D). The first two principal components explained 55.65% (PC1: 42.3% and PC2: 13.4%; Fig. 5B) 48.11% (PC1: 34.6% and PC2: 13.5%; Fig. 5C), and 39.58% (PC1: 27.7% and PC2: 11.8%; Fig. 5D) of the total variance within the datasets for endothelial, inflammatory, and combined endothelial and inflammatory biomarkers, respectively.

In the annotated plot, there was minimal overlapping between the control group and the patient clusters. Within the patient subset, we found that a subgroup of long COVID patients could be distinguished from ME/CFS patients (Fig. 5B, C and D). Specifically, the analysis suggested that incorporating endothelial biomarker data could correctly distinguish 51% of patients with ME/CFS from those with long COVID, while under the same conditions the inclusion of inflammatory biomarkers only classified 35% of patients (Fig. 5C). Finally, the incorporation of both endothelial and inflammatory biomarkers improved the classification results to 59% (Fig. 5D).

In an attempt to identify the optimal panel that could serve as markers to differentiate ME/CFS patients from the long COVID and healthy control groups, we performed a discriminant function analysis including only those biomarkers that varied significantly among study groups (ET-1, TSP-1, VCAM-1, serpin E1 (PAI-1), E-selectin, NOx, IL-1β, IL-4, IL-6, IL-10, and TNF-α), as well as scores on the patient-reported outcome measures (FIS-40, HADS, COMPASS-31, PSQI, OGS, and SF-36).

Discriminant analysis revealed two canonical discriminant functions: the first explained 88.9% of the variance, canonical R² = 0.91, whereas the second only explained 11.2%, canonical R² = 0.57. In combination these discriminant functions significantly differentiate the groups (Wilks' lambda = 0.38; χ² = 239; p < 0.001). Overall, 95.3% of the original grouped cases (87.1% of ME/CFS patients, 100% of those with long COVID, and 100% of matched sedentary healthy controls) were correctly classified. If the same procedure was carried out with stepwise variable selection, the final model included serpin-E1 (PAI-1), NOx, IL-1β, IL-6, and FIS-40 (Wilks' lambda = 0.068; χ² = 215; p < 0.001) and 85.9% of participants were correctly classified (67.7% of ME/CFS patients, 91.3% of those with long COVID, and 100% of matched sedentary healthy controls).

The current study has some limitations. The small sample size and the fact that both ME/CFS and long COVID are multifactorial conditions characterized by their great heterogeneity and phenotypic complexity among participants may possibly have reduced the statistical power and the chance of detecting true consequences. This might especially concern the analyzes concerning the impact of the preliminary panel of putative circulating diagnostic biomarkers to distinguish ME/CFS from long COVID. Future studies are required to confirm the association between outcomes measures and circulating putative biomarkers in larger populations to increase the confidence and the statistical power of the analysis.

Furthermore, the limited diversity of potential covariates (confounders) in the available data reduced the number of possible factors of interest to adjust for in the correlation analyzes. In addition, the inability to report separate sex-specific association due to a low number of males in both the long COVID and ME/CFS groups. This is relevant because the disease burden of ME/CFS is higher in females, and recent studies have uncovered sex differences in its pathophysiology. The inclusion of both sexes in our sample nevertheless provides a comprehensive picture of disease characteristics, which may be regarded as a strength of the study. A further limitation is that we only report correlations, underscoring the need for further research to establish causal relationships between the variables considered. Similarly, our use of an observational cross-sectional design and the differing duration of illness in the two patient groups (7 years for ME/CFS and 2 years for long COVID) mean that more comprehensive longitudinal research is required to achieve a better understanding of disease progression in the two scenarios. Indeed, further research is crucial for unraveling the pathomechanism of endothelial dysfunction and inflammation and their role in classifying these two distinct conditions.

Another limitation of the study is that illness severity were assessed through self-report questionnaires and reflect the week prior to follow-up assessments. Moreover, the use of self-reported outcome measures to assess the natural course of illness (symptoms) is prone to recall bias and over-reporting. Furthermore, we measured outcomes that were reported in structured medical coding in electronic health records and we had no access to diagnoses and outcomes reported in free text format. These data may not, therefore, completely reflect diagnoses and reported outcomes from study population. Although we cannot rule out diagnostic errors relating to conditions with similar health outcomes, we believe that these are equally likely in both groups, in other words, misclassification of outcomes is non-differential. Patient-reported outcomes such as weakness, cognitive impairment, anosmia, and disgeusia are also less objective than are clinical diagnoses by physicians and might not be uniform and accurate.

Importantly, we cannot rule out potential behavioral and environmental factor differences between infected and uninfected people, which might cause overestimation of incidence among the infected population. Neither can we exclude the possibility of additional confounders affecting long-term outcomes of acute SARS-CoV-2 infection that were unavailable to us in this study. Furthermore, some outcomes were reported at low frequency and larger populations might be necessary to gain a more reliable picture.