Single-cell study shows immune changes including monocyte and natural killer cell shifts, T cell exhaustion, and loss of specific T cells linked to Galectin-9 in Long COVID with ME/CFS

A total of 20 human participants were enrolled for this study. Given that our LC cohort consists of over 75% females, we recruited ten female participants with LC/ME/CFS, with a median age of 48.25 ± 8.55 years. The control group included ten female individuals who had recovered (R) from acute COVID-19 disease without any persistent symptoms or complications, with a median age of 47.10 ± 10.3 years. All participants had their SARS-CoV-2 infection confirmed by PCR testing conducted at the University of Alberta Hospital (Edmonton, Canada). Both the LC and R cohorts were recruited approximately 12 months after the initial onset of infection. Freshly isolated PBMCs from five participants from each group were subjected to scRNAseq. In addition, PBMCs and plasma from these participants, together with samples from an additional five participants per group, were used for mechanistic studies. LC patients with ME/CFS were selected based on clinical evaluation, laboratory tests, and the administration of seven validated questionnaires consistent with the Canadian Consensus Criteria (CCC) for ME/CFS and WHO guidelines, as we have reported before (3, 11). Especially, we used the DePaul Symptom Questionnaire (DSQ) to identify patients meeting ME/CFS criteria, followed by assessments of fatigue severity using the FACIT Fatigue Scale (Version 4) and the Multidimensional Fatigue Inventory. A DSQ score of at least 5 out of 6 on key symptoms, such as PEM was required for ME/CFS diagnosis.

We excluded patients with confounding health conditions or those with the history of ICU admission at the time of acute SARS-CoV-2 infection or co-morbidities. Additionally, as detailed in our previously published work (3, 20), LC patients and R were carefully matched for age, sex, BMI, and time from infection, and individuals with underlying medical conditions were excluded. Based on this cohort design and the limited sample size, additional covariate regression using a DESeq2 design formula was not performed.

