This cross-sectional study (IMMUNOFATIGUE) was carried out at Erasmus Medical Center (MC) Rotterdam, the Netherlands. The study was performed within a prospective cohort study (CO-FLOW) on long-term outcomes in adult patients who had been hospitalized for COVID-19 in the Netherlands (36). Patients were eligible to participate in the IMMUNOFATIGUE study if they visited the outpatient clinic at Erasmus MC between 3 and 6 months after hospital discharge for persistent COVID-19 sequelae. In the Netherlands, it is standard practice to offer post-discharge follow-up to COVID-19 patients by the discharging hospital. Study blood samples were taken as part of this follow-up. Demographics and clinical characteristics at hospital admission were retrospectively collected from electronic medical records and via a questionnaire. We reported comorbidity as this was either reported in the medical records or self-reported.

A group of age- and sex-matched healthy controls (HCs) without a history of SARS-CoV-2 infection (self-reported) was recruited among hospital visitors. HCs were asked whether they had been vaccinated against COVID-19 and were screened for fatigue and depression using the Fatigue Assessment Scale (FAS) and Hospital Anxiety and Depression Scale (HADS) questionnaires. HCs who showed indications of chronic depression or fatigue, based on established cut-off scores that are described below, were excluded. All other HCs were included in the study, as they were appropriately age- and sex-matched to the fatigued long COVID group (no further selection necessary). Information on the sample size calculation is provided in the Supplementary Methods . All patients and HCs provided written informed consent before the start of study measurements. The Medical Ethics Committee of Erasmus MC approved the study (MEC-2020-1893).

PROMs were collected in all patients as part of the CO-FLOW study at 3, 6, and 12 months post-hospital discharge (36). For the current IMMUNOFATIGUE study, the Beck’s Depression Inventory (BDI-21) was added.

Fatigue was assessed with the FAS questionnaire, assessing physical and mental fatigue, which was initially validated in patients with sarcoidosis but also used in other diseases (37, 38). The FAS consists of 10-items rated on a 5-points Likert scale. The total FAS score ranges from 10 to 50, with higher scores indicating increased severity of fatigue. We used a ≥22 cut-off value to indicate substantial fatigue, categorizing patients into fatigued and non-fatigued groups. This value is derived from studies examining sarcoidosis-related fatigue and is also used to indicate fatigue in patients with other conditions, including long COVID (4, 39, 40). Patients filled out a symptom questionnaire (Corona Symptom Checklist [CSC], Supplementary Methods ) to assess the presence of new or worsened symptoms following acute SARS-CoV-2 infection on a binary scale (yes or no). In this study, we included a selection of 12 typical long COVID symptoms. Several symptoms (such as chest pain, sensory overload, and headache) could not be taken into account as they were added to this questionnaire in a later stage and contain incomplete data (5). Dyspnea was assessed with the Modified Medical Research Council (mMRC) Dyspnea Scale (41), the questionnaire scales the severity of dyspnea from 0 (no dyspnea) to 4 (severe dyspnea); scores ≥2 were considered representative for the presence of dyspnea symptom. Anxiety and depression were assessed with the HADS and BDI-21. The HADS consists of the subscales anxiety and depression, a subscale score (range 0-21) ≥11 is considered clinically significant (42). Depressive symptoms were also assessed with the BDI-21, scores (range 0-63) 0-13 denote no/minimal signs of depression, 14-19 mild depression, 20-28 moderate depression, and 29-63 severe depression (43, 44). Post- Traumatic Stress Disorder (PTSD) was measured with the Impact of Event Scale-Revised (IES-R) (45), with scores (range 0-88) ≥33 representing clinically significant PTSD. Cognitive failures in everyday life were assessed with the Cognitive Failure Questionnaire, with scores (range 0-100) >43 indicating cognitive failure (46, 47). Health-Related Quality of Life (HR-QoL) was assessed with the 36-item Short Form survey (SF-36) (48). The SF-36 measures general health status and consists of eight domains: physical functioning, social functioning, physical role functioning, emotional role functioning, mental health, vitality, bodily pain, and general health perception. Each domain score ranges 0-100, with lower scores indicating more disability.

