Direct and indirect RANK and CD40 signaling regulate the maintenance of thymic epithelial cell frequency and properties in the adult thymus

Medullary thymic epithelial cells (mTECs) play a crucial role in suppressing the onset of autoimmunity by eliminating autoreactive T cells and promoting the development of regulatory T cells in the thymus. Although mTECs undergo turnover in adults, the molecular mechanisms behind this process remain unclear. This study describes the direct and indirect roles of receptor activator of NF-κB (RANK) and CD40 signaling in TECs in the adult thymus. Flow cytometric and single-cell RNA-seq (scRNA-seq) analyses suggest that the depletion of both RANK and CD40 signaling inhibits mTEC differentiation from CCL21⁺ mTEC progenitors to transit-amplifying TECs in the adult thymus. Unexpectedly, this depletion also exerts indirect effects on the gene expression of TEC progenitors and cortical TECs. Additionally, the expression levels of AP-1 genes, which enable the further subdivision of TEC progenitors, are up-regulated following the depletion of RANK and CD40 signaling. Overall, our data propose that RANK and CD40 signaling cooperatively maintain mature mTEC frequency in the adult thymus and sustain the characteristics of TEC progenitors through an indirect mechanism.

Introduction

Thymic epithelial cells (TECs) are required for the differentiation of self-tolerant T cells and regulatory T cells in the thymus. TECs are separated into cortical TECs (cTECs) and medullary TECs (mTECs) depending on their localization in the thymus (1). In addition, each TEC subset has distinct properties and functions in T cell selection and differentiation. cTECs are critical for early T cell development and positive selection of thymocytes expressing both surface makers CD4 and CD8. In contrast, mTECs ectopically express tissue-restricted self-antigens (TSAs) to filter out a wide range of self-antigen reactive T cells by apoptosis or to convert them into regulatory T cells. The TSA expression in mTECs is regulated by transcriptional regulator AIRE, which is highly expressed in mTECs expressing high levels of MHC class II and co-stimulatory molecules.

During embryonic development, both mTEC and cTEC differentiate from common bipotent progenitor cells (2, 3). For mTEC development, claudine 3 and 4-positive TECs (4), Krt19-positive mTECs (5), Tnfrsf11a-positive TECs (6), Ccl21a-positive TECs (7), and Pdpn-expressing TECs (8) were reported as mTEC progenitors giving rise to mTECs expressing AIRE and TSAs. In the adult thymus, Aire⁺ mTECs undergo a turnover of approximately 2 weeks (9), indicating the presence of mTEC progenitor maintaining the cellularity of mature mTECs. Some studies propose the progenitor of TECs in the adult thymus (10, 11). However, the phenotypes of the proposed progenitors seem to be inconsistent, implying that multiple fractions of TECs may have the potential as TEC progenitors.

Single-cell RNA-sequencing (scRNA-seq) analysis is a powerful tool for distinguishing cell types with high resolution. Recent studies utilizing scRNA-seq on TECs have highlighted their significant heterogeneity. Beyond identifying AIRE⁺ mTECs and CCL21⁺ mTECs, data analysis has revealed the presence of transit-amplifying TECs (TA-TECs) (12–15), which are proliferative progenitors for AIRE⁺ mTECs, as well as post-AIRE mTECs, including tuft-like mTECs and mimetic TECs (16, 17). Additionally, a recent study suggested the existence of TEC progenitors expressing a wide variety of keratin molecules in the human thymus (18). Moreover, a combination of scRNA-seq analysis and barcode cell labeling has proposed the presence of the early and late types of TEC progenitors in postnatal mice (19) although these TEC progenitors have not been isolated and fully characterized yet.

Mechanistically, several studies have revealed the roles of TNF family cytokine signaling in mTEC differentiation. Receptor activator of NF-κB (RANK) and CD40 play partially redundant roles in mTEC differentiation during early thymic development by activating signal transducer TRAF6- and NF-κB inducing kinase-dependent activation of transcription factor NF-κB (20–22). Additionally, lymphotoxin signaling is involved in early mTEC differentiation by inducing the expression of RANK on embryonic mTEC progenitors (6, 23), postnatal development of CCL21⁺ mTECs (24), and differentiation of post-Aire mTECs (25). The administration of a RANK ligand (RANKL) neutralizing antibody (RANKL-Ab) results in a reduction of mature mTECs (14, 26), suggesting that RANK signaling is involved in the homeostatic maintenance of AIRE⁺ mTEC frequency in the adult thymus.

Activator Protein 1 (AP-1) is a family of dimeric transcription factors including JUN, FOS, ATF, and MAF family members. AP-1 is activated by various stimuli, including cytokines and growth factors, through mitogen-activated protein kinase (MAPK) cascades, and thereby regulates numerous cellular and physiological functions (27). In a study of TEC development, FOS expression driven by the H2-Kb promoter was shown to cause thymic hyperplasia by expanding TECs (28). Additionally, RANK and CD40 signaling can activate the MAPK cascade via TRAF6 (29), a signal transducer critical for mTEC differentiation (30), implying a possible role for AP-1 in this process.

In this study, we describe how RANK and CD40 signaling cooperatively support the differentiation of CCL21+ mTECs into TA-TECs, thereby maintaining the frequencies of Aire⁺ mTECs and Post-Aire mTECs in the postnatal thymus. Unexpectedly, depletion of both RANK and CD40 signaling also has indirect effects on the gene expression profiles of TEC progenitors and cortical TECs. Additionally, after the depletion of RANK and CD40 signaling, the expression levels of AP-1 genes, which facilitate further subdivision of TEC progenitors, are up-regulated. Overall, our data suggest that these TNF family cytokine signals directly and indirectly regulate TEC frequency and properties.

Results

RANK and CD40 signaling cooperatively maintain mature mTEC cellularity

Aligned with the reported role of RANK and CD40 signaling in maintaining mature mTECs in adult mice (26), flow cytometric analysis confirmed that blocking RANKL-RANK signaling with an anti-RANKL antibody (RANKL-Ab) significantly reduces the number of mTECs expressing high MHC class II (MHCII^hiUEA-1⁺Ly51^–TECs; mTEC^hi) two weeks after the administration in mice (WT-RANKL Ab mice) compared to control IgG administration (WT-Control) whereas total thymic cells were not significantly altered (Figures 1A–C). However, approximately 10% of the mTEC^hi population persisted in the thymus of WT-RANKL Ab mice (Figure 1D). We speculated that CD40 signaling might compensate for the absence of RANK and CD40 signaling in maintaining adult mTECs, similar to its role during mTEC development in embryonic and neonatal stages (20). To test this hypothesis, we administered RANKL-Ab to Cd40-deficient (Cd40^–/–RANKL Ab) mice. Indeed, neutralizing RANK and CD40 signaling in Cd40^–/– mice resulted in a reduction of mTEC^hi cell numbers to just a few percent of those in WT-RANKL Ab mice and Cd40-deficient mice receiving control IgG (Cd40^–/–Control mice) (Figure 1D). In contrast to mTEC fractions, cell numbers of Ly51⁺UEA^– TECs (cTECs) and Ly51^–UEA^–TECs were unaffected by the RANKL-Ab administration and the Cd40-deficiency. These results suggest that RANK and CD40 signaling contribute to maintaining the frequency of mature mTECs in the adult thymus in a partially redundant manner, but not the frequency of other TECs.

Figure 1

Figure 1. Flow cytometric analysis of TECs from wild-type and Cd40-deficient mice treated with neutralizing RANKL antibody. (A) Experimental scheme for depleting RANK and CD40 signaling by the administration of RANKL antibody in mice. (B) The total thymic cell number from wild-type (WT) treated with control IgG (WT-Control), WT treated with neutralizing RANKL antibody (WT-RANKL-Ab), Cd40-deficient (Cd40^–/–) mice treated with control-IgG (Cd40^–/–Control), and Cd40^–/– mice treated with RANKL-Ab (Cd40^–/–RANKL-Ab) at 6-week-old age are summarized in graphs (n = 5 each). Bars indicate the mean value. (C) Flow cytometric analysis of UEA-1 ligand and Ly51 expressions in TECs (CD45^–Ter-119^–EpCAM⁺) from WT-Control, WT-RANKL-Ab, Cd40^–/–Control, Cd40^–/–RANKL-Ab at 6-week-old age (n = 5 each). The percentages and numbers of UEA-1⁺Ly51^– (mTEC), UEA-1^–Ly51⁺ (cTEC), and UEA-1^–Ly51^– in TECs are summarized in graphs. RANKL-Ab or control IgG was subcutaneously injected in mice at 4-week-old age. Bars indicate the mean value. Data were statistically analyzed using one-way ANOVA followed with multiple comparisons by Tukey’s test. Significant differences are indicated by *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. (D) Flow cytometric analysis of MHC class II (MHCII) and UEA-1 ligand expressions in mTECs (UEA-1⁺Ly51^– TECs) from WT-Control, WT-RANKL-Ab, Cd40^–/–Control, Cd40^–/–RANKL-Ab at 6-week-old age (n = 5). The percentages and numbers of MHCII^hiUEA-1⁺ cells and MHCII^loUEA-1⁺ mTEC in mTECs are summarized in graphs. Bars indicate the mean value. Data were statistically analyzed using one-way ANOVA followed with multiple comparisons by Tukey’s test. Significant differences are indicated by*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. (E) Flow cytometric analysis of UEA-1 ligand and Ly51 expressions in TECs from Cd40^–/– mice 2, 4, and 6 weeks after the treatment with RANKL-Ab (Cd40^–/–RANKL-Ab) or control IgG (Cd40^–/–Control), and no treatment (no injection). N = 3 each. The percentages and numbers of UEA-1⁺Ly51^– (mTEC), UEA-1^–Ly51⁺ (cTEC), and UEA-1^–Ly51^– in TECs are summarized in graphs. Bars indicate the mean value. Significant differences are indicated by *p < 0.05, **p < 0.01, ***p < 0.001. (F) Flow cytometric analysis of MHC class II (MHCII) and UEA-1 ligand expressions in mTECs from Cd40^–/– mice 2, 4, and 6 weeks after the treatment with RANKL-Ab (Cd40^–/–RANKL-Ab) or control IgG (Cd40^–/–Control mice), and no treatment (no injection). N = 3 each. MHCII^hiUEA-1⁺ cells and MHCII^loUEA-1⁺ mTEC in mTECs are summarized in graphs. Bars indicate the mean value. Data are expressed as the mean ± SD. Data were statistically analyzed using unpaired t-test. Significant differences are indicated by *p < 0.05, **p < 0.01.

