Immunometabolic changes and potential biomarkers in CFS peripheral immune cells revealed by single-cell RNA sequencing

Abstract

The pathogenesis of Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) remains unclear, though increasing evidence suggests inflammatory processes play key roles. In this study, single-cell RNA sequencing (scRNA-seq) of peripheral blood mononuclear cells (PBMCs) was used to decipher the immunometabolic profile in 4 ME/CFS patients and 4 heathy controls. We analyzed changes in the composition of major PBMC subpopulations and observed an increased frequency of total T cells and a significant reduction in NKs, monocytes, cDCs and pDCs. Further investigation revealed even more complex changes in the proportions of cell subpopulations within each subpopulation. Gene expression patterns revealed upregulated transcription factors related to immune regulation, as well as genes associated with viral infections and neurodegenerative diseases.CD4⁺ and CD8⁺ T cells in ME/CFS patients show different differentiation states and altered trajectories, indicating a possible suppression of differentiation. Memory B cells in ME/CFS patients are found early in the pseudotime, indicating a unique subtype specific to ME/CFS, with increased differentiation to plasma cells suggesting B cell overactivity. NK cells in ME/CFS patients exhibit reduced cytotoxicity and impaired responses, with reduced expression of perforin and CD107a upon stimulation. Pseudotime analysis showed abnormal development of adaptive immune cells and an enhanced cell-cell communication network converging on monocytes in particular. Our analysis also identified the estrogen-related receptor alpha (ESRRA)-APP-CD74 signaling pathway as a potential biomarker for ME/CFS in peripheral blood. In addition, data from the GSE214284 database confirmed higher ESRRA expression in the monocyte cell types of male ME/CFS patients. These results suggest a link between immune and neurological symptoms. The results support a disease model of immune dysfunction ranging from autoimmunity to immunodeficiency and point to amyloidotic neurodegenerative signaling pathways in the pathogenesis of ME/CFS. While the study provides important insights, limitations include the modest sample size and the evaluation of peripheral blood only. These findings highlight potential targets for diagnostic biomarkers and therapeutic interventions. Further research is needed to validate these biomarkers and explore their clinical applications in managing ME/CFS.

Introduction

Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) is a complex, disabling disease that affects approximately 65 million people and is primarily due to a previous viral infection or sometimes other stressful events [1, 2]. ME/CFS is generally recognized as a neurological syndrome with a list of prominent multisystem symptoms: persistent and intolerable fatigue, post-exertional malaise, orthostatic intolerance and cognitive impairments such as memory or attention deficits and complex information processing disorders [3, 4]. As there is no established molecular diagnostic test, the diagnosis of ME/CFS is based on several clinical criteria, often involving a combination of specific symptoms over a prolonged period of time, usually six months, and the exclusion of other known disabling diseases or conditions [5,6,7,8,9].

Driven by the great efforts of communities and colleagues, significant progress has been made in the etiology and epidemiology of the disease over the last decade. About 80% of people with ME/CFS often began after a viral infection and research shows that previous coronaviruses triggered ME/CFS in about 11% of patients, followed by repeated or even persistent influenza-like conditions [10,11,12]. Viral infections have been associated with the onset of ME/CFS in most cases and are generally considered to be the most common trigger of the complex [12,13,14]. There are in fact two possible explanations for this link. A simple possibility is that people with ME/CFS are susceptible to infections and suffer from persistent or recurrent infections due to dysregulated immune system failure to defend against viruses [15, 16]. In a subset of ME/CFS patients, antiviral drugs can reduce fatigue and increase natural CD4⁺ T cells and natural killer(NK) cells [17]. However, there is no consistent link between persistent or chronic viral infection and ME/CFS. And so far, no specific viral etiology has been confirmed. However, the other possibility is the disruptive autoimmune-like response of our immune system to certain pathogens, which is associated with a purposeless, delayed, persistent activation and a legacy of inflammation. There is growing evidence that one or more of these viruses can cause long-term abnormalities in the immune system and lead to ME/CFS [18,19,20].

These immune system abnormalities, including altered cytokine profiles and autoantibody levels in plasma, an abnormal phenotype of T or B lymphocytes, and impaired cytotoxic responses such as a decrease in NK cell cytotoxicity, have been frequently described in individuals with ME/CFS [12, 21, 22]. Autoantibodies against various antigens, including neurotransmitter receptors, have recently been identified by several groups in some ME/CFS patients, which alone may be sufficient to create a susceptibility to transition to a chronic immune response that eventually leads to the neurological symptoms of ME/CFS [23, 24]. And clinical trials from Norway have shown that depletion of B cells with rituximab leads to clinical benefits in a subset of ME/CFS patients [25]. Although the exact pathogenesis is not yet known, the most plausible hypothesis is that dysregulation of the immune system together with the autonomic nervous system and metabolic disturbances contribute to this complex syndrome in which severe fatigue and cognitive impairment are a central manifestation. An even greater challenge is therefore to answer the question of why these patients have lost the ability to respond properly to viral or other stressful perturbations and then return to homeostasis as planned.

However, recent reviews and observational studies have yet to identify robust and reproducible immune biomarkers of disease or disease severity, as data remain inconsistent and inconclusive [26,27,28]. Current evidence has still not been able to pinpoint the possible mechanisms behind the scenes influencing the shift of an immune and/or inflammatory response from acute stress to chronic turbulence in ME/CFS. The development of molecular diagnostic tools, intervention targets or even objective biomarkers has been hampered by the lack of understanding of the cellular and molecular basis. In addition, it is crucial to understand how the brain and autonomic nervous system (ANS) manifest the neurological symptoms, likely with activation of the specific immune system and resulting neuroinflammation. Due to the great phenotypic and functional plasticity of leukocytes, molecular signaling within or between individual cell could underlie the differences and possibly shed light on these changes. ScRNA-seq allows for a more granular analysis of individual cell types and their interactions within the immune system, which is particularly crucial in a heterogeneous condition like ME/CFS [29]. In-depth studies of cellular and molecular interactions at the single cell level in peripheral leukocytes are essential for the identification of potential biomarkers and the development of immunomodulatory treatments for ME/CFS.

In this study, we aim to elucidate the immunometabolic alterations in ME/CFS by analyzing the leukocyte landscape in peripheral blood using single-cell RNA sequencing (scRNA-seq). This approach allows for a detailed characterization of individual cell types and their interactions, providing insights into the molecular pathomechanisms of ME/CFS. By defining these cellular changes, we hope to identify potential biomarkers and contribute to the development of targeted treatments for this challenging condition.

Materials and methods

Study participants

The study participants were recruited from the Jinan Biolnoo Institute of Immunometabolism. Enrolled ME/CFS participants were required to meet the research criteria of the International Chronic Fatigue Syndrome Study Group [30], the Canadian consensus criteria [31], and the diagnostic tool issued by the Institute of Medicine in 2015 [32]. All participants completed an approved questionnaire, including the validated ME/CFS Symptom Inventory and additional questions on symptoms, comorbidities and family health history. To ensure the accuracy and reliability of our findings, participants were carefully screened to exclude those with confounding factors that could directly impact immune cell populations. This includes individuals who had taken immunosuppressant drugs, antiviral therapy, or any other medication known to affect immune function within the 6 months prior to enrollment. Subjects under 18 years of age, those with progressive malignancies, acute progressive illnesses within a 2-week period, clinical or laboratory evidence of current infections, established or suspected mental or neurological disorders and autoimmunity, and those in other circumstances deemed unsuitable for study participation by the investigator were excluded. Finally, 4 patients (ME/CFS) vs. 4 age-paired heathy volunteers (HC) were enrolled for scRNAseq analysis. Lately, 10 vs. 10 subjects were enrolled for further validation of the primary findings of the scRNAseq analysis. Moreover, to strengthen our findings, we conducted external validation using publicly available scRNAseq datasets from an independent cohort (GSE214284, n = 58).

All procedures performed in the study with human participants were in accordance with the ethical standards of the institutional and/or national research committee and the 1964 Declaration of Helsinki and its subsequent amendments or comparable ethical standards. And this study was approved by the Ethics Committee of Qilu Hospital of Shandong College (No.: KYLL-2020(KS)-069). All subjects in our study voluntarily signed an informed consent form after review by the ethics committee.

Isolation of peripheral blood mononuclear cells (PBMCs)

