Abstract

   Myasthenia gravis (MG) is a neurological disease caused by
   autoantibodies against neuromuscular-associated proteins. While MG
   frequently develops in thymoma patients, the etiologic factors for MG
   are not well understood. Here, by constructing a comprehensive atlas of
   thymoma using bulk and single-cell RNA-sequencing, we identify ectopic
   expression of neuromuscular molecules in MG-type thymoma. These
   molecules are found within a distinct subpopulation of medullary thymic
   epithelial cells (mTECs), which we name neuromuscular mTECs (nmTECs).
   MG-thymoma also exhibits microenvironments dedicated to autoantibody
   production, including ectopic germinal center formation, T follicular
   helper cell accumulation, and type 2 conventional dendritic cell
   migration. Cell–cell interaction analysis also predicts the interaction
   between nmTECs and T/B cells via CXCL12-CXCR4. The enrichment of nmTECs
   presenting neuromuscular molecules within MG-thymoma is further
   confirmed immunohistochemically and by cellular composition estimation
   from the MG-thymoma transcriptome. Altogether, this study suggests that
   nmTECs have a significant function in MG pathogenesis via ectopic
   expression of neuromuscular molecules.

   Subject terms: Autoimmune diseases, Data integration, Neuroimmunology,
   Thymus
     __________________________________________________________________

   Myasthenia Gravis and thymoma are frequently associated with patients
   suffering from both diseases. Here the authors perform single cell
   sequencing of thymoma and find that there are autoimmune antigens such
   as neuromuscular proteins expressed aberrantly in neuromuscular mTECs
   in patients with both diseases.

Introduction

   Myasthenia gravis (MG) is the most common disorder of neuromuscular
   transmission caused by autoantibodies against the motor endplate, such
   as anti-acetylcholine receptor (AChR) antibodies. MG is often
   accompanied by thymoma, and thymoma-associated MG (TAMG) is more
   difficult to manage than other forms of MG because of its frequent
   crisis, the need for surgery, the difficulty of perioperative
   management and the need for intense immunotherapies^[66]1. As the
   epidemiology, 21% of MG patients experienced thymoma^[67]2, and 25% of
   thymoma patients experienced MG^[68]3, indicating that MG and
   abnormalities in the thymus are closely related to each other. This is
   also exemplified by that thymectomy is a well-established treatment for
   TAMG in addition to immunosuppressive treatments^[69]4.

   Abnormalities of the thymus, in which immature thymocytes differentiate
   into matured CD4^+ or CD8^+ T cells, are frequently associated with a
   variety of autoimmune diseases, such as pure red cell aplasia and Good
   syndrome^[70]5. In the thymus, T-cell maturation and selection are
   conducted by the interaction with antigen-presenting cells, including
   thymic epithelial cells (TECs), myeloid cells, and B cells^[71]6,[72]7.
   The positive selection of functional T cells is mediated by cortical
   TECs (cTECs), while the negative selection of auto-reactive T cells is
   mediated by medullary TECs (mTECs) presenting self-antigens on MHCs. In
   line with the critical role of mTECs, a loss of function mutation of
   AIRE, which is an essential transcription factor for producing
   self-antigens in mTECs^[73]8, causes systemic autoimmunity called
   Autoimmune Polyglandular Syndrome Type 1 (APS-1)^[74]9. In addition,
   dysregulation of the thymus, especially thymoma, is frequently
   associated not only with MG but also neurological disorders, including
   encephalitis, which is caused by a wide range of
   autoantibodies^[75]10–[76]15 Thus, abnormalities of the thymus are
   closely associated with the generation of self-reactive autoantibodies,
   which result in the development of autoimmune diseases.

   Given that MG is caused by self-reactive autoantibodies, MG-specific
   changes within thymoma may be a clue for understanding the pathogenesis
   of MG. It has been reported so far that the accumulation of
   neurofilaments, which is expressed in neurons under the normal
   condition, are highly detected in MG-thymoma^[77]16. In addition,
   germinal centers (GCs) and T follicular helper (T[FH]) cells, both of
   which play critical roles in antibody production, are also enriched in
   MG-thymoma^[78]17,[79]18. Despite the possible contribution of these
   changes in MG pathology, the complete picture of MG pathogenesis from
   abnormalities in thymoma to auto-reactive B cell maturation is still
   poorly understood due to intra- and inter-individual heterogeneity of
   the thymus.

   Here, by integratively analyzing bulk and single-cell transcriptomes of
   MG-thymoma to identify the complex pathogenicity of MG in thymoma, we
   show that a distinct subpopulation of mTECs ectopically express
   neuromuscular-associated molecules and could contribute to the
   pathogenesis of MG via presenting neuromuscular molecules to
   self-reactive immune cells developed in thymoma.

