Gene activity in common blood cells shows lower inflammation in classical monocytes of long COVID patients

A total of 44 patients with long COVID were recruited from the outpatient clinic of Rostock University Medical Center. Selection was based on the order of their visit, and diagnosis followed the WHO definition of "Post COVID-19 condition (long COVID)" (5, 29). The symptoms, disease progression, comorbidities, previous treatments, and quality of life are summarized in Figure 1A. Health-related quality of life (HRQoL) was assessed using the Short-Form 12 (SF-12) (30) and the Bell Disability Score (31). Routine laboratory tests included complete blood counts, coagulation parameters, C-reactive protein (CRP), liver enzymes, renal parameters, vitamin B12, folate, 25-OH-vitamin D, and thyroid-stimulating hormone (TSH), as well as basic autoimmune diagnostics. Differential diagnoses were excluded by the consulting cardiologists and pneumologists. A total of 44 age- and sex-matched controls were recruited from a monitored healthcare worker cohort (32 –34). Selection was based on primary messenger RNA (mRNA) vaccination and on optimal matching with patients. Of the controls, 43 were convalescent from SARS-CoV-2 infection and one was uninfected.

This study was approved by the Ethics Committee of the Rostock University Medical Center under file no. A 2020-0086 and no. A 2023-0081. Written informed consent was provided by all participants.

Venous blood was collected into EDTA vacutainers. PBMCs were isolated using SepMates (StemCell, Vancouver, BC, Canada) with Ficoll-Paque PLUS (Cytiva, Marlborough, MA, USA). Aliquots of 1–3 × 10⁶ cells each were stored at −80°C until further processing.

After thawing, 180,000 PBMCs were labeled with BD Human Multiplexing Sample Tags (Becton Dickinson, Franklin Lakes, NJ, USA) for the multiplexing of 12 samples per run. The viability of PBMCs was assessed via staining with Calcein and DRAQ7 (Thermo Fisher Scientific, Waltham, MA, USA). RNA was harvested using the Single-Cell Capture and cDNA Synthesis protocol with the BD Rhapsody Single-Cell Analysis System, and a total of eight runs were performed to enable single-cell capture from all patients and controls. Each run contained a balanced mixture of patient and control samples plus a spike control (isolated CD3⁺ T cells from an independent donor) to monitor potential batch effects. A total of 40,000–60,000 cells were pooled and loaded onto the Rhapsody cartridge. The incubation times for settlement of the cells varied between 15 and 25 min. The RNA libraries for scRNA-seq were prepared using the Library Preparation Protocol for mRNA Whole Transcriptome Analysis and Sample Tag according to the manufacturer's instructions and were purified using AMPure beads (Beckman Coulter, Brea, CA, USA). The libraries were indexed using a run-specific primer and were quantified and quality-checked using a Qubit 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 4200 TapeStation System with the Agilent HS D5000 screen tape assay (Agilent Technologies, Santa Clara, CA, USA). Sequencing was performed by Novogene Europe (Cambridge, UK, Munich subsidiary) using Illumina technology.

Raw scRNA-seq reads were aligned to the human genome and annotated by cell type using the BD Rhapsody™ Sequence Analysis Pipeline 2.2.1, complemented by manual assignment as previously described in Hillman et al. (35). Subsequent analysis was conducted using the Seurat 5.1.0 package in R, version 4.4.1 (R Core Team 2024, Vienna, Austria) (36, 37).

Quality control procedures involved the exclusion of cells exhibiting either an excessive or an insufficient number of unique RNA products or more than 25% of mitochondrial genes. A total of 132,230 PBMCs could successfully be recovered (44 patients and 44 controls). Of these, 118,843 passed quality control. Cartridge batch effects were assessed using a cell mix score calculated using the CellMixS 1.7.1 package for R for the CD3⁺ spike cells, the aim of which was to statistically validate the visual impression of homogeneous distribution provided by the UMAP (uniform manifold approximation and projection). This score is based on the calculation of a p-value for each cell, which represents the probability of belonging to a common underlying population (38). Subsequent analysis involved the normalization and scaling of the reads.

Cytokines were measured using a custom-designed Luminex^® Discovery Assay panel (R&D Systems, Minneapolis, MN, USA), with selection based on the identified differentially expressed genes (DEGs). The assay was performed according to the manufacturer's protocol using undiluted plasma samples. Data acquisition was carried out using the Luminex^® 100/200™ System, with calculations performed using xPONENT^® 3.1 software. The lowest manufacturer standard was considered as the lower limit of detection (LLOD). No extrapolation was performed.

