Single-cell RNA-sequencing of myasthenia gravis reveals transcriptional heterogeneity and dysfunction of immune cell populations

Patients

This study recruited two MG patients from the department of thoracic surgery of Beijing Hospital (Beijing, China), and two age- and sex-matched HCs. The diagnostic criteria of MG were as follows: (i) clinical symptoms; (ii) neostigmine test result; (iii) repetitive nerve electrical stimulation; and (iv) serum autoantibody levels.

Participants were excluded if they had (i) any other neuromuscular conduction disorders, such as muscular dystrophy syndrome; (ii) autoimmune diseases, such as systemic lupus erythematosus, rheumatoid arthritis, etc.; (iii) neuromuscular conduction disorders caused by long-term use of certain anticholinesterase drugs or antiarrhythmic drugs (such as β-blockers); (iv) other malignant diseases, acute and chronic inflammation, and/or hepatorenal dysfunction; or (v) MG caused by other causes, such as congenital muscular dystrophy.

The data of HC thymus donor were obtained from: https://developmental.cellatlas.io/studies/thymus-development#datasets.

Sample collection and single-cell suspension preparation

Thymic tissue and venous blood from the same donors were collected from two patients with MG within 24 h after admission. Thymic tissue samples were minced into approximately 1 mm³ fragments on ice using sterile surgical blades. The minced tissue was then subjected to enzymatic digestion by incubation in 0.25% (w/v) trypsin-EDTA solution at 37°C for 20 min with constant agitation. Following digestion, the resulting cell suspension was filtered through a 70 µm cell strainer to remove undigested debris. Both the thymic tissue and venous blood of the above two patients with MG were centrifuged at 400 g for 5 min. Then 10× red blood cell lysis buffer (Solarbio, Beijing, China) was added to the cell precipitation, which was incubated on ice for 10 min to remove red blood cells. The cell pellets of thymic tissue and venous blood were resuspended in Hanks' buffer before the mononuclear cells were isolated using Ficoll-Paque (TBD, Tianjin, China) according to the manufacturer’s instructions.

ScRNA-seq library preparation and sequencing

The scRNA-seq library was prepared with the 10×Genomics Chromium single-cell protocol using the v3 reagent kit (10×Genomics, Pleasanton, CA, USA), according to the manufacturer’s instructions. 90 μl cellular suspension, 40 μl barcoded Single Cell 5′ Gel Beads, and 270 μl Partitioning Oil were loaded onto a Chromium Chip A to generate single-cell gel bead-in-emulsions (GEMs). GEM-reverse transcriptions (GEM-RTs) were performed in a Veriti 96-well thermal cycler (Thermo Fisher Scientific, Waltham, MA, USA) together with the RT master mix. After RT, cDNAs were amplified (98°C for 45 s; 13–18 cycles of 98°C for 20 s, 67°C for 30 s, 72°C for 1 min; and 72°C for 1 min) and then used for 5′ gene expression library construction with the Chromium Single Cell 5′ Library Construction Kit (10×Genomics, Pleasanton, CA, USA). The single-cell RNA library was sequenced on an Illumina NovaSeq 6000 with a read length of 150 bp paired-end reads.

ScRNA-seq data analysis

Pseudotime trajectory analysis

The pseudotime developmental trajectories for CD4⁺ T cells and CD8⁺ T cells were constructed using R library Monocle3 (18) according to the standard workflow. After dimensionality reduction, cells were ordered in a temporal differentiation sequence based on the significant marker genes identified by FindAllMarkers function in Seurat.

RNA velocity analysis

The velocyto.py command-line tool (19) was used to annotate the spliced and unspliced reads from binary alignment/map files generated by Cell Ranger (v 6.0.1) software from 10× Genomics, with outputs stored as loom files. Specifically, GRCh38 reference genome annotation and expressed repeats annotation files were required to align reads to and mask expressed repetitive elements. Loom files were then processed using velocyto.R (v 0.6) to estimate RNA velocity based on gene structure parameters. The RNA velocity was visualized on the uniform manifold approximation and projection (UMAP) plot using Monocle DDRTree embedding method.

Interactions between cell types

Cell-cell interactions and key signaling pathways were investigated using CellChat (20). Quantitative contrasts and manifold learning were employed to identify signaling changes across distinct sample groups. After selecting potential signaling molecules that were perturbed or induced in MG patients, the official workflow in CellChat was followed with default parameters. Putative interaction pairs between crucial cell clusters in each group were then identified, and the results were visualized as circle plots.

Statistical analysis