Results

Ectopic expression of neuron-related molecules in MG-thymoma

   To characterize MG-specific changes in thymoma comprehensively, we
   first investigated gene expression profiles from surgically dissected
   thymoma samples enrolled by The Cancer Genome Atlas (TCGA)^[80]19
   (Fig. [81]1a). Of the 116 thymoma samples with RNA-seq data, 34 were
   complicated with MG. In the WHO classifications, which are commonly
   used for the thymoma staging based on histology, there are six
   classifications: Type A, AB, B1-B3, and C (i.e., A, spindle cells; AB,
   mixed spindle cells and lymphocytes; B1, lymphocytes > epithelial
   cells; B2, mixed lymphocytes and epithelial cells; B3, predominant
   epithelial cells; and C, carcinoma). In the present dataset, MG was
   associated with multiple types except for type C (Supplementary
   Fig. [82]4a), with its peak at type B3 (MG complication rate: 54.5%)
   and B2 (53.8%), whereas a previous report observed the peak at type B2
   (71.1%)^[83]20. When we investigated differentially expressed genes
   between thymoma with and without MG (Supplementary Data [84]1,
   Supplementary Fig. [85]4b), 93 and 91 genes were identified as
   upregulated and downregulated genes in MG, respectively. The
   upregulated genes contained neuromuscular-related molecules; NEFM,
   RYR3, GABRA5, and immunoreceptors; PLXNB3, IL13RA. We also observed a
   slight increase in the acetylcholine receptor CHRNA1, which is the main
   target of autoantibodies in TAMG (log[2] fold change = 1.07,
   P[adj] = 0.87 with DESeq2^[86]21; P = 0.0051 with two-sided
   Mann–Whitney U-test, Supplementary Fig. [87]4c).