This study was approved by the Human Research Ethics Board (HREB) at the University of Alberta (Protocol# Pro00099502), and written informed consent was obtained from all participants prior to enrollment.

Freshly collected whole blood was processed within 30 min following venipuncture. Briefly, blood was centrifuged at 300 × g for 10 min to collect plasma, which was aliquoted into cryogenic tubes. The remaining blood was diluted 1:1 with sterile PBS and layered onto 15 mL Ficoll Paque gradients (GE Healthcare, Chicago, IL, USA) in a 50 mL conical tube, followed by centrifugation at 300 × g for 20 min at room temperature (RT) with no brake. The buffy coat was collected and washed twice in RPMI 1640 (Sigma‐Aldrich, Burlington, MA, USA) supplemented with 10% FBS (Sigma‐Aldrich), and 1% penicillin/streptomycin (Sigma‐Aldrich) at 300 × g at RT. Immediately after washing, PBMCs were either fixed for scRNAseq or subjected to flow cytometric analysis or functional assays according to our standard laboratory protocols (21, 22). In some experiments, PBMCs were cultured without or with human recombinant Gal-9 (rGal-9; 0.5 μg/ml, Gal Pharma) overnight for apoptotic assay.

Fluorophore-labeled antibodies with specificity to human cell antigens were purchased from BD Biosciences and Thermo Fisher Scientific. The following Abs were used in our study: anti-CD3 (SK7), anti-CD4 (RPA-T4), anti-CD8 (RPA-T8), anti-TIM-3 (7D3), anti-CD26 (M-A261), anti-CD161(HP-3G10), anti-TV α7.2 (3C10), anti-IL-18R α (H44), and anti TCRγδ (B1). To exclude dead cells, live/dead fixable Aqua staining (Thermo Fisher Scientific cat # L34965↗) was used. For defining the gating strategy and antibody specificity, appropriate FMO (fluorescence minus one) and isotype control antibodies were used per the supplier’s recommendation. Cells were stained in 96-well round-bottom plates at 4°C for 20 min using a live/dead stain, followed by 30 min incubation at 4 C with antibodies against other surface markers. Annexin V (BD Biosciences) was used to assess apoptosis in immune cells. To prevent variations, we used consistent flow cytometry panels with the same gating strategy for all study subjects, and a minimum of 100,000 events was acquired for each cell subset according to our protocols (23–25). Also, BD Biosciences CompBeads were used for control compensation. Following staining, cells were fixed and acquired on an LSRFortessa-SORP or Fortessa- X20 (BD Biosciences) and analyzed using the FlowJo software (version 10).

Frozen plasma samples at -20°C were thawed and centrifuged for 15 min at 1500g followed by dilution for Gal-9 concentration using the ELISA kit (R&D, DY 2045).

PBMCs were fixed using the Fixation Kit (Parse Biosciences, Seattle, WA, USA) according to the manufacturer’s protocol. After fixation, the cell count was determined using the Countess 3 Automated Cell Counter (Thermo Fischer Scientific, Waltham, MA, USA). A total of approximately 322,000 cells were used from 10 donors, resulting in an average of approximately 32,200 cells per donor. These cells were subjected to the Evercode WT V3 single cell whole transcriptome kit (Parse Biosciences, Seattle, WA, USA) for library construction. The scRNA-seq libraries were sequenced using the Illumina NovaSeq X Plus platform, at a sequencing depth of approximately 35,000 reads per cell. Sequencing data was aligned and quantified against the GRCh38 human reference genome using TRAILMAKER pipeline v1.4.1. The resulting filtered count matrix was subsequently used for downstream analysis. The scRNA-seq data generated in this study are made publicly available under SRA accession number PRJNA1257454.

All quality control procedures were performed using Seurat v5.1.0 within the R statistical environment (v4.4.2). A Seurat object was created using the CreateSeuratObject function, including genes expressed in at least 3 cells and cells expressing a minimum of 200 genes. The object was then filtered using the subset function to remove potential doublets or multiplets (cells with > 4000 detected genes), as well as low-quality or dying cells with a high proportion of mitochondrial gene expression (> 5%). Individual datasets were merged using the merge function. To correct for batch effects across samples, dataset integration was performed using the Harmony package via the RunHarmony function (26) and the donor identity was retained. Data normalization and variance stabilization were carried out using the SCTransform function, specifying 3000 highly variable features. Dimensionality reduction was carried out using the RunPCA, FindNeighbours, and FindClusters functions in succession, retaining 25 principal components (PCs) as determined by inspection of the PCA elbow plot. Clustering was performed by varying the resolution parameter from 0 to 1 in increments of 0.1 using the FindClusters function. The optimal clustering resolution of 0.5 was selected based on the clustering hierarchy visualized with the Clustree package (27).

Differential gene expression analysis was performed by comparing the LC and R libraries. Upregulated and downregulated differentially expressed genes (DEGs) were identified using the FindAllMarkers function in Seurat, with parameters set to min.pct = 0.3 and logFC threshold = 0.3. A total of 227 DEGs with adjusted p-values< 0.05 were retained for downstream analysis. Heatmaps were generated using the DoHeatmap function to visualize DEG expression patterns. Additional visualization methods included density plots (DensityPlot), feature plots (FeaturePlot), violin plots (VlnPlot), ridge plots (RidgePlot), and dot plots (DotPlot), each used to display gene expression distributions and marker profiles across clusters.

Common exhaustion genes (CTLA4, TOX, EOMES, TIGIT, PDCD1, LAG3, HAVCR2, CD160, ENTPD1, and BATF) (28, 29), MAIT activation genes (KLRB1, IL7R, STAT1, STAT3, STAT4, CD69, RORA, JUN, JUNB, and SELL) (30), and NK activation genes (NKG7, KLRD1, NCAM1, PRF1, IFNG, IL2RB, GZMB, GNLY, NCR1, KLRF1, BCL2, FCGR3A, CD69, GZMA, LAMP1, KLRB1, CASP8, and RIPK1) (31, 32) were used to calculate T cell exhaustion, MAIT, and NK activation scores, respectively.

Single-cell gene set enrichment analysis was performed to identify differentially enriched pathways between LC and R groups within each cell cluster. For each qualifying cluster, pseudo bulk expression profiles were generated by aggregating raw gene counts across biological replicates within each condition. Differential expression between LC and R was assessed using edgeR with a quasi-likelihood framework. Genes were ranked by a signed significance metric, the sign of the log fold change multiplied by the negative log10 of the raw p-value. This approach prioritizes genes with large directional changes and high statistical confidence. GSEA was carried out using the fgsea package in R. Pathways were sourced from the MSigDB Hallmark collection and Reactome pathways. Only pathways with an adjusted p-value (padj)< 0.25 were retained. Normalized Enrichment Scores (NES), nominal p-values, adjusted p-values, and leading-edge genes were computed and reported for each cluster. Positive NES values indicate pathway upregulation in LC, while negative NES values indicate downregulation. For visualization at the compartment level, immune cell clusters were grouped into adaptive and innate compartments based on canonical lineage identity. NES for selected Hallmark pathways were summarized within each compartment using significance-weighted averaging across constituent cell types. Positive NES values indicate pathways enriched in LC relative to the R controls.

Data distribution was evaluated using the Shapiro–Wilk test. For comparisons involving non-normally distributed data, statistical significance was assessed using the Mann–Whitney U test or the Wilcoxon matched-pairs signed-rank test, as appropriate. Results are presented as mean ± SEM, and P values< 0.05 were considered statistically significant. In violin plots, the central line denotes the median, while the lower and upper boundaries correspond to the first and third quartiles, respectively. No randomization procedures were applied, and no data points were excluded from the analyses.

Differential gene expression analysis was conducted on count-based data using the DESeq2 R package (R version 4.2.0). Genes were classified as differentially expressed if they met the criteria of log₂ fold change greater than +1 or less than −1, together with an adjusted P value (Padj)< 0.05. Bubble plots were generated using custom Python scripts that integrate log₂ fold change and adjusted P values, as previously described (11, 33).

Freshly isolated PBMCs from ten individuals, five LC patients with ME/CFS and five recovered individuals (R) following acute SARS-CoV-2 infection without any LC symptoms, were processed for scRNAseq. LC patients with ME/CFS were diagnosed according to the criteria as outlined in the Canadian Consensus Criteria (CCC) for ME/CFS and WHO (34, 35). A total of 50,856 cells passed quality control and were retained for downstream analysis. We clustered cells using Seurat (33) based on their transcriptomic profiles and visualized them in two-dimensional space. This analysis identified 24 distinct cell clusters based on the expression of canonical marker genes (Figure 1A). Among T cells, CD4 T cells (clusters 0 and 2) were characterized by CD3E and CD4 expression, while CD8 T cells (clusters 3 and 6) expressed CD3E, CD8A, and CD8B (Supplementary Figure 1A). Tregs (cluster 10) were defined by FOXP3 (Supplementary Figure 1B), whereas cytotoxic T lymphocytes (CTLs; cluster 9) exhibited TBX21 expression (Supplementary Figure 1C). MAIT cells (cluster 13) were annotated based on KLRB1, IL18RAP, RORC, and ZBTB16 genes (Supplementary Figures 1D, E), while gamma delta (γδ) T cells (cluster 17) were identified by TRDC, TRGC1, and TRGC2 genes (Supplementary Figure 1F). Additionally, NKG7, GNLY, and NCAM1 genes were used to identify NK cells (cluster 4) (Supplementary Figure 1G) and CD79A gene to recognize B cells (clusters 5, 8, 15) (Supplementary Figure 1H).

Monocytes were subdivided into CD14+ monocytes (clusters 1, 12), characterized by the expression of S100A8 and CD14 as well as non-classical monocytes (cluster 7), distinguished by CDKN1C gene (Supplementary Figure 1I). Intriguingly, we identified a B-monocyte (B-Mono) cell (cluster 22) exhibited markers of both B cells and monocytes. Other identified blood cell types included platelets (PLTs; cluster 21), defined by TUBB1 expression (Supplementary Figure 1J), dendritic cells (DCs; cluster 16), characterized by ITGAX and CLEC10A (Supplementary Figure 2A), and basophils (Baso; cluster 14), identified by HDC expression (Supplementary Figure 2B). Plasma cells (PCs; cluster 19) were marked by JCHAIN and IGHA1 (Supplementary Figure 2C), while low-density neutrophils (LDNs; cluster 18) were defined by LTF (Supplementary Figure 2D). Furthermore, hematopoietic stem and progenitor cells (HSPCs; cluster 11) and innate lymphoid cells (ILCs; cluster 23) were classified as lineage-negative populations (36) (Supplementary Figure 2E). HSPCs were identified by HNRNPA1 and TMSB4X, whereas ILCs were marked by KIT expression (Supplementary Figure 2E). We additionally validated our cell type annotations through computational reference-based classification using SingleR (Supplementary Figure 2F). These findings provide a comprehensive overview of the major immune cell subsets identified through clustering analysis of PBMCs.

We next compared the frequency of each immune cell cluster between R and LC patients by calculating its percentage relative to the total immune cell population within each group, as determined by scRNAseq analysis (Figure 1B, Supplementary Figure 3A). These analyses revealed the enrichment of clusters effector 2 (effector CD4 T cells), 3 (effector CD8 T cells), 5 (a subcluster of B cells), 15 (a subcluster of B cells), 16 (DCs), 18 (LDNs), and 21 (PLTs) and depletion of clusters 0 (naïve CD4 T cells), 6 (naïve CD8 T cells),13 (MAIT cells), and 17 (γδ T cells) in LC patients, while other subsets show negligible changes (Figure 1B). When we compared the clusters based on cell numbers, we observed a similar pattern (Supplementary Figure 3B). These observations suggest differential effects of LC on the frequency of various immune cell subsets.

To further assess the biological impact of LC with ME/CFS across immune cell populations, we quantified the number of significantly differentially expressed genes per cell cluster. This analysis revealed that classical monocytes (cluster 1) were the most transcriptionally affected subset, followed by CD4 T cell subsets (clusters 0 and 2), NK cells (cluster 4), CD8 T cell subsets (clusters 3 and 6), and non-classical monocytes (Figure 2A). Collectively, these findings indicate that LC exerts differential transcriptional effects across distinct immune cell types.

We performed GSEA on cluster-specific differentially expressed genes derived from scRNAseq data. Across immune compartments, LC patients exhibited distinct yet partially overlapping transcriptional pathway alterations, with enrichment of TNF-α signaling via NF-κB, interferon-associated, inflammatory, and metabolic pathways, indicating widespread immune remodeling in LC (,). 1 1

CD4⁺ T cells (clusters 0 and 2) showed significant enrichment of TNFα signaling via NFκB and IFN-γ response pathways in LC patients compared with R individuals. In addition, pathways related to RNA processing and mRNA splicing were upregulated, consistent with increased transcriptional activity (Supplementary Table 1, Supplementary Figure 4A). CD8⁺ T cells (clusters 3 and 6) similarly demonstrated enrichment of interferon-associated and inflammatory pathways, accompanied by modulation of metabolic programs, including oxidative phosphorylation in select subsets. These transcriptional changes reflected altered activation states within CD8⁺ T cells in LC (Supplementary Table 1, Supplementary Figure 4A). Tregs (cluster 10) exhibited enrichment of interferon-associated and activation-related pathways in LC patients, suggesting transcriptional remodeling and instability within immunoregulatory compartments. CTLs (cluster 9) also showed enrichment of IFN-γ–responsive and inflammatory pathways, indicating persistent cytokine-associated transcriptional activity (Supplementary Table 1, Supplementary Figure 4A). Similarly, MAIT cells (cluster 13) and γδ T cells (cluster 17) displayed prominent enrichment of interferon-related pathways in LC patients (Supplementary Table 1, Supplementary Figure 4A). NK cells (cluster 4) showed enrichment of interferon-responsive and inflammatory pathways in LC patients. B-cell subsets (clusters 5, 8, and 15) demonstrated enrichment of inflammatory and complement-associated pathways in LC patients, together with modulation of metabolic programs in selected clusters (Supplementary Table 1, Supplementary Figure 4A). In contrast, plasma cells (cluster 19) exhibited downregulation of oxidative phosphorylation and metabolic pathways, indicating altered metabolic activity within antibody-secreting cell populations in LC (Supplementary Table 1, Supplementary Figure 4A). Notably, CD14⁺ monocytes (clusters 1 and 12) displayed some of the most extensive transcriptional alterations in LC patients. These cells exhibited strong enrichment of TNFα/NFκB signaling, type I and type II interferon responses, and complement activation pathways, along with increased transcriptional and RNA-processing programs (Supplementary Table 1, Supplementary Figure 4A). Non-classical monocytes (cluster 7) also demonstrated enrichment of activation- and inflammation-associated pathways, although to a lesser extent than classical monocytes and B–monocyte hybrid population (cluster 22) showed mixed innate and adaptive transcriptional signatures, consistent with altered immune-state plasticity (Supplementary Table 1, Supplementary Figure 4A). LDNs (cluster 18) exhibited strong enrichment of both type I and type II interferon pathways, together with inflammatory related gene sets and DCs (cluster 16) demonstrated altered activation- and antigen-processing–associated pathways, suggesting transcriptional remodeling of antigen-presenting functions (Supplementary Table 1, Supplementary Figure 4A). Finaly, Basophils (cluster 14), despite low abundance, also displayed activation-associated transcriptional signatures (Supplementary Table 1, Supplementary Figure 4A).

Notably, heatmap visualization of TNFα signaling via NFκB as the dominant enriched pathway across immune cell clusters revealed the highest pathway enrichment in classical monocytes, MAIT cells, CTLs, CD4 T cells, and NK cells followed by intermediate enrichment in multiple lymphoid subsets, while platelets exhibited minimal pathway enrichment (Figure 2B).

To provide a higher-level view of pathway alterations across immune compartments, Hallmark pathway enrichment scores were summarized at the compartment level by grouping immune cell subsets into adaptive and innate populations. Radar plot visualization revealed broadly shared enrichment of inflammatory and interferon-associated pathways across both compartments in LC patients with ME/CFS (). Specifically, TNFα signaling via NFκB, interferon-γ response, and inflammatory response pathways were enriched in both adaptive and innate immune cells. In contrast, type I interferon (IFN-α) signaling showed relatively stronger enrichment within innate immune compartments. IL-6–JAK–STAT3 signaling and oxidative phosphorylation pathways displayed more moderate and heterogeneous enrichment patterns across compartments (). Together, this analysis highlights coordinated yet non-identical transcriptional pathway remodeling between adaptive and innate immune systems in LC. 1 1

In our previous study when bulk RNAseq on the whole blood was performed (11), we described the top 10 upregulated and top 10 downregulated genes in LC patients compared to the R group. Given that whole blood contains RNA from diverse tissue sources, we questioned whether these genes were detectable in PBMCs of LC patients. Among the upregulated genes, RELN, RPL17-C18orf32, MEIS2, ZNF469, and ZFHX3 were detected in our scRNAseq dataset (Figures 3A–G). Notably, we found that RELN and RPL17-C18orf32 were broadly expressed across various immune cell types (Figure 3A, B). While MEIS2 expression was primarily observed in basophils (Figure 3C), ZNF469 was predominantly expressed in DCs and non-classic monocytes (Figure 3D), and ZFHX3 was highly expressed in both types of monocytes (Figures 3E, G). Consistent with our bulk RNAseq findings, while RELN, RPL17-C18orf32, ZFHX3, and ZNF469 were generally upregulated, MEIS2 expression was paradoxically downregulated in LC patients (Figure 3F).

Next, we examined the expression of the top 10 downregulated genes in bulk RNAseq in scRNAseq dataset. Among these, PTPRU, GMPR, WDR62, MUC12, TMEM191B, PLEKHN1, PCDHGC5, and AIRE were detected and found to be broadly distributed across all cell clusters (Figures 4A–H). We found the downregulation of TMEM191B, PLEKHN1, and PCDHGC5 in LC patients, which is consistent with our bulk RNAseq data. In contrast, the expression of PTPRU, GMPR, WDR62, MUC12, and AIRE were paradoxically upregulated in LC individuals (Figure 4I).

These observations highlight differential sources of specific genes in immune and non-immune cells.

Our whole blood RNAseq analysis revealed the upregulation of 25 olfactory receptor (OR) genes and downregulation of 2 olfactory-related genes (OLFM1 and OLFM2) in LC patients (11). Intriguingly, none of these 25 OR genes were detected in our scRNAseq dataset. However, consistent with our bulk dataset, we observed downregulation of OLFM1 and OLFM2 in LC patients (Figures 4I–K). Interestingly, OLFM1 was primarily expressed in monocytes, as well as DCs, while OLFM2 showed diffuse distribution across PBMCs (Figures 4J, K). These findings indicate non-immune cell sources for detected 25 ORs in bulk RNAseq analysis (11).

Additionally, our bulk RNAseq data had revealed the upregulation of hormone-related receptors, including estrogen receptors 1 and 2 (ESR1, ESR2), androgen receptor (AR), and the glucocorticoid receptor (NR3C1), in female LC patients compared to R females (11). Our scRNAseq analysis confirmed the same pattern in PBMCs of LC patients (Figure 4L). Notably, NR3C1, AR, and ESR2 showed the highest expression in the pre-B cell population (cluster 20), while ESR1 was highly expressed in plasma cells (cluster 19) (S Figure 5A). The expression levels of these genes were significantly upregulated in LC patients compared to the R group (Figure 4M). These findings highlight distinct transcriptional signatures in PBMCs of LC patients, including dysregulation of immune-related genes, olfactory signaling components, and hormone receptor pathways, suggesting complex immune and hormonal involvement in the pathophysiology of LC.

We then compared changes in T cell populations between R and LC patients. As shown in Figure 1A, CD4⁺ T cells were represented by clusters 0 and 2, while CD8⁺ T cells corresponded to clusters 3 and 6. Closer examination revealed that cluster 0 expressed naïve CD4⁺ T cell marker genes such as CCR7 and SELL, whereas cluster 2 expressed genes typically associated with effector CD4⁺ T cells, including CD69, KLRG1, and TOX (Figure 5A). Similarly, within the CD8⁺ T cell population, cluster 6 exhibited markers of naïve CD8⁺ T cells, while cluster 3 was enriched for genes linked to effector CD8⁺ T cells (Figure 5B).

Notably, clusters 2 and 3, representing activated CD4⁺ and CD8⁺ T cells, respectively, were expanded in LC patients, whereas clusters 0 and 6—representing naïve T cells—were depleted (Figure 1B, Supplementary Figure 3B). These findings may reflect a transition from naïve to effector states, potentially driven by increased immune activation in LC patients, as we previously reported by flow cytometry analysis (3).

Next, we subsetted cells with a CD3E expression level greater than 0.5 (Figure 5C) and compared the expression of cytokine genes typically secreted by T cells, as well as the genes encoding perforin (PRF1) and granzyme B (GZMB), two key cytotoxic molecules, between R and LC patients. In this subset, we observed higher expression of T cell-associated cytokine genes including TNF, IFNG, IL2, IL12A, IL6, IL10, IL1B, CXCL8 (IL-8), IL13, IL27, PRF1, and GZMB compared to the R group (Figure 5D). Given that persistent immune activation may contribute to T cell exhaustion (37, 38), we observed a significantly higher exhaustion score in LC patients (Figure 5E), consistent with our previous findings (3). Additional GSEA analysis revealed significant enrichment of a CD8 T cell exhaustion-associated gene signature in clusters 3, 6, and 9 in LC compared with R controls, with exhaustion-related genes preferentially ranked among transcripts upregulated in LC (Figure 5F). Among the most highly ranked genes were CTLA4, ENTPD1, and TOX, respectively.

While our previous bulk RNAseq analysis revealed a reduction in Tregs subset in LC patients (11), the proportion of Tregs (cluster 10) appeared unchanged between LC and R groups in the current dataset (Figure 1B, Supplementary Figure 3B). However, further analysis revealed a substantial reduction in the expression of FOXP3 gene—both across the entire PBMC population and specifically within the Treg cluster in LC patients (Figures 5G, H). These findings support our previous observations, highlighting increased T cell exhaustion and possibly impaired Treg function in LC patients.

Finally, to gain deeper insight into the biological impact of LC on different T cell subsets, we quantified the number of DEGs within each cluster. This analysis revealed that cluster 3 (effector CD8 T cells) exhibited the highest number of DE genes, followed by cluster 0 (naïve CD4 T cells) and naïve CD8 T cells (). In contrast, LC was associated with a smaller number of transcriptional changes in cluster 9 (CTLs) and cluster 10 (Tregs) (). These results demonstrate that effector CD8 T cells are disproportionately affected by LC at the transcriptional level. 1 1

We then evaluated unconventional T cell populations, including MAIT and γδ T cells in patients with LC. We observed a reduction in MAIT cells, represented by cluster 13, in LC patients (Figure 1B, Supplementary Figure 3B). Further subsetting of this cluster revealed five MAIT cell subpopulations (Figure 6A). Interestingly, while the proportions of subclusters 0 and 4 as most dominant subsets were decreased, there was an increase in the proportion of subclusters 2 and 3 in LC patients (Figure 6B). Gene expression analysis revealed that subclusters 2 and 3 were enriched for genes associated with activated MAIT cells, such as STAT1, STAT3, CD69, JUN, JUNB, RORA, IL23R, SLC4A10, TFRC, and KLRB1, and exhibited lower expression of genes typically associated with a naïve MAIT phenotype, including SELL, KLRG1, and IL7R (Figure 6C). We also observed downregulation of IL7R and KLRG1 and upregulation of STAT1, STAT3, JUN, JUNB, RORA, and CD69 in the overall MAIT cell cluster in LC patients (Figure 6D). Consistently, MAIT cells also exhibited a higher activation score in LC patients (Figure 6E). These results suggest that the overall decrease in MAIT cells in LC patients is largely driven by a loss of naïve MAIT cells, while the proportion of activated MAIT cells relatively increased.

We also observed a reduction in γδ T cells in LC patients. Specifically, cluster 17, representing γδ T cells, was reduced in LC patients (Figure 1B, Supplementary Figure 3B). Further sub-setting of this cluster revealed three distinct γδ T cell subpopulations, all of which were depleted in LC patients (Figures 6F, G) with differential pattern of gene expression profile in these subsets (Supplementary Figure 6A). These findings suggest that LC may be associated with a reduction in multiple γδ T cell subsets.

Evaluation of the NK cell population (cluster 4) (Figure 1A) revealed a reduced frequency in LC patients, consistent with our previous bulk RNAseq data (11) (Figure 1B, Supplementary Figure 3B). Sub-clustering of the NK cell population within cluster 4 revealed three distinct NK cell subsets (Figure 7A). Examination of CD56 (NCAM1) and CD16 (FCGR3A) gene expression within these subclusters showed that FCGR3A was expressed in subclusters 0 and 1, while NCAM1 was predominantly expressed in subcluster 2 (Figure 7B). This pattern may suggest a shift in NK cell composition in LC patients, with a relative decrease in the NCAM1-enriched subcluster (Figure 7C). Traditionally, CD56⁺ CD16⁻ NK cells are considered prominent cytokine producers with limited cytotoxicity, whereas CD56^dim CD16⁺ NK cells are known for their potent cytolytic activity but limited cytokine production (9, 39). Based on this observation, we evaluated the expression of various cytokine genes, activation markers, and cytotoxic molecules within the NK cell cluster in LC patients. This analysis revealed a decrease in the expression of IFNG, NKG7, KLRD1, KLRF1, PRF1, IL2RB, GZMA, GNLY, NCR1, and BCL2 in NK cells from LC patients. Conversely, we observed an increase expression of CD69, GZMB, LAMP1, KLRB1, RIPK1, and CASP8 in NK cells from the same group (Figure 7D). Given the mixed pattern of gene expression, some indicative of enhanced NK cell function and others of impaired activity, we calculated an activation score based on the expression of these genes. This score was significantly lower in NK cells from LC patients (Figure 7E). These results suggest that NK cells in LC patients not only occur at reduced frequencies but also show differences in the expression of transcripts associated with NK cell function.

Identification of B cells based on CD79A expression revealed that clusters 5, 8, and 15 correspond to B cells (Figure 1A, Supplementary Figure 1H), with clusters 5 and 15 showing expansion in LC patients (Figure 1B, Supplementary Figure 3B). To further characterize these populations, we examined the expression of genes associated with different B cell states. Genes typically associated with naïve B cells, such as SELL, PAX5, BCL2, CD19, CD22, MS4A1, IGHD, and IGHM, were upregulated in R individuals. In contrast, genes associated with effector or activated B cells, including CR2, FCERs, CD69, FOXP1, CD40, JUN, and FOS, were more expressed in LC patients (Figure 8A). These findings may suggest an increased activation state of B cells in LC patients compared to the R group.

Plasma cells (PCs) were primarily represented by cluster 19 (Figure 1A, Supplementary Figure 2C). Upon subsetting cluster 19, we identified two distinct subclusters (Figure 8B). The proportion of subcluster 0 was lower, while subcluster 1 was more prevalent in LC individuals (Figure 8C). Interestingly, subcluster 0 expressed canonical plasma cell markers such as IGHA1, IGHA2, IGLC2, and IGLC3, confirming its identity as a conventional plasma cell population. In contrast, subcluster 1 (19-1) displayed expression of several genes typically associated with NK cells (Figure 8D).

Given that subcluster 1 was positioned near the NK cell cluster (cluster 4) in the UMAP space (Figure 1A), we next questioned whether this population might represent misclassified NK cells. However, heatmap analysis revealed distinct gene expression patterns between subcluster 1 and the NK cell cluster, suggesting that despite the shared markers, subcluster 1 likely represents a unique plasma cell subset or a transitional population rather than true NK cells (Figure 8E).

When the frequency of basophils (cluster 14) in LC patients versus the R group was assessed (Figure 1A), no significant difference in the frequency was observed (Figure 1B, Supplementary Figure 3B). Further analysis of this population revealed the presence of two distinct subclusters (Figure 9A). When examining genes linked to basophil activation, such as CD9, HDC, GATA2, FCER1A, CR1, IL1RL1, IL5RA, MS4A2, CLC, AKAP12, and C5AR1 (40, 41), we found that these genes were enriched in subcluster 1 (Figure 9B). Although the proportion of subcluster 0 was higher in LC patients, subcluster 1 was proportionally reduced (Figure 9C). We also observed a downregulation of these basophil-related activation genes in the overall basophil cluster from LC patients (Figure 9D). In line with these findings, the activation score derived from the expression of these genes was significantly lower in LC patients (Figure 9E). These findings suggest reduced activation-associated transcriptional profiles in basophils from LC patients.

We next characterized neutrophils from R and LC patients based on transcriptional profiles in our scRNA-seq dataset. Although mature neutrophils are typically excluded during Ficoll-Paque separation due to their higher density, LDNs can be found among PBMCs following leukocyte separation by density gradient centrifugation (42, 43). Analysis of Cluster 18, which corresponds to LDNs, revealed their expansion in LC patients (Figures 1A, B, Supplementary Figure 3B). Further examination of this cluster identified three subclusters, with subcluster 1 displaying a distinct gene expression profile, consistent with its separation in UMAP space (Figure 9F). Interestingly, only subcluster 1 was expanded in LC individuals (Figure 9F). A heatmap of the most highly expressed genes in each subcluster showed that subcluster 1 was enriched in genes such as MME and CXCL8, both of which are associated with neutrophil activation (Figure 9G). Together, these results suggest the expansion of a transcriptionally distinct LDN subset with elevated expression of activation-associated genes in LC patients.

Monocytes were classified into three distinct clusters (1, 12, and 7), corresponding to classical, intermediate, and non-classical monocytes, respectively (Figure 1A). This classification was based on the expression levels of key marker genes—CD14, FCGR3A, and CDKN1C (Figure 10A, Supplementary Figure 1I). Cluster 1 exhibited the highest expression of CD14, while Cluster 7 showed the highest expression of FCGR3A and CDKN1C (Figure 10A). Cluster 12 displayed intermediate expression of these genes and was therefore identified as the intermediate monocyte population (Figure 10A). Notably, there were no significant differences in the frequency of these monocyte subsets between LC patients and the R group (Figure 1B, Supplementary Figure 3B). To evaluate monocyte function, we subsetted cells with CD14 expression levels greater than 0.5 and assessed the expression of genes associated with phagocytosis (TLR1, TLR2, TLR3, TLR4, TLR8, CD80, and CD86) and cytokine signaling (IL1B, TNF, CXCL8, IL23A, IL12A, IL27, and IL15). The phagocytosis-associated gene score was reduced (Figure 10B), whereas a cytokine-associated gene score was elevated in monocytes from LC patients (Figure 10C). These findings reflect a shift in transcriptional profile of monocytes in LC patients, characterized by reduced expression of phagocytic-related genes but enhanced expression of pro-inflammatory cytokine genes.

Finally, we analyzed platelets assigned to cluster 21 and compared their distribution between recovered (R) individuals and LC patients. As shown in Figure 1A and Supplementary Figure 3B, platelets were markedly expanded in LC patients. Further subclustering of this population revealed two distinct platelet subclusters, both of which were also expanded in the LC group (Supplementary Figure 6B). In addition, genes associated with platelet activation were upregulated in platelets from LC patients (Figure 10D).

To determine whether the immune cell changes observed in ME/CFS associated with LC are also present in idiopathic ME/CFS, we re-analyzed an independent, publicly available scRNA-seq dataset comparing PBMC gene expression profiles from individuals with ME/CFS and healthy controls (19). To ensure consistency with our dataset, we included only female participants, comprising 20 healthy controls and 19 ME/CFS patients. All samples were collected prior to the COVID-19 pandemic, ruling out the possibility of LC misclassification. After quality control, 45,047 cells were retained for downstream analysis. UMAP clustering identified 13 distinct cell clusters based on the expression of canonical marker genes (Figure 10E, Supplementary Figure 7A). In addition to the previously described markers, we used CLEC4C and SCT to define plasmacytoid dendritic cells (pDCs; cluster 9), and HBA1 and HBA2 to identify red blood cells (RBCs; cluster 11) (Figure 10E, Supplementary Figure 7A).

We merged the cells from two conditions and evaluated the function of different immune cells (Supplementary Figure 7B). Cells with CD3E expression greater than 0.5 were selected as T cells (Figure 10F), and their activation status was assessed using a module score based on the following genes: IFNG, GZMB, PRF1, CTLA4, FOS, TNF, IL2RA, EOMES, TOX, IL2, CD44, MYC, TBX21, ICOS, CD69, and JUN. This analysis revealed higher T cell activation in ME/CFS patients compared to healthy controls (Figure 10G). In contrast, no significant difference in T cell exhaustion status was observed between groups when assessed using expression of exhaustion-associated genes: CTLA4, EOMES, TOX, TIGIT, PDCD1, LAG3, HAVCR2, CD160, ENTPD1, and BATF (Figure 10H). Notably, unlike our cohort, FOXP3 expression was upregulated in idiopathic ME/CFS patients (Figure 10I).

We were unable to identify a distinct cluster corresponding specifically to MAIT cells or NK cells using canonical markers in this dataset. However, cluster 2 exhibited high expression of genes associated with MAIT cells, γδ T cells, and NK cells (Supplementary Figure 7C). We therefore subsetted cluster 2 and found that subclusters 0 and 3 showed elevated expression of MAIT cell–related genes, including ZBTB16 and IL18RAP, and also included NK cells, as indicated by low expression of CD4, CD8A, and CD8B, and high expression of NK cell markers such as KLRD1, NCAM1, NCR1, FCGR3A, KLRK1, PRF1, GZMB, and GNLY (Supplementary Figures 7D, E). Based on these findings, we calculated activation scores for both MAIT and NK cells within these subclusters. This analysis revealed no significant difference in NK cell activation between healthy individuals and ME/CFS patients (Supplementary Figure 7F). However, we observed a higher activation score in MAIT cells from ME/CFS patients compared to healthy controls (Supplementary Figure 7G) without any changes in the proportion of subclusters 0 and 3 between groups (Supplementary Figure 8A).

We next assessed the function of B cells by analyzing the expression of genes associated with naïve and activated B cell subsets, using the same gene set applied in our dataset. We observed both upregulation and downregulation of genes corresponding to different B cell subsets, suggesting a dysregulated distribution of B cell types between ME/CFS patients and healthy individuals (Supplementary Figure 8B). Additionally, platelet-related gene expression in ME/CFS patients showed mixed patterns of upregulation and downregulation, indicating alterations in platelet activation status (Supplementary Figure 8C). In contrast to our LC/ME/CFS cohort, we found no significant differences in activation or cytokine scores in monocytes (cluster 1) between the two groups (Supplementary Figures 8D, E). We also evaluated the activation status of pDCs by examining the expression of genes associated with their activation, including IRF7, IRF3, IRF4, TLR7, TLR9, ISG15, MX1, IFI6, IFI44L, IFIT1, IFIT3, OAS1, OAS2, OASL, TNFSF10, IL6, CD86, CD83, HLA-DRA, HLA-DRB1, LAMP3, BST2, LY6E, and CLEC4C. This analysis revealed upregulation of certain genes, such as CLEC4C and LY6E, and downregulation of others, including HLA-DRA (Supplementary Figure 8F). However, the overall activation score of pDCs did not differ significantly between ME/CFS patients and healthy controls (Supplementary Figure 8G).

In agreement with our scRNAseq data in the present study, we observed a significant reduction in the frequency of CD8⁺CD26^hi T cells (Figures 11A, B), which express MAIT-associated markers CD161, TCR Vα7.2, and IL-18Rα (Figure 11C, Supplementary Figure 8H) (44). Consistent with our previous findings (2), plasma Gal-9 levels were elevated in LC patients compared to the R group (Figure 11D). Given the activation state of MAIT cells (Figure 6E), we found that these cells exhibited significantly higher levels of TIM-3 expression than their counterparts in the R group (Figures 11E, F). Since Gal-9 is known to induce apoptosis and suppress TIM-3+ T cells (7, 45), we examined its effect on MAIT cells using total PBMCs. These experiments revealed that recombinant Gal-9 (rGal-9) enhanced MAIT cell apoptosis (Figures 10G, H), which subsequently resulted in a reduction in MAIT cells (Figures 11I, J). Notably, the addition of lactose (30 mM) (46) partially but significantly reversed rGal-9-induced MAIT cell apoptosis (Figure 11K), consistent with previous reports in other models (3, 44).

We also found a significant reduction in the frequency of γδ T cells in LC patients (Figures 12A, B). Similar to MAIT cells, γδ T cells in LC patients exhibited elevated TIM-3 expression (Figure 12C). Notably, increased apoptotic potential of γδ T cells in LC patients was evident in freshly isolated cells (Figure 12D). This was further confirmed by in vitro experiments, where overnight culture of PBMCs with rGal-9 promoted γδ T cells depletion via enhanced apoptosis (Figures 12E–G). However, this effect was not observed in γδ T cells from the R group (Figures 12F, H). Collectively, these observations provide a potential mechanistic explanation for MAIT and γδ T cell activation and apoptosis in LC patients.