Although mTEC^hi was severely reduced in the thymus of Cd40^–/–RANKL-Ab mice, mTECs expressing low levels of MHC class II (mTEC^lo) were less affected. The number of mTEC^lo cells was reduced by approximately half, with a substantial number remaining in the thymus of Cd40^–/–RANKL-Ab mice. Given that the mTEC^lo fraction includes immature mTECs in addition to post-Aire mTECs (16), it is likely that precursors for mature mTECs persist in the mTEC^lo fraction in these mice. Indeed, mature mTECs were restored in Cd40^–/–RANKL-Ab mice 4 weeks after RANKL administration, likely due to the homeostatic clearance of the injected antibody (Figures 1E, F). Moreover, 6 weeks after administration, the ratio of mTEC^lo to mTEC^hi shifted; the proportion of mTEC^lo decreased while the proportion of mTEC^hi increased in total mTECs compared to age-matched controls. This observation supports the idea that the mTEC^lo pool serves as a precursor for mTEC^hi during the rapid recovery, leading to a reduction in the relative proportion of mTEC^lo. Overall, these data suggest that immature mTECs remain in the mTEC^lo fraction in Cd40^–/–RANKL-Ab mice 2 weeks after antibody administration and differentiate into mTEC^hi following the clearance of RANKL-Ab.

RANK and CD40 signaling up-regulate cell-cycle related genes and down-regulates Ccl21a expression in mTEC^lo fraction

Given that the mTEC^lo fraction remaining after the depletion of RANK and CD40 signaling might represent the phenotype of mTEC progenitors prior to receiving these cytokine signals, we aimed to investigate the gene expression profile of a specific subfraction of mTEC^lo cells in Cd40^–/–RANKL-Ab mice. To minimize contamination from post-Aire mTECs, we selectively sorted cells within the mTEC^lo fraction that were negative for Ly6d (a marker for post-Aire mTECs) and L1CAM (a marker for tuft-like TECs) (31) (Supplementary Figure 1A). These sorted cells were then subjected to RNA sequencing (RNA-seq) analysis to elucidate their gene expression profiles. Principal component analysis (PCA) of the RNA-seq data showed that the gene expression profiles of the mTEC^lo subfraction differed significantly among wild-type (WT), WT-RANKL-Ab, Cd40^–/–, and Cd40^–/–RANKL-Ab mice (Supplementary Figure 1B). Differentially expressed genes (DEGs) were identified using a threshold of a 2-fold change with an FDR P-value < 0.05 (Figure 2A, Supplementary Table 1). The administration of RANKL-Ab to wild-type (WT) mice led to the up-regulation of 111 genes and down-regulation of 274 genes. The deletion of CD40 resulted in the up-regulation of 255 genes and down-regulation of 277 genes. Notably, administering RANKL-Ab to Cd40^–/– mice induced the up-regulation of 313 genes and down-regulation of 819 genes compared to WT-RANKL-Ab mice. Venn diagram analysis of the DEG sets revealed a significant reduction of 492 genes specifically in Cd40^–/–RANKL-Ab mice (Figure 2B). These results highlight the redundant and additive effects of RANK and CD40 signaling in regulating gene expression in immature mTECs.

Figure 2

Figure 2. RNA-seq analysis of mTEC^lo fraction from wild-type and Cd40-deficient mice receiving neutralizing RANKL antibody. (A) Volcano plots of differentially expressed genes from bulk RNA-seq data of wild-type (WT) treated with control IgG (WT-Control), WT treated with neutralizing RANKL antibody (WT-RANKL-Ab), Cd40-deficient (Cd40^–/–) mice treated with control-IgG (Cd40^–/–Control), and Cd40^–/– mice treated with RANKL-Ab (Cd40^–/–RANKL-Ab) at 6-week-old age. Red dots in volcano plots indicate genes for which expression differed significantly between the two samples (FDR P-value < 0.05, Fold change > 2). Numbers of differentially expressed genes are shown in the panels. The log2 fold change is plotted on the x-axis, and the log10 P-value is plotted on the y-axis. P-values were determined by Baggerley’s test (32) (B) The Venn diagram illustrates the overlap of down-regulated genes among three samples compared to WT mice treated with control-IgG. (C) Gene ontology enrichment analysis of the down-regulated genes and up-regulated genes in Cd40^–/– mice treated with RANKL-Ab compared to WT mice treated with control-IgG. (D) Volcano plots of differential expression of cell cycle-related gene sets. Red dots in volcano plots indicate genes for which expression differed significantly between the two samples (FDR P-value < 0.05, Fold change > 2). For cell cycle-related gene sets, mouse orthologues of the previously reported human cell cycle-related gene sets (33) were used. Numbers of differentially expressed genes are shown in the panels. The log2 fold change is plotted on the x-axis, and the log10 P-value is plotted on the y-axis. P-values were determined by Baggerley’s test (32). (E) Dot plot showing normalized gene expression value of Ccl21a. The horizontal lines show the mean. Data were statistically analyzed using one-way ANOVA followed with multiple comparisons by Tukey’s test. Significant differences are indicated by *p < 0.05, **p < 0.01, ***p < 0.001. (F) Dot plots showing normalized gene expression value of some cTEC-associated genes. Data were statistically analyzed using one-way ANOVA followed with multiple comparisons by Tukey’s test. The horizontal lines show the mean. Significant differences are indicated by *p < 0.05, **p < 0.01.

Gene Ontology (GO) analysis of down-regulated gene sets in the mTEC^lo subfraction from Cd40^–/–RANKL-Ab mice, compared to WT-control mice, revealed significant enrichment in GO terms associated with ion transport, cytoskeletal organization, and microtubule motor activity (Figure 2C, Supplementary Table 2). This suggests that RANK and CD40 signaling promote the expression of these gene sets in the mTEC^lo subfraction. Alternatively, there may be a reduction in the frequency of post-Aire mimetic mTECs that express these gene sets but not L1CAM and Ly6d, potentially influenced by RANK and CD40 signaling. Additionally, GO analysis of up-regulated gene sets indicated an increase in certain genes encoding extracellular matrix proteins, such as those in the collagen family, and cell adhesion molecules, including integrins (Figure 2C, Supplementary Table 2). Interestingly, genes coding for specific growth factor families and frizzled-binding molecules were also up-regulated (Figure 2C, Supplementary Table 2).

In addition to the GO analysis, we found that cell cycle-related gene sets were down-regulated following the disruption of RANK and CD40 signaling (Figure 2D). Furthermore, Ccl21a expression was up-regulated in the mTEC^lo subfraction of Cd40^–/–RANKL-Ab mice (Figure 2E). These findings suggest that RANK and CD40 signaling may initiate the differentiation of CCL21⁺ mTECs into transit-amplifying TECs, which serve as precursor cells for Aire⁺ mTECs (15). Alternatively, the increased expression of Ccl21a in the mTEC^lo subfraction might be due to a higher proportion of Ccl21a-expressing cells, resulting from a reduction in the frequency of post-Aire mimetic mTECs in this subfraction.

Interestingly, genes typically associated with cTECs, including Psmb11, Prss16, Tbata, and Ccl25, were up-regulated in the mTEC^lo subfraction of Cd40^–/–RANKL-Ab mice (Figure 2F). The expression level of mRNA coding FOXN1, which regulates the expression of these genes, showed a similar trend, although the change was not statistically significant (Figure 2F). This observation suggests that RANK and CD40 signaling may help suppress the aberrant expression of certain cTEC-associated genes in mTECs, potentially through the downregulation of FOXN1.

Single-cell RNA-seq analysis suggested that RANK and CD40 signaling indirectly regulate gene expressions in cTECs and TEC progenitors in the thymus

Given the high heterogeneity of TECs, the mTEC^lo subfraction identified by flow cytometric analysis may encompass multiple TEC subsets including various types of post-Aire mTECs as well as immature mTECs. To gain a more comprehensive understanding of the changes in frequency and gene expression profiles of TECs following the depletion of these cytokine signals, we conducted single-cell RNA sequencing (scRNA-seq) analysis. Droplet-based scRNA-seq was performed on the TEC fraction (EpCAM⁺CD45^–TER119^–) isolated from WT-control, WT-RANKL-Ab, Cd40^–/–control, and Cd40^–/–RANKL-Ab mice. After quality control (Supplementary Figure 2) and integration of these scRNA-seq data (Figure 3A), TEC clusters were defined based on the expression of marker genes (Figure 3A, Supplementary Figure 3). In WT-control mice, percentages of CCL21⁺ mTECs, Aire⁺ mTECs, TA-TECs, post-Aire mimetic cells, tuft-like TECs, and cTECs were 34.2%, 41.9%, 8.1%, 6.5%, 5.6%, and 1.7%, respectively (Supplementary Table 3). Given that CCL21⁺ mTECs, tuft-like mTECs, and a portion of Post-Aire mimetic cells belong to the mTEC^lo population, while Aire⁺ mTECs and the majority of TA-TECs fall within the mTEC^hi population, these data are relatively consistent with the flow cytometric data shown in Figure 1. However, cTECs appear at lower frequencies than in the flow cytometry data, which may be due to cell loss during preparation and quality control or because the cTEC fraction (Ly51⁺UEA-1^–) in flow cytometric analysis includes other cell types. In addition to the relatively well-characterized TEC subsets, we determined a cluster (Cluster 9 in Figure 3A) that appears to correspond to the early TEC progenitor population previously described (19), characterized by expression of Psmb11, Prss16, Pdpn and Krt5 (Figure 3B, Supplementary Figure 3).