Fig. 1. Transcriptome profiling of thymoma and MG-specific expression of
neuro-related genes.

   [88]Fig. 1
   [89]Open in a new tab

   a PCA plots for transcription profile of thymomas from 116 patients.
   The left panel shows the disease status, MG or non-MG, and the right
   panel shows WHO classification based on histology. b Gene modules
   defined using WGCNA and the association with MG, WHO classification,
   gender, and age at diagnosis. Numbers in colored boxes on the left are
   the number of genes included in each module. The numbers in the heatmap
   show the correlation (upper) and the P-value (lower). c Eigengenes of
   each module for each patient on the PCA plot. d Heatmap of the gene
   expression of keratins in the yellow and blue modules. The color
   represents the Z-score of normalized expression by DESeq2. WHO
   classification and MG status were shown at the top of the heatmap. e
   Immunohistochemical (IHC) staining of KRT6 and KRT17 in MG and
   non-MG-thymoma. The scale bar: 100 µm. f Protein levels of KRT6 and
   KR17, and KRT6 normalized by KRT17 in MG and non-MG-thymoma quantified
   using microscopic images (MG n = 8, nonMG n = 9, Details are provided
   in Supplementary Fig. [90]2, [91]3 and Methods. Source data are
   provided as a Source Data file.). The signals were analyzed using a
   two-sided Mann–Whitney U-test. g Significantly enriched REACTOME
   pathways in the yellow module. The node size represents the number of
   genes included in each pathway, and the color represents the adjusted
   P-value of the enrichment. The pathways were sorted by the ratio of
   genes included in the yellow module. h Genes in enriched REACTOME
   pathway in the yellow module. Genes with log[2] fold change > 1 in
   comparison of MG and non-MG were selected. i Venn diagram showing
   overlap of targets of autoantibodies associated with thymoma with genes
   in the yellow module and upregulated genes in MG. Data were analyzed
   using a two-sided Fisher’s exact test.

   To dissect transcriptome changes unbiasedly, we next adopted the
   unsupervised gene clustering approach; Weighted Gene Co-expression
   Network Analysis (WGCNA)^[92]22 on the large-scale thymoma samples
   (Supplementary Fig. [93]4d). We initially constructed gene “modules”,
   each of which were composed of a set of genes showing correlated gene
   expression. In the thymoma samples, seven modules consisting of 30
   -1102 genes were obtained and represented by colors (Supplementary
   Data [94]2). We next investigated the association of clinical
   information with a representative gene expression of each module,
   calculated as eigengene (Fig. [95]1b). MG had a most significant
   correlation with the yellow module (ρ = 0.55, P = 6 × 10^−10,
   Supplementary Fig. [96]4e) among modules. Each type of the WHO
   classification corresponded to different modules, respectively; i.e.,
   A, AB—black, turquoise; B1, B2 - blue; B2, B3 - yellow; C—green, red.
   The gray module was strongly associated with gender, and the black and
   turquoise modules with age at diagnosis. On the PCA plot based on the
   transcriptome, the enrichment of each module was well-coordinated with
   the profile of WHO types, suggesting that the heterogeneity of thymoma
   can be represented by the gene modules (Fig. [97]1c). Next, we
   investigated the detailed gene profiles of the modules associated with
   MG. In accordance with that MG was particularly enriched in the WHO
   type B epidemiologically^[98]20 and in TCGA samples (Supplementary
   Fig. [99]4a), the yellow module linked to MG was associated with type
   B2 and B3 (Fig. [100]1c). In contrast, the blue module, which was
   independent of MG, was also associated with types B1 and B2
   (Fig. [101]1c). To distinguish the yellow module from the blue one, we
   selected cytokeratins, which have various isotypes specific to tissues.
   Within cytokeratin isotypes, KRT6A, KRT6C, KRT15 were specific to the
   yellow module, whereas KRT7, KRT17, KRT18 were to the blue module
   (Fig. [102]1d). To confirm the difference histologically, we stained
   KRT6 and KRT17 proteins in thymoma tissue sections and observed the
   corresponding staining patterns to MG and non-MG-thymoma, respectively
   (Fig. [103]1e, f). We next examined the enriched pathways in the yellow
   module. Intriguingly, when we examined the enriched pathways in the
   yellow module, the most significantly enriched pathway was neuronal
   systems which included GABA receptors (GABRA5, GABRB3), neurofilaments
   (NEFM, NEFL), voltage-gated potassium channels (KCNC1, KCNH2, KCNH5,
   KCND3), and an NMDA receptor (GRIN2A) (Fig. [104]1g, h, Supplementary
   Data [105]3). Formation of the cornified envelope including
   cytokeratins and ion channel transport were also enriched in the yellow
   module. On the other hand, the other modules did not show the
   enrichment in neuronal systems but instead showed the enrichment of
   different types of pathways, such as interleukin-4 and interleukin-13
   signaling in the blue module, extracellular matrix organization in the
   green module, and interleukin-10 signaling in the red module
   (Supplementary Fig. [106]4f, Supplementary Data [107]4).

   Thymoma has been shown to associate with paraneoplastic neurological
   diseases, such as encephalitis and myositis, besides myasthenia
   gravis^[108]3,[109]23. We observed the significant overlap of candidate
   target antigens of thymoma-relating autoantibodies (Supplementary
   Data [110]6) with the yellow module genes (odds ratio = 7.87,
   P = 1.65 × 10^−6) and more weakly with the differentially expressed
   genes between MG and non-MG (odds ratio = 7.42, P = 2.67 × 10^−3). The
   overlap included an NMDA receptor (GRIN2), voltage-gated potassium
   channels (KCNH2, KCNC1, KCNA1, KCNC2, KCND3, KCNH5), a glycine receptor
   (GLRA4), a GABA receptor (GABRA5), and a ryanodine receptor (RYR3)
   (Fig. [111]1i).

   We also examined immune profiles of MG, such as T-cell receptor (TCR)/B
   cell receptor (BCR) diversity and viral infections. The diversity of
   immunoglobulins in MG-thymoma was lower than that of non-MG-thymoma,
   suggesting that B cell maturation and expansion were occurred in
   MG-thymoma (Supplementary Fig. [112]5a). The diversity of TCR
   represented by the Gini index was mostly unchanged between them
   (Supplementary Fig. [113]5a), but the composition rate of a TCR alpha
   chain J, TRAJ24, was high in MG (P[adj] = 9.6 × 10^−4; Supplementary
   Fig. [114]5b), and especially the TRAJ24-TRAV13-2 combination was
   7.50-fold more frequent in MG-thymoma (Supplementary Fig. [115]5c). To
   assess the effect of HLA on MG susceptibility, we determined major
   alleles of the HLA class I and II in MG using the same TCGA bulk
   RNA-seq dataset. The strongest association was observed in DQA1*01:04
   (odds ratio = 4.43, P = 0.050), followed by DQB1*05:03 (odds
   ratio = 4.25, p = 0.056), A*24:02 (odds ratio = 2.84, P = 0.058, also
   reported by Machens et al.^[116]24) though all associations were below
   the significance level (Supplementary Fig. [117]5d). MG development has
   also been shown to associate with viral infections, including
   SARS-CoV2^[118]25 and Epstein-Barr virus^[119]26,[120]27. Therefore, we
   next examined the presence of viral transcripts in thymoma using the
   TCGA dataset. Although various viral transcripts such as Epstein-Barr
   virus and herpesvirus 6A were detected in MG-thymoma, no significant
   association with viruses was observed (False Discovery Rate
   (FDR) < 0.1, Supplementary Fig. [121]5e). In addition, we could not
   find any significant somatic mutations associated with MG, whereas
   missense mutations in GTF2I were observed in 49% of thymoma patients as
   previously reported^[122]28 (Supplementary Fig. [123]6). Altogether,
   the unbiased large-scale omics analysis made visible MG-specific
   expression of neuromuscular molecules and the distortion in the
   diversity of TCRs and BCRs.