The marker genes for transcriptional subpopulations were determined by the FindAllMarkers function in Seurat with the Wilcoxon rank-sum test. One-way analysis of variance followed by the Newman–Keuls post hoc test was employed to analyze the differences among multiple groups. The differences between paired samples were analyzed by the Wilcoxon signed-rank test, while the nonparametric Spearman test was used for correlation analysis. Statistical analyses and graphic processing were carried out by R. Unless otherwise stated, data are considered statistically significant in terms of *P < .05, **P < .01, and ***P < .001.

Composition and classification of immune cells in thymus tissue and peripheral blood of MG patients

To explore the immune pathogenesis and the immune microenvironment of MG, we performed scRNA-seq to study the transcriptomic profiles of thymus tissue and PBMCs (MG-T and MG-B) in two patients with MG (Table 1) as well as PBMCs from two gender- and age-matched healthy controls (HC-B). After quality control and cell filtering, 23 846 and 9701 cells, respectively, originated from MG-T and MG-B, and 6930 cells from HC-B were used for analysis. According to the classical genetic markers, all immune cells were divided into four categories: T cells, B cells, myeloid cells, and NK cells (Fig. 1a and c), of which T cells accounted for the vast majority. Then, we compared the proportion of four immune cell groups in different samples. Because of the difficulty of obtaining absolutely normal thymus tissue intraoperatively, we cited data from other studies (a male donor aged 35–40) in the literature on thymus tissue from healthy controls (HC-T) as a control group (21). Obviously, each sample has a different frequency of cell distribution (Fig. 1b). The proportion of T cells in MG-T group (98.24%) was much higher than that in other groups including HC-T group (58.23%). Compared with HC-B group, the proportions of T cells and B cells were higher (64.93% vs. 55.90%, 11.62% vs. 4.60%) and the proportions of myeloid cells and NK cells were lower than those in MG-B group (8.61% vs. 12.82%, 14.85% vs. 26.68%). It can be seen that most of the DEGs were involved in T cell activation and differentiation signal pathways in both MG-T vs. MG-B group and MG-B vs. HC-B group, providing further evidence that the role of T cells in MG is primarily mediated by the activation of their functions, rather than merely by their numerical dominance in lesioned tissues (Fig. 1d). We further investigated cell-specific expression patterns of these MG-associated genes and found that the human leukocyte antigen (HLA) genes demonstrated the highest expression levels in MG-T cells, consistent with previous studies (22). Moreover, the predominant expression of HLA-DQB1, HLA-DRA, HLA-DRB1, and HLA-DQA1 in B cells and myeloid cells suggests that these cellular populations may form the cellular basis for antigen presentation and immune activation. Genes with the function of regulating T cell activation were most highly expressed in T cells (e.g. CTLA4 gene), while the expression levels of genes such as AP2B1, HSPA8, and AP1B1 were apparently different in the MG-T, MG-B, and HC-B groups (Fig. 1e).

Characterization of T cells derived from scRNA-seq data

We regrouped the T cells from MG-T, MG-B, and HC-B into four populations, including CD4⁺ T cells, CD8⁺ T cells, double-negative T cells (DNT), and double-positive T cells (DPT). On this basis, we focus on the analysis of CD4⁺ T cells and CD8⁺ T cells, in which CD4⁺ T cells can be further divided into naïve, effector CD4⁺ T/Tfh, Treg, and Tfh cells, CD8⁺ T cells can be divided into naïve, cytotoxic, and exhausted CD8⁺ T cells (Fig. 2a).

The proportions of T cell subsets in MG-B and HC-B were very similar, but the proportions of T cell subsets in MG-T were apparently different from those of MG-B and HC-T groups. In particular, the proportions of naïve CD8⁺ T cells (40.62% vs. 18.19%, 40.62% vs. 0.12%) and DPT (41.14% vs. 1.17%, 41.14% vs. 0.11%) in MG-T were much higher than those in MG-B and HC-T groups, indicating that there were a large number of immature T cells in the thymus of patients with MG (Fig. 2b).

Figure 2c demonstrated that most of the marker genes were highly expressed only in their corresponding cluster. For instance, FOXP3 was a recognized marker gene in Treg cells, which was only highly expressed in CD4⁺ Treg cells, but had low expression in other clusters, which further proves the accuracy of cell clustering.

Then, we used the Monocle3 method to construct the differentiation trajectory of CD4⁺ T cells and CD8⁺ T cells. With the development of the differentiation process, CD4⁺CD8⁻ T cells gradually developed from a naïve state to effector CD4⁺ T/Tfh, Treg, and Tfh states. However, almost all the CD4⁺CD8⁻ T cells in MG-T remain in a naïve state, resulting in disruption of immune cell activation and differentiation processes (Fig. 2d, Supplementary Figure S1).