Figure 3

Figure 3. Single RNA-seq analysis of TECs from wild-type and Cd40-deficient mice receiving neutralizing RANKL antibody. (A) Uniform Manifold Approximation and Projection (UMAP) plot of droplet-based scRNA-seq data of TECs (CD45^–Ter-119^–EpCAM⁺). The scRNA-seq data from wild-type (WT) mice treated with control IgG, WT mice treated with a neutralizing RANKL antibody (RANKL-Ab), Cd40-deficient (Cd40^–/–) mice treated with control IgG, and Cd40^–/– mice treated with RANKL-Ab were integrated using the Seurat package. These mice were sacrificed at 6-week-old age, 2 weeks after antibody treatment. In the plot, cell clusters are distinguished by colors and numbers and are assigned based on marker gene expression (Supplementary Figure 3). T cell cluster was removed, and cluster 12 showed expression of Lck (Supplementary Figure 3) and were assigned as nurse cTECs. (B) Violin plots showing the expression levels of Prss16, Psmb11, Pdpn, and Ccl21a in each cluster. The expression levels of these genes in WT treated with control IgG were exhibited. Each dot represents expression levels in individual cells. (C) Monocle trajectory and pseudotime analyses of TEC scRNA-seq data from WT treated with control IgG. The green circle indicates the root node when cluster 9 is considered as TEC progenitors. (D) Plot of the area under the recovery curve (AUC), which reflects the enrichment of each regulon, versus pseudotime predicted from Monocle. The regulon of ATF3 and JUN were determined by the SCENIC program. The number of regulons for each gene is shown in parentheses. ATF3-”extended” includes both high-confidence annotations (which are based on direct annotation and those inferred by orthology) as well as lower-confidence annotations that are inferred by motif similarity.

Monocle trajectory analysis (34) using WT-Control TEC clusters indicated that cluster 9 is situated between the CCL21⁺ mTEC and cTEC clusters (Figure 3C), supporting the idea that this cluster likely represents a progenitor TEC population. When the root node was set in cluster 9 (indicated by the green circle in Figure 3C), a combined analysis using Monocle and SCENIC tools (35) demonstrated an increased activity of ATF3- and JUN-associated regulons—groups of genes regulated by shared transcription factors—during the differentiation of progenitor clusters into CCL21⁺ mTECs (Figure 3D, Supplementary Table 2). These findings suggest that the activities of these AP-1 transcription factors may play a role in driving the differentiation of progenitor cells into the mTEC lineage.

Comparison of the scRNA-seq data clusters among WT-control, WT-RANKL-Ab, Cd40^–/–control, and Cd40^–/–RANKL-Ab mice revealed a marked reduction in the frequencies of Aire⁺ mTECs, TA-TECs, post-Aire mTECs, and tuft-like mTECs in Cd40^–/–RANKL-Ab mice (Figures 4A, B). In contrast, RANKL-Ab administration in wild-type mice resulted in a milder reduction of these mTEC subsets (Figures 4A, B). These findings are consistent with those from flow cytometric analysis, further supporting the functional overlapping of RANK and CD40 signaling in mTEC maintenance in the adult thymus. In contrast, the CCL21⁺ mTEC, cTEC, and TEC progenitor clusters appeared to remain in the thymus of Cd40^–/–RANKL-Ab mice. Within the Aire⁺ mTEC clusters, Cluster 7, which likely represents a non-proliferative transition stage between CCL21⁺ mTECs and Aire⁺ TA-TECs, was less affected in Cd40^–/–RANKL-Ab mice. This observation suggests that depletion of both RANKL and CD40 signaling may lead to a differentiation arrest of mTECs at this stage.

Figure 4

Figure 4. Differentially expressed gene and trajectory analyses of single-cell RNA-seq data from wild-type and Cd40-deficient mice receiving neutralizing RANKL antibody. (A) UMAP plot of droplet-based scRNA-seq data of TECs (CD45^–Ter-119^–EpCAM⁺) from wild-type (WT) mice treated with control IgG (WT-Control), WT mice treated with a neutralizing RANKL antibody (WT-RANKL-Ab), Cd40-deficient (Cd40^–/–) mice treated with control IgG (Cd40^–/– Control), and Cd40^–/– mice treated with RANKL-Ab (Cd40^–/– RANKL-Ab). The integrated UMAP plot in Figure 3 was separated into each data set. The number of sequenced cells after the quality control is indicated in parentheses. (B) Percentages of cell subsets in total TECs were compared among the scRNA-seq data from WT-Control, WT-RANKL-Ab, Cd40^–/– Control, and Cd40^–/– RANKL-Ab mice. (C) MA plots show differentially expressed genes between WT-control and Cd40^–/– RANKL-Ab mice from scRNA-seq data for each cell cluster subset. The log2 average expression level (CPM) is plotted on the x-axis, and the log2 fold change is plotted on the y-axis. Genes with log2 fold change greater than 0.15 or less than –0.15 (FDR P <0.05) are represented by blue dots. Red dots indicate cTEC-associated genes, violet dots indicate AP-1 transcription factor genes, and green dots indicate interferon-stimulated genes. The orange dot represents Ccl21a. (D) Total expression of ATF3 and JUN regulon genes in the TEC progenitors, CCL21⁺ mTECs, and cTECs subclusters of WT-Control, WT-RANKL-Ab, Cd40^–/– Control, and Cd40^–/– RANKL-Ab mice. * *p < 0.05, **p < 0.01, and ****p < 0.0001.

As previously reported (15), the TA-TEC subcluster was divided into Aire⁺ TA-TECs and Ccl21⁺ TA-TECs (Supplementary Figure 4). To address the influence of the depletion of RANK and CD40 signaling on proliferative activity of CCL21⁺ mTECs and AIRE⁺ mTECs, we estimated the proportion of Ccl21⁺ TA-TECs within the total Ccl21⁺ mTECs, and similarly for Aire⁺ mTECs. Data suggested that the proliferative activity of CCL21⁺ mTECs seemed to be influenced by the depletion of these signals (Supplementary Figure 4), which is consistent with bulk RNA-seq analysis (Figure 2D). Overall, scRNA-seq analysis suggested that RANK and CD40 signaling maintain the frequency of Aire⁺ mTECs and post-Aire mimetic mTECs by promoting the differentiation of CCL21⁺ mTEC into TA-TECs in the adult thymus.

In line with findings from bulk RNA-seq analysis, differential gene expression analysis of scRNA-seq data revealed that the depletion of RANK and CD40 signaling leads to the up-regulation of cTEC-associated genes and Ccl21a in CCL21⁺ mTEC clusters (Figure 4C). Additionally, several interferon-stimulated genes (ISGs) were notably down-regulated in these clusters (Figure 4C). Interestingly, these changes in gene expression were observed not only in CCL21⁺ mTECs, which are the primary recipients of RANK and CD40 signaling, the primary recipients of RANK and CD40 signaling, but also in TEC progenitors and cTEC clusters (Figure 4C), with subcluster composition confirmed as unchanged in cTECs. (Supplementary Figure 5). Given that the RANK expression of these cell types is virtually absent (Supplementary Figure 5), this unexpected finding suggests an indirect regulatory mechanism of gene expression driven by RANK and CD40 signaling. Furthermore, the loss of RANK and CD40 signaling resulted in the up-regulation of some AP-1 transcription factor genes within the progenitor cell subset (Figure 4C). Consistently, the SCENIC analysis suggested an increase in the activity of ATF3- and JUN-inducing regulons in the progenitor TECs (Figure 4D). Collectively, these results imply that under normal conditions, RANK and CD40 signaling may act to indirectly suppress gene regulatory networks governed by AP-1 transcription factors in progenitor cells, highlighting a complex interplay of direct and indirect signaling pathways in maintaining TEC homeostasis.

TEC progenitors are classified into subpopulations with unique gene expression profiles

Given that our data suggest TEC progenitors are indirectly influenced by RANK and CD40 signaling, we focused our analysis on these cells. Subclustering of the progenitor cluster from the scRNA-seq data revealed four distinct subclusters with unique gene expression profiles (Figure 5A, Supplementary Table 4), all exhibiting similar levels of Pdpn expression (Figure 5B). Cluster S1 showed high expression of cTEC-associated genes, such as Prss16 and Psmb11 (Figure 5B), suggesting a bias toward the cTEC lineage. In contrast, cluster S2 displayed high levels of Ccl21a expression in a part of the cells, indicating a bias toward the mTEC lineage. Clusters S0 and S3 exhibited low expression of both cTEC-associated genes and Ccl21a (Figure 5B). Notably, cluster S3 was characterized by the elevated expression of AP-1 transcription factor family genes (Figure 5B, Supplementary Figure 6). These findings underscore the heterogeneous composition of TEC progenitors, categorized by their expression levels of cTEC-associated genes and AP-1 family genes.

Figure 5