Single-cell profiling of thymoma and PBMCs from anti-AChR antibody positive
patients

   To clarify the source of neurological molecules and the surrounding
   immune environments in MG-type thymoma, we conducted single-cell RNA
   sequencing (scRNAseq) experiments of thymoma and peripheral blood
   mononuclear cells (PBMCs) derived from four anti-AChR antibody positive
   patients (Fig. [124]2a, Supplementary Fig. [125]7). The patients
   consisted of three females and one male, had not received
   immunosuppressive therapy preoperatively except for one patient, ranged
   in age from 35 to 55 years, and had thymoma type AB-B2 (Supplementary
   Data [126]7). Using a droplet-based single-cell isolation method, we
   profiled 33,839 cells from thymoma and 30,810 cells from PBMCs and
   identified 49 clusters upon them (Fig. [127]2b,c, Supplementary
   Fig. [128]8a, b). To analyze the similarity and differences between
   thymic cells and PBMCs, clustering was performed for the pooled cells.
   The cell annotation of PBMCs and the thymus was well-concordant with
   the previously reported scRNAseq experiments for healthy PBMCs^[129]29
   and the thymus^[130]30 (Supplementary Fig. [131]8d, e), and each
   cluster was well-separated by the specifically expressed genes.
   (Fig. [132]2d, Supplementary Fig. [133]8c, Supplementary Data [134]9,
   [135]10). In the latter parts, we analyzed the detailed expression
   profiles of the major clusters; stromal cells, T cells, and B cells, of
   MG-type thymoma.

Fig. 2. Overview of scRNAseq of thymoma and blood from anti-AChR receptor
antibody positive patients.

   [136]Fig. 2
   [137]Open in a new tab

   a The experimental design of scRNAseq. Immune cells and non-immune
   cells from MG-type thymoma and immune cells from the blood of
   corresponding patients were collected for scRNAseq. b UMAP plot for
   65,935 cells displaying the 49 clusters from thymoma and blood of MG
   patients. c UMAP plot of marker genes, inferred cell cycle, and tissue
   origins. d Dot plot depicting signature genes’ mean expression levels
   and percentage of cells expressing them across clusters. The detailed
   dot plot is shown in Supplementary Fig. [138]8c.