Cryopreserved plasma samples from both patients and controls were analyzed using an enzyme-linked immunosorbent assay (ELISA) specific for immunoglobulin G (IgG) antibodies against the nucleocapsid protein of SARS-CoV-2 [anti-SARS-CoV-2-NCP-ELISA (IgG)] (EUROIMMUN Medizinische Labordiagnostika AG, Lübeck, Germany) following the manufacturer's protocol. The samples were centrifuged at 10,000 × g for 5 min prior to 101-fold dilution. Photometric measurements were obtained at 450 nm with a reference wavelength of 620–650 nm using the Infinite M200 spectrophotometer (Tecan, Männeheim, Switzerland). The results were calculated as the ratio of the sample extinction to that of the calibrator.

DEGs were determined using a non-parametric Wilcoxon rank-sum test and were classified based on a fold change greater than 3 and a (Bonferroni-corrected) p-value less than 10⁻¹⁸. UMAPs were constructed using nearest-neighbor clustering of all the genes expressed. Spearman's rank correlation coefficient was employed for correlation analyses. Heatmaps were constructed using a column/row-wise clustering approach with complete linkage and Euclidean distance metrics. Gaussian distribution was assessed using Kolmogorov–Smirnov tests. The Mann–Whitney U test was employed for the comparison of data that did not adhere to a normal distribution. Statistical analysis was performed using SPSS version 29.0.1.1 (IBM, Armonk, NY, USA) and GraphPad Instat version 3.1 (GraphPad Software, Boston, MA, USA).

Estimation of transcription factor (TF) activity was performed using the DoRothEA R package version 1.16.0, a gene regulatory network comprising signed interactions between TFs and their target genes (39 –41). The VIPER algorithm (42) was employed to estimate the activity of the TFs from the expression data of all monocytes of all genes. Group comparisons were made based on the mean differences.

The raw reads of all genes were divided into spliced and unspliced read matrices before estimation of the RNA velocities using gene-relative slopes. The velocity graph was visualized on the previously calculated UMAP of all monocytes embedding, employing a correlation-based transition probability matrix within the k-nearest neighbor graph. All analyses were performed using velocyto.R and velocyto.py (43).

Pathway analysis was conducted with ShinyGO software (44), which compared all downregulated genes in all monocytes of the long COVID cohort with the Reactome 2022 (45) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) 2021 (46 –48) databases. Gene set enrichment analyses (GSEA) were performed by alignment of the normalized counts of all monocytes with the Hallmark Gene Sets of Molecular Signatures Database (49 –51). Enriched pathways were identified using the normalized enrichment score (NES) and the false discovery rate (FDR).

We analyzed the PBMCs from 44 individuals diagnosed with long COVID using scRNA-seq, with the detailed clinical characteristics, comorbidities, and treatments summarized in Table 1. Majority of the patients reported mild (n = 27) or moderate (n = 13) initial SARS-CoV2 infections (29). Quality of life, assessed in 21 and 16 long COVID patients using the Short-Form 12 (30) and the Bell score (31), respectively, was significantly reduced compared with that of standardized populations [MCS (mental component summary): t (20) = −13.17, p < 0.001; PCS (physical component summary): t (20) = −15.33, p < 0.001]. Patients were grouped according to symptoms such as dizziness, pain, cognitive or cardiorespiratory constraints, and fatigue. The number of patients in each symptom combination group is presented in Figure 1A. The severity of the symptoms and post-exertional malaise, classified as extreme exhaustion after minimal stress, were mentioned, when available. Controls (n = 44) were age- and sex-matched, mostly convalescent (n = 43) plus one naive to SARS-CoV-2 (n = 1), as determined based on the anti-nucleocapsid antibodies. Both groups received comparable numbers of SARS-CoV2 vaccinations and had received mRNA vaccines only (Table 1). The controls, mainly healthcare workers, were primed and boosted at the turnaround 2020/2021, followed by the patients in the spring/summer of 2021. First SARS-CoV-2 infections occurred mainly in 2022, suggesting Omicron as the predominant variant in all study participants (52). Figure 1B summarizes, for both cohorts, the time intervals between blood sampling and completion of the primary immunization, as well as between blood sampling and first infections.