Figure 5. Subclustering analysis of the TEC progenitor subset in scRNA-seq data from wild-type and Cd40-deficient mice receiving neutralizing RANKL antibody. (A) UMAP plot showing TEC progenitor clusters in the droplet-based scRNA-seq data. (B) Violin plots depicting expression levels of Pdpn, Prss16, Psmb11, Ccl21a, Atf3, and Jun across each subcluster of TEC progenitors. Each dot represents the expression level in individual cells. (C) The left panel shows a trajectory analysis of subclusters predicted by Monocle 3, visualized by UMAP. Cells are color-coded by pseudotime, transitioning from purple to yellow as pseudotime progresses. The right panel displays the dynamics of gene expression along pseudotime for Prss16, Psmb11, Ccl21a, and Pdpn. (D) The integrated UMAP plot from panel (A) is separated by dataset. The percentages of each cell cluster across the four scRNA-seq datasets are shown in the graph. (E) Box plots illustrating the expression levels of Atf3, Fos, Junb, Egr1, and Jun in each subcluster for WT-control and Cd40^–/–RANKL-Ab. Each dot represents the expression level in an individual cell.

To understand the lineage connections among these clusters, we applied Monocle trajectory analysis to the progenitor cluster. Interestingly, the trajectory analysis suggested an ordering of the clusters in the sequence S3, S1, S0, and S2 (Figure 5C). Assuming that the S3 cluster represents the root node, the pseudotime analysis indicated that the cTEC-biased cluster S1 may differentiate into the mTEC-biased cluster S2 through the non-biased cluster S0.

Comparing the frequencies of subcluster subsets across WT-control, WT-RANKL-Ab, Cd40^–/–control, and Cd40^–/–RANKL-Ab mice suggested that the frequency of the S3 cluster increase additively with the elimination of RANK and CD40 signaling (Figure 5D). In addition, expression of Atf3, Fos, Junb, Egr1 in other subclusters including cTEC-biased cluster S1 was increased in Cd40^–/–RANKL-Ab mice (Figure 5E). Thus, the depletion of RANK and CD40 signaling increases the expression level of AP-1 family genes and the frequency of subsets expressing AP-1 family genes in TEC progenitors.

Integrative analysis of droplet-based scRNA-seq, well-based scRNA-seq and flow cytometric analyses suggested that the TEC progenitors are present in Ly51^-UEA-1^-TEC and cTEC fractions

Our data indicated that TEC progenitor cells, as identified in scRNA-seq analysis, are divided into four clusters depending on gene expression profile. To further validate these subpopulations, we aimed to correlate the scRNA-seq clusters with mTEC and cTEC surface markers in flow cytometric analysis. To this end, we first performed single-cell sorting of UEA-1⁺Ly51^– TECs (mTEC-enriched), UEA-1^–Ly51⁺ TECs (cTEC-enriched), and UEA-1^–Ly51^– TECs (other TECs) from wild-type mouse thymus, followed by RNA-seq of the individual sorted cells (Figure 6A). We then integrated these well-based scRNA-seq data with the droplet-based scRNA-seq data. After assigning each cluster to typical TEC subsets (Supplementary Figure 7), we determined the cell types of the individual cells sorted by flow cytometric analysis (Figure 6B).

Figure 6

Figure 6. Integration of droplet-based and well-based scRNA-seq data of TECs. (A) Schematic diagram illustrating the integration analysis of well-based scRNA-seq data with droplet-based scRNA-seq data. Single cells were sorted from UEA-1⁺Ly51^- TECs (mTECs), UEA-1^-Ly51⁺ TECs (cTECs), and UEA-1^-Ly51^- TECs (DN) fractions of wild-type 6-week-old mice, and subjected to well-based scRNA-seq analysis. The well-based scRNA-seq data was integrated with the droplet-based scRNA-seq data shown in Figure 4. (B) UMAP plot and clustering after the integration of droplet-based and well-based scRNA-seq data. Cell clusters were reassigned based on marker gene expression following the integration (Supplementary Figure 4). UMAP plots of individual cells sorted from UEA-1⁺Ly51^- TECs (mTECs), UEA-1^-Ly51⁺ TECs (cTECs), and UEA-1^-Ly51^- TECs (DN) were overlaid on the droplet-based scRNA-seq data (in gray). (C) Percentages of TEC subsets among total single cells in UEA-1⁺Ly51^- TECs (mTECs), UEA-1^-Ly51⁺ TECs (cTECs), and UEA-1^-Ly51^- TECs (DN).

As expected, our analysis revealed that the UEA-1⁺Ly51^– mTEC fraction includes both CCL21⁺ mTECs and AIRE⁺ mTECs (Figure 6B). Notably, the UEA-1^–Ly51^– TEC fraction contains approximately 30% of cells classified as TEC progenitors (Figure 6C), along with some contamination from various mTEC subsets, likely due to the loss of UEA-1 binding ligands during TEC sample preparation using collagenase digestion. Additionally, the UEA-1^–Ly51⁺ cTEC fraction also contains TEC progenitors alongside mature cTECs. Overall, our data suggest that the TEC progenitors identified in scRNA-seq analysis are negative for UEA-1 binding ligands and are further distinguished based on Ly51 expression levels in flow cytometric analysis.

We next assigned sorted individual progenitor cells to the subpopulations of TEC progenitors identified through droplet-based scRNA-seq analysis. Data analysis revealed that the UEA-1^–Ly51^– TEC fraction contains all types of the TEC progenitor subpopulation (Figure 7A). In contrast, with one exception, progenitor cells sorted from the UEA-1^–Ly51⁺ TECs predominantly belong to the cluster S1 (Figure 7A), which showed high expressions of cTEC genes (Figure 5B).

Figure 7

Figure 7. Subclustering analysis of TEC progenitors in integrated scRNA-seq data. (A) UMAP and subclustering of TEC progenitors in the integrated data of the droplet-based scRNA-seq and the well-based scRNA-seq. (B) Dot plots of the normalized expression values of Prss16, Psmb11, Pdpn, Atf3, and Jun for each cell cluster. Horizontal lines show the median. (C) Scatter plot showing the fluorescence intensity of UEA-1 ligand and Ly51 for individual sorted cells assigned as TEC progenitor cells. One dot represents one cell, color-coded for each cluster. Total TECs are indicated as gray dots. Gates used for cell sorting are indicated.

We further investigated the expression level of the cTEC-associated genes and others in the sorted individual cells. In consistent with the droplet-based scRNA-seq data, expression levels of Prss16 and Psmb11 were highest in individually sorted cells assigned as the S1 subpopulation (Figure 7B). However, their expression levels were remarkably lower as compared to those in mature cTECs. Expression of Pdpn was detected in all subpopulations with almost the same level and may be slightly higher than that in CCL21⁺ mTECs (Figure 7B). Sorted single cells assigned as the cluster S3 exhibited high levels of Atf3 and Jun expressions. In contrast, Atf3 and Jun expression levels were lower in other subpopulations in addition to CCL21⁺ mTECs (Figure 7B). These data further confirmed that TEC progenitors are subdivided by expression levels of some AP-1 genes and cTEC-associated genes.

We confirmed the distribution of these progenitor subpopulations in the flow cytometric profile. As expected, subpopulations S0, S2, S3, and part of S1 were derived from the UEA-1^–Ly51^– fraction and could not be distinguished based on the expression level of these markers (Figure 7C). Interestingly, cells in the S1 cluster within the UEA-1^–Ly51⁺ TEC fraction exhibited lower surface Ly51 expression compared to cells classified as mature cTECs. This finding suggests that part of the cTEC-biased subpopulation of TEC progenitors is present within the Ly51^loUEA-1^– fraction in flow cytometric analysis.

Depletion of RANK and CD40 signaling in adult thymus leads to the reduction in frequencies of regulator T cells, natural killer T cells and eosinophils

Flow cytometric analysis on thymocyte fractions showed the depletion of both CD40 and RANK signaling increased in the ratio of CD4SP cells, whereas their cell number was not significantly increased (Figure 8A). Interestingly, a significant change in double negative fractions was observed, which may be due to the decrement in early thymocyte progenitors or minor cell subsets (Figure 8A). Consistently, the frequencies of eosinophil, which is critical for thymus regeneration (36), and natural killer T cells were reduced in the thymus of Cd40^–/– mice receiving RANKL-Ab (Figures 8B, C). In addition, the frequency of CCR6^–Foxp3⁺ regulatory T cells was reduced (Figure 8D), implying the mTEC-mediated function of Treg selection might be impaired. These data suggested that the depletion of both RANK and CD40 signaling in adult thymus may impact thymic self-tolerance and recovery from the thymic injury.

Figure 8

Figure 8. Flow cytometric analysis of Thymocytes from wild-type and Cd40-deficient mice treated with neutralizing RANKL antibody. (A) Flow cytometric analysis of CD8a and CD4 expressions in thymocytes from wild-type (WT) treated with control IgG (WT-Control), Cd40-deficient (Cd40^–/–) mice treated with control-IgG (Cd40^–/–Control), and Cd40^–/– mice treated with RANKL-Ab (Cd40^–/–RANKL-Ab) at 6-week-old age (n = 5 each). The percentages and numbers of CD4^–CD8^–, CD4⁺CD8⁺, CD4⁺CD8^–, and CD4^–CD8⁺ in thymocytes are summarized in graphs. RANKL-Ab or control IgG was subcutaneously injected in mice at 4-week-old age. (B) Flow cytometric analysis of Foxp3 and CD25 expressions in CCR6^–CD4SP thymocytes (CD4^–CD8^–CCR6^–) from WT-Control, WT-RANKL-Ab, Cd40^–/–Control, Cd40^–/–RANKL-Ab at 6-week-old age (n = 5 each). The percentages and numbers of CCR6^–Foxp3⁺ regulatory T cells in CCR6^–CD4SP thymocytes are summarized in graphs. RANKL-Ab or control IgG was subcutaneously injected in mice at 4-week-old age. (C) Flow cytometric analysis of Siglec-F and SSC-A expressions in thymocytes from WT-Control, WT-RANKL-Ab, Cd40^–/–Control, Cd40^–/–RANKL-Ab at 6-week-old age (n = 5 each). The percentages and numbers of eosinophil in thymocytes are summarized in graphs. RANKL-Ab or control IgG was subcutaneously injected in mice at 4-week-old age. (D) Flow cytometric analysis of TCRβ and CD1d expressions in thymocytes from WT-Control, WT-RANKL-Ab, Cd40^–/–Control, Cd40^–/–RANKL-Ab at 6-week-old age (n = 5 each). The percentages and numbers of natural killer T cells in thymocytes are summarized in graphs. RANKL-Ab or control IgG was subcutaneously injected in mice at 4-week-old age. (A-D) Bars indicate the mean value. Data were statistically analyzed using one-way ANOVA followed with multiple comparisons by Tukey’s test. Significant differences are indicated by *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Discussion

