Single-Cell Transcriptomics Identifies a Prominent Role for the MIF-CD74 Axis in Myasthenia Gravis Thymus

We recruited 12 immunotherapy-naive and corticosteroid-naïve patients with AChR-Ab⁺ EOMG from the National Hospital for Nervous Diseases, London, United Kingdom, and Oxford University Hospitals, United Kingdom (eTable 1), for scRNA-seq analysis; 32 MG serum samples for MIF ELISA analysis (eTable 2) and 3 others from the Charité University Hospital Berlin for multiplex histology analysis of thymic tissue (eTable 3); and 10 matched healthy controls from the Münster University Hospital for MIF ELISA analysis.

Thymic samples were dispersed mechanically within 1.5 hours of removal.^5-7,16 The resulting cell suspensions were then washed and frozen within further 6 hours.⁵ In brief, thymic cells were frozen in 5% dimethyl sulfoxide/95% fetal calf serum at 100–250 × 10⁶ cells in vials held in methanol that was cooled at 5°C per min and subsequently stored over liquid nitrogen. On thawing, they were diluted by doubling the volume every 45 seconds with Roswell Park Memorial Institute medium (RPMI) plus staphylococcal nuclease (to prevent trapping of viable cells in clumps of congealed DNA) and immediately washed. All samples were carefully thawed and processed in parallel. As expected, we recovered approximately 40% of the input cells.