Fresh whole blood was collected in the morning in 5-ml tubes treated with EDTA immediately at room temperature. PBMCs were isolated by Ficoll-Hypaque density gradient centrifugation (Haoyang Biological Manufacture, Tianjin, China). In brief, the blood sample was diluted 1:1 with PBS (Gino Biological Pharmaceutical Co., Ltd, Hangzhou, China) and then Ficoll was added. The mixture was centrifuged at 400 g for 20 min at room temperature. The enriched PBMCs were washed with D-Hanks containing 2% BSA (Sigma). After each wash, they were centrifuged at 100 g for 15 min and this procedure was repeated twice. The cell number and viability were then determined.

Single-cell RNA sequencing

The PBMC samples were immediately placed in GEXSCOPETM tissue preservation solution (Singleron) and transported on ice to the Singleron laboratory. The samples were rinsed three times with Hank’s balanced salt solution (HBSS) and then cut into 2 mm pieces. After digestion, the materials were filtered through sterile 40-micron sieves and centrifuged at 1,000 rpm for 5 minutes. The supernatant was discarded and the debris was resuspended in 1 ml PBS (HyClone). At 25 °C for 10 min, 2 ml of GEXSCOPETM red blood cell lysis buffer (Singleron) was added to eliminate the red blood cells. After centrifugation at 500 g for 5 minutes, the solution was suspended in PBS. Samples were stained with trypan blue (Sigma) and examined under a microscope. Single cell suspensions were prepared at a concentration of 1 × 10⁵ cells/mL in PBS (HyClone). Subsequently, these single-cell suspensions were loaded onto microfluidic devices and single-cell RNA sequencing libraries (scRNA-Seq) were prepared using the GEXSCOPER Single-Cell RNA Library Kit (Singleton Biotechnologies) according to the Singleron GEXSCOPER protocol. The verified scRNA-seq data of the samples from the studies were obtained as fastq raw files from the Gene Expression Omnibus database (GSE137804) [29].

ScRNA-seq data analysis

In our study, we employed the Seurat R package to transform scRNA-seq data into Seurat objects. Cell-level quality control was conducted, filtering cells based on the following criteria: (1) total UMI counts not exceeding 1,000; (2) gene numbers no higher than 200; or (3) mitochondrial gene percentage of greater than 20. Normalization of gene expression levels within each cell was carried out using the Normalize Data function, utilizing the Log Normalize method with a scale factor of 10,000. This step aimed to mitigate the impact of sequencing library size, converting expression values from UMI counts to ln [10,000× UMI counts/total UMI counts in cell + 1].

To address batch effects, integration of individual samples was performed using the canonical correlation analysis (CCA) pipeline within Seurat. For feature selection in integration, the ‘Select Integration Features’ function was applied, prioritizing features based on their detection in multiple datasets. Subsequently, the ‘Find Integration Anchors’ function identified 2,000 anchors across different samples using the top 50 dimensions from CCA, defining the neighborhood search space. Integration of datasets was achieved through the ‘Integrate Data’ function, using pre-computed anchors, and the integrated dataset was scaled using ‘Scale Data’.

Principal component analysis (PCA) and UMAP dimension reduction were carried out based on the top 20 principal components. The resulting clusters were projected onto 2D maps using t-distributed stochastic neighbor embedding (t-SNE) or UMAP methods. Clustering of cells was accomplished using Seurat’s ‘Find Neighbors’ with dimensions 1–20, followed by ‘Find Clusters’ with a resolution of 0.5. For the detection of gene expression markers, the ‘Find All Markers’ function was employed. Subsequently, cell type annotation in our study was performed using the Single R package and the Cell Marker dataset. To facilitate downstream analysis, the ‘Subset Data’ function was utilized to extract subclusters. After identifying clusters and gene expression markers within subclusters using ‘Find Clusters’ and ‘Find All Markers’ functions, UMAP analysis was also conducted using the ‘RUNUMAP’ function. The process of annotating subclusters followed the aforementioned approach.

Pseudo-time analysis

Developmental pseudo-time analysis was conducted using the Monocle3 package to investigate cell differentiation and state transitions. Initially, we used the normalized gene expression matrix to infer potential lineage differentiation trajectories. A new Monocle3 object was created using the new_cell_data_set function, followed by processing with the preprocess_cds function using default parameters.

To facilitate trajectory learning, we reduced data dimensionality using the reduce_dimension function with parameters set to umap.min_dist = 0.5 and umap.n_neighbors = 35. This dimensionality reduction step enabled clearer visualization and interpretation of the cell states and transitions.

In the context of ME/CFS, pseudo-time analysis helps to map out how cells progress through different states of differentiation and how these transitions might be disrupted in disease conditions. By modeling the trajectory of cell states over pseudo-time, we gain insights into the dynamic changes in gene expression associated with cell maturation or dysfunction. This approach allows us to identify key regulatory genes and pathways that are altered during the progression of ME/CFS, providing a deeper understanding of the disease mechanism and potential targets for therapeutic intervention.

SCENIC analysis