Identification of a unique thymic epithelial cell cluster in MG-type thymoma

   We profiled stromal cells of thymoma firstly. Clustered stromal cells
   corresponded to endothelial cells (positive for PECAM1/CD31, VWF),
   normal fibroblasts (FN1, EGFL6), tumor-associated fibroblasts (TAFs;
   PDGFRA, ADH1B), and thymic epithelial cells (TECs; KRT19, S100A14)
   (Fig. [139]3a, Supplementary Fig. [140]9a). We then extracted the TEC
   cluster and re-clustered them into cTEC (CCL25, PSMB11) and mTEC
   (CCL19, KRT7) clusters (Fig. [141]3b, c). The mTECs further fell into 3
   clusters; mTEC(I) specifically expressing KRT15 and IFI27; mTEC(II)
   expressing CLDN4 and KRT7; and the unique mTECs expressing
   neuromuscular-related molecules (Fig. [142]3d). Cells in this unique
   mTEC cluster also expressed brain-specific genes included in the yellow
   module, such as GABRA5, MAP2, NEFL, NEFM, SOX15, TF. Their ectopic
   expression was also confirmed immunohistochemically in MG-thymoma
   tissue sections (Fig. [143]3e, Supplementary Fig. [144]9c, d). Among
   the yellow module genes, we selected KRT6 and GABRA5 as marker genes of
   nmTECs in the following criteria; (1) the expression was increased in
   MG patients in TCGA bulk RNA-seq dataset, (2) in the scRNAseq dataset,
   stably and preferentially expressed in nmTECs (Supplementary
   Fig. [145]9b). By immunohistochemistry (IHC), GABRA5 as one of the
   neuronal molecules expressed in the unique cluster and the cytokeratin
   KRT6, which belongs to the yellow module, were detected in identical
   cells (odds ratio = 50.6, P < 10^−16), with the cytoplasm and the
   pericellular localization of the cells, respectively (Fig. [146]3g, h,
   Supplementary Fig. [147]10). Due to the atypical expression profile of
   the cluster, we named the population neuromuscular-mTECs or nmTECs.
   nmTECs also expressed some of the targets of autoantibodies in
   thymoma-associated neuromuscular disorders highlighted in TCGA bulk
   RNA-seq analysis in Fig. [148]1i (Supplementary Fig. [149]9e). To
   assess the counterparts in the normal thymus, we compared scRNAseq data
   of thymoma TECs with that of the normal thymus previously
   published^[150]30. Thymoma nmTECs were partially correlated with an
   immature TEC cluster (mcTECs) and not with myoid cells (TEC(myo)) and
   neuroendocrine cells (TEC(neuro)) in the normal thymus (Fig. [151]3f).
   Next, to clarify the biological characteristics of nmTECs, gene set
   enrichment analysis of nmTECs was performed using the REACTOME gene
   sets. nmTECs showed the enhancement of pathways such as TP53 activation
   and pathways in cancer (Fig. [152]3i, k, Supplementary Data [153]11).
   nmTECs showed the highest number of detected reads per cell, while
   markers for other cell-types such as T cells and B cells were not
   detected in their expression (Fig. [154]3m, Supplementary
   Fig. [155]9f), suggesting that nmTECs are tumorous cells and not
   doublets. nmTECs also showed the enrichment in E3 ubiquitin ligases,
   IFNγ signaling, IFNα signaling, and class I MHC mediated antigen
   processing and presentation (Fig. [156]3i, j). Although MHC class II
   antigen presentation was not significantly enriched in the pathway
   analysis, nmTECs showed upregulation of HLA class II molecules and IFNγ
   together with downstream molecules of IFNγ signaling; STAT1, IRF1,
   CIITA, which have been shown to activate MHC class II
   regulations^[157]31,[158]32. It suggests that nmTECs may possess a high
   capability of antigen-presentation via HLA class II (Fig. [159]3n).
   AIRE and FEZF2, which have been reported to be involved in the
   production of self-antigens in mTECs, and tissue-restricted antigens
   (TRAs) were expressed in a few cells of mTEC(I) cells, but not in
   nmTECs (Supplementary Fig. [160]9g). These observations thus suggest
   that nmTECs is a unique population producing neuromuscular-related
   molecules with active antigen presentation via MHC class I and II
   molecules.