Previous studies showed that the administration of RANKL-Ab causes a reduction in AIRE⁺ and MHCII^hi mTECs (14, 26). Consistently, our data also indicate the reduction in these mTEC subsets by the RANKL-Ab administration. The extent of reduction in AIRE⁺ and MHCII^hi mTECs appears milder in our study than in previous studies. Moreover, the recovery kinetics was faster in our study compared to that in the previous study (14, 26). These differences may be attributed to the repeated injections of RANKL-Ab in the prior research (14, 26) compared to the single injection in our study. Our data also show that the reduction in AIRE⁺ mTECs caused by RANKL-Ab injection was considerably more pronounced on the CD40-deficient background, suggesting a partial compensation for the loss of RANK and CD40 signaling by CD40 signaling, which plays a role in early mTEC development (20). Thus, cytokine signaling that promotes mTEC differentiation may be similar during the developmental processes in the embryonic and neonatal periods and the turnover process in the adult thymus.

Whereas the depletion of both RANK and CD40 signaling led to a marked reduction in AIRE⁺ mTECs, Post-Aire⁺ mTECs, and TA-TECs, CCL21⁺ mTECs persisted, though their gene expression profile was altered. These findings suggest that RANK and CD40 signaling may facilitate the differentiation of CCL21⁺ mTECs into TA-TECs. Additionally, our scRNA-seq analysis indicates that these signals may also promote the proliferation of CCL21⁺ mTECs, consistent with previous research showing CD40 signaling’s role in supporting CD80^–MHCII^– mTEC proliferation (37)). Notably, RANK and CD40 expression levels were relatively higher in AIRE⁺ mTECs compared to CCL21⁺ mTECs, suggesting that these signals may further enhance both proliferation and gene expression in AIRE+ mTECs. This idea aligns with previous findings that NF-κB signaling, activated by these pathways, might enhance AIRE expression in mTECs (38). Future studies are needed to clarify the specific roles of RANK and CD40 signaling in AIRE+ mTECs.

One limitation of our study is the use of Cd40^–/– mice, in which CD40 signaling is eliminated in mTECs postnatally because CD40L is virtually undetectable in the fetal thymus (20). Consequently, mTECs in Cd40^–/– mice are differentiated and maintained solely by RANKL signaling from embryonic stages through adulthood. This exclusive reliance on RANKL may lead mTECs to develop an unusual dependency on RANKL signaling compared to normal conditions. To address this concern, using CD40L-neutralizing antibodies instead of Cd40^–/– mice might provide a more refined approach. Although there are currently no studies demonstrating that CD40L-neutralizing antibodies can transiently reduce mTEC populations, future research should investigate the effects of a simultaneous and temporary reduction in both RANK and CD40 signaling, given the evidence supporting the effectiveness of this approach.

A previous study suggested that Pdpn-expressing TECs are localized at the cortico-medullary junction of the thymus, where they are referred to as junctional TECs (jTECs) (8). Additionally, RNA-seq analysis in that study revealed that jTECs express Ccl21a (8). Consistent with these findings, our scRNA-seq analysis showed that CCL21⁺ mTEC clusters express Pdpn. Notably, both our scRNA-seq analysis and previous studies have demonstrated that Pdpn is expressed in TEC progenitors (19). Sub-clustering analysis further revealed that Pdpn is present in all subpopulations of TEC progenitors. Together, these results suggest that Pdpn marks not only the CCL21⁺ mTEC precursor pool but also TEC progenitors.

DEG analysis of scRNA-seq data suggested that the depleting RANK and CD40 signaling affedts gene expression profiles not only in RANK-expressing CCL21⁺ mTECs, but also in cTECs and TEC progenitors. First, several interferon-stimulated genes were down-regulated in these cell types in Cd40^–/–RANKL-Ab mice. Type I and III interferons were reportedly expressed in a part of AIRE⁺ mTECs, thereby influencing phenotypes of thymic antigen-presenting cells such as conventional dendritic cells (39). Consequently, interferon signaling could also impact gene expression profiles in cTECs and TEC progenitor cells, suggesting intercellular communications between mTEC and both cTEC and TEC progenitors. which may affect TEC development (40). Second, AP-1 family gene expression was upregulated in progenitors. This observation suggests that, beyond interferon signaling, another intercellular communication between mTECs and progenitors may indirectly regulate AP-1 expression levels in progenitor cells (Supplementary Figure 8). Finally, an indirect signaling mechanism appears to suppress the upregulation of cTEC genes in TEC progenitors and cTECs. The upregulation of cTEC-associated genes following RANK and CD40 signaling depletion was minimal, suggesting that this mechanism likely has a limited impact on TEC functions. Overall, further research is required to clarify the mechanisms underlying these indirect regulatory pathways and their influence on TEC phenotypes.

A previous study reported that enhancing RANK and CD40 signaling after thymic injury induces lymphotoxin α expression in lymphoid tissue inducer (Lti) cells, thereby promoting TEC regeneration, including cTECs. It is possible that a severe reduction in RANK receptors due to the loss of mature mTECs leads to an increase in free RANKL concentration, potentially activating Lti cells to upregulate lymphotoxin α, which in turn supports cTEC gene expression. This hypothesis warrants further investigation in future studies.

Deletion of RANK and CD40 signaling causes up-regulation of AP-1 transcription factor genes selectively in TEC progenitors. Sub-clustering analysis suggested that the increment of TEC progenitor subpopulation expressing high levels of AP-1 transcription factor genes and up-regulation of these genes in some TEC progenitor subpopulations. This finding suggests that, besides the interferon signaling, another indirect signaling between TEC progenitors and mTECs may suppress expression of these genes. In a previous study, H2-Kb promoter-driven Fos expression leads to thymic hyperplasia via the expansion of TECs (28). Thus, the increment of TEC progenitors expressing AP-1 genes may result in enhancing the development of TECs. Thus, RANKL and CD40 signaling homeostatically can suppress AP-1 transcription factor genes in TEC progenitor by a negative feedback loop. Abolishing this negative feedback system may cause the increment of the Fos-expressing TEC progenitors that give rise to mTECs, which could contribute to the recovery of mTECs.

TEC progenitors were separated into four subpopulations. TEC progenitors expressing a high level of some AP-1 genes may differentiate into a cTEC-biased subpopulation and subsequently into an mTEC-biased subpopulation. A previous study suggested that adult mTECs differentiate from mTEC lineage progenitors derived from Psmb11-positive cells (41). Furthermore, a recent study defined this TEC progenitor as an early TEC progenitor based on this observation (19). Our single-cell study revealed that the expression level of Psmb11 in TEC progenitors was approximately ten times lower than that in cTECs. Consequently, this suggests that Psmb11^loPdpn⁺ TEC progenitor subpopulation could contribute to maintaining the frequency of adult mTECs. A fate-mapping study using specific marker genes in this subpopulation would be crucial for addressing this issue. Ultimately, our findings illuminate the crucial roles of RANK and CD40 signaling in maintaining mTEC frequency and TEC progenitor properties in the postnatal thymus, offering promising avenues for developing strategies to address thymic hypofunction associated with aging and various stressors.

Materials and methods

Mice and antibody treatment

Female wild-type C57BL/6 mice, aged 3-4-weeks-old, were purchased from CLEA Japan. Cd40-deficient mice were established on a C57BL/6 background. All mice were maintained in standard controlled conditions with a 12-h lighting cycle and access to chow and water ad libitum, housed under specific pathogen-free conditions and handled in accordance with Guidelines of the Institutional Animal Care and Use Committee of RIKEN, Yokohama Branch (2018-075). Rat IgG-Isotype Control antibody (abcam, R&D Systems) or Anti-mouse RANK ligand neutralizing antibody (Anti-RANKL antibody: Mab clone OYC1, Oriental enzyme) (42) was injected subcutaneously at 5 mg/kg into C57BL/6J background wild-type mice or Cd40-deficient mice.

Isolation and flow cytometric analysis of TECs from mice