Serum and sodium heparinized peripheral blood samples were collected from patients during post-discharge follow-up (March to October 2021) and from HCs between July and November 2021. Peripheral blood mononuclear cells (PBMC) were isolated by low-density gradient centrifugation using Ficoll-density separation shortly after blood draw to avoid erythrophagy-related activation of the monocytes. Isolated PBMC were frozen in 10%-dimetylsulfoxide and stored in liquid nitrogen.

Staining A, the absolute counts of total leukocytes (CD45⁺), Natural Killer (NK) cells (CD3^-CD16⁺CD56⁺), B-lymphocytes (CD19⁺), T-lymphocytes (CD3⁺), CD4⁺ T-lymphocytes, and CD8⁺ T-lymphocytes were determined with a clinical laboratory ISO 15189 accredited method. Staining was performed on whole blood using multitest 6-color T-B-NK reagent (BD Biosciences). Flowcytometric analyses were performed using a FACS Canto II instrument (BD Biosciences).

Staining B, the percentages of CD4⁺ T helper lymphocyte subsets were determined using a 8-color (membrane and intracellular) staining. PBMC were thawed and the recovery and viability of cells were 68% and 97% respectively, as assessed by Trypan blue staining. A total of 1 x 10⁶ of defrosted PBMCs were stimulated for 4 h at 37 °C in RPMI-1640 culture medium with 50 ng/ml phorbol 12-myristrate 13-acetate (PMA; Sigma Aldrich, St. Louis, MO, USA) and 1.0 μg/ml ionomycin (Sigma) in the presence of Golgistop (BD Biosciences). CD4⁺ T-lymphocyte subsets were identified by their secreting cytokines: T helper (Th)1 (CD3⁺CD4⁺IFNγ⁺), Th2 (CD3⁺CD4⁺IL4⁺), Th17 (CD3⁺CD4⁺IL17A⁺). Regulatory T-lymphocytes (T_reg) were identified by their transcription factor FOXP3 (CD3⁺CD4⁺CD25highFOXP3⁺). The percentages of these CD4⁺ T-lymphocyte subsets are expressed as percentage of total lymphocytes.

Staining C, a more in-depth analysis of CD4⁺ and CD8⁺ T-lymphocyte subsets was performed using a second vial of PBMCs that had been thawed. After thawing, the recovery and viability of cells were 69% and 97% respectively, as assessed by Trypan blue staining. A total of 1,5 - 2 x 10⁶ PBMCs were stained with a cocktail of CD45-V500, CD45RA-BB515, CD3-Alexa Fluor700, CD4-BUV805, CD8-BUV395, CD197-BV421, CD279-RB780 (BD Biosciences), CD28-BV711, CD27-APC and CD57-PE (BioLegend) for 15 minutes at room temperature, washed twice with PBS, pH 7.8 and subsequently stained with viability dye Zombie NIR (BioLegend). 500.000 events in a live/CD45 stopping gate were collected on an Aurora-5 laser instrument. The analysis was performed using FlowJo software. Quadrant gating strategy on CD45RA and CD197 (CCR7) was used to define the CD4⁺ and CD8⁺ T-lymphocyte subsets: naïve-like (CD45RA⁺CD197⁺), central memory (T_CM, CD45RA^-CD197⁺), effector memory (T_EM, CD45RA^-CD197^-) and effector memory RA (T_EMRA, CD45RA⁺CD197^-) ( Supplementary Figure S1 ). Within each indicated T-lymphocyte subset, the expressions of CD27 and CD28 as well as CD279 and CD57 were assessed. The percentages of the T-lymphocyte subsets in staining C are presented as percentages of total T-lymphocytes. The monoclonal antibodies used in each staining are in more detail described in Supplementary Table S1 .