Single thymic cell suspensions were stained for 20 minutes at 4°C. After washing, all samples were analyzed using a CytoFLEX flow cytometer (Beckman Coulter, Brea, CA) and FlowJoTM Software V10.8.1; the gating strategy is shown in eFigure 1. The antibodies used were purchased from Biolegend, CD4 (SK3; #344606), CD8 (SK1; #344722), and CD45 (2D1; #368508), and from Becton, Dickinson and Company (BD) Bioscience, CD3 (UCHT1; #557943) and CD19 (SJ25C1; #562653).

We enriched CD45⁺ cells using CD45 (tumor-infiltrating lymphocytes) microbeads (Miltenyi Biotec., 130-118-780) according to the manufacturer's protocol. We checked cell numbers and viabilities before and afterward to check for significant cell death.

CD45⁺ magnetic-activated cell sorting-enriched single-cell suspensions were loaded onto the Chromium Single Cell Controller using the Chromium Next Gel Bead-in-Emulsion (GEM) Single Cell 5′ Library and Gel Bead Kit v1.1 chemistry (10x Genomics). The derived barcoded complementary DNA was partly used in Chromium Single Cell V(D)J Enrichment Kits (10x Genomics) for further scBCR-seq. We processed samples and libraries according to the manufacturer's instructions using solid phase reversible immobilization select beads (Beckman Coulter). Sequencing was performed on a local Illumina Nextseq 2000 using the P3 Reagents 100-Cycle Kit with a 26-8-0-91 read setup. One of the samples was then excluded because of the low quality of the data recovered (MG07, eTable 1). We processed sequenced data with the cellranger pipeline v.6.1.1 (10X Genomics) according to the manufacturer's instructions. CellBender software was applied to remove the counts because of ambient RNA molecules and random barcode swapping from the count matrices.

Downstream analysis was performed with the R package Seurat v.4.2.0¹⁷ using R v.4.2.1. Only genes present in more than 3 cells and cells with more than 500 genes were retained. Low-quality cells and cell doublets were removed by filtering based on gene count (<500 or >3,000) and high mitochondrial percentages (>10%). After quality control, a total of 33,613 cells remained for further analyses including data normalization, scaling, and principal component analysis (PCA) where cell-cycle gene sets defined in a previous study¹⁸ were removed from the list of variable genes to exclude cell cycle–associated variation and clustering. We removed 2 clusters from the data set that expressed mainly ribosomal genes and most likely represented low-quality cells (11.68% of the cells). Cluster cell identity was assigned by manual annotation using known marker genes and by applying anchor-based label transfer, implemented in Seurat, using as reference a thymi dataset.¹⁸ To achieve a high-resolution annotation in the B lineage, subclustering was performed and batch effects were removed using Harmony.¹⁹ Because B-cell types are so rare in the healthy thymus, the annotation was performed using the anchor-based label transfer Azimuth taking as a reference the tonsil data set.²⁰ Finally, pathway enrichment was conducted using ReactomeGSA.²¹

Spliced/unspliced counts were generated with velocyto v0.17.²² We used the run10x pipeline with the cellranger output, the cellranger gene annotation file, and the human repeat masker file from the University of California, Santa Cruz (UCSC) Genome Browser. The resulting loom files were imported into Seurat with the ReadVelocity wrapper. The spliced/unspliced assays were added to the existing Seurat object. We extracted B-lineage cell clusters from the data set and re-ran the Seurat pipeline. Further downstream analysis was performed with scVelo v.0.25²³ according to the official tutorial. The data were preprocessed by normalization, log transformation (top 2,000 genes), PCA, neighborhood graph, and moments estimation. We then used the dynamical model from scVelo to estimate the RNA velocities. Pseudotime analysis was performed using slingshot v2.6 with default parameters, with the naïve B-cell cluster defined as the starting cluster.

Single-cell B cell receptor (BCR) data were analyzed using scRepertoire following the official vignette.²⁴ The Ig heavy and light chain sequences were combined, and those cells missing 1 or with more than two of the immune receptor chains were excluded. scBCR data were merged with the scRNA-seq object containing the B-lineage cell clusters. BCR clonotypes were classified using the combinations of variable–diversity–joining (VDJ) genes and their CDR3 nucleotide sequences. Relative antibody subclass frequency analyses were performed with Alakazam.²⁵

We analyzed cellular communications with CellChat.²⁶ Processed and clustered scRNA-seq data from thymic samples were used as an input. We used a truncated mean of 0.1; cell-cell communications expressed by less than 10 cells were removed. The resulting ligand-receptor pairs are based on the CellChat repository.

Characteristics of the 3 donors are presented in eTable 3 (ethical approval EA1/391/16). Their thymic structural pathology was checked by a board-certified neuropathologist, excluding infarcts, abscesses, and hemorrhages. Tissue preparation and multiplex microscopy image acquisition were performed as previously published.^27,28 Thymus samples were subjected to multiplexed immunofluorescence histology using the antibodies listed in eTables 4 and 5, stained in the indicated order. Images contain 2048 × 2048 pixels and are generated using an inverted wide-field fluorescence microscope with a ×20 objective, a lateral resolution of 325 nm, and an axial resolution above 5 µm. The fluorescence images were subsequently analyzed using Fiji. All used antibodies were validated in human tonsil or lymphnode tissue control tissue.

Plasma macrophage migration inhibitory factor (MIF) levels (eTable 1) were determined using the human MIF Quantikine ELISA kit (DMF00B, R&D systems) according to the manufacturer's instructions. Serum MIF levels (eTable 2) were assayed with ELISA kits (ProteinTech Group) according to the manufacturer's instructions. The concentration was determined according to the standard curve.

We checked for correlations between different cell type frequencies and between severity scores and MIF expression levels in serum by calculating the Spearman correlation coefficient (R). The Mann-Whitney test was performed on the plasma MIF ELISA results using GraphPad Prism 10. A p value < 0.05 was considered statistically significant.

The study was conducted according to the Declaration of Helsinki and approved by the Ethics Committees of the Universities of Münster (registration nos. 2010-262-f-S, 2011-665-f-S, 2013-350-f-S, 2014-068-f-S, and 2016-053-f-S) and Charité-Universitätsmedizin Berlin (EA1/281/10), London, and Oxford; written informed consents were given by all participants.

Single-cell RNA-seq data have been deposited in the Gene Expression Omnibus repository under reference number GSE233180 and are publicly available as of the date of publication. This article does not report original code. Any additional information is available from the lead contact on request.

To characterize the immune cellular landscape in EOMG thymi, we first performed scRNA-seq analysis in cryopreserved thymic cell suspensions derived from 11 immunotherapy-naïve patients with AChR-Ab⁺ EOMG with a representative range of demographics (eTable 1). After enrichment for CD45⁺ hematopoietic cells, we profiled 29,688 individual cells and annotated 15 distinct populations (Figure 1A), defined by their expression of cell type–specific transcripts (Figure 1B) and alignment with a reference data set.¹⁸ Each population was represented in each patient, with broadly overlapping distributions and no obvious outliers (Figure 1C).

As expected, T cells were most abundant in each thymic sample. Substantial numbers of B-lineage cells were found in all patients except MG01 (Figure 2A), constituting >5% of the total TFH hematopoietic cells (Figure 2, A–D). This is around 10-fold higher than the previously reported values (0.1%–0.5%) for B cells in non-EOMG thymi.^29,30

The EOMG thymic leukocytes contained 8 different T-lineage cell clusters (numbered in Figure 1A): (1) double CD4⁻CD8⁻ T-cell precursors, (2) double-positive T cells, (3) single CD4⁺ T cells, (4) CD4⁺ T memory cells, (5) T regulatory cells, (6) CD8⁺ T cells, (7) CD8⁺ T memory cells, and (8) CD8αα⁺. Non–T-cell and non–B-cell populations present in EOMG thymi included natural killer cells, type 3 innate lymphoid cells, and dendritic cells (DCs) including conventional DC1 and DC2 subsets and plasmacytoid DCs.

The lymphocyte frequencies from scRNA-seq correlated strongly with those determined by flow cytometry (eFigure 2, A–C). Frequencies of B cells determined by scRNA-seq among lymphocytes correlated positively with plasma anti-AChR titers (rs = 0.66, p = 0.027; eFigure 2D).

We identified 6 distinct subpopulations of B cells as defined by their transcript combinations (Figure 3, A and C), and comparison was made with a reference data set including tonsils as a control for the lymph node–like infiltrates.²⁰ Each subpopulation was observed in each patient (Figure 3B); the most abundant were (1) memory cells (CD44, KLF2, TNFRSF13B), followed by (2) naïve B cells (FCER2, FCMR, SELL), (3) plasmablasts/plasma cells ([PCs]: JCHAIN, PRDM1, XBP1, MZB1), (4) GC B cells (CD27, CD38, BCL6), and (5) a cycling population expressing cell cycle and proliferation genes (absence of BCL6 and expression of HMGB2, TUBA1B, MKI67, UBE2C). Cycling B cells are separate in their own cluster at higher dimensions because their profile is dominated by these general cycling-specific genes, but they cluster with GC B cells at lower dimensional analysis, likely because of shared proliferation-related gene expression. We also identified (6) a distinct cluster of activated B cells (CD69, MYC), which expressed features of naïve (FCER2), memory (CD44 and TNFRSF13B), and pre-GC stages (MIR155HG) (Figure 3C).

To profile the proliferative capacity of thymic B-lineage cells, we performed a comparative reactome analysis of pathways indicative of cell proliferation and apoptosis (Figure 3D). Apoptosis-related transcripts were enriched in GC, memory, and activated B cells relative to other B-lineage subsets. As expected, cell cycle–related transcripts were much more prominent in the cycling and GC B cells than in the naïve, memory, and activated B cells and, especially, the PCs (Figure 3D).

To assess directionality and developmental relationships within the B-cell lineage, we performed RNA velocity and slingshot pseudotime analyses (Figure 4, A–C) and identified 2 developmental trajectories: (1) one starting from naïve going through memory and activated B cells to GC cells (Figure 4B); (2) an additional trajectory going from naïve through activated and cycling B cells and culminating in PCs (Figure 4C). The activated cluster, expressing features from naïve, memory, and pre-GC cells, occupied an intermediate position in both trajectories (Figure 4, B and C). GC cells and PCs have diverged significantly, as they acquire unique traits such as activation-induced cytidine deaminase (AID)/somatic hypermutation and terminal differentiation. The low number of cells in their clusters suggests that early intermediate stages may have been missed. Thus, our analyses highlight distinct transcriptional signatures but may underrepresent intermediate states, which could be less apparent due to uniform manifold approximation and projection (UMAP)'s data organization.

We next used single-cell B-cell receptor (scBCR)-seq combined with the scRNA-seq to assess the BCR repertoire of EOMG thymic B-lineage cells and identified 4,536 cells with matching BCR information (eFigure 3). Clonally expanded sequences were most enriched in the GC and cycling B-cell subsets (Figure 4D, eFigure 4A). GCs are critical for the formation of long-lived plasma cells and memory B cells. Within these structures, mature B cells undergo somatic hypermutation (SHM), Ig class switch DNA recombination (CSR), and clonal selection to optimize antigen recognition. Activation-induced cytidine deaminase, encoded by the human activation-induced cytidine deaminase gene (AICDA), is critical for SHM, CSR, and the maturation of Ab responses to foreign and self-antigens. Hence, AICDA-related transcripts were abundantly expressed and enriched in the thymic GC cluster (eFigure 4B). CSR-related proteins and transcripts were predominantly observed in GC and cycling B cells, which exhibited a pre-GC phenotype (eFigure 4B). Analysis of scVDJ-derived antibody class and subclass frequencies across the B lineage identified IgG1 as the most frequently observed immunoglobulin produced by class-switched EOMG thymic B-lineage cells, followed by IgA1, IgG2, and IgG3 (eFigure 4C). IgA⁺ PCs were also prominent in a previous study on EOMG thymi.⁷

To profile cell-cell communications in EOMG thymi, we integrated single-cell data with a curated ligand-receptor pair database through a bioinformatics application, CellChat (Figure 5A). We included mature T cells and the B-lineage cells in these network analyses. They showed strong interactions between B-lineage cells and T-cell subsets, notably including CD8⁺ T cells (Figure 5, B and C). To better define which signals might support or sustain B-lineage cells, we profiled individual ligand↔receptor interactions from T-cell to B-cell subsets and identified 21 significant pairs, all of them functionally associated with immune cell migration and signaling (Figure 5C).

The strongest interaction was predicted for macrophage migration inhibitory factor (MIF; Figure 5C), so we next focused on interaction strengths between different cell types. Surprisingly, in these EOMG thymic samples, we observed that CD8⁺ and regulatory T cells displayed the highest levels of outgoing interactions within the MIF pathway, interacting with its receptors and co-receptors CD74, CXCR4, and CD44 on B-lineage cells, in which the interactions were mainly incoming (Figure 5D). However, GC cells and PCs are strongly human leukocyte antigen (HLA)-class II positive, supposed to express CD74, and so are capable of delivering signals by presenting antigens; malignant plasma cells can also express CD74.³¹ Thus, it would make biological sense for GC B cells and plasma B cells to have low incoming interaction strengths when their role is primarily to present antigen (GC B cells) or secrete antibodies (PCs, which have very few surface receptors and are known to receive minimal incoming interactions). This suggests a central role for these cell types in the MIF pathway. CD74, also referred to as invariant chain, is a type II transmembrane protein expressed on antigen-presenting cells, including the thymus,¹⁸ and also B cells, in which it mediates assembly and trafficking of peptide:major histocompatibility complex class II complexes. CD74 further functions as a survival receptor on B cells, essential for transcriptional regulation of genes that control their proliferation, differentiation, and survival.²⁰ To better assess spatial CD74 expression at the protein level, we used multiepitope ligand cartography in 3 EOMG thymic tissue samples, which confirmed the expected strong and abundant CD74 and CD44 expression on thymic B cells. The CD74⁺ B-cell accumulations were located mainly in the perivascular space but also in the medullary areas. T cells (major source of MIF on gene expression level) were in proximity with B cells (Figure 6, eFigures 5–7). In addition to B cells, there were some cells with myeloid morphology and strong CD74 expression in frequencies as expected from previous studies.³²

MIF is a proinflammatory cytokine mediating activation in addition to survival of B cells through CD74 ligation. Its expression and secretion are also elevated in most solid and hematogenous cancers, and serum levels of soluble MIF are substantially elevated in patients with B cell–driven autoimmune diseases such as systemic lupus erythematosus.³³

Notably, serum MIF levels reportedly correlate with higher Myasthenia Gravis Foundation of America (MGFA) status and disease severity, as also assessed using the quantitative MG disease score (QMG), in patients with generalized MG.³⁴ Because modern assessments of disease severity are not available for the exploratory cohort, we quantified serum MIF levels in an independent cohort of patients with AChR-Ab–positive EOMG. They were all receiving immunotherapies at the time of sampling (eTable 2). In this study, serum MIF levels correlated with higher MGFA status (Figure 7A) although not with QMG or activities of daily living scores (Figure 7, B and C) assessed at sampling.