At the time of blood collection between July 2023 and February 2024, the controls, as part of a vaccination study (32), had been regularly assessed for anti-nucleocapsid antibodies, revealing a median of two infections each. Regarding the patients, the infection that triggered long COVID had been documented; thereafter, any reliable information on additional SARS-CoV-2 infections was sketchy. We therefore quantified the patients' anti-nucleocapsid IgG antibodies at the time of blood collection and found evidence of recent infections in some patients, but also confirmed fewer and/or longer ago infections for the majority, resulting in a significantly lower median for the anti-nucleocapsid titers in patients compared with the controls (Table 1).

Quality control of the scRNA-seq included the exclusion of cells with aberrant counts of unique RNA products (<200 and >4,000) and cells expressing more than 25% of mitochondrial genes. Assessment of batch effects and calculation of the cell mixing scores of the spike controls from each experimental run confirmed a high degree of comparability between all experimental runs, justifying the combined processing of all data (). 1

Using the BD Rhapsody™ Sequence Analysis Pipeline on the basis of the highly expressed genes, in conjunction with a manual assignment of monocytes according to the expression of CD14 and CD16, as previously described in Hillman et al. (35), allowed for the identification of eight distinct cell types (Figure 2A). In order to mitigate potential biases due to sample-to-sample variations in cell yields, we here compared the relative cell counts. Figure 2A shows comparable immune cell compositions between patients and controls, with only minor differences in NK cells (Figure 2B).

UMAP clustering of cells revealed distinct expression profiles between patients and controls in classical and intermediate, but not non-classical, monocytes (Figures 3A, B). Differential gene expression analysis, based on a fold change threshold greater than 3 and a p-value threshold below 10⁻¹⁸, identified 40 DEGs in classical, 13 in intermediate, and 6 in non-classical monocytes, while other cell types showed minimal changes (Figure 3C).

The volcano plot in Figure 4A details these findings for the classical monocytes. In total, there were 37 downregulated genes in patients with long COVID, among them IL1B, several chemokines (CCL3, CCL4, CXCL1, CXCL2, CXCL3, and CXCL8), and the TNFα-induced protein 3 (TNFAIP3). Further downregulated genes included DUSP2, encoding a phosphatase that downregulates members of the mitogen-activated protein (MAP) kinase superfamily (53), and MIR155HG, a microRNA involved in the regulation of MHCII antigen presentation (54) (Figure 4A). In order to eliminate the possibility that the observed differences between patients and controls were merely a consequence of the varying intervals between the blood sampling and the most recent SARS-CoV-2 infection, we employed the anti-nucleocapsid titers as a surrogate and correlated them with the respective raw read counts for each of the 40 DEGs. Supplementary Table S2 presents the correlation coefficients, with range between 0.0133 and 0.185 or −0.0003 and −0.217, thus ruling out any simple linkage. To further evaluate the use of convalescents as a control group, we compared our data to those of publicly available datasets generated via scRNA-seq of healthy controls from two different studies (55, 56). These external control groups were free of immune disorders and, when compared with our long COVID group, confirmed a consistent downregulation of the majority of the DEGs (Supplementary Figure S2).

As monocytes displayed the most DEGs, GSEA was performed for this subset using NES > 0 and FDR < 0.01. Patients with long COVID showed a significantly reduced expression of the inflammation-related pathways, including TNFα signaling via NF-κB, consistently through all subsets of monocytes (Figure 4B). Additional pathways, such as inflammatory response and apoptosis, were also downregulated, while heme metabolism and hallmark genes of oxidative phosphorylation (although failing to reach significance) demonstrated a higher NES in patients compared to the controls (Figure 4C). Of note is that the downregulated genes in the monocytes of patients formed a homogeneous pattern, clearly separating patients and controls, and most prominently in classical and intermediate monocytes, although to a lower extent also visible in non-classical monocytes (Figures 5A–C).

The three most significantly upregulated genes in the classical monocytes of patients with long COVID were ENSG00000287256, GRIP1, and TNFRSF12A, the latter a weak inducer of apoptosis linked to fibrosis-related pathways (57) (Figure 4A). However, application of Euclidean distance metrics revealed for the former two an upregulation in only very few patients, while the majority remained inconspicuous (Figure 5A). Supplementary Figure S3 summarizes the differences in the expression patterns between patients and controls based on median counts per million (CPM) per DEG. Consistent with the GSEA results, intermediate monocytes showed significant expression of the genes related to inflammation, such as IL-1b, CXCL2, and CXCL8, with uniformly lower results in patients.