Staining D, the percentages of classical monocytes, intermediate monocytes, and non-classical monocytes were determined using a 5-color membrane staining. Thawed PBMC used for staining C (see above) were also used for further monocyte subset analysis. A total of 1 x 10⁶ PBMCs were stained with a cocktail of CD45-PO (Life Technologies), CD64-APC (BD Pharmingen), CD66b-BV421 (BD Horizon), CD14-APC-H7 (BD Biosciences), CD16-PE-Cy7 (Invitrogen) for 30 minutes at RT, lysed with BD Lysing solution (BD Biosciences) for 10 min at RT and washed with PBS pH 7.8. Monocytes were measured using a BD FACS Lyric flowcytometer. The analysis was performed using BD FACSSuite software. Monocytes (CD45⁺CD64⁺CD66^-) were defined as: classical monocytes (CD14⁺⁺CD16^-), intermediate monocytes (CD14⁺⁺CD16⁺), and non-classical monocytes (CD14⁺CD16⁺⁺) (49). For the monoclonal antibodies used in staining D, see Supplementary Table S1 .

We assessed 34 genes of a previously established inflammation-related gene signature found in mood-disorder patients (15, 16, 18, 32, 33) and patients with ME/CFS (unpublished data). Expression of inflammation-related genes in monocytes were assessed using procedures that have been described in previous publications (18). Briefly, CD14⁺ monocytes were isolated from thawed PBMCs by magnetic cell sorting (Automacs Pro Miltenyi Biotec, Bergisch Gladbach, Germany) and RNA was isolated (Qiagen RNeasy mini kit). The average viability and the purity of monocytes were both 95% as determined by Trypan blue staining and flow cytometry. Subsequently, RNA (0,5 µg) was reverse transcribed (High Capacity cDNA Reverse Transcription Kit; Applied Biosystems, Thermo Fisher Scientific) to obtain cDNA for quantitative-polymerase chain reaction (q-PCR). qPCR was performed using TaqMan Gene expression assays (Applied Biosystems). The expression levels of genes were determined using the comparative cycle (CT) method. All values were normalized to the housekeeping gene ABL1 (ΔCT values). Gene expression values of patients were also expressed relative to the average ΔCT value of HCs (ΔΔCT values). The following genes were evaluated: ABCA1, ABCG1, ADM, BAX, BCL10, BCL2A1, CCL2, CCL20, CCL7, CXCL2, DUSP2, EGR1, EGR2, EMP1, HMOX1, IFI44, IFI44L, IFIT1, IFIT3, IL1A, IL1B, IL1R1, IL6, MAFF, MAPK6, MRC1, MVK, MX1, MXD1, NR1H3, PTX3, SERPINB2, TNF, and TNFAIP3 ( Supplementary Table S2 ).

The following cytokines and soluble cell surface molecules were measured in serum with a Luminex multiplex bead-based assay (R&D Systems Europe, Abingdon, United Kingdom): brain-derived neurotrophic factor (BDNF), C-C motif chemokine ligand (CCL)2, CCL7, C-X-C motif chemokine ligand (CXCL)9 and CXCL10, cluster of differentiation 163 (CD163), Galectin-9, granulocyte macrophage-colony stimulating factor (GM-CSF), interferon (IFN)-α, IFN-β, IFN-γ, interleukin (IL)-6, IL-7, IL-10, IL-12, P-selectin, T-cell immunoglobulin and mucin domain 1 (TIM-1), and tumor necrosis factor-alpha (TNF-α). IL-6 levels were determined by a high sensitivity ELISA (apDia, Turnhout, Belgium) and serine protease inhibitor B2 (SERPINB2) levels by ELISA (R&D Systems Europe). All assays were performed following the manufacturer’s protocol.