Fig. 3. Neuromuscular thymic epithelial cells (nmTECs) expressed
neuromuscular genes, IFN gamma signaling pathway genes, and HLA molecules.

   [161]Fig. 3
   [162]Open in a new tab

   a, b UMAP embedding for stromal clusters (a) and thymic epithelial
   cells (TECs) clusters (b) in thymoma. c Gene expression of marker genes
   on UMAP embedding. d, Violin plots of mean expression of the REACTOME
   gene sets; Neuronal System (left) and Muscle contraction (right) in TEC
   clusters. e Dot plot of the yellow module genes. Corresponding protein
   expressions were also confirmed using IHC (Supplementary Fig. [163]9b).
   f Heatmap showing correlation of transcriptional profile with TEC cells
   in thymoma (this publication) and a normal thymus (Park et
   al.^[164]30). g Immunofluorescence staining for confirming the presence
   of nmTECs positive for GABRA5 (red), KRT6 (green), and DAPI (blue).
   Scale bars: 20 µm. h Cross table showing cell numbers of GABRA5
   positive/negative and KRT6 positive/negative cells in a representative
   IHC slide. Data were analyzed using a two-sided Fisher’s exact test. i
   Volcano plot showing REACTOME gene sets enriched in nmTECs. j–l Violin
   plots of mean expression of the gene sets (j–k) and the number of
   detected reads per cell (l). j Represents the significantly enriched
   REACTOME gene sets and k represents the KEGG gene set, Pathways in
   cancer. m Dot plot of gene expression of HLA class II-related molecules
   in TECs.

Dynamics of myeloid cells in MG-type thymoma

   To explore the MG-specific immune environment in thymoma, we next
   profiled myeloid cells. We identified six myeloid clusters in thymoma
   and PBMCs (Supplementary Fig. [165]11a, b). Monocytes were dominated in
   PBMCs, while macrophages and dendritic cells were populated mostly in
   thymoma (Supplementary Fig. [166]11e). Among clusters, type 2
   conventional dendritic cells or cDC2s (CLEC10A, FCER1A, ITGAX/CD11c),
   which preferentially polarize toward T[H]2, T[H]17, and T[FH]
   responses^[167]33,[168]34, were inferred to migrate from the periphery
   into thymoma from RNA velocity^[169]35 (Supplementary Fig. [170]11c–f).

B cell maturation with ectopic germinal center formation in MG-type thymoma

   Since B cells are the source of the autoantibodies causative for MG, we
   next assessed B cell dynamics in MG-type thymoma. To determine the
   subpopulations of B cells, we categorized them into eight distinct B
   cell clusters. Notably, we found a population forming a germinal center
   (GC; positive for BCL6, MEF2B) in MG-type thymoma (Fig. [171]4a), while
   GC B cells were not detected in the normal thymus (Supplementary
   Fig. [172]12d). The formation of ectopic germinal centers in MG-thymoma
   was also histologically confirmed (Fig. [173]4f). Based on the
   expression of immunoglobulins, B cells were divided into three groups;
   (1) Naive, GC, pre-GC (IGHM, IGHD, IGHG3 high);(2) memory B cells
   (IGHA1, IGHA2, IGHG2 high); and (3) plasmablasts (IGHG1, IGHG3, IGHG4
   high) (Fig. [174]4b). We also observed that a pre-GC B cell population
   (STMN1, TCL1A) was preferentially enriched in thymoma (Fig. [175]4c).
   The RNA velocity analysis showed that pre-GC cells were directed from
   naive B cells toward GC B cells, memory B cells, and plasmablasts in
   MG-type thymoma (Supplementary Fig. [176]12a), suggesting that the B
   cell maturation progresses normally in the MG-type thymoma. In
   addition, Pre-GC, GC, thymic memory B cells, and plasmablasts were
   enriched in thymoma compared to PBMCs (Fig. [177]4d, e).

Fig. 4. Immune cell landscape elucidates GC formation, T[H]0-T[FH]
enhancement, and Treg recirculation in MG-type thymoma.

   [178]Fig. 4
   [179]Open in a new tab

   a UMAP embedding for B cell clusters of thymoma and peripheral blood. b
   Heatmap of immunoglobulin expressions in each B cell cluster. Mean
   expressions in each group are shown as a heatmap. c Dot plot of gene
   expression of marker genes of each B cell cluster. d Density plots
   showing B cell accumulation in the periphery (left) and thymus (right).
   e Cell proportion of each B cell cluster in thymoma and peripheral
   blood. Periphery n = 2, Thymoma n = 4. Error bars show 98% highest
   density interval. *FDR < 0.05. (Methods) f Representative IHC image of
   germinal center in MG-thymoma stained for CD79A. Scale bar: 100 µm. The
   presence of GC was also checked in Fig. [180]6j and Source Data. g UMAP
   embedding for CD4^+ T-cell clusters of thymoma and peripheral blood. h
   Dot plot of gene expression of marker genes of each T-cell cluster. i
   Density plots showing T-cell accumulation in the periphery (left) and
   thymus (right). j Cell proportion of each T-cell cluster in thymoma and
   peripheral blood. Periphery n = 2, Thymoma n = 4. Error bars show 98%
   highest density interval. Center points represent the mean of the
   posterior distributions. *FDR < 0.05. (Methods) k Bar plot of thymus
   specific genes across CD4^+ T-cell clusters ranked by the number of
   cell types where each gene was upregulated (P[adj] < 0.05 and log[2]
   fold change > 1) in mature CD4^+ T cells. l TCR similarity between
   peripheral blood and thymoma. The thicknesses of edges represents TCR
   similarity. *FDR < 0.05 in e and j. Statical procedures in e and j are
   described in Methods. GC germinal center.