The plasma concentrations of CXCL2 confirmed the RNA sequencing (RNAseq) data, while other chemokines, such as CCL20, showed a trend that did not quite reach statistical significance (Figure 6). The plasma concentrations of IL-6, TNFα, IFNγ, IL-1β, CXCL1, CXCL8, and CCL4 remained below their respective detection thresholds (Supplementary Table S3).

As there was a minor quantitative difference observed for the NK cell populations, differential expression of the genes was also analyzed. However, only two genes allowed for the clustering of patients and controls. These genes were identified as TNFα-induced protein 8-like protein 2 (TNFAIP8L2) and CISH, both of which have been implicated in the regulation of immune processes. Both were downregulated in patients with long COVID (Supplementary Figure S4).

In summary, we here showed that the differences in the immune landscape between patients with long COVID and controls were i) restricted to classical and intermediate monocytes; ii) reflected in the DEGs that are predominantly downregulated; and iii) confirmed as reduced protein concentrations in the plasma for CXCL2 and CCL20.

Subsequent subclustering analysis yielded 10 distinguishable monocyte populations (Figure 7A). Subcluster 3 was predominantly associated with non-classical monocytes, as well as a portion of subcluster 7. All the other subclusters were associated with classical and intermediate monocytes (Figure 7B). The bubble heatmap in Figure 7C provides a detailed visualization of all 40 DEGs, highlighting the subclusters in which each gene was expressed and the proportion of cells within those subclusters expressing the respective genes. Figure 7D attributes subclusters 0, 5, and 7 to monocytes from the patients and subclusters 4 and 6 to the controls.

Moreover, the RNA velocity trajectories suggested an origin of differentiation and thereby majority of the immature monocytes in subcluster 0 and neighboring areas that comprise classical and intermediate monocytes from patients in subclusters 1 and 7 (Figure 7D). From this origin, the trajectories extended in two directions, either via subcluster 1 toward subcluster 4 or toward subcluster 2. However, interestingly, the trajectories pointed from subcluster 6 toward subcluster 5, suggesting another origin of differentiation. Similarly, the trajectories from subcluster 3, consisting of non-classical monocytes, pointed to various directions, except toward the patients' monocytes. However, subclusters 3, 5, 7, 8, and 9 harbored the least DEGs (see Figure 7C).

Gene Ontology (GO) enrichment analysis of the 37 downregulated genes in all monocytes from patients, ranked based on the FDR, highlighted functional alterations in long COVID. Cross-referencing with the KEGG and Reactome databases revealed that the infection-related pathways (IL-17, IL-10, and TNFα signaling) were significantly enriched among the downregulated genes, indicating the reduced activity of these pathways in the classical and intermediate monocytes of patients. Similarly, pathways linked to pro-inflammatory conditions such as rheumatoid arthritis and legionellosis were found to be suppressed in patients (Figure 8A). To assess transcriptional regulation in patients with long COVID, the TF activities were imputed from the expression data of all the genes expressed in all monocytes using DoRothEA. The heatmap in Figure 8B shows that the monocytes from patients appeared to be primarily regulated by SNAI1, which is known to suppress inflammasome activity (58). Conversely, ATF2, a TF associated with monocyte differentiation (59), displayed a markedly reduced activity in patients. A more extensive spectrum of TFs exhibited weaker differences, encompassing WT1, ATF2, ETV4, DDIT3, and HMGA1.

In summary, the transcriptional profile of classical monocytes implied not only a lower expression of pro-inflammatory cytokines and various chemokines but also an involvement of TFs that suppress inflammasome activity in patients. A conspicuous deficiency of transcriptional regulation promoting monocyte differentiation was observed, indicating an immature phenotype in the classical monocytes of patients. The RNA velocity data corroborated these findings, as the patient-derived subcluster 0 not only contained the largest number of DEGs but also the least differentiated monocytes.

In order to explore potential diagnostic indicators for long COVID, we set out to investigate the correlations between routine pathology parameters, quality of life, and the DEGs in classical monocytes. We here concentrated on classical monocytes as they not only expressed the most DEGs but also constituted by far the largest fraction of peripheral monocytes. As demonstrated in, there was no linkage, and neither was there any correlation between the main clinical symptoms and the DEGs. 1