Active EBV and CMV were determined in randomly selected patients from the fatigued and non-fatigued long COVID groups to assess whether symptoms could be attributed to viral reactivation, as suggested in the literature (50). EBV and CMV DNA load was measured using internally controlled quantitative real-time Taqman PCR based on assays performed as published previously (51, 52). For EBV a value over the lower limit of quantification of >100 IU/mL indicated the presence of active virus, and for CMV this was >50 IU/mL.

Continuous variables are presented as the median and interquartile range (IQR) and categorical variables as a number and percentage. Differences in demographics and clinical characteristics at hospital admission between groups of fatigued and non-fatigued long COVID patients were assessed using the Mann-Whitney U test for continuous variables and a Fisher’s Exact test for dichotomous categorical variables. The number of clinical symptoms was calculated using 14 typical long COVID symptoms, 12 symptoms from the CSC and the symptoms fatigue (FAS) and dyspnea (mMRC Dyspnea Scale).

For gene expression levels in monocytes, we first identified clusters of mutually correlating genes. Hierarchical cluster analysis of gene expression levels in monocytes was performed using Spearman’s rank correlation coefficient matrix. Missing gene values (0.7% of COVID-19 patients and 1.4% of HCs) were imputed with the median of patient or HC value. For single gene expressions we used ΔΔCT values and p-values were calculated with Wilcoxon signed rank test using the Benjamini-Hochberg-method for multiple testing. In the serum analysis, cytokines and soluble cell surface molecules positive in >20% of patients were used in the analysis ( Supplementary Table S3 ). A chi-square test was performed to assess differences in the number of patients with and without detectable cytokine and soluble cell surface molecule levels across groups of fatigued and non-fatigued long COVID patients and HCs ( Supplementary Table S3 ). Values below the lower limit of detection (LOD) were imputed by half of the lowest value observed of a given cytokine and extreme outliers (TNF-α n=1, IFN-β n=1, IFN-γ n=1, CXCL9 n=2, SERPINB2 n=2) were removed due to potential erratic measurements. We performed a Mann-Whitney U test to compare immune characteristics between the groups of long COVID patients and HCs. To compare immune features of fatigued and non-fatigued long COVID patients and HCs, we performed a Kruskal-Wallis test followed by a post-hoc test using Bonferroni correction for multiple group comparisons. In addition, for cytokines with imputed data, data were also categorized (<LOD vs. ≥LOD) and compared across groups using a chi-square test.

We also analyzed the severity of fatigue (total FAS score) as a continuous outcome in relation to immune characteristics. In preliminary analyses, we considered age, sex, BMI, pre-existing diabetes, the number of days in the hospital, and the number of days between SARS-CoV-2 infection and follow-up as factors potentially associated with fatigue severity. The association between continuous variables and fatigue severity was assessed using Spearman’s rank correlation coefficient, while differences in fatigue severity for categorical variables were assessed with the Mann-Whitney U test; none of these variables were significantly associated with the total FAS score and were therefore not included in further analyses as potential confounders. The associations between fatigue severity and immune parameters were assessed using Spearman’s rank correlation coefficient. A p-value <0.05 was considered statistically significant. Statistical analyses were performed with IBM SPSS Statistics version 28 (SPSS Inc., Chicago, IL, USA) and R software version 1.4.2.

We included 37 long COVID patients who experienced substantial fatigue (total FAS score ≥22) at the time of collecting blood samples, from here on referred to as fatigued long COVID patients. As a contrast group, we included 36 long COVID patients who did not experience substantial fatigue (total FAS score <22, non-fatigued long COVID), representing patients with mild or no symptoms of fatigue. We included a group of 42 age- and sex-matched HCs; the characteristics of HCs’ are shown in Supplementary Table S4 . FAS items scores in fatigued and non-fatigued long COVID groups and HCs are presented in Supplementary Table S5 . All patients had been discharged from hospital between October 2020 and May 2021, representing patients with SARS-CoV-2 alpha variant. The median follow-up time (day of blood sampling) since hospital discharge was 107.0 (IQR 92.5-138.5) days.