T-cell polarization in MG-type thymoma

   In the thymus, T cells are characteristically educated by
   antigen-presenting cells, including mTECs, and abnormalities of the
   antigen-presentation frequently associate with autoimmune diseases via
   T-cell dysfunction^[181]36,[182]37. Therefore, we next investigated
   whether T cells in thymoma were engaged in MG pathogenesis
   (Supplementary Fig. [183]12d). We observed the existence of immature
   and mature T cells in thymoma, suggesting that the physiological T-cell
   development was maintained even in MG-type thymoma. Among populations,
   we identified a thymoma-specific mature T-cell population, CD8^+
   tissue-resident memory T cell (CD8 T[RM]) expressing CXCR6 as seen in
   other tissues such as the lung^[184]38,[185]39 and skin^[186]40
   (Supplementary Fig. [187]12d). We next focused on CD4^+ T-cell
   clusters, which are essential for B cell activation. In thymoma and
   PBMCs, we identified 13 specific clusters, which corresponded to cells
   in the process of differentiation, i.e., immature thymic CD4^+ T cells
   to terminally differentiated effector memory CD4^+ T cells (CD4
   T[EMRA]). T cells after thymic selection contained CD4^+ naive T cells
   (CD4 T[NAIVE]; CCR7^+ FAS^-), CD4^+ central memory T cells (CD4 T[CM];
   CCR7^+ FAS^+), effector memory T cells (CD4 T[EM]; CCR7^- FAS^+), and
   terminally differentiated effector memory CD4^+ T cells (CD4 T[EMRA];
   FAS^+ CD28^-) (Fig. [188]4g, h). We also identified T-cell
   polarizations using characteristic transcription factors and chemokine
   receptors such as T[H]1 (TBX21/Tbet) in T[EM] and T[EMRA], T[H]2
   (GATA3, CCR4) in T[CM], T[H]17 (RORC, CCR6) in T[CM], T follicular
   helper cells (T[FH]; CXCR5, PDCD1) in T[CM] (Fig. [189]4h). These cell
   annotations were also concordant with the bulk RNA-seq dataset of
   purified T cells^[190]41 (Supplementary Fig. [191]12e). When we
   assessed the tissue localization of these cells, CD4 T[CM] (T[H]0) was
   more abundant in the thymus, and CD4 T[CM] (T[FH]) was equally abundant
   in the thymus, whereas other memory T cells such as CD4 T[CM] (T[H]2),
   CD4 T[CM] (T[H]17) were significantly more abundant in the periphery,
   suggesting that T[H]0-T[FH] axis are prominent in the thymus
   (Fig. [192]4i, j). Next, to infer T-cell dynamics between thymoma and
   periphery, we investigated the commonalities of T-cell receptor (TCR)
   repertoires of these cell populations. Strong clonal expansions were
   observed in T[H]1 prone clusters, including CD4 T[EM] (T[H]1/17), CD4
   T[EM] (T[H]1) and T[EMRA] (T[H]1), and also slight clonal expansions in
   CD4 T[CM] (T[H]2), CD4 T[CM] (T[H]17), and activated T[reg] cells
   (Supplementary Fig. [193]12f). By examining TCR similarity between the
   thymus and periphery for each cluster, T[reg] cells showed higher
   levels of TCR similarity between the thymus and the periphery, compared
   to the other cell populations (Fig. [194]4l). This suggests that naive
   T[reg] cells are activated in thymoma aberrantly and circulated into
   the periphery. In addition, a chemokine receptor CXCR4 was
   preferentially expressed in thymic mature T cells (Fig. [195]4k,
   Supplementary Fig. [196]12g). A couple of thymic T cell-specific genes
   such as CD69, SOCS1 (STAT-Induced STAT Inhibitor 1) and RGS1 (Regulator
   Of G-Protein Signaling 1) were also expressed in B cells in thymoma,
   but not in periphery (Supplementary Fig. [197]12b, c), suggesting that
   these genes were regulated by a shared tissue-specific program between
   T and B cells. Taken together, a detailed analysis of B and T cells
   revealed that MG-type thymoma kept primary lymphoid tissue
   characteristics for T-cell education and gained abnormal inflammatory
   profiles with ectopic GC formations.