In CD4⁻CD8⁺ T cells, trajectory analysis revealed that naïve CD8⁺ T cells differentiated into cytotoxic CD8⁺ T cells in MG-B, while in MG-T, a small number of naïve CD8⁺ T cells differentiated into exhausted CD8⁺ T cells. This may contribute to the imbalance of thymic immune function in MG patients. In addition, the KEGG enrichment analysis of up-regulated and down-regulated genes of MG-B vs. HC-B showed CD4⁺ was highly correlated with NOD-like receptor signaling pathway and primary immunodeficiency, and CD8⁺ was highly correlated with NOD-like receptor signaling pathway and Fc gamma R-mediated phagocytosis, indicating that candidate signals are significantly enriched within MG peripheral blood T cells that may be involved in disease pathogenesis (Fig. 2e and f).

Characterization of B cells derived from scRNA-seq data

Similarly, we subdivided B cells into nine subgroups and named them according to the reported marker genes. These included four memory B cell clusters (0, 1, 2, 3), three naïve B cell clusters (4, 5, 6), and two plasma cell clusters (7, 8) (Fig. 3a and c).

Compared with HC-T, the proportion of memory B cells in MG-T increased (74.51% vs. 57.33%), while the number of naïve B cells decreased (14.52% vs. 23.55%). Similarly, compared with HC-B, the proportion of memory B cells in MG-B was observed to increase (51.46% vs. 39.57%), while the number of naïve B cells was observed to decrease (44.81% vs. 57.47%) (Fig. 3b). Meanwhile, compared with HC-B, we identified 460 up-regulated and 399 down-regulated genes of B cells in MG-B (Fig. 3d). Furthermore, RNA velocity analysis of B cells across various samples showed that the arrows in cluster 8 of MG-T and MG-B were longer than those of HC-B, indicating a higher presence of unspliced RNAs in this cell cluster of MG patients. It revealed that B cells in MG patients exhibited a more active state of differentiation (Fig. 3e).

Myeloid cell clustering analysis of MG patients and HCs

Next, myeloid cells were sub-grouped into seven clusters, including three classical monocytes (clusters 0, 1, and 3), nonclassical monocytes (clusters 2), intermediate monocytes (clusters 5), conventional dendritic cells (cDCs, clusters 4), and plasmacytoid dendritic cells (pDCs, clusters 6) (Fig. 4a). And, myeloid cells were apparently fewer in MG-T. The proportion of myeloid cell subtypes showed considerable heterogeneity (Fig. 4b). Compared to HC-B, MG-B showed an apparent decrease in classical monocytes (60.81% vs. 67.61%) and an apparent increase in non-classical monocytes (22.24% vs. 14.89%). DCs were predominant in HC-T (nearly 80%), whereas monocytes and DCs each accounted for 50% of the total in MG-T.

Figure 4c demonstrated that the marker genes of monocytes and DCs were specifically highly expressed in the corresponding clusters. Compared with HC-B, we identified 854 up-regulated and 640 down-regulated genes of myeloid cells in MG-B. Moreover, the dys-regulated genes were mainly involved in the secretion of cytokines (Fig. 4d). The dys-regulated genes of MG-T vs. MG-B were mainly involved in antigen presentation and monocyte differentiation (Fig. 4e, Supplementary Figure S2). Furthermore, we analyzed the expression level of the MG-related genes in myeloid cell subpopulations and found that there was also cellular subpopulation heterogeneity (Fig. 4f). Genes such as HLA-DQA1, HLA-DQB1, HLA-DRA, and HLA-DRB1, whose expression was appreciably up-regulated in MG-T, were predominantly expressed in intermediate monocytes, and genes such as AP1S3 and AP1G1, which were appreciably over-expressed in the patients compared to controls, were predominantly distributed in pDCs.

Cell-to-cell communication network analysis among samples

In addition to intracellular information, scRNA-seq can also explore interactions between cells by integrating ligand and receptor information. We analyzed the possible interactions between the major immune cells in the MG patient group and the control group and found that monocytes interacted most strongly with other cells in the PBMCs of MG patients (Fig. 5a). Whereas in HC-B samples, plasma cells exhibited the strongest interactions with other cells, but this predominant effect was no longer observed in MG-B (Fig. 5b). Compared to MG-B, COMPLEMENT, XCR, and IFN-II signaling pathways were more enriched within MG-T (Fig. 5c). While CXCL, GAS, and CD30 signaling pathways were more enriched within MG-B compared with HC-B (Fig. 5d). Finally, we focused on analyzing the role of MG-related signaling pathways CD40 and the interleukin-16 (IL-16) pathway among various cell subpopulations. Compared with HC-B, Effector CD4⁺ T/Tfh cells and Treg cells interacted more strongly with NK cells and monocytes through CD40 signaling pathway, and this interaction of Treg cells was also reflected in MG-T. In MG-B and HC-B, the IL-16 signaling pathway is ubiquitous across cell populations, while in MG-T, this pathway communicates more strongly between T cells and B cells (Fig. 5e).