The patient’s demographic and clinical characteristics are shown in Supplementary Table S6 . Among the fatigued long COVID patients, the median age was 58.0 (55.0-66.0) years, 24 (64.9%) were male, 18 (48.6%) had been treated in the intensive care unit (ICU) for COVID-19, and the length of stay (LOS) in hospital was 17.0 (9.0-26.0) days. In non-fatigued long COVID patients, the median age was 61.0 (52.3-67.0) years, 26 (72.2%) were male, 22 (61.1%) had been treated in the ICU for COVID-19, and the LOS in hospital was 15.0 (10.0-26.8) days. Demographic and clinical characteristics at hospital admission did not differ significantly between fatigued and non-fatigued long COVID patients ( Supplementary Table S6 ). The median age of HCs was 62.0 (51.8-68.3) years (p=0.95, comparison with fatigued long COVID group) and 26 (61.9%) were male (p=0.82).

The PROMs in the fatigued and non-fatigued long COVID patients at the time of collecting blood samples are presented in Table 1 . In the fatigued long COVID group, low energy levels throughout the day and difficulties with concentration were prominent features of the fatigue ( Supplementary Table S5 ). Regarding the main symptoms, collectively, all fatigued long COVID patients reported ≥3 symptoms. Other PROMs showed that 40.5% of these patients experienced cognitive failure and signs of depression were found via HADS-D score in 24.3% and via BDI-21 score in 37.2% of the patients. The fatigued long COVID patients experienced significantly more symptoms and reduced HRQoL outcomes as compared to non-fatigued long COVID patients ( Table 1 ). In the non-fatigued long COVID group, 21 (65.6%) patients reported ≥3 symptoms and 3 (9.1%) patients experienced cognitive failure, signs of depression were not found.

Figure 1 presents the recovery profile of the main symptoms during the first year after hospital discharge. The fatigued patients showed hardly any clinical improvement over time, as 95.0% of the patients reported ≥3 symptoms at one year follow-up. Symptoms were overall less prevalent in non-fatigued patients as compared to the fatigued patients. However, still 55.6% of non-fatigued patients reported ≥3 symptoms at one year follow-up.

We performed enumerations of circulating leukocytes, NK cells, B-lymphocytes, T-lymphocytes, and CD4⁺ and CD8⁺ T-lymphocytes; Supplementary Table S7 shows the outcomes per group. Counts of these sets of cells did not show significant differences between fatigued and non-fatigued long COVID groups, only in comparison to healthy controls (HCs). When analyzed as a single group, the group of long COVID patients showed significantly increased counts of leukocytes and total T-lymphocytes compared to HCs; these increases were due to increases in CD8⁺ T-lymphocyte counts ( Figure 2A ). Considering fatigue as a graded outcome, we found that increased counts of CD8⁺ T-lymphocytes significantly correlated with increased fatigue severity (total FAS score, r=0.24, p=0.043) ( Figure 2A ).

In a more in-depth analysis of the CD4⁺ and CD8⁺ T-lymphocyte subsets (see Supplementary Table S8 for the outcomes per group), the percentages of these subsets did not differ significantly between fatigued and non-fatigued long COVID groups. We found a significant increase in the percentage of CD8⁺ T_EMRA-lymphocytes and the subsets of late stage/exhausted CD279⁺CD57⁺ CD8⁺ T_EMRA-lymphocytes and CD27^-CD28^- CD8⁺ T_EMRA-lymphocytes in the entire group of long COVID patients compared to HCs ( Figure 2B ). Percentages of naïve CD4⁺ T-lymphocytes (CD45RA⁺CCR7^-) and CD4⁺ T_reg-lymphocytes (CD25⁺FOXP3⁺) were significantly decreased in long COVID patients compared to HCs, while the percentage of exhausted/senescent CD4⁺ T_EM-lymphocytes (CD279⁺CD57⁺) were increased ( Figure 2C ). There were no significant correlations found between the percentages of the various CD4⁺ and CD8⁺ T-lymphocyte subsets and the fatigue severity in long COVID patients.