Cell–cell interaction inference

   To analyze the communications among cells, we next inferred cell–cell
   interaction by integrating single-cell data with a curated
   ligand-receptor pair database through a bioinformatics application,
   CellPhoneDB^[198]42 (Fig. [199]5a). The cell fraction possessing the
   highest number of intercellular interactions was nmTECs (Fig. [200]5b).
   The cell–cell interaction network analysis showed that nmTECs acted as
   a hub in the network and interacted with myeloid cells, T cells, B
   cells, tumor-associated fibroblasts, and endothelial cells
   (Fig. [201]5c). nmTECs and tumor-associated fibroblasts preferentially
   expressed CXCL12, and thymic B cells and helper T cells, including
   T[FH]and T[reg] cells, expressed its receptor, CXCR4. Given that the
   CXCR4-CXCL12 axis has been shown to be important in T-cell homing in
   synovial tissues of rheumatoid arthritis^[202]43, neurogenesis^[203]44,
   and maintenance of hematopoietic stem cells^[204]45, the interaction
   may be important for nmTEC-mediated T-cell regulation (Fig. [205]5d).
   CXCR5 was expressed in B cells, T[FH], and CD8 T[RM], while its ligand,
   CXCL13, was expressed in T[FH] and CD8 T[RM], suggesting that the
   putative role of CXCR5-CXCL13 for T-B interaction in thymoma. This
   predicted interaction was consistent with the previous
   findings^[206]46. We also predicted the interactions of nmTECs with
   vascular endothelial cells via VEGFA and VEGFE and with
   tumor-associated fibroblasts via PDGFA-PDGFRA. To verify the predicted
   interaction, we performed immunostaining of CD31 on MG-type thymoma
   sections and observed that GABRA5^+ nmTECs were in proximity to CD31^+
   vascular endothelial cells (Fig. [207]5e–g, Supplementary
   Fig. [208]13). These observations suggest that nmTECs may promote
   angiogenesis via the interaction of vascular endothelial cells, which
   is consistent with previous reports^[209]17.

Fig. 5. nmTECs strongly associated with epithelial cells, myeloid cells, and
T cells with characteristic ligand-receptor pairs.

   [210]Fig. 5
   [211]Open in a new tab

   a Schematic view of the cell–cell interaction analysis. b Bar chart
   showing the number of significant interactions with other cell types in
   each cell type. c Cell–cell interaction network inferred from scRNAseq
   data. Each node represents a cell type, and the thickness of each edge
   represents the number of significant interactions. Edges with <75
   significant interactions were removed. d Dot plot of gene expression of
   ligand-receptor pairs involved in trafficking, immunomodulation, and
   angiogenesis in CD4^+ T cells, B cells, and stromal cells. e, f
   Representative images of the colocalization of nmTECs (GABRA5; DAB) and
   endothelial cells (CD31; purple) in the vicinity of GABRA5^+ cells (e)
   and not in the vicinity of GABRA5^+ cells (f). Binarized signals are
   shown in yellow. Scale bar: 100 µm. g Protein levels of CD31 near and
   not near from GABRA5+ cells in MG-thymoma quantified using microscopic
   images (Near n = 7, Non-Near n = 7. Details are provided in
   Supplementary Fig. [212]13 and Methods. Source data are also provided
   as a Source Data file.). For each group, seven areas from four MG
   patients were quantified. The signals were analyzed using a two-sided
   Mann–Whitney U-test.

Integrative analysis of MG pathology across cell types

   Recently, a computational method has been developed to infer the cell
   proportions from bulk RNA-seq datasets using references constituted by