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.

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.

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.

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.

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 ).

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.

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.

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.

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.

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 ).

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 ).

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.

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.

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.

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.

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).

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.

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.

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.

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.