Figure 3A presents the percentages of classical (CD14⁺CD16^-), intermediate (CD14⁺⁺CD16⁺), and non-classical (CD14⁺CD16⁺⁺) monocytes across groups. The percentages of these populations did not differ significantly between fatigued and non-fatigued long COVID groups ( Supplementary Table S9 ). The entire group of long COVID patients showed a significantly reduced percentage of classical monocytes and an increased percentage of non-classical monocytes compared to HCs ( Figure 3A ). We found a significant negative correlation between the percentage of classical monocytes and the fatigue severity (r=-0.28, p=0.02), as well as a significant positive correlation between the percentages of intermediate (r=0.28, p=0.02) and non-classical (r=0.31, p=0.009) monocytes and the fatigue severity ( Figure 3B ).

A gene expression analysis was performed in the total population of CD14⁺ monocytes. Hierarchical clustering of gene expression levels in monocytes revealed three main gene clusters ( Figure 4A ), similar to the gene clusters found in previous studies in major depressive disorder (MDD), bipolar disorder, and various autoimmune disorders (thyroid autoimmune disease, type 1 diabetes, Sjögren disease and SLE) (16, 18, 53, 54). These clusters represent strong mutually correlating genes within each cluster; only cluster C genes correlated weaker amongst themselves. Gene cluster A was composed of inflammation-regulation genes and genes related to adhesion, chemotaxis, apoptosis, and pyroptotic mechanisms. Cluster B consisted of type 1 IFN driven inflammation-related genes. Cluster C consisted of genes involved in mitochondrial anti-inflammatory action and cholesterol pump genes.

Figure 4B shows the gene expression pattern in monocytes for the fatigued and non-fatigued patients relatively to the expression levels of HCs; differences between fatigued and non-fatigued patients did not reach significance and are therefore not shown. Both long COVID groups were characterized by a significantly overexpression of many cluster A inflammation-regulating genes (e.g. CCL7, CCL20, IL-6) as well as some cluster C genes (BAX, HMOX1) as compared to HCs. The cholesterol pump genes (ABCA1, ABCG1, NR1H3) and the M2 macrophage marker MRC1 were significantly downregulated. Normal expression levels were found for the type 1 IFN induced genes (ISGs) in cluster B. This profile represents a strong pro-inflammatory pyrogenic state of the monocytes.

Upregulated levels of the cluster A inflammatory genes CCL2 (r=0.29, p=0.016), CCL7 (r=0.24, p=0.048), and SERPINB2 (r=0.28, p=0.022) in monocytes were significantly correlated to increased fatigue severity ( Figure 4C ); significant correlations were not found for the other genes.

The level of the various tested inflammation regulating cytokines and soluble cell surface molecules in serum was evaluated to further investigate the inflammatory state of long COVID patients. These levels did not differ significantly between the fatigued and non-fatigued long COVID groups, both groups did show altered levels in comparison to HCs ( Supplementary Table S10 ). The entire long COVID group showed significantly increased serum levels of Galectin-9, IL-6, TNF-α, CXCL10, CD163, and CCL2 compared to HCs ( Figure 5A ). Levels of CXCL9, SERPINB2, IFN-β, and IFN-γ were significantly reduced in long COVID patients compared to HCs ( Figure 5A ).

In terms of fatigue severity, a significant positive correlation was found between serum Galectin-9 levels and the fatigue severity (r=0.24, p=0.039) ( Figure 5B ), a trend toward significance was found for CD163 levels (r=0.21, p=0.078), but not for the other measured inflammatory mediators.

We randomly selected 10 (27.0%) fatigued and 19 (52.8%) non-fatigued long COVID patients to test for the viral load of EBV and CMV. None of the tested patients showed viral loads that exceeded the limit of quantification, and we therefore did not perform further tests on the remaining patients.