Mice were sacrificed using CO2, and thymi were dissected and placed into cold 1× PBS. Adhering non-thymus tissue was carefully cleared off using sharp tweezers under a fluorescence stereomicroscope. Thymi were minced with a razor blade and pipetted up and down in 1 mL of RPMI 1640 (Wako) to remove lymphocytes. Then, thymic fragments were digested in RPMI 1640 containing Liberase (Roche, 0.05U/mL) and DNase I (Sigma-Aldrich, 0.01% w/v) by incubating three times at 37°C for 12 min each. The supernatant was collected, added to 2 mL of FACS buffer (D-PBS (-) with 2% FBS) containing 1 mM EDTA, and centrifuged at 1500 rpm for 5 min. The supernatant was removed and suspended in FACS buffer. After filtering through a 67-µm nylon monofilament mesh, the resulting cell suspension was incubated with anti-mouse CD16/32 (BioLegend, Cat#101302) in FACS buffer to block nonspecific binding. For flow cytometric analysis and bulk RNA-seq, cells were stained with primary antibodies (APCCy7-labeled anti-CD45; BioLegend Cat#103116, APCCy7-labeled anti-TER119; BioLegend Cat#116223, BV510-labeled anti-EpCAM; BioLegend Cat#118231, PerCPCy5.5-labeled anti-Ly51; BioLegend Cat#108316, Alexa647-labeled anti-L1CAM; R&D Cat#FAB5674R, PECy7-labeled anti-I-A/I-E; BioLegend Cat#107630, BV711-labeled anti-CD104; BD Cat#123609, FITC-labeled anti-Ly-6D; BioLegend Cat#138606, biotinylated UEA-1; Vector Laboratories Cat#B-1065) in FACS buffer and sequentially incubated with secondary reagent (Alexa700-labeled Streptavidin; Invitrogen Cat#S2183) in FACS buffer. Dead cells were excluded by staining with SYTOX™ Blue. For droplet-based scRNA-seq, cells were stained with antibodies (APCCy7-labeled anti-CD45, APCCy7-labeled anti-TER119, FITC-labeled anti-EpCAM; BioLegend Cat #118208) in FACS buffer. For well-based scRNA-seq, cells were stained with primary antibodies (APCCy7-labeled anti-CD45, APCCy7-labeled anti-TER119, FITC-labeled anti-EpCAM, Alexa647-labeled anti-Ly51; BioLegend Cat#108312, biotinylated UEA-1) in FACS buffer and depleted of hematopoietic cells and erythrocytes by Magnetic-Activated Cell Sorting (MACS) using APC-MicroBeads (Miltenyi Biotec). Then, cells were stained with secondary reagent (PECy7-labeled Streptavidin; Invitrogen Cat# 25-4317-82) in FACS buffer. Dead cells were excluded by staining with 7-Aminoactinomycin D. Cells were sorted using a FACS Aria instrument (BD). Data were analyzed using Flowjo 10.

Bulk RNA-seq analysis

Cells were sorted using a cell sorter (Aria; BD) into 1.5 ml tube with 20 μL of cell lysis solution (2xTCL, 2-Mercaptoethanol). Cell lysis solution or RNase-free water was added to the sorted sample to achieve the final 1x TCL, mixed using a vortex, and the mixture was kept on ice for 5 min. The mixture was centrifuged at 13,000 rpm for 1 min and then stored at -80°C. Cell lysate was dissolved on ice and then purified using total x2.2 volumes of RNAClean XP Beads using Magna Stand. The mixtures were eluted with 40U of RNasin^® Plus Ribonuclease Inhibitor in RNase-free water. The supernatant was collected using Magna Stand and denatured at 65°C for 5 min. The mixture was rapidly cooled on ice for 2 min, and then added 10 μL of DNase I solution (PrimeScript Buffer and 2U of DNase I, Amplification Grade in RNase-free water). The mixtures were incubated in a thermal cycler at 30°C for 15 min. Ten μL of first strand cDNA synthesis solution (PrimeScript Buffer, PrimeScript RT Enzyme Mix I, 1 μg of T4 Gene 32 Protein, 6 pmol Oligo(dT)18 Primer and 100 pmol 1st-NSR primer) was added to the DNase I-treated mixture. The mixtures were incubated in a thermal cycler at 25°C for 10 min, 30°C for 10 min, 37°C for 30 min, 50°C for 5 min and 94°C for 5 min. Twenty μL of second strand cDNA synthesis solution (NEBuffer™ 2, 0.625mM dNTP Solution Mix, 500 pmol 2nd-NSR primer and 6.5 U of Klenow Fragment (3’→5’ exo-) in RNase-free water) was added to the first strand cDNA lysate. The mixtures were incubated in a thermal cycler at 16°C for 60 min, 70°C for 10 min. The mixtures were purified with 100 μL of AMPure XP SPRI beads (Beckman Coulter) using Magna Stand, and the concentration were then quantified using Qubit™ dsDNA Quantification Assay Kits. Of the purified dsDNA, 1 ng was used for library preparation, and the rest was stored at -80°C. Thirty μL of tagmentation solution (10 mM Tris-HClpH 8.5, 5 mM MgCl2, 10% N, N-Dimethylformamide and Tn5-linker complex in RNase-free water) and incubated at 55°C for 10 min. Zero-point two percent SDS were added to the mixture and incubated at room temperature for 5 min. Then, the mixtures were purified using Monarch^® PCR & DNA Cleanup Kit and eluted into 15 μL of buffer EB. After ligation of adapters using PCR on 25 µl of the purified mixture, sequencing library DNA was purified with x1.2 volumes of AMPure XP SPRI beads and eluted into 15 µL of buffer EB. The sequencing library was sequenced in multiplex on the HiSeqX_Ten platform. FASTQ files were processed using Fastp (43) and then quantified for annotated genes using CLC Genomics Workbench (Version 21.0.6, QIAGEN). Differential expression analysis was performed using Proportion-based Statical Analysis on CLC (Version 23.0.4).

Droplet-based scRNA-seq analysis

For scRNA-seq analysis, cell suspensions of thymi from three mice were prepared and pooled for each individual scRNA-seq experiment. Cellular suspensions were loaded onto a Chromium instrument (10× Genomics) to generate a single cell emulsion. scRNA-seq libraries were prepared using Chromium Next GEM Single Cell 3' GEM, Library & Gel Bead Kit v3.1 and sequenced in multiplex on the HiSeqX Ten platform. FASTQ files were processed using Fastp. Reads were demultiplexed and mapped to the mm10 reference genome using Cell Ranger (v 5.0.1). Processing of data with the Cell Ranger pipeline was performed using the HOKUSAI supercomputer at RIKEN and the NIG supercomputer at ROIS National Institute of Genetics. Expression count matrices were prepared by counting unique molecule identifiers. Downstream single-cell analyses (integration of datasets, correction of dataset-specific batch effects, UMAP dimensional reduction, cell cluster identification, conserved marker identification, and regressing out cell cycle genes) were performed using Seurat v4. Genes that were expressed in more than five cells and cells expressing at least 200 genes were selected for analysis. Cells that contained a percentage of mitochondrial transcripts greater than 13% to 25% were filtered out. Four scRNA-seq datasets were integrated with a combination of Find Integration Anchors and Integrate Data functions. Resolution was set as 0.43 for the FindClusters function. Murine cell cycle genes equivalent to human cell cycle genes listed in Seurat were used for assigning cell cycle scores. Trajectory analysis was performed using Monocle 3.

Well-based scRNA-seq analysis

Single cells were sorted using a cell sorter (Aria; BD) into 96-well PCR plates with 1 μL of cell lysis solution (1:10 Cell Lysis buffer [Roche], 10 U/μL Rnasin plus Ribonuclease inhibitor [Promega]) in each well, shaken at 1400 rpm for 1 min using a thermo mixer and then stored at -80°C. Cell lysate was dissolved on ice and then denatured at 70°C for 90 sec. To eliminate genomic DNA contamination, 1 μL of genomic DNA digestion solution (PrimeScript Buffer and 0.2 U of DNase I Amplification Grade in RNase-free water) was added to each denatured sample. The mixtures were shaken at 1400 rpm for 1 min using a thermo mixer, and then incubated in a thermal cycler at 30°C for 5 min and held on ice until the next step. One μL of first strand cDNA synthesis solution (PrimeScript Buffer, 8 pmol 1st-NSR primer, 0.6 pmol Oligo(dT)18 Primer, 100 ng of T4 gene 32 protein and PrimeScript RT Enzyme Mix I in RNase-free water) was added to each digested lysate. The mixtures were shaken at 1400 rpm for 1 min using a thermo mixer, and then incubated in a thermal cycler at 25°C for 10 min, 30°C for 10 min, 37°C for 30 min, 50°C for 5 min and 94°C for 5 min. Two μL of second strand synthesis solution (NEBuffer™ 2, 0.625 mM dNTP Solution Mix, 40 pmol 2nd-NSR primer and 0.75U Klenow Fragment (3’→5’ exo-) in RNase-free water) was added to each first strand cDNA lysate. The mixtures were shaken at 1400 rpm for 1 min using a thermo mixer, and then incubated in a thermal cycler at 16°C for 60 min, 70°C for 10 min. The mixtures were purified 15 μL of AMPure XP SPRI beads (Beckman Coulter) diluted two-fold with Pooling buffer (20% PEG8000, 2.5 M NaCl, 10 mM Tris-HClpH 8.0, 1 mM EDTA, 0.01% NP40) using Magna Stand. The mixtures were eluted with 3.75 μL of tagmentation solution (10 mM Tris-HClpH 8.5, 5 mM MgCl2 and 10% N, N-Dimethylformamide in RNase-free water). 1.25 μL of diluted Tn5-linker complex was added to the eluate and the mixtures were incubated at 55°C for 10 min, and then One point two five μL of 0.2% SDS was added and incubated at room temperature for 5 min. After PCR for adaptor ligation, sequencing library DNA was purified using AMPure XP SPRI beads and eluted into 25 μL of buffer EB. Reads were demultiplexed and mapped to the mm10 reference genome with STAR. Cells with less than or more than half the average count of reads detected were excluded from the analysis. Integration of well-based scRNA-seq data with droplet-based scRNA-seq data and UMAP dimension were performed using Seurat. Genes that were expressed in more than five cells and cells expressing at least 200 genes were selected for analysis.

Isolation and flow cytometric analysis of thymocytes from mice

Mice were sacrificed using CO2, and thymi were dissected and placed into cold 1× PBS. Adhering non-thymus tissue was carefully cleared off using sharp tweezers under a fluorescence stereomicroscope. The thymus was grinded using glass slides (MATSUNAMI) in 3 mL of RPMI 1640 (Wako). The suspension was collected and centrifuged at 1500 rpm for 5 min. The supernatant was removed and suspended in FACS buffer. After filtering through a 67-µm nylon monofilament mesh, the resulting cell suspension was incubated with anti-mouse CD16/32 in FACS buffer to block nonspecific binding. For flow cytometric analysis of eosinophil and natural killer T cells, cells were stained with antibodies (APCCy7-labeled anti-CD4; BioLegend Cat#100526, PECy7-labeled anti-CD8a; BioLegend Cat#100722, FITC-labeled anti-CD69; BioLegend Cat#104506, BV510-labeled anti-TCRb; BioLegend Cat#109233, Alexa647-labeled anti-Siglec-F; BD Cat#562680, PE-labeled CD1d-PBS-57 tetramer or CD1d-unloaded tetramer) in FACS buffer. Dead cells were excluded by staining with SYTOX™ Blue. For flow cytometric analysis of CCR6–Foxp3+ regulatory T cells, cells were stained with primary antibody (PECy7-labeled anti-CCR6; BioLegend Cat #129816) in FACS buffer, then sequentially incubated with secondary antibodies (APCCy7-labeled anti-CD4, APC-labeled anti-CD8a; BioLegend Cat #100712, FITC-labeled anti-CD25; BioLegend Cat #102006) in FACS buffer and tertiary antibody (PE-labeled anti-Foxp3; eBioscience Cat #12-5773-82) in Permeabilization buffer. Dead cells were excluded by staining with Zombie Aqua. Cells were analyzed using a FACS Aria instrument and a FACS CantoII instrument (BD). Data were analyzed using Flowjo 10.

Statistical analysis

Statistically significant differences between mean values were determined using unpaired t-test or one-way ANOVA followed with multiple comparisons by Tukey’s test in GraphPad Prism (* p≦0.05, ** p≦0.01, *** p≦0.001, **** p≦0.0001). Principle component analysis was performed using edgeR package. P-value correction for differential gene expression analysis in bulk RNA-seq was performed using Baggerley’s test on CLC Genomics Workbench.

Data availability statement

The bulk RNA-Seq data have been deposited at NCBI Sequence Read Archive under series accession number SUB14671957. The droplet-based scRNA-seq data were deposited at NCBI GEO under accession number GSE276809. The well-based scRNA-seq data have been deposited at NCBI Sequence Read Archive under series accession number SUB14826438, SUB14826492, SUB14826498, and SUB14826963.

Ethics statement

The animal study was approved by Institutional Animal Care and Use Committee of RIKEN, Yokohama Branch (2018-075). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

MH: Writing – original draft, Investigation, Validation, Data curation, Formal analysis. HI: Data curation, Formal analysis, Investigation, Validation, Writing – review & editing. MM: Investigation, Validation, Writing – review & editing. MY: Investigation, Writing – review & editing. NH: Writing – review & editing, Data curation. WM: Writing – review & editing, Investigation, Validation. KN: Investigation, Validation, Writing – review & editing. RE: Investigation, Validation, Writing – review & editing. TM: Funding acquisition, Investigation, Supervision, Validation, Writing – review & editing. NA: Funding acquisition, Investigation, Supervision, Writing – review & editing. TA: Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by Grants-in-Aid for Scientific Research from JSPS (21K19391 and 23K27399 to TA, 23K06385 to NA, and 24K18386 to TM) and by CREST from the Japan Science and Technology Agency (JPMJCR2011 to TA). This work was supported by RIKEN Junior Research Associate Program to MH.

Acknowledgments

Computations were performed on supercomputers at the National Institute of Genetics and at ISD, RIKEN (HOKUSAI).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1500908/full#supplementary-material

Supplementary Figure 1 | Gating strategy and PCA analysis for bulk RNA-seq analysis 1 sis of 2 mTEClo fraction 3 (A) Gating strategy for sorting the mTEC^lo fraction expressing low levels of L1CAM and Ly6d. 4 (B) PCA analysis of bulk RNA-seq data for the mTEC^lo fraction.

Supplementary Figure 2 | Quality control data of single cell RNA-seq data.

Supplementary Figure 3 | Volcano plot of marker gene expression in cell clusters from scRNA9 seq analysis of wild-type mice injected with control IgG.

Supplementary Figure 4 | Subclustering analysis of TA-TECs. 12 (A) UMAP projection for subclustering of total TA-TECs. 13 (B) Expression of Aire and Ccl21a in each subcluster of TA-TECs. 14 (C) Separation of UMAP projection into individual data sets. 15 (D) Percentages of proliferating cells (TA-TEC subset) in total Ccl21+ mTECs and total Aire+ 16 mTECs.

Supplementary Figure 5 | Frequency of subclusters in cTEC clusters (A) and expression of 19 RANK and CD40 in TEC clusters from scRNA-seq analysis of wild-type mice injected with 20 control IgG (B).

Supplementary Figure 6 | Expression of Fos, JunB, and Egr1 in subclusters of TEC progenitors 23 in droplet-based scRNA-seq analysis.

Supplementary Figure 7 | UMAP and clustering of integrated scRNA-seq data combining 26 droplet-based scRNA-seq and well-based scRNA-seq data.

Supplementary Figure 8 | A hypothesis for direct and indirect RANK and CD40 signaling in 29 TECs.

References

1. Kadouri N, Nevo S, Goldfarb Y, Abramson J. Thymic epithelial cell heterogeneity: TEC by TEC. Nat Rev Immunol. (2020) 20:239–53. doi: 10.1038/s41577-019-0238-0

PubMed Abstract | Crossref Full Text | Google Scholar

2. Bleul CC, Corbeaux T, Reuter A, Fisch P, Monting JS, Boehm T. Formation of a functional thymus initiated by a postnatal epithelial progenitor cell. Nature. (2006) 441:992–6. doi: 10.1038/nature04850

PubMed Abstract | Crossref Full Text | Google Scholar

3. Rossi SW, Jenkinson WE, Anderson G, Jenkinson EJ. Clonal analysis reveals a common progenitor for thymic cortical and medullary epithelium. Nature. (2006) 441:988–91. doi: 10.1038/nature04813

PubMed Abstract | Crossref Full Text | Google Scholar

4. Sekai M, Hamazaki Y, Minato N. Medullary thymic epithelial stem cells maintain a functional thymus to ensure lifelong central T cell tolerance. Immunity. (2014) 41:753–61. doi: 10.1016/j.immuni.2014.10.011

PubMed Abstract | Crossref Full Text | Google Scholar

5. Lucas B, White AJ, Klein F, Veiga-Villauriz C, Handel A, Bacon A, et al. Embryonic keratin19(+) progenitors generate multiple functionally distinct progeny to maintain epithelial diversity in the adult thymus medulla. Nat Commun. (2023) 14:2066. doi: 10.1038/s41467-023-37589-4

PubMed Abstract | Crossref Full Text | Google Scholar

6. Akiyama N, Takizawa N, Miyauchi M, Yanai H, Tateishi R, Shinzawa M, et al. Identification of embryonic precursor cells that differentiate into thymic epithelial cells expressing autoimmune regulator. J Exp Med. (2016) 213:1441–58. doi: 10.1084/jem.20151780

PubMed Abstract | Crossref Full Text | Google Scholar

7. Ohigashi I, White AJ, Yang MT, Fujimori S, Tanaka Y, Jacques A, et al. Developmental conversion of thymocyte-attracting cells into self-antigen-displaying cells in embryonic thymus medulla epithelium. Elife. (2024) 12:RP92552. doi: 10.7554/eLife.92552

PubMed Abstract | Crossref Full Text | Google Scholar

8. Onder L, Nindl V, Scandella E, Chai Q, Cheng HW, Caviezel-Firner S, et al. Alternative NF-kappa B signaling regulates mTEC differentiation from podoplanin-expressing presursors in the cortico-medullary junction. Eur J Immunol. (2015) 45:2218–31. doi: 10.1002/eji.201545677

PubMed Abstract | Crossref Full Text | Google Scholar

9. Gray D, Abramson J, Benoist C, Mathis D. Proliferative arrest and rapid turnover of thymic epithelial cells expressing Aire. J Exp Med. (2007) 204:2521–8. doi: 10.1084/jem.20070795

PubMed Abstract | Crossref Full Text | Google Scholar

10. Wong K, Lister NL, Barsanti M, Lim JM, Hammett MV, Khong DM, et al. Multilineage potential and self-renewal define an epithelial progenitor cell population in the adult thymus. Cell Rep. (2014) 8:1198–209. doi: 10.1016/j.celrep.2014.07.029

PubMed Abstract | Crossref Full Text | Google Scholar

11. Ulyanchenko S, O’Neill KE, Medley T, Farley AM, Vaidya HJ, Cook AM, et al. Identification of a bipotent epithelial progenitor population in the adult thymus. Cell Rep. (2016) 14:2819–32. doi: 10.1016/j.celrep.2016.02.080

PubMed Abstract | Crossref Full Text | Google Scholar

12. Baran-Gale J, Morgan MD, Maio S, Dhalla F, Calvo-Asensio I, Deadman ME, et al. Ageing compromises mouse thymus function and remodels epithelial cell differentiation. Elife. (2020) 9:e56221. doi: 10.7554/eLife.56221

PubMed Abstract | Crossref Full Text | Google Scholar

13. Dhalla F, Baran-Gale J, Maio S, Chappell L, Hollander GA, Ponting CP. Biologically indeterminate yet ordered promiscuous gene expression in single medullary thymic epithelial cells. EMBO J. (2020) 39:e101828. doi: 10.15252/embj.2019101828

PubMed Abstract | Crossref Full Text | Google Scholar

14. Wells KL, Miller CN, Gschwind AR, Wei W, Phipps JD, Anderson MS, et al. Combined transient ablation and single cell RNA sequencing reveals the development of medullary thymic epithelial cells. Elife. (2020) 9:e60188. doi: 10.7554/eLife.60188

PubMed Abstract | Crossref Full Text | Google Scholar

15. Miyao T, Miyauchi M, Kelly ST, Terooatea TW, Ishikawa T, Oh E, et al. Integrative analysis of scRNA-seq and scATAC-seq revealed transit-amplifying thymic epithelial cells expressing autoimmune regulator. Elife. (2022) 11:e73998. doi: 10.7554/eLife.73998

PubMed Abstract | Crossref Full Text | Google Scholar

16. Michelson DA, Hase K, Kaisho T, Benoist C, Mathis D. Thymic epithelial cells co-opt lineage-defining transcription factors to eliminate autoreactive T cells. Cell. (2022) 185:2542–2558 e2518. doi: 10.1016/j.cell.2022.05.018

PubMed Abstract | Crossref Full Text | Google Scholar

17. Givony T, Leshkowitz D, Del Castillo D, Nevo S, Kadouri N, Dassa B, et al. Thymic mimetic cells function beyond self-tolerance. Nature. (2023) 622:164–72. doi: 10.1038/s41586-023-06512-8

PubMed Abstract | Crossref Full Text | Google Scholar

18. Ragazzini R, Boeing S, Zanieri L, Green M, D’Agostino G, Bartolovic K, et al. Defining the identity and the niches of epithelial stem cells with highly pleiotropic multilineage potency in the human thymus. Dev Cell. (2023) 58:2428–2446 e2429. doi: 10.1016/j.devcel.2023.08.017

PubMed Abstract | Crossref Full Text | Google Scholar

19. Nusser A, Sagar, Swann JB, Krauth B, Diekhoff D, Calderon L, et al. Developmental dynamics of two bipotent thymic epithelial progenitor types. Nature. (2022) 606:165–71. doi: 10.1038/s41586-022-04752-8

PubMed Abstract | Crossref Full Text | Google Scholar

20. Akiyama T, Shimo Y, Yanai H, Qin J, Ohshima D, Maruyama Y, et al. The tumor necrosis factor family receptors RANK and CD40 cooperatively establish the thymic medullary microenvironment and self-tolerance. Immunity. (2008) 29:423–37. doi: 10.1016/j.immuni.2008.06.015

PubMed Abstract | Crossref Full Text | Google Scholar

21. Hikosaka Y, Nitta T, Ohigashi I, Yano K, Ishimaru N, Hayashi Y, et al. The cytokine RANKL produced by positively selected thymocytes fosters medullary thymic epithelial cells that express autoimmune regulator. Immunity. (2008) 29:438–50. doi: 10.1016/j.immuni.2008.06.018

PubMed Abstract | Crossref Full Text | Google Scholar

22. Irla M, Hugues S, Gill J, Nitta T, Hikosaka Y, Williams IR, et al. Autoantigen-specific interactions with CD4+ thymocytes control mature medullary thymic epithelial cell cellularity. Immunity. (2008) 29:451–63. doi: 10.1016/j.immuni.2008.08.007

PubMed Abstract | Crossref Full Text | Google Scholar

23. Mouri Y, Yano M, Shinzawa M, Shimo Y, Hirota F, Nishikawa Y, et al. Lymphotoxin signal promotes thymic organogenesis by eliciting RANK expression in the embryonic thymic stroma. J Immunol. (2011) 186:5047–57. doi: 10.4049/jimmunol.1003533

PubMed Abstract | Crossref Full Text | Google Scholar

24. Lkhagvasuren E, Sakata M, Ohigashi I, Takahama Y. Lymphotoxin beta receptor regulates the development of CCL21-expressing subset of postnatal medullary thymic epithelial cells. J Immunol. (2013) 190:5110–7. doi: 10.4049/jimmunol.1203203

PubMed Abstract | Crossref Full Text | Google Scholar

25. White AJ, Nakamura K, Jenkinson WE, Saini M, Sinclair C, Seddon B, et al. Lymphotoxin signals from positively selected thymocytes regulate the terminal differentiation of medullary thymic epithelial cells. J Immunol. (2010) 185:4769–76. doi: 10.4049/jimmunol.1002151

PubMed Abstract | Crossref Full Text | Google Scholar

26. Khan IS, Mouchess ML, Zhu ML, Conley B, Fasano KJ, Hou Y, et al. Enhancement of an anti-tumor immune response by transient blockade of central T cell tolerance. J Exp Med. (2014) 211:761–8. doi: 10.1084/jem.20131889

PubMed Abstract | Crossref Full Text | Google Scholar

27. Hess J, Angel P, Schorpp-Kistner M. AP-1 subunits: quarrel and harmony among siblings. J Cell Sci. (2004) 117:5965–73. doi: 10.1242/jcs.01589

PubMed Abstract | Crossref Full Text | Google Scholar

28. Ruther U, Muller W, Sumida T, Tokuhisa T, Rajewsky K, Wagner EF. c-fos expression interferes with thymus development in transgenic mice. Cell. (1988) 53:847–56. doi: 10.1016/s0092-8674(88)90289-9

PubMed Abstract | Crossref Full Text | Google Scholar

29. Wagner EF, Eferl R. Fos/AP-1 proteins in bone and the immune system. Immunol Rev. (2005) 208:126–40. doi: 10.1111/j.0105-2896.2005.00332.x

PubMed Abstract | Crossref Full Text | Google Scholar

30. Akiyama T, Maeda S, Yamane S, Ogino K, Kasai M, Kajiura F, et al. Dependence of self-tolerance on TRAF6-directed development of thymic stroma. Science. (2005) 308:248–51. doi: 10.1126/science.1105677

PubMed Abstract | Crossref Full Text | Google Scholar

31. Bornstein C, Nevo S, Giladi A, Kadouri N, Pouzolles M, Gerbe F, et al. Single-cell mapping of the thymic stroma identifies IL-25-producing tuft epithelial cells. Nature. (2018) 559:622–6. doi: 10.1038/s41586-018-0346-1

PubMed Abstract | Crossref Full Text | Google Scholar

32. Baggerly KA, Deng L, Morris JS, Aldaz CM. Differential expression in SAGE: accounting for normal between-library variation. Bioinformatics. (2003) 19:1477–83. doi: 10.1093/bioinformatics/btg173

PubMed Abstract | Crossref Full Text | Google Scholar

33. Tirosh I, Izar B, Prakadan SM, Wadsworth MH 2nd, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. (2016) 352:189–96. doi: 10.1126/science.aad0501

PubMed Abstract | Crossref Full Text | Google Scholar

34. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. (2014) 32:381–6. doi: 10.1038/nbt.2859

PubMed Abstract | Crossref Full Text | Google Scholar

35. Aibar S, Gonzalez-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. (2017) 14:1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | Crossref Full Text | Google Scholar

36. Cosway EJ, White AJ, Parnell SM, Schweighoffer E, Jolin HE, Bacon A, et al. Eosinophils are an essential element of a type 2 immune axis that controls thymus regeneration. Sci Immunol. (2022) 7:eabn3286. doi: 10.1126/sciimmunol.abn3286

PubMed Abstract | Crossref Full Text | Google Scholar

37. Desanti GE, Cowan JE, Baik S, Parnell SM, White AJ, Penninger JM, et al. Developmentally regulated availability of RANKL and CD40 ligand reveals distinct mechanisms of fetal and adult cross-talk in the thymus medulla. J Immunol. (2012) 189:5519–26. doi: 10.4049/jimmunol.1201815

PubMed Abstract | Crossref Full Text | Google Scholar

38. LaFlam TN, Seumois G, Miller CN, Lwin W, Fasano KJ, Waterfield M, et al. Identification of a novel cis-regulatory element essential for immune tolerance. J Exp Med. (2015) 212:1993–2002. doi: 10.1084/jem.20151069

PubMed Abstract | Crossref Full Text | Google Scholar

39. Ashby KM, Voboril M, Salgado OC, Lee ST, Martinez RJ, O’Connor CH, et al. Sterile production of interferons in the thymus affects T cell repertoire selection. Sci Immunol. (2024) 9:eadp1139. doi: 10.1126/sciimmunol.adp1139

PubMed Abstract | Crossref Full Text | Google Scholar

40. Otero DC, Baker DP, David M. IRF7-dependent IFN-beta production in response to RANKL promotes medullary thymic epithelial cell development. J Immunol. (2013) 190:3289–98. doi: 10.4049/jimmunol.1203086

PubMed Abstract | Crossref Full Text | Google Scholar

41. Ohigashi I, Zuklys S, Sakata M, Mayer CE, Hamazaki Y, Minato N, et al. Adult thymic medullary epithelium is maintained and regenerated by lineage-restricted cells rather than bipotent progenitors. Cell Rep. (2015) 13:1432–43. doi: 10.1016/j.celrep.2015.10.012

PubMed Abstract | Crossref Full Text | Google Scholar

42. Furuya Y, Mori K, Ninomiya T, Tomimori Y, Tanaka S, Takahashi N, et al. Increased bone mass in mice after single injection of anti-receptor activator of nuclear factor-kappaB ligand-neutralizing antibody: evidence for bone anabolic effect of parathyroid hormone in mice with few osteoclasts. J Biol Chem. (2011) 286:37023–31. doi: 10.1074/jbc.M111.246280

PubMed Abstract | Crossref Full Text | Google Scholar

43. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessorBioinfomatics. J Biol Chem. (2011) 34:i884–i890. doi: 10.1093/bioinformatics/bty560

PubMed Abstract | Crossref Full Text | Google Scholar