The SCENIC analysis was conducted using the SCENIC pipeline (version 1.1.2.2) with the RcisTarget (version 1.2.1) and AUCell (version 1.4.1) packages, applying default parameters. We used the motifs database available for RcisTarget to identify over-represented transcription factor (TF) binding motifs within our gene list. This step allowed us to determine which TFs are potentially regulating gene expression in our dataset. To evaluate the activity of each regulon (a group of co-expressed genes regulated by a specific TF) in individual cells, we utilized the AUCell package. This tool assesses the activity score of each regulon based on the gene expression data from single cells. For assessing cell type specificity of the predicted regulons, we calculated the regulon specificity score (RSS). The RSS was computed using the Jensen-Shannon divergence (JSD), which measures the similarity between the probability distributions of regulon activity across different cell types. Specifically, we calculated the JSD between each binary regulon activity vector and the cell type assignments to quantify how distinctly each regulon is associated with specific cell types. Additionally, the connection specificity index (CSI) for all regulons was determined using the scFunctions package (https://github.com/FloWuenne/scFunctions/). The CSI provides an additional measure of how specific a regulon is to individual cell types based on its interaction patterns.

Cell chat analysis

The R package CellChat (http://www.cellchat.org/) was employed to conduct cell-cell interaction analysis. In essence, CellChat utilizes manually curated databases that encompass established structural configurations of ligand-receptor interactions to deduce and assess intercellular communication networks from scRNA-seq data, employing network analysis and pattern recognition approaches. The Seurat object, comprising the count matrix and clustering outcomes from each dataset, was integrated into CellChat. For the analysis of single-cell datasets, the default human database (KEGG, STRING, and others that compile interaction information) was utilized. These databases were instrumental in defining the interaction networks and assessing their biological significance.

Statistical analysis

The data are expressed as the mean ± SEM. All statistical analyses were conducted using R (http://www.r-project.org), while GraphPad Prism 8 (GraphPad Software, Inc.) was employed for specific statistical procedures. To generate box plots, we employed the R Base package with default settings. The unpaired t-test was applied to compare between the two groups. Assumptions for the t-test, including normality and equal variance, were assessed using the Shapiro-Wilk test and Levene’s test, respectively. For comparisons involving multiple groups, one-way ANOVA was employed, followed by Tukey’s post-hoc test to adjust for multiple comparisons. A significance level of p < 0.05 was considered statistical significant. All statistical tests were two-tailed.

Results

Alterations in peripheral immune cell composition in ME/CFS patients

PBMCs were collected from healthy donors (HC1-4) and patients with ME/CFS (CFS1-4). Prior to blood collection, all HC and ME/CFS participants underwent basic clinical evaluation (Supplementary Materials 1-2). Peripheral blood NGS was performed on each subject to rule out acute or persistent infection with pathogens (both viral and bacterial). First, routine normalization and quality control is performed for single-cell RNA sequencing data (Fig. S1A-S1H). Using the scRNA-seq platform, we analyzed a total of 28,180 cells from all patients after filtering the data with strict quality control. This resulted in a mean value of 6,390 nUMIs per cell and the detection of an average of 1,680 genes per cell. Using t-distributed stochastic neighbor embedding (tSNE) for highly variable genes, we identified 18 unbiased clusters that were free of sampling or experimental batch bias (Fig. 1A, Fig. S2A). These clusters were categorized into 9 cell types using established marker genes (Fig. 1B and C, Fig. S2B).

Fig. 1

Single-cell transcriptomes of PBMCs from ME/CFS Patients. A. Schematic of single-cell transcriptome profiling of PBMCs from Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) patients (CFS) and healthy controls (HC). The flowchart outlines the process as follows: single-cell RNA sequencing (ScRNAseq) is performed to identify molecular pathways. This data is then subjected to internal verification and supplemented with external data verification to ensure robustness and accuracy of the findings. B. Heatmap displaying the top 10 marker gene expressions in PBMC cell subpopulations. C. Uniform Manifold Approximation and Projection (UMAP) projections of 45,633 PBMCs from the HC group (4 samples, 28,180 cells) and ME/CFS group (4 samples, 17,453 cells), depicting single-cell transcriptomes from ME/CFS and HC groups. D. Cell density plots comparing ME/CFS group and HC group segregated by ancestry. E-I. Proportion of cell types in each individual, encompassing T cells (E, p = 0.0045), NK cells (F, p = 0.0044), Monocyte cells (G, p = 0.047), B cells (H, p = 0.5), cDCs (I, p = 0.029), and pDCs cells (J, p = 0.035). Group information is represented by colors. ns, not significant; *p < 0.05, **p < 0.01; ME/CFS group compared with HC group

In the subsequent analysis, we focused exclusively on the primary immune cell types - T cells, NK cells, B cells, monocytes, conventional dendritic cells (cDCs) and plasmacytoid dendritic cells (pDCs)— - to exclusion of platelets, erythrocytes and hematopoietic stem/progenitor cells. A remarkable increase in the proportion of T cells was observed in the ME/CFS group, accompanied by a relative decrease in the proportion of NK and monocyte cells (Fig. 1D). As a preliminary illustration of the primary immunologic changes, we have shown the relative proportions of immune cell subsets among PBMCs in ME/CFS compared to healthy controls (HCs) (Fig. 1E and J). A remarkable increase in the proportion of T cells was observed in the ME/CFS group, accompanied by significant decreases in the proportions of NK, monocytes, cDCs and pDCs (Fig. 1E-G, I and J). No significant changes were observed in the proportions of B cells (Fig. 1H).

Characterization of CD4+ and CD8+ T cell subsets

To investigate the variations of CD4⁺ T cell subtypes in the ME/CFS group, we categorized CD4⁺ T cells into different groups based on the distribution and expression of classical subtype markers (Fig. S3A). We successfully identified naïve CD4⁺ subclusters (CD4⁺ naïve), CD4⁺ central memory T cells (CD4⁺ TCM), CD4⁺ effector memory T cells (CD4⁺ TEM) and CD4⁺ cytolytic T cells (CD4⁺ CTL) (Fig. 2A and S3B-S3C). Density plots of CD4⁺ T cell subsets in the ME/CFS and Con groups, categorized by lineage, revealed observable effects of different pools and processing batches on cell distribution (Fig. 2B). In the ME/CFS group, CD4⁺ naïve cells showed a significant increase (2.18-fold, p = 0.035), while CD4⁺ TEM cells showed a significant decrease (2.23-fold, p = 0.045) in total CD4⁺ T cells. However, no significant changes were observed in the relative proportions of CD4⁺ TCM and CD4⁺ CTL to total CD4⁺ T cells between the two groups (Fig. 2C and F).

Fig. 2

Cluster analysis of phenotypic and functional differences between CD4⁺T cell subsets. A. UMAP plot illustrating the distribution of CD4⁺T cell subsets, with each color representing a distinct cell subset. B. Density plots of CD4⁺ T cell subsets for ME/CFS and HC groups, segregated by ancestry. C-F. Cell numbers of CD4⁺ T cell types in each individual, including CD4⁺ Naive (C, p = 0.032), CD4⁺ TCM (D, p = 0.21), CD4⁺ TEM (E, p = 0.047), and CD4⁺ CTL cells (F, p = 0.99). Colors denote cell type information. ns, no significant. *p < 0.05. ME/CFS group compared to HC group. G-J. Monocle3-based pseudotemporal analysis of CD4⁺ T cell subsets. Cells are color-coded by CD4⁺ T cell subsets and the pseudotime trajectory in HC group (G-H) or ME/CFS groups (I-J). K-L. Regulon modules based on the correlation matrix of shared regulons of HC(K) and ME/CFS(L) of CD4⁺ T cells. The network shows 515 orthologous transcription factor regulons grouped into 5 major modules (red boxes) with representative transcription factor (TF) regulons. M. Bubble diagram show the top 10 enriched regulons for the CD4⁺ T cells. N-O. STRING plot depicting confidence of gene-gene interactions among phase-specific regulons of HC(N) and ME/CFS(O) based on the STRING database and gene ontology (GO) terms highlighted. P. Heatmap of up and down regulated genes from the HC group enriched in GO Biological Process (BP) terms in CD4⁺ Naive, CD4⁺ TCM, CD4⁺ TEM, and CD4⁺CTL cells. *p < 0.05, **p < 0.01, ***p < 0.001. note: p means the significance of all paths in the GO Biological Process (BP) pathway. Q. Heatmap of up and down regulated genes from the ME/CFS group enriched in KEGG pathway in CD4⁺Naive, CD4⁺ TCM, CD4⁺ TEM, and CD4⁺ CTL cells. *p < 0.05, **p < 0.01, ***p < 0.001. note: p means the significance of all paths in the KEGG pathway

We then used Monocle3 to visualize the trajectories of CD4⁺ T cells, which revealed different cell state transitions in ME/CFS. The CD4⁺ T cells in ME/CFS were in the mid to late pseudotime, while the HC cells were in the early pseudotime, which is similar to naïve CD4⁺ T cells (Fig. 2H and J). A unique branch of the CD4⁺ TCM cluster, which was absent in HCs, suggests that the differentiation capacity of CD4⁺ TEM is increased in individuals with ME/CFS (Fig. 2G and I).

Based on the differences in trajection patterns, CD4⁺ TCM cells from ME/CFS exhibited a more diverse differentiation status compared to HC cells. Furthermore, we analyzed the regulons (transcription factors, TFs) of the Con and ME/CFS groups separately. The common TFs were categorized into five modules corresponding to specific T cell clusters of both the ME/CFS and Con groups (Fig. 2K and L). In the ME/CFS group, the highly activated TFs were mainly enriched in the regulation of primary metabolic processes (Fig. 2N and O). However, we identified a unique group of highly activated TFs that become disorganized in ME/CFS, including KLF2, JUND, ZBTB7A, USF2, RELA, and KDM5A (Fig. 2M). Importantly, the ME/CFS regulons of IRF1, KLF2, FOS, ZBTB7A, FLT1 and RELA show significant enrichment in cell differentiation, IL-18 and IL-17 signaling pathways.

To investigate transcriptomic variations in ME/CFS, we compared the expression profiles of CD4⁺ T cell. We then performed GO BP and KEGG enrichment analyzes for upregulated and downregulated genes to identify significant signaling pathways. GO BP analysis revealed significant upregulation in T cell activation, including the T cell receptor signaling pathway and positive regulation of T cell activation (Fig. 2P and Q). KEGG analysis revealed significant upregulation in cellular senescence and neurodegenerative signaling pathways, including amyotrophic lateral sclerosis (ALS), prion disease, Parkinson’s disease (PD) and Huntington’s disease (HD). Further classification within the category “human diseases” revealed a significant enrichment of pathways related to neurodegenerative diseases and infectious diseases (Fig. S4A).

We then used a Venn map analysis to determine the major coregulated genes in the GO BP and KEGG pathway (Fig. S4A). CCR7, CD81, PTPRC, LCK, TESPA1, CSK, FYN, ZAP70 and CD28 were upregulated in the pathways responsible for regulating immune responses, including T cell activation and TCR signaling (Fig. S4B). Genes such as ATXN2L, FUS, OPTN, PSMC5, MAP2K2, SDHA, SIGMAR1 and VCP are actively involved in neurodegeneration (multiple diseases) signaling pathways. Genes such as SDHD, NDUFB2 and ATP5PO are integral components of the mitochondrial respiratory chain and enable electron transfer and ATP synthesis.

In addition, we successfully identified subclusters including naïve CD8⁺ cells, central memory T cells (CD8⁺ TCM), effector memory T cells (CD8⁺ TEM) and gamma-delta (γδ) T cells for analysis (Fig. 3A and B and Fig. S5C-5D). Similarly, CD8⁺ T cell also show activation of T cell signaling and upregulation of the neurodegenerative disease pathway. Density plots of CD8⁺ T cell subsets sorted by lineage showed distinct effects of different pools and processing batches on cell distribution in ME/CFS (Fig. 3B). A significant increase (2.27-fold, p = 0.028) in CD8⁺ naïve cells was observed among total CD8⁺ T cells in ME/CFS, with no significant changes in CD8 + TCM cells, CD8⁺ TEM cells and γδ T cells (Fig. 3C and F).

Fig. 3

Cluster analysis of variations in cellular proportion and functional differences among CD8⁺T cell subsets. A. UMAP plot displaying the distribution of CD8⁺ T cell subsets, with each color indicating a specific cell subset. B. Density plots depicting CD8⁺ T cell subset distributions for ME/CFS and HC groups, categorized by ancestry. C-F. Proportion of cell types in each individual, including CD8⁺ Naïve (C, p = 0.039), CD8⁺ TCM (D, p = 0.42), CD8⁺ TEM (E, p = 0.95), and γδT cells (F, p = 0.042). Colors signify cell type information. ns, no significant. ME/CFS group compared to the HC group. H-I. Pseudotime analysis revealing distinct trajectories of CD8⁺T cell subset differentiation in ME/CFS. Cells are color-coded based on CD8⁺T cell subsets, along with pseudotime trajectories for HC group (G-H) or ME/CFS group (I-J). L-L. Regulon modules based on the correlation matrix of shared regulons of HC(K) and ME/CFS(L). The network shows 499 orthologous transcription factor regulons grouped into 4 major modules (red boxes) with representative transcription factor (TF) regulons. M. Bubble diagram show the top 10 enriched regulons for the CD8⁺ T cells. N-O. STRING plot depicting confidence of gene-gene interactions among phase-specific regulons of HC(N) and ME/CFS(O) based on the STRING database and GO terms highlighted. P. Heatmap displaying upregulated genes enriched in Gene Ontology Biological Processes (GO BP) terms for CD8⁺ Naive, CD8⁺ TCM, CD8⁺ TEM, and γδT cells in the ME/CFS group. *p < 0.05, **p < 0.01, ***p < 0.001, p means the significance of all paths in the GO Biological Process (BP) pathway. Q. Vertical axis depicting CD8⁺ naive upregulated genes of the ME/CFS group enriched in annotated KEGG pathways, with the graph indicating the number of metabolites annotated within each secondary classification under the primary pathway classification. R. Bubble chart analysis of upregulated genes enriched in GO BP terms for CD8⁺ Naive cells in the ME/CFS group

The analysis of pseudotime patterns revealed dynamic changes over time in CD8⁺ T cell states (Fig. 3H and J). A more distinct trajectory branch was observed in CD8⁺ naïve cells of the ME/CFS group, which was absent in the HCs. This indicates a potential suppression of differentiation from CD8⁺ naïve to CD8⁺ TCM-like cells in ME/CFS (Fig. 3G and I). Furthermore, a less distinct trajectory branch was noted in CD8⁺ TEM cells of the ME/CFS group, contrasting with a more pronounced branch in the Con group. This suggests a potential suppression of differentiation from CD8⁺ TCM to CD8⁺ TEM cells in ME/CFS (Fig. 3G and I).

To investigate the regulation of gene expression in ME/CFS, we also examined the transcription factors (TFs) of CD8⁺ T cells. The TFs that were common between the two groups were grouped into 4 modules that corresponded to specific T cell clusters (Fig. 3K and M). In healthy controls (HCs), the highly activated TFs were mainly involved in regulating integrated stress response signaling and cytokine signaling in the immune system (Fig. 3N). Notably, we observed higher activation of TFs in ME/CFS, including ZBTB7A, JUND, ESRRA, SP3, and KDM5A, among others (Fig. 3O, Fig. S5F-S5G). Importantly, the highly activated TFs in ME/CFS, such as PSMD8, NFKBIA, PSMD7, NFKB1, REST, TAF7, etc. and MYC, showed significant enrichment in processes related to cell proliferation and differentiation, cellular senescence, infectious diseases, and neurodegenerative diseases like AD and HD, shedding light on potential molecular mechanisms underlying the observed differences (Fig. 3O).

Furthermore, as in CD4⁺ T cells, genes involved in Amyotrophic lateral sclerosis of T cell activation were upregulated in all CD8⁺ T cells in ME/CFS (Fig. 3P). Our focus turned to CD8⁺ naïve cell subgroups. The 662 elevated genes found in this study were linked to 290 KEGG pathways. Further classification of these KEGG pathways indicated a significant enrichment in pathways relevant to human disorders, specifically neurological diseases and infectious diseases caused by both viral and bacterial pathogens (Fig. 3Q). In addition, we used GO BP enrichment analysis on these elevated genes to identify significantly enriched pathways in CD8⁺ naïve cells. Differentially expressed genes (DEGs) such as CCR7, CD28, LCK, and SMAD3 play important roles in immune cell activation, signaling, viral processes, and cytokine synthesis. DEGs such as EZR, CSK, and ACTB, as well as ITGB7, encode proteins that contribute to cytoskeleton dynamics, cell adhesion and motility. Furthermore, proteins encoded by GATA3, ZFPM1, MYC, LEF1, and XBP1 play critical roles in gene regulation and transcription factor function. (Fig. 3R).

Comprehensive insight into B cell dynamics in ME/CFS

To better understand the dynamic changes in various B cell subtypes, B cells were divided into discrete subsets based on the location and expression of traditional subtype markers. Transcriptome analysis allowed us to successfully identify naive B cells, memory B cells, and plasma cells (Fig. 4A and B, S6A-S6C). A trend towards higher plasma cell proportions (of total B cells) was noted, but statistical significance was not achieved due to the small sample size (Fig. 4C and E). To acquire a more complete knowledge of B cell dynamics, we examined DEGs within each subcluster (Fig. S6D-S6F).

Fig. 4

Cluster analysis of the phenotypic and functional differences between B cell subsets. A. UMAP plot depicting the distribution of B cell subsets, with each color representing a distinct cell subset. B. Density plots of B cell subsets for ME/CFS and HC groups, segregated by ancestry. C-E. Proportion of cell types in each individual, including naive B cells (C, p = 0.81), memory B cells (D, p = 0.38), and plasma cells (E, p = 0.25). The different colors indicate cell type information.“ns” indicates no significant difference. The ME/CFS group compared to HC group. F-I. Pseudotime analysis revealing distinct trajectories of B cell subset differentiation in the ME/CFS group. Cells are color-coded by B cell subset and pseudotime trajectory in the HC group (F-G) or ME/CFS group (H-I). J-K. Regulon modules based on the correlation matrix of shared regulons of HC(J) and ME/CFS(K) of NK cells. The network shows 500 orthologous transcription factor regulons grouped into 5 major modules (red boxes) with representative transcription factor (TF) regulons. L. Bubble diagram show the top 10 enriched regulons for the B cells. M. Heatmap of up-regulated and down regulated genes in the ME/CFS group, enriched in GO BP terms in naive B cells, memory B cells, and plasma cells. *p < 0.05,**p < 0.01,***p < 0.001. note: p mean the significance of all paths in the GO BP pathway. N. Venn diagram of the intersection genes of the Figure M red marked signaling pathway. O. Heatmap of up-regulated and down regulated genes in the ME/CFS group, enriched in KEGG terms in naive B cells, memory B cells, and plasma cells. *p < 0.05,**p < 0.01,***p < 0.001. note: p mean the significance of all paths in the KEGG pathway. P. Venn diagram of the intersection genes of the Figure O red marked signaling pathway

Furthermore, B cell state transitions were represented by a trajectory graph based on changes in B cell subset gene signatures. Memory B cells from the ME/CFS group were positioned in early pseudotime, whereas those from the HCs were positioned in mid-to-late pseudotime (Fig. 4G). Notably, a distinct branch appeared in the trajectory map for both the memory B cluster and the plasma B cluster within the ME/CFS group, but not in the HCs (Fig. 4I). These findings indicate that this plasma cell cluster may represent a distinct subset of circulating B lymphocytes specific to ME/CFS (Fig. 4F and H). The ME/CFS group’s naïve and memory B clusters had lower differentiation status compared to those from HCs, based on trajectory inference differences. In contrast to T cells, mature effector B cells (plasma cells producing antibodies) appeared to be more abundant in ME/CFS. This observation could indicate an excessive ongoing activated posture in B lymphocytes, with a higher proportion of plasma cells possibly originating from memory subsets rather than effector subsets in response to the present pathogen or stressor.

Regulon analysis revealed that shared regulons were arranged into five modules that corresponded to specific B cell clusters in all individuals (Fig. 4J and K). Representational difference analysis identified putative important regulons that control B-cell growth in ME/CFS, including KLF2, TFDP1, MBD4, and TFDP1 (Fig. 5L, S6G-S6H).

Fig. 5

Analysis of phenotypic and functional differences among NK cell subsets. A. UMAP plot displaying the distribution of NK cell subsets, with each color denoting a distinct cell subset. B. Density plots illustrating NK cell subsets in the ME/CFS and HC groups, differentiated by ancestry. C-E. Proportion of cell types in each individual, encompassing NK (C, p = 0.16), NK-CD56^bright (D, p = 0.28) and NKT (E, p = 0.031) cells. Colors indicate cell type information. *p < 0.05, ns, no significant. D-I. Pseudotime analysis unveiling distinctive trajectories of NK cell subset differentiation in the ME/CFS context. Cells are color-coded based on NK cell subsets, along with their pseudotime trajectory in HC groups (F-G) or ME/CFS groups (H-I). J-K. Regulon modules based on the correlation matrix of shared regulons of HC(J) and ME/CFS(K) of NK cells. The network shows 495 orthologous transcription factor regulons grouped into 5 major modules (red boxes) with representative transcription factor (TF) regulons. L. Bubble diagram show the top 10 enriched regulons for the NK cells. M-N. STRING plot depicting confidence of gene-gene interactions among phase-specific regulons of HC(M) and ME/CFS(N) based on the STRING database and gene ontology (GO) terms highlighted. O. Heatmap analysis depicting NK cells upregulated and downregulated genes enrichment of GO BP terms pathways of ME/CFS group. *p < 0.05,**p < 0.01,***p < 0.001. note: p mean the significance of all paths in the GO BP pathway. P. Heatmap analysis depicting NK cells upregulated and downregulated genes enrichment of KEGG pathways of ME/CFS group. *p < 0.05,**p < 0.01,***p < 0.001. note: p mean the significance of all paths in the KEGG pathway. Q-S. The frequency of Granzyme B (Q), CD107a (R) and Perforin (S) in NK cells from PBMCs of ME/CFS patients. NK cells from ME/CFS patients were cocultured with target cells for 6 h. Flow cytometry assay(FACS, n = 4 vs. 4) was used to detect cell function in target cells. T. NK cells from ME/CFS patients were cocultured with target cells for 6 h. Flow cytometry assay was used to detect 7-AAD staining in target cells

We next used GO BP and KEGG pathway enrichment to look for shared DEGs between the three B cell subsets. As expected, the GO BP analysis revealed significant upregulation in signaling pathways related to B cell activation, positive regulation of cell activation, immune response regulation, immune response-activating cell surface receptor signaling, receptor signaling, and regulation of B cell activation in ME/CFS B cell subsets (Fig. 4M). A Venn map indicated common regulated genes such as CD81, PTPN6, INPP5D, PTPRC, IGHG4, and IGHA2 (Fig. 4N). These genes jointly contribute to several areas of B cell function, such as signaling, immune response control, and antibody synthesis, demonstrating their significance in sustaining a strong and effective B cell-mediated immune response. Furthermore, KEGG pathway analysis revealed increase in signaling pathways linked with Parkinson’s disease, Prion disease, ALS, and Epstein-Barr virus (EBV) infection (Fig. 4O). The Venn map revealed that the specific genes BAX, ADRM1, PSMC5, and PSMD1 play critical roles in regulating crucial biological processes within B cells. Together, these genes highlight the complex regulatory processes that keep B cells functioning and homeostasis. (Fig. 4P).

NK cell dynamics and functionality in ME/CFS

In addition to T cells and B cells, natural killer cells (NK cells) are the third main class of lymphocytes. They are involved in the regulation of the immune system, the prevention of viral infections, and the treatment of tumors. Additionally, they are involved in the development of autoimmune diseases and hypersensitivity. We were able to effectively identify classical NK cells (NK-CD56^dim), NK-CD56^birght, and included NKT cells in this subpopulation analysis (Fig. 5A and Fig. S7A-S7D). Density plots representing NK cell subsets, segregated by ancestry, displayed observable effects differences from various pools and processing batches on cell distribution (Fig. 5B). Although this trend did not achieve statistical significance, we observed a discernible decrease in the proportion of CD56^bright NK cells (relative to total NK cells) in the ME/CFS group. Interestingly, we observed a statistically significant and noteworthy increase in the frequency of NKT cells (3.15-fold, p = 0.017) (Fig. 5C and E).

In contrast to HCs, which occupied a mid-to-late pseudotime trajectory (Fig. 5G), NKT and NK-CD56bright cells from the ME/CFS patients were positioned along the early pseudotime trajectory (Fig. 5E) in single-cell trajectories. This observation suggests that the NK-CD56^bright and NKT clusters exhibited suppressed differentiation in ME/CFS (Fig. 5F and H). We also conducted a thorough examination of the highly activated TFs in ME/CFS. The shared TFs were arranged into five modules that corresponded to particular NK cell clusters (Fig. 5J and K). In ME/CFS, we identified TFs that were highly activated, including key factors such as ZBTB7A, KLF2, KLF13, RELA, CEBPB, and CEBPD (Fig. 5L, Fig. S7H-S7I). The ME/CFS-specific highly activated TFs, which include ZBTB7A/7B, RELA, KLF2, NFKB1, ASCL2, and POLE3, exhibited a notable enrichment in cellular processes such as population proliferation, cell differentiation, nucleic acid metabolic processes, cellular response to stress, and pathways such as IL-17, IL-18 signaling, and IL-6 production (Fig. 5M). In contrast, the genes FOS, JUN, JUNB, IRF1, and BCLAF1 were the most frequently associated with highly activated TFs in HCs. This enrichment was limited to processes related to the regulation of integrated stress response signaling, leukocyte differentiation, and cytokine signaling in the immune system (Fig. 5N).

Furthermore, the 603 upregulated and 213 downregulated DEGs in ME/CFS demonstrated enrichment in biological processes that were predominantly associated with the positive regulation of the response to external stimuli, immunoglobulin-mediated immune response, antigen processing, and presentation of exogenous antigens, as well as the regulation of adaptive immune response based on somatic recombination of immune receptors (Fig. 5O). The KEGG enrichment analysis of upregulated genes in the ME/CFS group revealed associations with disease pathways related to infectious diseases (HIV, EBV infection, COVID-19, etc.), Amyloidotic neurodegenerative diseases including ALS and PD, and autoimmunity. Conversely, downregulated genes were involved in IL-17 signaling and NK cell-mediated cytotoxicity (Fig. 5P).

Afterward, we implemented flow cytometry to evaluate the functional status of NK cells in ME/CFS. Granzyme B expression was at the lower limit of the range in NK cells and unaffected in CTLs (ΔGranzyme B NK = 77.09%, Normal ≥ 77%; ΔGranzyme B CTL = 19.86%, Normal ≥ 6%) (Fig. 5Q). The CD107a surface expression of NK cells from ME/CFS was significantly impaired (ΔCD107a = 11.05%, Normal ≥ 10%) in comparison to CTLs (ΔCD107a = 11.05%, Normal ≥ 10%) upon PBMC stimulation with PMA/Ionomycin (Fig. 5R). Furthermore, the expression of perforin was significantly diminished in NK cells (ΔPerforin = 49.58%, Normal ≥ 81%) (Fig. 5S). Furthermore, NK cells from ME/CFS patients exhibited significantly lower cytotoxicity against target cells (CD4⁺ T cells) compared to normal (13.48%, Normal ≥ 15.11%), which provides evidence that NK functions were impaired in ME/CFS patients (Fig. 5T).

Upregulated genes enriched in phagocytosis related function in monocytes

Monocytes, which are essential to the immune system, are capable of fulfilling a variety of functions. Recognized for their exceptional phagocytosis, they protect themselves by ingesting pathogens and facilitating specific immune responses through antigen presentation. Multiple combinations of receptors are utilized by these cells to coordinate immune communication through the production of cytokines. They facilitate tissue repair by transforming into macrophages or dendritic cells. Inflammation is crucially balanced by monocytes, which are responsible for the initiation and regulation of immune responses. The distribution and expression of classical subtype markers were used to classify monocyte clusters into three subsets in order to reveal the fractional abundances of monocyte cell states (Fig. S8A-S8E). According to this classification, classical (CD14) monocyte cells, non-classical (CD16) monocyte cells, and conventional DC cells (cDC) were identified (Fig. 6A and B). Although trends suggested a reduction in CD16 and cDC cells in the ME/CFS group, statistical significance was not achieved (Fig. 6C and E).

Fig. 6

Cluster analysis of phenotypic differences in monocyte cell subsets. A. UMAP plot displaying the distribution of monocyte cell subsets, with each color representing a specific cell subset. B. Density plots of monocyte cell subsets in ME/CFS and HC groups, categorized by ancestry. C-E. Proportion of cell types in each individual, including CD14 monocytes (C, p = 0.22), CD16 monocytes (D, p = 0.35), and cDC cells (E, p = 0.18). The colors indicate cell type information. ns, indicates no significance. ME/CFS group compared to HC group. F-I. Pseudotime analysis revealing distinct trajectories of monocyte cell subset differentiation in ME/CFS. Cells are color-coded by monocyte cell subsets and pseudotime trajectories in HC group (F-G) or ME/CFS group (H-I). J-K. Regulon modules based on the correlation matrix of shared regulons of HC(J) and ME/CFS(K) of monocyte cells. The network shows 530 orthologous transcription factor regulons grouped into 5 major modules (red boxes) with representative transcription factor (TF) regulons. L. Bubble diagram show the top 10 enriched regulons for the monocyte cells. M-N. STRING plot depicting confidence of gene-gene interactions among phase-specific regulons of HC(M) and ME/CFS(N) based on the STRING database and gene ontology (GO) terms highlighted. O-P. Heatmap analysis illustrating the upregulated and downregulated genes enriched in GO BP (O) and KEGG(P) terms in monocyte cells of ME/CFS group. *p < 0.05,**p < 0.01,***p < 0.001. note: p mean the significance of all paths in the GO BP or KEGG pathway datebase

In order to investigate developmental connections among monocyte cell states, we implemented Monocle3 to generate single-cell trajectories. In the trajectory map, the CD14 monocyte cells in the ME/CFS group exhibited a reduced number of branches. This finding suggests that the differentiation diversity of CD14 monocyte cells in the ME/CFS group has decreased (Fig. 6F and I). They are distinguished from HCs by the reduced variety in their differentiation status, as demonstrated by their less distinct trajectory patterns.

In order to elucidate the complexities of gene expression regulation in ME/CFS, we conducted a comprehensive examination of regulons. The highly activated TFs were categorized into five modules that corresponded to particular monocyte cell clusters (Fig. 6J and K). One notable observation was the increased activity of transcription factors (TFs) in the ME/CFS group, which included pivotal factors such as ZBTB7A, KLF13, RXRA, IRF2, and SPI1, among others (Fig. 6L, Fig. S5F-S5G). The ME/CFS highly activated TFs, which include STAT3, IRF7, MYC, CEBPD POU2F2, ESRRA, CEBPB, CEBPD, and SPI1, exhibited a notable enrichment in cellular processes such as the regulation of monocyte differentiation, EBV infection, prion disease pathway, overlap between signal transduction pathways contributing to LMNA laminopathies, nuclear receptors, and pathways like IL-17, IL-2, and IL-5 signaling (Fig. 6M). Conversely, regulons in HCs were primarily associated with the genes EGR1, EGR2, FOSB, EGR3, and DDIT3, indicating an enrichment in processes that were exclusively associated with the response to endogenous stimulus (Fig. 6N).

Additionally, we investigated the 581 upregulated DEGs in monocytes by utilizing GO BP and KEGG pathway enrichment. We observed that these DEGs were more actively involved in critical processes, including phagocytosis, endocytosis, and Fc gamma R-mediated phagocytosis. The signaling share the specified genes, FCGR1A, RAC1, PTPRC, PRKCD, LYN, GSN, and PAK1, collectively contribute to the crucial biological processes in monocytes (Fig. 6O and P). Collectively, these DEGs and regulatory patterns demonstrated substantial shifts in ME/CFS based on monocyte cell behavior, such as migration, phagocytic activities, and immune response modulation, underscoring their importance in the innate immune system.

Cell-cell communications in ME/CFS indicates that APP is core signaling of the observed interactions between monocytes and other cells

The coordination of cell-cell communications and interplays among the immune cell subpopulations is essential for the maintenance of a normal immune response to extrinsic antigens and homeostasis to intrinsic stresses. For instance, the orchestration between T cells and antigen presenting cells (APCs) is a prime example. The intricate interplay of signaling pathways in cell-cell communication in ME/CFS is still not fully comprehended. We endeavored to depict the cell-to-cell communication relationships among the 12 main cell types in a circle plot network (Fig. 7A and B) in this study. The cell communication networks within monocyte clusters were significantly enhanced in comparison to other cell clusters in the ME/CFS group. This included “CD4 Native,” “CD4 TCM,” “CD4 TEM,” “CD4 CTL,” “CD8 Naïve,” “CD8 TCM,” “CD8 TEM,” “γδT,” and “B cell.” The inferred incoming and outgoing interaction strengths are also elaborated in Fig. 7C and D, where monocytes are the central cell type that express ligands and receptors that are actively altered and involved in the cellular interactions in ME/CFS. These findings indicate that monocytes may play a critical role in the unprofessional interactions observed in ME/CFS.

Fig. 7

Cell-cell communication networks. A-B. Circle plot depicting cell-to-cell communication relationships between 12 main cell types and monocyte cluster cells of the HC(A) and ME/CFS(B) group. The quantitative network diagram uses nodes to represent different cell types, arrows to indicate interaction signals from ligand cells to recipient cells, and line thickness to denote the number of significant ligand-receptor interaction pairs detected between different cell types. C-D. Illustration of the incoming and outgoing interaction strengths for each of the cell types of the HC(C) and ME/CFS(D) group. E-F. The outgoing signaling pathways of each cell type of the HC(E) and ME/CFS(F) group. G-H. The incoming signaling pathways of each cell type of the HC(G) and ME/CFS(H) group. I. Significant signaling pathways ranked based on differences in overall information flow within inferred networks between ME/CFS and HC. Enriched signaling pathways more prominent in ME/CFSS are shown in blue. J. Chord diagrams of APP signals pathway from monocyte cells to other cell types in the ME/CFS group. monocyte cells are color matched to the origin, with each arc representing one pathway and arc length depicting the strength. K. The heatmap shows the communication probability of the APP signaling pathway. L-M. UMAP plot displaying the APP distribution of PBMC(L) and monocyte(M) cell subsets(L). N-O. Violin plots of APP expression in each ME/CFS patients(N) and group(O) of monocyte cells. P. Immunohistochemical staining reveal the distribution and expression of the APP. Q. Quantification of protein expression shown in figure P. ***p < 0.001, ns p>0.05, ME/CFS group(n = 3) compare with HC group (n = 3). R. Representative western blot images showing the expression of APP in the PBMC cells of the ME/CFS and HC group. β-actin was used as a loading control. S. Quantification of Western blot bands with ImageJ software. The bar graph shows the intensities of the immunoblot bands in Figure R, which were quantified using ImageJ software. ***p < 0.001, ME/CFS group (n = 3) compare with HC group (n = 3). T. The gene regulatory networks for monocyte cells enriched regulons of APP. U. Violin plots of ESRRA expression in total group of monocyte cell types (n = 4). V. Bar graph of ESRRA expression in male patients’ monocyte cell types from the GSE214284 database. *p < 0.05, ME/CFS group(n = 10) compare with HC group(n = 10)

In order to identify the signaling pathways that are responsible for the intricate intercellular communication, we conducted additional investigations into the intensity of each signaling pathway’s cellular outgoing and incoming interactions. The potential molecular contributors to the exaggerated cell communications in monocytes were identified as signaling pathways involving CD99, ITGB2, ICAM, CD23, MHC-II, MHC-I, CD45, CD22, SELPLG, LCK, IL16, VISFATIN, BAG, BAFF, GRN, MIF, CD86, SEMA4, CD40, PARs, CD48, RESISTIN, CCL, and APP (Fig. 7E and H). In addition, we investigated the specific ligand-receptor interactions between monocyte cell clusters and other immune cell clusters. In comparison to the Con group, the ME/CFS group exhibited a significant upregulation of ligand-receptor pairs, including APP-CD74, HLA-DRB5-CD4, HLA-DMA-CD4, HLA-DMB-CD4, HLA-DPA1-CD4, HLA-DPB1-CD4, HLA-DQA1-CD4, HLA-DQA2-CD4, HLA-DQB1-CD4, HLA-DRA-CD4, HLA-DRB1-CD4, IL16-CD4, and RETN-CAP1 (sFig.9 A).

In monocyte cells, APP signaling is the predominant pathway, as evidenced by its substantial outgoing and incoming interaction strengths in comparison to other signaling pathways (Fig. 7I). This is a significant finding. The APP signaling was also one of the distinctive features of ME/CFS, but it did not occur in HCs. Furthermore, it appears that the Amyloid-beta precursor protein (APP)/CD74 pathway is primarily involved in the interactions between monocytes and other cells in ME/CFS (Fig. 7J and K, sFig.11 A). CD74 is a membrane protein that functions as an MHC class II chaperone, thereby facilitating the antigen presentation process. CD74’s function may be linked to neuroinflammatory mechanisms or aberrant antigen presentation in the context of conditions such as autoimmune disorders or neurodegenerative diseases. These data collectively suggest that monocytes may play a critical role in the expansion and malfunction of lymphocytes (T and B cells), which would contribute to the unexpected yet critical discoveries regarding amyloidotic neurodegenerative signaling in immune cells that extend beyond immune disturbance.

In order to independently verify our findings, we conducted an analysis of single-cell data for APP expression, which demonstrated its prevalent presence in monocyte cells among all monocyte subtypes (Fig. 7L and M). The violin diagrams emphasize the increased expression of APP in monocyte cells of the ME/CFS group (Fig. 7O and P). Following this, we conducted immunohistochemistry (IHC) on APP, Aβ (Amyloid-beta), and Tau, and discovered a significant increase in APP protein accumulation, which is an unprecedented discovery in ME/CFS (Fig. 7P and Q, SFig. 11B). Western blot analysis further verified the elevated APP expression in PBMCs from ME/CFS patients (Fig. 7R and S). In addition, we identified 22 regulons that have the potential to affect APP expression through additional regulon analyses (Fig. 7T). The significant correlation between ESRRA and APP expression is indicated by the modifications to higher regulons that are unique to ME/CFS. This suggests that the production of Amyloid β peptides in ME/CFS may be influenced by ESRRA’s influence on APP gene expression. Estrogen-related receptor-α (ESRRA) is a nuclear receptor of transcription factor that binds to estrogen and estrogen-related responsive elements. It was recently identified as a molecular immune-metabolic antitumor target. ESRRA encodes this receptor. ESRRA has been demonstrated to be one of the most highly enriched DEGs in the aforementioned cell types. Our scRNA-seq analysis also corroborated these findings by verifying that the expression of ESRRA was higher in each ME/CFS patient than in HCs (Fig. 7U, SFig.11 C). Furthermore, our observations were further validated by the ESRRA expression in the monocyte cell types of male patients from the GSE214284 database, which yielded the same results (Fig. 7V). In the context of ME/CFS, these compelling discoveries establish a preliminary but novel link between cellular components, immune disturbance, and neuro symptoms.

Discussion

A significant amount of evidence indicates that the pathogenesis of ME/CFS is the result of perturbative immune and inflammatory responses [17, 33]. Nevertheless, the current data is inconclusive and discrepant in relation to the peripheral immune cell components. Numerous studies have reported conflicting results, which may be due to the limited resolution of low-throughput, mixed-in bulk-analysis of the screening assays, or the significant heterogeneity and unpredictable dynamics of peripheral immune cells [34, 35]. In addition, there are limited studies that have investigated the circulating immune landscape in patients with a single-cell resolution [36]. This ignorance has impeded the development of precise diagnostic tests and treatment targets. In this investigation, we employed scRNA-seq to analyze the intricate immune responses of PBMCs in the constant state of ME/CFS. It is crucial to note that we excluded individuals with a history of taking immunomodulatory medications, those with evident comorbidities, and those with evidence of ongoing acute or chronic infection from this study in order to prevent confounding variables that could potentially affect immune cell profiles. As an illustration, our investigation was conducted prior to the participants’ exposure to COVID-19, which transpired to have the potential to influence their immune profiles and induce symptoms similar to those of chronic fatigue syndrome, or “long-Covid.” The chronic nature of the disease is underscored by the presence of immunological defects in ME/CFS patients, as evidenced by the data previously mentioned [29].

The classifications of peripheral immune cell subsets in ME/CFS individuals were characterized by distinct characteristics in this study. Our preliminary investigation has revealed that the remarkable alterations in the populations of lymphocytes, monocytes, and NK cells, which are indicative of altered developmental processes and functional status. The most notable changes were observed in the increased total numbers of T cells and T cell subcluster composition in ME/CFS in this study. In ME/CFS, we observed unusually elevated numbers and percentages of naïve T cells (both CD4⁺ and CD8⁺), absolute numbers of central memory T cells, and decreased percentages of CD4⁺ effector memory T cells. This may indicate that T cells in ME/CFS exhibit a non-antigen-activated early expansion, but they are relatively insensitive to pathogens for the development of effector T cells. The other main components of the adaptive immune system, the B cells, also underwent similar changes, resulting in an increase in the production of plasma cells that produce antibodies. And we observed substantial decreases in the proportions of monocytes and NK cells. Nevertheless, these discoveries have not yet been able to produce a definitive conclusion in the field of immunopathology. It is imperative that we reevaluate these unique phenotypes in the context of pathogen-immune interaction, cellular processes or pathway regulation, and cell-cell dynamics. The dysregulation of classical monocytes in ME/CFS patients is consistent with the findings of the previous study, suggesting that there may be a malfunction in the differentiation and migration of these cells to tissues [29].

It is also crucial to note that we discovered significant discrepancies in gene expression patterns and cellular biological processes in ME/CFS that were rarely examined in previous research. Three aspects of the main findings are particularly noteworthy. First and foremost, we identified numerous transcription factors and the interactive gene expression that drives the differential processes of immune cells. One example is the substantial upregulation of genes such as CD81, H1FX, and KLF2, and the downregulation of CXCL8 and NFKBIA across all T cell subclusters. The preposition of unspecific T cell activation and chronic inflammatory status in ME/CFS is further justified by the critical functions of these genes in immune regulation and signaling. In CD8 TEM, naïve CD4, CD4⁺ CTL, γδT, and NKT cells, we also observed down-regulated gene expressions of HLA-A, -B, -C, and FCGR3B, among others. NK cells are recruited and play a dual function in the limitation of tissue autoimmunity and the clearance of viral infections. This function is essential for the regulation of inflammation, the restoration of immune equilibrium, and the preservation of tissue health, and it is contingent upon NK targeting T cells [37].

Secondly, we were taken aback by the discovery of numerous cellular pathways that are implicated in certain well-known diseases. Our attention was drawn to two classes of diseases within the “Human Disease” KEGG category. One is the virus infection-related pathways, which include non-specific viral processes and specific virus-specific pathways (such as HSV, HPV, CMV, EBV, or COVID-19). It is important to acknowledge that the current viral infections have been excluded from the participants through NGS sequencing and serological screening. Consequently, the clinical manifestations of ME/CFS are more likely to be consequential than etiological or pathogenic, as the activation of these disease-related pathways. The second unexpected category of diseases consisted of amyloid neurodegenerative disorders, such as AD, ALS, PD, HD, or Prion diseases. It is intriguing that the aforementioned pathways, which involve both infectious and noninfectious amyloidosis, were “coincidentally” upregulated in specific peripheral immune cell subsets in ME/CFS. It is crucial to note that we revised the data and identified genes that are similar but distinct (e.g., ADRM1 and PSMD2) and are implicated in these pathways. ADRM1 and PSMD2 are both linked to the 26 S proteasome, a substantial protein complex that is responsible for the degradation of ubiquitinated proteins in eukaryotic cells [38].

Finally, the gene expression patterns substantially increased the signal pathways involved in cell senescence and exhaustion, particularly in T cells. The state of cellular senescence is characterized by the irreversible arrest of the cell cycle and is linked to chronic inflammatory conditions and aging [39, 40]. The notion that immune dysfunction may be associated with chronic inflammation and aging-related processes is emphasized by the presence of senescent T cells in ME/CFS patients. These data partially elucidated the anergy of immune cells to specific pathogens, particularly viruses that are frequently infected, despite the aberrant expansion of cell deposits and the upregulation of viral-related pathways previously described. The complexity of immune perturbations in ME/CFS is underscored by the involvement of the cellular senescence, T cell receptor signaling, NF-κB pathway, and pathways related to infectious and neurodegenerative diseases. In this condition, immunodeficiency and autoimmunity appear to be the two sides of the same coin. Similar to the previous study, both CD4⁺ and CD8⁺ T cells from ME/CFS patients exhibit reduced glycolysis at rest, with CD8⁺ T cells also exhibiting diminished glycolysis following activation. This implies that the impaired immune function of ME/CFS patients may be attributed to a fundamental alteration in energy metabolism in T cells. Consistent with previous researches, ME/CFS patients often exhibit an abundance of exhausted T cells—a reversible state associated with chronic viral infections that impairs cell proliferation, survival, and cytotoxicity [17, 41]. In general, our analyses enhance our understanding of the immunological changes in ME/CFS, which will be beneficial for future research in the development of potential therapeutic strategies and the identification of the underlying causes of the debilitating symptoms.

Another major intriguing finding in our study is clues of the distinguished trajectory of development, activation, expansion and differentiation of immune cell subsets in ME/CFS. These processes were first visualized in Pseudotime analysis. We found a unique branch in both naïve CD4⁺ and CD8⁺ T cells and CD4⁺ TCM cells. More importantly, we unveiled specific lineage development of CD4⁺ CTL originated directly from TCM-TEM, and trajectory connecting CD8⁺ naïve T cells to TCM cells. Interestingly, we also identified a unique group of expanded plasma (mature antibody-producing B) cell population that was originated directly from memory B cells. These implied possibly excrescent developmental links in the adaptive immune cell subpopulations. Overactive B cells can secrete pro-inflammatory cytokines, such as IL-6 and TNF-α, which contribute to the chronic inflammatory state. If B cell-driven autoimmunity and inflammation are confirmed as key mechanisms in ME/CFS, treatments that modulate B cell activity, such as B cell depletion therapy or immunomodulatory drugs, could offer new therapeutic avenues. Additionally, identifying specific autoantibodies in ME/CFS patients could lead to better diagnostic markers and personalized treatment approaches. Further research is necessary to fully understand these connections and to translate this knowledge into effective clinical interventions.

Comparatively, in health control, CD4⁺ TEM, CD4⁺ CTL and plasma cells were more dependent from TCM and memory B cell subpopulation, respectively. However, lack of developmental divergencies were seen in the B cells, and the effector and innate immune populations, NK and monocytes. These may again reflect autoimmunity-like recognition (stress-sensing or antigen processing and presenting) or recurrent activation, but yet deficiency in normal cytotoxicity or effective clearance of pathogens. Moreover, we further analyzed expressions of dominant TFs and their corresponding in activation but gene-gene interactions among phase-specific regulons. Apart from several increased expressions of TFs, like AP-1(JUND, FOS), NFκB (RELA), and KLF2 above-mentioned, we also unveiled surprisingly more intricated “input signals” in ME/CFS. In regulons of health controls, it was solely primary metabolic processes that exerted a dominant role in governing the development status of CD4⁺ T cells. However, several interactive genes involving inflammation, cytokine or stress pathways babelized with these processes in ME/CFS.

Similar or even more chaotic situations also occurred in CTL cells, B cells, NK cells, and monocytes. NK cells are critical for the early control of viral infections. Reduced NK cell function in ME/CFS patients could lead to prolonged or chronic viral infections, which might exacerbate or perpetuate the disease. In ME/CFS, several studies have consistently reported impaired NK cell activity, characterized by reduced cytotoxicity and altered cytokine production [14, 42]. Monocytes have the ability to cross the blood-brain barrier, particularly when activated [43]. Once in the central nervous system (CNS), they can differentiate into macrophages and contribute to neuroinflammation [44]. In ME/CFS, neuroinflammation has been proposed as a significant contributing factor to cognitive dysfunction, pain, and sleep disturbances. Activated monocytes and the cytokines they release can interact with microglia, the resident immune cells of the brain, promoting a neuroinflammatory response that may underlie some of the neurological symptoms of ME/CFS. Studies have also identified alterations in monocyte subsets in ME/CFS patients [27, 29]. The shift in monocyte subsets may therefore play a central role in sustaining the chronic immune activation observed in ME/CFS. The identification of these targets provides a foundation for the development of targeted therapies aimed at addressing the underlying immune dysfunction in ME/CFS. Continued exploration of these areas may lead to significant advancements in the management of ME/CFS.

Finally, we have described a specific signaling pathway that is particularly appealing in cell-to-cell communication. The intricate interplay of signaling pathways in cell-cell communication in ME/CFS is still not fully comprehended. Initially, we observed heightened communication networks between monocyte clusters and other cell populations, which confirmed that the monocyte was the central cell type involved in the cellular interactions in ME/CFS. CD74 is also known as the invariant chain of the MHC class II complex, playing a critical role in antigen presentation and immune response modulation [45]. APP-CD74 was the most specific and predominant contributor to exaggerated cell communication in monocytes among the various ligand-receptor pairings. Recent research has further supported the relevance of this pathway: Matsuda et al. demonstrated that CD74 interacts with APP and suppresses the production of Abeta [46]. Xu et al. highlighted the role of APP in blood vessels and its implications for neurodegenerative diseases [47], while Ma et al. showed that therapeutic modulation of the APP-CD74 axis can activate phagocytosis in glioblastoma [48]. Its interaction with APP in monocytes could be a mechanism driving chronic immune activation in ME/CFS, contributing to the systemic inflammation observed in the disease. This pathway could be responsible for maintaining a pro-inflammatory state, which might lead to the persistent fatigue and muscle pain experienced by patients.

We further confirmed our discovery of APP aggravation through in situ histology and WB analysis. The Amyloid-beta precursor protein was the most plausible explanation linking immune and metabolic disturbance with neurologic manifestations, in light of the aforementioned unexpected yet crucial findings regarding amyloidotic neurodegenerative signaling in immune cells. Furthermore, subsequent regulon analyses revealed an orphan nuclear receptor, ERRRA, as the primary regulator of APP transcription in monocyte cells, providing evidence for its potential involvement in the pathogenesis of ME/CFS. Furthermore, our observations were further validated by the ESRRA expression in the monocyte cell types of male patients from the GSE214284 database, which corroborated our findings. These data, when combined, revealed a previously unknown elevated ESRRA-APP-CD74 signaling in ME/CFS. Targeting the APP-CD74 signaling pathway presents a novel therapeutic opportunity in ME/CFS. If this pathway is indeed contributing to both immune dysregulation and neurodegeneration, modulating its activity could potentially alleviate symptoms and slow disease progression. For instance, small molecules or antibodies that inhibit the APP-CD74 interaction could be explored as potential treatments. Furthermore, reducing APP expression or blocking its downstream effects might mitigate the neuroinflammatory aspects of ME/CFS. This finding has the potential to provide insight into the transition from pathophysiologic alterations to clinical manifestations and to be further developed into a novel peripheral blood-based biomarker of ME/CFS .

Our investigation is subject to numerous constraints. This limited sample size compromises the reliability and generalizability of the biomarker findings, constituting an exploratory nature of our findings. To strengthen our findings, we conducted cross- validation using publicly available scRNAseq datasets from an independent cohort (GSE214284). We were unable to incorporate female data into the analysis due to the GSE214284 dataset’s lack of comprehensive information and the influence of hormone levels on ESRRA. It is also important to mention that our study did not categorize ME/CFS patients according to specific stages or subtypes, such as post-infectious ME/CFS or ME/CFS with comorbidities. In addition, the cross-sectional nature of this investigation results in a static perspective that fails to provide any insights into the changing immune cell profiles over time. This heterogeneity has the potential to obscure immune cell alterations that are associated with distinct subtypes and diverse periods. While scRNA-seq provides invaluable insights into gene expression profiles, it does not directly assess the functional activity of immune cells. In order to overcome these limitations, it is essential to conduct additional research that incorporates functional assays, longitudinal methodologies, and a broader and more diverse patient cohort to substantiate and broaden our findings.

The study has established an intricate atlas of the immune landscape, which encompasses the development, activation, and differentiation of immune cells, as well as the dynamics of TF or function genes. Additionally, we have examined the cellular constituents of subsets and, most importantly, cell-cell communications. The results presented here have the potential to provide new perspectives on the immunological complexities of ME/CFS. In particular, we identified an excessive communication that was initiated by monocytes and transmitted to other immune cell components through ESRRA-APP-CD74. This communication could serve as a potential biomarker and facilitate the development of molecular diagnostic tools.