Integrative analysis of scRNA-seq and scATAC-seq revealed transit-amplifying thymic epithelial cells expressing autoimmune regulator

  1. Takahisa Miyao
  2. Maki Miyauchi
  3. S Thomas Kelly
  4. Tommy W Terooatea
  5. Tatsuya Ishikawa
  6. Eugene Oh
  7. Sotaro Hirai
  8. Kenta Horie
  9. Yuki Takakura
  10. Houko Ohki
  11. Mio Hayama
  12. Yuya Maruyama
  13. Takao Seki
  14. Hiroto Ishii
  15. Haruka Yabukami
  16. Masaki Yoshida
  17. Azusa Inoue
  18. Asako Sakaue-Sawano
  19. Atsushi Miyawaki
  20. Masafumi Muratani
  21. Aki Minoda
  22. Nobuko Akiyama  Is a corresponding author
  23. Taishin Akiyama  Is a corresponding author
  1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Japan
  2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Japan
  3. Laboratory for Cellular Epigenomics, RIKEN Center for Integrative Medical Sciences, Japan
  4. YCI Laboratory for Immunological Transcriptomics, RIKEN Center for Integrative Medical Sciences, Japan
  5. YCI Laboratory for Metabolic Epigenetics, RIKEN Center for Integrative Medical Sciences, Japan
  6. Laboratory for Cell Function Dynamics, RIKEN Center for Brain Science, Japan
  7. Transborder Medical Research Center, and Department of Genome Biology, Faculty of Medicine, University of Tsukuba, Japan

Abstract

Medullary thymic epithelial cells (mTECs) are critical for self-tolerance induction in T cells via promiscuous expression of tissue-specific antigens (TSAs), which are controlled by the transcriptional regulator, AIRE. Whereas AIRE-expressing (Aire+) mTECs undergo constant turnover in the adult thymus, mechanisms underlying differentiation of postnatal mTECs remain to be discovered. Integrative analysis of single-cell assays for transposase-accessible chromatin (scATAC-seq) and single-cell RNA sequencing (scRNA-seq) suggested the presence of proliferating mTECs with a specific chromatin structure, which express high levels of Aire and co-stimulatory molecules, CD80 (Aire+CD80hi). Proliferating Aire+CD80hi mTECs detected using Fucci technology express a minimal number of Aire-dependent TSAs and are converted into quiescent Aire+CD80hi mTECs expressing high levels of TSAs after a transit amplification. These data provide evidence for the existence of transit-amplifying Aire+mTEC precursors during the Aire+mTEC differentiation process of the postnatal thymus.

Editor's evaluation

This report shows by scRNAs-seq and scATAC-seq the presence of a population of proliferating medullary thymic epithelial cells (mTECs) with a specific chromatin structure and high expression of Aire and CD80. Such Aire-expressing transit-amplifying mTECs may play a key role in establishing immunological self-tolerance.

https://doi.org/10.7554/eLife.73998.sa0

Introduction

Medullary thymic epithelial cells (mTECs) are essential for induction of T cell self-tolerance in the thymus (Abramson and Anderson, 2017; Inglesfield et al., 2019). mTECs ectopically express thousands of tissue-specific antigens (TSAs), and this expression is regulated by transcription factors, AIRE and FEZF2 (Anderson et al., 2002, Takaba et al., 2015). TSAs are directly or indirectly presented to developing T cells, and T cells that recognize TSAs with high affinity undergo apoptosis or are converted into regulatory T cells, thereby suppressing the onset of autoimmune diseases (Abramson and Anderson, 2017; Inglesfield et al., 2019).

Several studies have suggested processes and underlying mechanisms of mTEC differentiation during thymic organogenesis (Abramson and Anderson, 2017; Inglesfield et al., 2019; Rossi et al., 2007; Akiyama et al., 2016; Akiyama et al., 2005; Akiyama et al., 2008; Hikosaka et al., 2008; Mouri et al., 2011; Kajiura et al., 2004). In addition, some previous studies have reported that mTEC turnover is homeostatic in the adult thymus, with a duration of approximately 2 weeks (Gäbler et al., 2007; Gray et al., 2007; Gray et al., 2006). Notably, however, cellular mechanisms underlying maintenance of adult mTECs remain unclear. mTEC subpopulations are largely classified based on their expression of cell surface markers (mainly CD80 and MHC class II) and Aire in the adult thymus (Abramson and Anderson, 2017). CD80lo and Aire-negative (Aire-) mTECs (mTEClo) are thought to be immature, and they differentiate into CD80hi Aire-expressing (Aire+) mTECs that are reportedly post-mitotic (Gray et al., 2007). Aire+ mTECs are further converted into Aire-negative mTECs (post-Aire mTECs) (Metzger et al., 2013; Michel et al., 2017; Nishikawa et al., 2014; Wang et al., 2012; White et al., 2010). Moreover, a previous study suggested that mTECs might be differentiated from stage-specific embryonic antigen-1+ (SSEA-1) claudine3/4+ mTEC stem cells (Sekai et al., 2014). These views are primarily based on fate mapping studies involving transfer and reaggregation of sorted cell populations with fetal thymus (Rossi et al., 2007; Gray et al., 2007; Sekai et al., 2014) and on experiments employing genetic marking (Metzger et al., 2013; Nishikawa et al., 2014).

Single-cell RNA sequencing (scRNA-seq) technology has yielded new insights into cell diversity and differentiation in various tissues. In TEC biology, previous scRNA-seq studies revealed the stochastic nature of TSA expression in mTECs (Sansom et al., 2014; Meredith et al., 2015) and high heterogeneity of TECs in mice (Bornstein et al., 2018; Miller et al., 2018; Dhalla et al., 2020; Baran-Gale et al., 2020). Bornstein et al. showed that mTECs in the postnatal thymus are separated into four subsets, mTEC I to IV (Bornstein et al., 2018). In addition to the classical mTEClo (mTEC I), Aire+ mTEC (mTEC II), and post-Aire mTEC (mTEC III) types, a tuft-like mTEC subset (mTEC IV) was identified (Bornstein et al., 2018; Miller et al., 2018). Subsequent scRNA-seq studies suggested further heterogeneity of TECs, such as cilium TECs (Dhalla et al., 2020), GP2+ TECs (Dhalla et al., 2020), intertypical TECs (Baran-Gale et al., 2020), neural TECs (Baran-Gale et al., 2020), and structural TECs (Baran-Gale et al., 2020), according to specific gene expression profiles. However, it has not yet been clarified whether this heterogeneity identified from gene expression profiles is correlated with differences in chromatin structure.

In general, transit-amplifying cells (TACs) are a proliferating cell population linking stem cells and differentiated cells (Lajtha, 1979). TACs are short-lived and undergo differentiation after a few cell divisions. To date, the presence of TACs has been confirmed in some tissues such as intestines (Clevers, 2013), hair follicles (Hsu et al., 2014), and neurons (Lui et al., 2011). Previous analyses of scRNA-seq data of murine adult TECs revealed a cell cluster expressing an abundance of cell cycle-regulated genes, which implies the presence of TACs for TECs (TA-TECs) (Dhalla et al., 2020; Wells et al., 2020). Computational trajectory analysis of scRNA-seq data suggested that this population might give rise to Aire-expressing mTECs (Dhalla et al., 2020; Baran-Gale et al., 2020). Intriguingly, another trajectory study predicted that this cell cluster might differentiate into Aire-expressing mTECs and an mTEC population expressing CCL21a (Wells et al., 2020). However, because the TA-TEC candidate has not been isolated and specific marker genes of TA-TECs have not been reported, exact properties of TA-TECs, in addition to their cellular fates, remain to be clarified.

In this study, droplet-based scRNA-seq and single-cell assays for transposase-accessible chromatin sequencing (scATAC-seq) of murine TECs were performed to characterize TEC heterogeneity and differentiation dynamics. Integrative analysis of these data showed that Aire+ mTECs are separated into at least two clusters with different gene expression profiles and chromatin accessibilities. One of these Aire+ mTEC clusters exhibited high expression of cell cycle-related genes, which accords with a previously proposed TAC population of mTECs (Dhalla et al., 2020; Wells et al., 2020). By using the Fucci technology (Mort et al., 2014), proliferating mTECs expressing Aire and maturation marker CD80 were isolated as TA-TEC candidates. This proliferating Aire+ CD80hi mTEC subpopulation showed minimal expression of TSAs regulated by AIRE, in contrast to quiescent Aire+ CD80hi mTECs. Moreover, in vivo BrdU pulse-labeling, and in vitro reaggregated thymic organ culture suggested that proliferating Aire+ CD80hi mTECs are short-lived and that they differentiate into quiescent Aire+ CD80hi mTECs, post-Aire mTECs, and tuft-like mTECs. Consequently, these data strongly suggest that proliferating Aire+ CD80hi mTECs are TACs for mTECs expressing TSAs.

Results

Droplet-based scATAC-seq reveals heterogeneity of TEC chromatin structure

Given that chromatin structures can be changed during cell differentiation, scATAC-seq analysis of TECs may offer some insights into TEC heterogeneity and differentiation dynamics. Droplet-based scATAC-seq analysis was carried out with EpCAM+ CD45 cells that were sorted and pooled from thymi of two mice, 4 weeks of age. scATAC-seq analysis was repeated twice and integrated with removal of batch effects via a combination of the Signac R package (Stuart et al., 2021) and the Harmony algorithm (Korsunsky et al., 2019). Unsupervised graph-based clustering and dimensional reduction via uniform manifold approximation and production (UMAP) revealed 11-cell clusters from 15,255 cells (7884 for Experiment #1 and 7371 for Experiment #2) (Figure 1A). Chromatin accessibility of previously known TEC marker genes (gene coordinates including their 2 kbp upstream region) was analyzed. Clusters 0, 4, 5, 6, 8, 9, and 11 contained relatively higher numbers of cells having the open chromatin structure of the Cd80 gene, a maturation marker of TECs (Figure 1B and C). Among these clusters, the cis-regulatory element of the Aire gene (LaFlam et al., 2015) (about 2 kbp upstream from the transcriptional start site) is opened in clusters 0 and 4 (Figure 1D), suggesting that these clusters may be concordant with Aire-expressing mTECs (Aire+ mTECs, also referred to as mTEC II Bornstein et al., 2018). In contrast, the cis element of Aire genes is closed in clusters 5, 6, 8, 9, and 11 (Figure 1D), suggesting that these clusters may correspond to post-Aire mTECs and other Aire-negative mature mTECs (Bornstein et al., 2018). Because the Irga2 gene (also called Lrmp) region is accessible in cluster 6 (Figure 1B and Figure 1—figure supplement 1), this cluster may be equivalent to tuft-like mTECs (mTEC IV) (Bornstein et al., 2018; Miller et al., 2018). CD80 and Aire gene regions in clusters 1, 2, and 3 are relatively closed, whereas the mTEC marker Tnfrsf11a is relatively accessible (Figure 1B and C, and Figure 1—figure supplement 1). Therefore, these clusters should be equivalent to mTECs expressing low levels of CD80 and Aire (mTEClo). Cluster 7 should be cTECs, because a cTEC marker Psmb11 gene region is opened (Figure 1B and Figure 1—figure supplement 1). Finally, cluster 10 was deemed thymocyte contamination because the Rag1 gene was opened (Figure 1—figure supplement 1). Comparison between two independent scATAC-seq experiments suggested relatively high reproducibility (Figure 1—figure supplement 2).

Figure 1 with 3 supplements see all
Droplet-based single-cell assays for transposase-accessible chromatin sequencing (scATAC-seq) analysis of thymic epithelial cells (TECs) in 4-week mice.

(A) Uniform manifold approximation and production (UMAP) plot of scATAC-seq data from TECs (EpCAM+ CD45 TER119) from 4-week mice. Cell clusters are separated by colors and numbers in the plot. The two datasets were integrated using the Seurat package. The graph on the right shows percentages of each cluster in the total number of cells detected (15,255 cells). (B) Chromatin accessibility of typical marker genes of TECs. Accessibility in each gene region is represented in red. (C) Violin plot depicting chromatin accessibility in Aire and Cd80 gene regions in each cluster. (D) Pseudo-bulk accessibility tracks of the Aire gene region in each cluster (upper panels) and frequency of sequenced fragments within the Aire gene region of individual cells in cluster 0, 2, and 4 (lower panels).

We next sought to correlate scATAC data with TEC scRNA-seq data. Droplet-based scRNA-seq analysis of EpCAM+CD45 cells from age- and gender-matched mice (4-week female mice) was performed. scRNA-seq analysis was repeated twice and integrated with the removal of batch effects via the Seurat package (Butler et al., 2018). Analysis of integrated data (total 11,475 cells) revealed 18-cell clusters in the UMAP dimension (Figure 2A, Figure 2—figure supplement 1, and Supplementary file 1). Comparison between two independent experiments suggested high reproducibility of the scRNA-seq analysis (Figure 2—figure supplement 2).

Figure 2 with 6 supplements see all
Droplet-based single-cell RNA sequencing (scRNA-seq) analysis of thymic epithelial cells (TECs) in 4-week mice.

(A) Uniform manifold approximation and production (UMAP) plot of scRNA-seq data from TECs (EpCAM+ CD45 TER119) from 4-week mice. Cell clusters (R0 to R17) are indicated by colors and numbers in the plot. The graph on the right shows the percentages of each cluster in the total number of cells detected (11,792 cells). (B) Violin plots depicting expression levels of typical TEC marker genes in each cluster.

These TEC clusters were assigned according to expression of TEC marker genes (Figure 2 and Figure 2—figure supplement 3). Clusters R0, R1, R3, and R9 showed high expression of Aire (Figure 2B), suggesting that these clusters are equivalent to Aire+ mTECs (also referred to as mTECs II). Clusters R2, R4, and R5 include cells showing relatively higher levels of Itga6 and Ccl21a expression and a very low level of Aire expression (Figure 2B), corresponding to mTEC I (Bornstein et al., 2018), CCL21-expressing mTECs (Lucas et al., 2020), and possibly intertypical TECs (Baran-Gale et al., 2020). Cluster R6 expresses Irga2 (Figure 2B) and should contain tuft-like mTECs (mTEC IV) (Bornstein et al., 2018). Clusters R7 and R10 were marked with Krt10 and Pigr genes, respectively (Figure 2B). Accordingly, these clusters should be categorized as post-Aire mTECs (also referred to as mTECs III Bornstein et al., 2018). Cluster R13 showed high expression of chemokines, Ccl6 and Gp2 (Figure 2B and Figure 2—figure supplement 3), which is concordant with Gp2+ TECs, as reported recently (Dhalla et al., 2020). Clusters R8 and R11 exhibited high expression of typical cTEC marker genes, Psmb11 and Dll4 (Figure 2B and Figure 2—figure supplement 3), and should be equivalent to cTECs. Given that thymocyte genes are highly expressed, cluster R11 was most likely thymic nurse cells enclosing thymocytes (Nakagawa et al., 2012). Cluster R12 showed relatively high expression of Pdpn (Figure 2—figure supplement 3), which may comprise junctional TECs (Onder et al., 2015). Cluster R14 was considered thymocyte contamination because thymocyte markers, but not TEC markers, were detected. Cluster R15 apparently corresponds to structural TECs, reported recently, because of their expression of Car8 and Cd177 (Baran-Gale et al., 2020; Figure 2—figure supplement 3). Cells in cluster R16 highly express Tppp3 and Fam183b (Figure 2—figure supplement 3). Since these genes are expressed in ciliated cells (Orosz and Ovádi, 2008; Beckers et al., 2018), this cluster may be equivalent to ciliated columnar TECs associated with thymic cystic structure (Dhalla et al., 2020; Khosla and Ovalle, 1986; Park et al., 2020). We failed to assign cluster R17, which may be contaminated with endothelial cells or macrophages, because they express Ly6c1 and Aqp1, but low levels of Epcam (Figure 2—figure supplement 3). Correlation of our scRNA-seq data with reported scRNA-seq data (Bornstein et al., 2018; Dhalla et al., 2020; Baran-Gale et al., 2020; Wells et al., 2020) was investigated by integrating datasets (Figure 2—figure supplement 4) or checking expression of differentially expressed genes in clusters of other datasets (Figure 2—figure supplements 5 and 6). Overall, our data and assignments were reasonably correlated with previous scRNA-seq data analyses.

In next, we bioinformatically integrated the scRNA-seq data with scATAC-seq data. Gene expression, predicted from accessible chromatin regions of scATAC-seq data, was correlated with scRNA-seq data using the Signac R package (Figure 3 and Figure 3—figure supplement 1). As described, clusters 0 and 4 in scATAC-seq analysis contain cells with the accessible cis-regulatory element of the Aire gene (Figure 1D). Consistently, cluster 0 in scATAC-seq was mostly transferred to clusters R0 (56.8 %) and R3 (12.9%) in scRNA-seq analysis (Figure 3B and C, and Supplementary file 2), which were assigned as Aire+ mTECs (Figure 2). Cells transferred to R0 and R3 appear to be separately embedded in cluster 0 in the UMAP dimension, implying that these Aire+ mTEC subsets have slightly different chromatin structures. Cluster 4 was mostly transferred to cluster R1 (57.1%) (Figure 3B and C), also designated as Aire+ mTECs, and in part R2 (23.5%), which belongs to mTEC I. Interestingly, cells transferred to cluster R9 are embedded around the junction between cluster 0 and 4 (Figure 3B), suggesting that cluster R9 may be a transitional stage between R1 and R0. Clusters 1, 2, and 3 are closely embedded in the UMAP dimension and principally assigned to clusters R2, R4, and R5 (Figure 3B and C), suggesting that these clusters are concordant with mTEC I or intertypical TECs assigned in the scRNA-seq data. Cluster 5 mainly contains cells transferred to cluster R3 (56.0%) and R10 (27.4%) (Figure 3B), which were assigned as late-Aire mTECs (mTEC IIc) and post-Aire mTECs (mTEC III), respectively. As expected, cluster 6 with an open Irga2 gene was transferred to cluster R6, a tuft-like mTEC subset (mTEC IV). Cluster 7 was transferred to cluster R8 and R12, assigned as cTECs and jTECs, respectively. Cluster 9 was assigned as cluster R7, which is Krt10+ mTEC III subset (Figure 3C). Cluster 8 contains clusters R15 (35.4%) and R16 (12.2%), which express markers of structural TECs and cilia TECs, respectively (Figure 3C and Figure 3—figure supplement 1), in addition to R7 and R10 assigned as mTEC III. Cluster 11 was concordant with R13, which was Gp2+ TECs (Figure 3B and C). Finally, cluster 10 was transferred to clusters R11 and R14, which are assigned as T cells and nurse TECs (Figure 3C and Figure 3—figure supplement 1). Although a few cells were transferred to R17 in scRNA-seq, these cells did not form cluster in this analysis. These assignments were practically the same in the two scATAC datasets (Figure 3—figure supplement 2). Overall, the integration analysis suggested that TEC heterogeneity predicted from scRNA-seq may be ascribed to differences in chromatin structure.

Figure 3 with 4 supplements see all
Integrative analysis of single-cell assays for transposase-accessible chromatin sequencing (scATAC-seq) data and single-cell RNA sequencing (scRNA-seq) data of thymic epithelial cells (TECs).

(A) Gene expression was predicted from scATAC-seq data using Signac. Individual cells in the cluster from scATAC data (clusters 0 to 11) were assigned and transferred to the uniform manifold approximation and production (UMAP) plot of scRNA-seq cluster (R0 to R17). (B) Correlation between clusters derived from scATAC-seq and scRNA-seq datasets of TECs. Cell types were annotated in scATAC dataset of TECs by transferring clusters from an scRNA-seq dataset. (C) Ratios of cells assigned to each scRNA-seq cluster in each scATAC cluster.

Aire-positive mTECs are divided into two subsets having distinct chromatin structures

Previous scRNA-seq studies proposed the existence of a TEC population showing high expression of cell cycle-regulated genes (Dhalla et al., 2020; Baran-Gale et al., 2020; Wells et al., 2020). In our scRNA-seq data, cluster R1 (mTEC IIb) appears equivalent to such a TEC subset, expressing cell cycle-related genes (Figure 2—figure supplement 3B). Integrative analysis of scRNA-seq and scATAC-seq suggested that cells in cluster 4 in scATAC-seq were mainly transferred to cluster R1 (Figure 3B). Although both clusters 4 and 0 have the accessible enhancer element of the Aire gene (Figure 1), 269 genomic regions were significantly opened, and 147 regions were closed in cluster 4, in contrast to cluster 0 (Supplementary file 3 and Figure 3—figure supplement 3). Thus, it is likely that the Aire+ mTECs are divided into two subsets based on expression of proliferation markers and chromatin accessibility.

RNA velocity, which recapitulates differentiation dynamics by comparing unspliced and spliced RNA in scRNA-seq data (La Manno et al., 2018), predicted that cluster R1 may differentiate into other Aire+ mTECs (clusters R0, R3, and R9) (Figure 3—figure supplement 4A), which is consistent previous analyses (Dhalla et al., 2020). Moreover, trajectory analysis of scATAC-seq data using Monocle3 suggested a possible transition between cluster 4 and cluster 0 (Figure 3—figure supplement 4B). Overall, integrative analysis of scATAC-seq and scRNA-seq data imply that the Aire+ mTEC subset expressing cell cycle-related genes (cluster 4 in scATAC-seq and cluster R1 in scRNA-seq) may be equivalent to transiently amplifying cells (TA cells) with a distinct chromatin structure. Although a previous study suggested a trajectory of proliferating TEC clusters to mTEC I31, RNA velocity of our scRNA-seq data did not clearly recapitulate it (Figure 3—figure supplement 4A).

A proliferating TEC cluster is sub-divided into an Aire-expressing subcluster and the Aire-negative Ccl21ahigh subcluster

Subclustering of cluster R1 showed its separation into seven subclusters (R1A to R1G in Figure 4A and B). Clusters R1A, R1B, R1C, R1D, and R1E showed expression of Aire and Cd80 (Figure 4C and D). In contrast, Ccl21a, but not Aire, is highly expressed in clusters R1F and R1G (Figure 4C and D), which may be consistent with a previous study (Wells et al., 2020). Interestingly, cell cycle scoring analysis suggested that R1A, R1B, and R1F contain mainly G2M phase cells. In contrast, R1D, R1E, and R1G contain S phase cells and R1C contains both G2M and S phase cells (Figure 4E and F). Thus, the Aire+ subclusters and the Aire-negative Ccl21ahigh subclusters may be further divided by cell cycle phase. These data suggested that proliferating TECs consist of two distinct subsets, distinguished by expression levels of Aire and Ccl21a. Correlation between subclusters of R1 (R1A to R1G) and scATAC-seq clusters was analyzed (Figure 4—figure supplement 1). Cells in cluster 4 were transferred to both the Aire+ (R1A to R1E) and the Aire-negative Ccl21ahigh (R1F and R1G) subclusters (Supplementary file 4), implying that subclusters in R1 may have similar chromatin accessibility, although this point should be explored further in the future.

Figure 4 with 1 supplement see all
Subcluster analysis of the thymic epithelial cell (TEC) subset expressing a high level of cell cycle-related genes.

(A) Uniform manifold approximation and production (UMAP) plot of single-cell RNA sequencing (scRNA-seq) data of each subcluster (R1A to R1G) in R1. Cell subclusters (R1A to R1G) are separated by colors and numbers in the plot. The graph on the right shows the percentages of each cluster in the parent R1 cluster. (B) Heatmap of the top five genes selectively expressed in each subcluster. Yellow color indicates high expression. (C) Expression levels of Aire and Ccl21a in the subcluster are exhibited as violin plots. (D) Expression levels of Aire, Ccl21a, and Cd80 in the subcluster are shown in dot plots. (E) Expression levels of marker genes for G2/M phase (upper, Mki67 and Hmgb2) and S phase (lower, Tyms and Slbp) in the subcluster are exhibited as violin plots. (F) Percentage of cells predicted as each cell cycle (G1, G2M, and S phases) in the subclusters.

A proliferative cell subset is present in Aire+ mTECs

TA cells were generally defined as a proliferative, short-lived cell subset generated from progenitor or stem cells and differentiating into mature quiescent cells (Lajtha, 1979; Zhang and Hsu, 2017). Analysis of scRNA-seq and scATAC-seq data suggested that proliferating Aire+ mTECs may be equivalent to TA cells with distinct chromatin accessibility. To search for evidence supporting the presence of TA cells of mTECs (TA-TECs), we first sought to isolate the proliferating Aire+ CD80hi mTEC subset as candidate TA-TECs. Fucci2a mice, in which cell cycle progression can be monitored with mCherry (G1 and G0 phases) and Venus (G2, M, and S phases) fluorescence, were used to isolate such proliferating cells (Figure 5A; Mort et al., 2014; Lazzeri et al., 2018; Wong et al., 2018; Antonica et al., 2019), and were crossed with Aire-GFP-reporter mice to facilitate detection of Aire expression (Yano et al., 2008). Flow cytometric analysis indicated that Venus+ cells are present among mTECs expressing high levels of CD80 (mTEChi) (Figure 5—figure supplement 1A), although the expression level of CD80 might be slightly lower than that of CD80+ mCherryhi mTECs (Figure 5—figure supplement 1B). Moreover, these Venus+ mTEChi cells expressed GFP (Figure 5B), indicating expression of AIRE. Thus, these data revealed the presence of dividing cells in the Aire+ CD80hi mTEC fraction. The fluorescence intensity of Aire-GFP in Venus+ CD80hi mTECs showed a broad peak and was slightly lower than that of Venus mTEChi cells, which may be due to the relatively lower expression of Aire in Venus+ CD80hi mTECs. However, the compensation between GFP and Venus proteins, which have very close fluorescence spectra, hampered an exact comparison of Aire expression levels between Venus+ mTEChi cells and Venus mTEChi cells. We next confirmed Aire protein expression in proliferating mature mTECs. Immunostaining with an anti-Aire-antibody revealed the presence of Aire protein localized in the nucleus of sorted Venus+ CD80hi mTECs (Figure 5C). Immunostaining of the thymic section from Foxn1-specific Fucci2a mice revealed that Venus+ cells are localized in the medulla, and some of the Aire+ mTECs were Venus+ (Figure 5D and Figure 5—figure supplement 2). Taken together, these data confirm the presence of proliferating Aire+ CD80hi mTECs in the thymic medulla.

Figure 5 with 3 supplements see all
A highly proliferative subset of Aire+ CD80hi medullary thymic epithelial cells (mTECs).

(A) Schematic depiction of cell cycles and Fucci fluorescence. (B) Flow cytometric analysis of TECs from Fucci2a mice crossed with Aire-GFP-reporter mice. The gating strategy is shown. Intensities of GFP to monitor Aire expression in each subset (Venus+ CD80hi mTEC, Venus CD80hi mTEC, and CD80lo mTEClo) are shown in the right panels. Left, Airegfp/+:: Fucci2a; right, control::Fucci2a. Typical figures of three independent experiments are exhibited. (C) Immunostaining of a sorted Venus+ CD80hi mTEC subset via anti-Aire antibody and DAPI (nucleus staining). Typical panels of three independent experiments are exhibited. Scale bars, 10 μm. (D) Immunostaining of thymic sections from Fucci2a mice with anti-Aire and anti-keratin-5 (Krt5) antibodies. Typical panels of three independent experiments are exhibited. Scale bars, 10 μm. (E) Scatter plots of RNA sequencing data from Venus+ CD80hi mTEC and Venus CD80hi mTEC subsets. The left panel shows a plot of all detected genes and the right panel shows tissue-specific antigen (TSA) genes detected. N = 3. (F) Atypical RNA sequencing tracks of Aire, typical Aire-dependent TSA genes (Ins and Sst), Fezf2, and Top2a (a marker of G2/M phase). (G) Scatter plots and volcano plots of RNA sequencing data from Venus+ CD80hi mTEC and Venus CD80hi mTEC subsets. Upper panels show Aire-dependent TSAs, lower panels show Aire-independent TSAs. Red dots in volcano plots indicate genes for which expression differed significantly (twofold change and FDR p < 0.05) in Venus+ and Venus CD80hi mTEC subsets. Numbers of differentially expressed genes are shown in the panels. N = 3. Y axis is log10 of FDR p-value.

Proliferating Aire+ mTECs express low levels of Aire-dependent TSAs

We next addressed whether the proliferating Aire+CD80hi mTECs subset has a molecular signature distinct from that of quiescent Aire+CD80hi mTECs. RNA-seq analysis of sorted cells from Fucci mice suggested that Venus+ CD80hi mTECs and Venus CD80hi mTECs subsets have considerably different gene expression profiles (Figure 5E). As expected, gene ontology analysis confirmed enrichment of cell cycle-related genes in Venus+ CD80hi mTECs compared with Venus CD80hi mTECs (Supplementary file 5). Notably, the Venus+ CD80hi mTEC subset expressed lower levels of Aire-dependent TSAs (Supplementary file 6) than the Venus CD80hi mTECs subset (Figure 5F and G and Figure 5—figure supplement 3). Whereas expression of Aire-independent TSAs was also low in the Venus+ CD80hi mTEC subset, the difference was smaller than in Aire-dependent TSAs (Figure 5G and Figure 5—figure supplement 3). These data suggested that proliferating Aire+CD80hi mTECs are phenotypically immature, compared to quiescent Aire+CD80hi mTECs.

Proliferating Aire+ mTECs are precursors of mature mTECs

Because TA cells are defined as short-lived cells differentiating into mature cells (Lajtha, 1979), we next addressed this issue regarding the proliferating Aire+ CD80hi mTECs. In vivo pulse labeling of TECs with 5-bromo-2’-deoxyuridine (BrdU) was performed. Because mCherryhi cells and mCherrylo were generally in G0 and G1 stages of the cell cycle, respectively (Tomura et al., 2013), each fraction in CD80hi mTECs was sorted separately after i.p. administration of BrdU, and thereafter stained with anti-BrdU antibody (Figure 6A). This procedure was necessary because mCherry fluorescence is lost after BrdU staining. Flow cytometric analysis showed that approximately 35% of mCherrylo CD80hi mTECs (hereafter referred to as mCherrylo) were labeled 12 hr (day 0.5) after the BrdU injection (Figure 6B). In contrast, about 3% of mCherryhi CD80hi mTECs (referred to as mCherryhi) were BrdU-positive (Figure 6B). Thus, as expected, cell cycle progression of mCherrylo is much faster than mCherryhi. Importantly, cell number and the ratio of BrdU-positive cells in the mCherrylo fraction were significantly decreased 3.5 days after the BrdU injection (Figure 6B and C). On the other hand, the frequency of BrdU-positive cells in mCherryhi was increased by day 3.5, and plateaued from day 3.5 to day 6.5 (Figure 6B and C). Notably, mean fluorescence intensity (MFI) of BrdU staining in mCherryhi at day 3.5 was about 50% lower than that in mCherrylo at day 0.5 (Figure 6D), suggesting that mCherryhi cells at day 3.5 were generated after cell division. Overall, these data suggest that mCherrylo cells are transiently proliferating, and after cell division, they are converted to mCherryhi, having low proliferative activity.

Fate mapping study with in vivo BrdU pulse-labeling of Fucci thymic epithelial cells (TECs).

(A) Schematic procedure of in vivo BrdU pulse labeling of Fucci mice, and analysis of BrdU staining in mCherryhiCD80hi and mCherryloCD80hi medullary TECs (mTECs) by flow cytometiric analysis. BrdU staining was done after sorting each cell fraction. (B) Typical flow cytometric profile of BrdU staining in mCherryloCD80hi mTECs (upper panels) and mCherryhiCD80hi mTECs (lower panels) at days 0.5, 3.5, and 6.5 after the BrdU injection. Data for the ratio of BrdU+ cells in each mTEC fraction are summarized in right-hand figures. N = 7 for 0.5 day after the BrdU injection, N = 3 for 3.5 and 6.5 days after the injection. Two-tailed Student’s t-tests. **p < 0.01 and *p < 0.05. NS, not significant (p > 0.05). p = 1.5 × 10–3 for the upper figure and p = 0.033 for the lower figure. Original data were shown in Figure 6—source data 1. (C) Cell number of BrdU+mCherryloCD80hi mTECs and BrdU+mCherryhiCD80hi mTECs at days 0.5, 3.5, and 6.5 after the BrdU injection. Two-tailed Student’s t-tests. **p < 0.01. NS, not significant (p > 0.05). p = 4.3 × 10–3 for the left figure and p = 5.1 × 10–3 for the right figure. Original data were shown in Figure 6—source data 1. (D) Mean fluorescence intensity (MFI) of BrdU staining in mCherryloCD80hi at day 0.5 and mCherryhiCD80hi at days 3.5 and 6.5. MFIs of other time points were difficult to evaluate because of very low cell numbers. Two-tailed Student’s t-tests. *p = 0.015 and **p = 6.5 × 10–3. NS, not significant (p > 0.05). Original data were shown in Figure 6—source data 1.

To verify that mCherrylo cells are precursors of mCherryhi, we performed an in vitro reaggregated thymic organ culture (RTOC) experiment (Figure 7A). The mCherrylo fraction was sorted (Figure 7—figure supplement 1A) and reaggregated with wild type embryonic thymic cells. After 5 days of culture, mCherryhi was detected in RTOC (Figure 7A). Because Venus+mCherrylo cells were practically absent in RTOC (Figure 7—figure supplement 1B), surviving mCherrylo cells were mostly converted into mCherryhi in RTOC. The possibility that mCherryhi contaminating cells during cell sorting survived in RTOC was ruled out by control reaggregation experiments using only mCherryhi (Figure 7—figure supplement 1C) in addition to the higher expression of pro-apoptotic genes (Liberzon et al., 2015) in mCherryhi than mCherrylo (Figure 7—figure supplement 1D). Interestingly, reaggregation with allogenic fetal thymus (Balb/cA background) was not sufficient for conversion to mCherryhi (Figure 7—figure supplement 1E), implying that high affinity interaction between TCR and MHC contributes to survival and maintenance of mCherrylo TECs, as described previously (Irla et al., 2008). Next, we sorted mCherryhi cells in the RTOC (referred to as mCherryhi-RTOC) in addition to mCherrylo and mCherryhi from the Fucci thymus, and analyzed gene expression by RNA-seq. As expected, the mCherrylo fraction expressed a lower level of Aire-dependent TSAs, compared to mCherryhi (Figure 7B and Figure 7—figure supplement 2A), although Aire and Mki67 were highly expressed (Figure 7C and Figure 7—figure supplement 2A). Importantly, in comparison to the mCherrylo fraction, the mCherryhi-RTOC fraction showed higher levels of Aire-dependent TSAs (Figure 7B). Moreover, beside cell cycle-related genes, some genes were highly expressed in all mCherrylo, Venus+ cells, and cluster R1 cells (Figure 7—figure supplement 2B and Supplementary file 7). Notably, these gene set were down-regulated in mCherryhi-RTOC (Figure 7D and Figure 7—figure supplement 2C). These data suggest that mCherrylo cells were converted into mCherryhi in RTOC.

Figure 7 with 3 supplements see all
Fate mapping study of proliferating Aire+ medullary thymic epithelial cells (mTECs) in in vitro reaggregated thymic organ culture (RTOC).

(A) RTOC experiment to test the differentiation capacity of proliferating Aire+ mTECs. Proliferating Aire+ mTECs (mCherrylo) and E15.5 embryonic thymic cells were reaggregated and subsequently cultured for 5 days. Reaggregated thymic organ (RTO) was analyzed by flow cytometry. Representative flow cytometric profiles of RTOC are shown. N = 5. The ratio of mCherryhi cells in TECs is summarized in the right-hand figure. *p < 0.05. p = 0.027 between CD80hi and CD80lo in mCherrylo and p = 0.024 between CD80hi mCherrylo and CD80hi RTOC control. (B) Volcano plots of RNA-seq data from mCherrylo CD80hi mTECs (mCherrylo), mCherryhi CD80hi mTECs (mCherryhi), and mCherryhi CD80hi mTECs in RTOC (mCherryhi in RTOC). Red dots in volcano plots indicate genes for which expression differed significantly between the two subsets. Numbers of differentially expressed genes are shown in the panels. N = 3. Y axis is log10 of FDR p-value. (C) Expression levels of Aire and Mki67 in mCherrylo, mCherryhi, and mCherryhi in RTOC. (D) Scatter plot of normalized expression values of TA-TEC marker candidates in mCherrylo and mCherryhi in RTOC. TA-TEC marker candidate genes were selected from bulk RNA-seq data and scRNA-seq data in Figure 7—figure supplement 2 (E) Integration of well-based single-cell random displacement amplification sequencing (scRamDA-seq) data (mCherrylo, mCherryhi, and mCherryhi in RTOC) with the droplet-based scRNA-seq data in Figure 2. (F) Frequency of each cell cluster in scRamDA-seq data of mCherrylo, mCherryhi, and mCherryhi -RTOC. (G) Volcano plot of tissue-specific antigen (TSA) expression in each cell cluster in scRamDa-seq data of mCherryhi-RTOC as compared to mCherrylo. Red dots indicate significantly changed TSA genes.

In order to detail phenotypes of mCherryhi-RTOC, we next performed well-based scRNA-seq. mCherryhi-RTOC, in addition to mCherryloCD80hi and mCherryhiCD80hi mTECs from the murine thymus, were single-cell sorted by flow cytometry, and then gene expression in individual cells was determined by random displacement amplification sequencing (RamDA-seq) technology (Hayashi et al., 2018). After quality control of the data, gene expression matrix data of single-cell RamDA-seq (scRamDA-seq) were integrated with the droplet-based scRNA-seq data (Figure 7E). Although this integration slightly changed the UMAP dimension and clustering compared to Figure 2, assignment of each cluster was successfully achieved in practically the same fashion (Figure 7—figure supplement 3), except that cluster R15 (s-TEC) in Figure 2 was incorporated into cluster R10 (mTEC IIIb) and one new cluster was separated from cluster R2 and R3.

Cells from the mCherryloCD80hi mTEC fraction (total 36 cells) were assigned mainly to clusters R1 (17 cells) and R9 (11 cells) (Figure 7E and F, and Supplementary file 8). Some cells were assigned to clusters R0 (3 cells) and R2 (2 cells). Although other cells were assigned to clusters R4, R7, and R14, the embedded position was separated from each parent cluster, which may be due to misclustering. In contrast, cells in the mCherryhiCD80hi mTEC fraction (total 35 cells) were more heterogenous and consisted of cells assigned mainly to clusters R0 (7 cells), R3 (9 cells), R5 (4 cells), R7 (3 cells), R10/15 (5 cells), and R13 (2 cells) (Figure 7E and F, and Supplementary file 8). Except for cluster R5, these clusters were concordant with Aire+ mTECs, post-Aire mTECs, and GP2+ TECs. Notably, after RTOC, heterogenous cell populations including clusters R0 (18 cells), R3 (13 cells), R5 (5 cells), R6 (3 cells), R7 (8 cells), and R10/15 (5 cells) were found in the mCherryhi-RTOC population (total 65 cells). Its composition was relatively similar to that of the mCherryhiCD80hi mTEC fraction (Figure 7F). Moreover, these mCherryhi-RTOC cells expressed high levels of TSAs (Figure 7G). Interestingly, 5 cells in mCherryhi-RTOC were assigned to cluster R5, which also reside in the mCherryhiCD80hi mTEC fraction from the adult thymus. This finding is consistent with the idea of an ‘intertypical’ mTEC cluster, which reportedly contains both CD80hi mTECs and CD80lomTECs (Baran-Gale et al., 2020). Overall, these data suggest that mCherryloCD80hi mTECs differentiate into quiescent mature mTECs expressing high levels of TSAs, including Aire+ mTECs (mTEC II), post-Aire mTECs (mTEC III), and tuft-like mTECs (mTEC IV).

Proliferating Aire+ mTECs are present after puberty in mice

We investigated whether proliferative Aire+ mTECs persisted in thymi of older mice. TECs were analyzed in 4-, 8-, and 19-week Fucci Aire-GFP mice. Flow cytometric analysis showed that a Venus+ mTEChi subset was present in 19-week mice as well as younger mice (Figure 8A). Moreover, Venus+ mTEChi cells expressed Aire genes (Figure 8A). These data strongly suggested that TA-TECs persist in the adult thymus as a source of mature mTECs.

Figure 8 with 2 supplements see all
Proliferating Aire+ CD80hi medullary thymic epithelial cells (mTECs) persist in older mice.

(A) Flow cytometry analysis of CD80hi mTEC subsets from Fucci2a mice aged 4, 8, and 19 weeks. Representative data are shown. Percentages of Venus+ cells in CD80hi mTEC subsets are summarized in the graph in the right panel. N = 4 each for Airegfp/+:: Fucci2a (circle) and control::Fucci2a (closed triangles). (B) Schematic depiction of the proposed process of Aire+ mTEC development in the adult thymus. Transit-amplifying TSAlo Aire+ TECs give rise to mature mTECs. Precursor cells to the transit-amplifying TECs have not been determined yet. Cluster numbers in Figure 1 are shown together with the model of mTEC subsets I to IV.

Integrative computational analysis of our scRNA-seq data with a previously reported dataset of fetal TECs (E12 to E18) (Kernfeld et al., 2018) showed considerably different cell embedding between adult TECs and fetal TECs (Figure 8—figure supplement 1). This trend was common when other scRNA-seq data of adult TECs (Bornstein et al., 2018; Dhalla et al., 2020; Wells et al., 2020) was integrated with the fetal TEC data (Figure 8—figure supplement 2). A TEC-expressing subset was present in the fetal thymus whereas Aire expression was low (clusters F3 and F12, Figure 8—figure supplement 1B). This implies that fetal proliferating mTECs may have a different gene expression profile than adult proliferating Aire+CD80hi mTECs.

Discussion

With regard to mTEC differentiation in the adult thymus (Figure 8B), we hypothesize that Aire+ TA-TECs were generated from their Aire-negative progenitors. Aire+ TA-TECs (cluster 3) undergo cell division and then differentiate into quiescent Aire+ mTECs (cluster 0) through a transition stage, which corresponds to cluster R9 in scRNA-seq data. This differentiation process is accompanied by a chromatin structure change. Post-mitotic Aire+ mTECs begin high-level TSA expression, and further differentiate into post-Aire mTECs (R7, R10, and R13) by closing the Aire enhancer region. Differentiation of mTECs expressing TSAs may have to coordinate differentiation with cell cycle regulation, as proposed in neural cells and muscle differentiation (Ruijtenberg and van den Heuvel, 2016).

Previous scRNA-seq of TEC and thymic stroma suggested the presence of a TEC subset expressing proliferative marker genes (Dhalla et al., 2020; Baran-Gale et al., 2020; Wells et al., 2020). Wells et al. proposed that TA-TECs belong to MHC class II (MHC II)loAire subsets (Wells et al., 2020). In contrast, our analysis of Fucci and Aire reporter mice and cell fate mapping indicated the presence of CD80hi AIRE+ TA-TECs. This inconsistency may be due to the fact that surface expression of MHC class II and CD80 in TA-TECs could be slightly lower than mature mTECs. Interestingly, analysis using the Fucci reporter suggested the existence of CD80lo Aire proliferating mTECs. It is important to address whether Aire proliferating TECs are TA-TECs or a proliferative population of CD80lo mTECs in the future.

Interestingly, although TA-TECs express AIRE protein, the expression level of Aire-dependent TSAs is much lower than that of quiescent AIRE+ TECs. This suggests that AIRE protein is required, but not sufficient for induction of this TSA expression. Several scenarios may explain this fact. First, AIRE requires other regulators for TSA expression. It was reported that AIRE binds to various proteins, possibly regulators of AIRE-dependent TSA expression. Some of these essential regulators may be absent in Aire+ TA-TECs. Another possibility is that cell proliferation inhibits TSA expression induced by AIRE mTECs. There may be mechanisms that suppress the AIRE function in TA-TECs, and thereby inhibit incidence of unfavorable cell states, such as tumor onset due to promiscuous gene expression.

Generally, in other tissues, TACs constitute a link between stem cells and mature cells (Zhang and Hsu, 2017). An important question is, ‘What cells differentiate into proliferating Aire+ mTECs?’ Previous studies have suggested that mTECslo expressing low levels of maturation markers, that is, CD80 or MHC II, are precursors (Abramson and Anderson, 2017; Gray et al., 2007). However, several recent studies have suggested that mTEClo contains several subsets, including CCL21a-positive mTECs, tuft-like mTECs, and others. One possible explanation for this is that a small number of mTEC stem cells or other precursor cells may be present in the mTEClo subset (Gray et al., 2007). Consistently, RNA velocity analysis also suggested that most mTEClo cells do not appear to differentiate into Aire-expressing mTECs. Given that transit-amplifying mTECs are present, a small number of stem/precursor cells would theoretically be sufficient for mTEC reconstitution. A previous study proposed that TECs expressing claudin 3/4 and SSEA-1 had characteristic features of mTEC stem cells in embryonic thymus (Sekai et al., 2014). We failed to detect a corresponding cluster of mTEC stem cells as a subset of adult scRNA clusters. This may be because corresponding mTEC stem cells in adult thymus are included in the ‘intertypical’ TEC cluster, which may be a mixture of various TECs (Baran-Gale et al., 2020). More detailed characterization of mTEC stem cells in the adult thymus is necessary to illuminate differentiation dynamics of mTECs.

In this study, we performed and combined two experiments in scRNA-seq and scATAC-seq of TECs prepared from age- and gender-matched mice. Integration and splitting analysis of datasets suggested substantial reproducibility in both single-cell analyses. Especially, UMAP and clustering of scRNA-seq data were similar between the two experiments, even without batch effect reduction (Figure 2—figure supplement 2). This good consistency in cell clustering may be due to the multiple dimensions of single-cell data. However, it should be noted that only two droplet-based, single-cell analyses may be less reliable in terms of differential gene expression, because only a few thousand genes are normally detected in individual cells. Thus, it may be important to confirm the results obtained from single-cell analysis by using other ways, such as cytometric analysis.

Overall, the scRNA-seq analysis in the present study suggested the presence of a novel differentiation process of TECs in adult thymus. Disturbance of adult TEC homeostasis may cause thymoma, autoimmunity, and other diseases. Further characterization of molecular mechanisms underlying differentiation and maintenance processes in TECs will aid development of novel therapeutic strategies against these thymus-related diseases.

Materials and methods

Mice

C57BL/6 mice were purchased from Clea Japan. Littermates or age-matched, wild-type mice from the same colonies as the mutant mice were used as controls. Aire-GFP mice (CDB0479K, http://www2.brc.riken.jp/lab/animal/detail.php?brc_no=03515) and B6;129-Gt(ROSA)26Sor < tm1(Fucci2aR)Jkn> (RBRC06511) (Fucci2a) (Mort et al., 2014) were provided by the RIKEN BRC through the National Bio-Resource Project of the MEXT, in Japan. CAG-Cre transgenic mice were kindly provided by Dr Jun-ichi Miyazaki (Sakai and Miyazaki, 1997). B6(Cg)-Foxn1tm3(cre)Nrm/J are from Jackson Laboratory (Gordon et al., 2007). Fucci2a mice were crossed with CAG-Cre or Foxn1-Cre mice to activate mCherry and Venus expression. Fucci mice crossed with CAG-Cre were used for all experiments except for immunostaining experiments (Figure 5). All mice were maintained 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). Almost all of available mutant and control mice were randomly used for experiments without any selection.

Preparation of TEC suspensions and flow cytometry analysis

Request a detailed protocol

Murine thymi were minced using razor blades. Thymic fragments were then pipetted up and down to remove lymphocytes. Then, fragments were digested in RPMI 1640 medium containing Liberase (Roche, 0.05 U/mL) plus DNase I (Sigma-Aldrich) via incubation at 37°C for 12 min three times. Single-cell suspensions were stained with anti-mouse antibodies. Dead cells were excluded via 7-aminoactinomycin D staining. Cells were sorted using a FACS Aria instrument (BD). Post-sorted cell subsets were determined to contain >95% of relevant cell types. Data were analyzed using Flowjo 10. No data points or mice were excluded from the study. Randomization and blinding were not used.

Droplet-based scRNA-seq analysis

Request a detailed protocol

For scRNA-seq analysis, cell suspensions of thymi from three mice were prepared and pooled for each individual scRNA-seq experiment. Two experiments were performed. Cellular suspensions were loaded onto a Chromium instrument (10× Genomics) to generate a single-cell emulsion. scRNA-seq libraries were prepared using Chromium Single Cell 3′ Reagent Kits v2 Chemistry and sequenced in multiplex on the Illumina HiSeq2000 platform (rapid mode). FASTQ files were processed using Fastp (Chen et al., 2018). Reads were demultiplexed and mapped to the mm10 reference genome using Cell Ranger (v3.0.0). 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 two datasets, correction of dataset-specific batch effects, UMAP dimension reduction, cell cluster identification, conserved marker identification, and regressing out cell cycle genes) were performed using Seurat (Butler et al., 2018). Briefly, cells that contained a percentage of mitochondrial transcripts > 15% were filtered out. Genes that were expressed in more than five cells and cells expressing at least 200 genes were selected for analysis. Two scRNA-seq datasets were integrated with a combination of Find Integration Anchors and Integrate Data functions (Stuart et al., 2019). Resolution was set as 0.6 for the FindCluster function. In subclustering analysis of cluster R1, the resolution of FindCluster function was set to 0.7. Murine cell cycle genes equivalent to human cell cycle genes listed in Seurat were used for assigning cell cycle scores.

For comparison with a previously reported RNA-seq dataset obtained via a well-based study (Bornstein et al., 2018), the expression matrix of unique molecule identifiers was used. Integration of the two datasets was performed using the Seurat package as described above. RNA velocity analysis was performed using velocyto. Bam/sam files obtained from the Cell Ranger pipeline were transformed to loom format on velocyto.py. RNA velocity was estimated and visualized using loom files by the velocyto R package and pagoda2.

Droplet-based scATAC-seq analysis

Request a detailed protocol

In scRNA-seq analysis, cell suspensions of thymi from two mice were prepared and pooled for each individual scRNA-seq experiment. The EpCAM+CD45TER119 fraction was collected using a cell sorter (BD Aria). After washing with PBS containing 0.04% BSA, sorted cells were suspended in lysis buffer containing 10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.1% Tween-20, 0.1% NP-40, 0.01% Digitonin, and 1% BSA on ice for 3 min. Wash buffer containing 10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.1% Tween-20, and 1% BSA was added to the lysed cells. After centrifuging the solution, a nuclear pellet was obtained by removing the supernatant and the pellet was re-suspended in wash buffer. The concentration of nuclei and their viability were determined by staining with acridine orange/propidium iodide, and 10,000 nuclear suspensions were loaded onto a Chromium instrument (10× Genomics) to generate a single-cell emulsion. scATAC-seq libraries were prepared using Chromium Next GEM Single-Cell ATAC Reagent Kits v1.1 and sequenced in multiplex on an Illumina HiSeq X Ten platform. Reads were demultiplexed and mapped to the mm10 reference genome with Cell Ranger ATAC. Processing data with the Cell Ranger pipeline was performed using the NIG supercomputer at ROIS National Institute of Genetics. Downstream analyses (integration of two datasets, correction of dataset-specific batch effects, UMAP dimension reduction, cluster identification, and identification of differentially accessible regions) were performed using Signac (v1.6) and the Harmony algorithm. Dimension reduction was done using latent semantic indexing. Cells were filtered according to the following parameters: peak_region_fragments 2000–20,000, percentage of fragments in peaks > 50%, blacklist_ratio < 0.03, nucleosome_signal < 0.8, and TSS enrichment > 2. To create a gene activity matrix from scATAC-seq data, the number of fragments in gene coordinates and their 2 kb upstream regions were counted. Integration of scATAC-seq and scRNA-seq data was performed with Signac (v1.6.) using the gene activity matrix in scATAC-seq. The gene activity matrix in scATAC-seq was transferred to scRNA-seq data.

Well-based scRNA-seq analysis

Request a detailed protocol

Single cells were sorted into PCR tubes containing 1 μL of cell lysis solution (1:10 Cell Lysis buffer [Roche], 10 U/μL Rnasin plus Ribonuclease inhibitor [Promega]) using a cell sorter, shaken at 1400 rpm for 1 min with a thermo mixer, and then stored at –80°C. Cell lysates were denatured at 70°C for 90 s and held at 4°C until the next step. To eliminate genomic DNA contamination, 1 μL of genomic DNA digestion mix (0.5 × PrimeScript Buffer, 0.2 U of DNase I Amplification Grade, in RNase-free water) was added to 1 μL of the denatured sample. The mixtures were mixed by gentle tapping, incubated in a T100 thermal cycler at 30°C for 5 min and held at 4°C until the next step. One microliter of RT-RamDA mix (2.5× PrimeScript Buffer, 0.6 pmol oligo(dT)18, 8 pmol 1st-NSRs, 100 ng of T4 gene 32 protein, and 3× PrimeScript enzyme mix in RNase-free water) was added to 2 μL of the digested lysates. The mixtures were mixed with gentle tapping, and incubated 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. Then, the mixtures were held at 4°C until the next step. After RT, the samples were added to 2 μL of second-strand synthesis mix containing 2.25 × NEB buffer 2 (NEB), 0.625 mM each dNTP Mixture (NEB), 40 pmol 2nd-NSRs, and 0.75 U of Klenow Fragment (NEB) in RNase-free water. Mixtures were again mixed by gentle tapping, and incubated at 16°C for 60 min, 70°C 10 min and then at 4°C until the next step. The above-described double-stranded cDNA was purified using 15 μL of AMPure XP SPRI beads (Beckman Coulter) diluted twofold with Pooling buffer (20% PEG8000, 2.5 M NaCl, 10 mM Tris-HCl pH 8.0, 1 mM EDTA, 0.01% NP40) and Magna Stand (Nippon Genetics). Washed AMPure XP beads attached to double-stranded cDNAs were directly eluted using 3.75 μL of 1× Tagment DNA Buffer (10 mM Tris-HCl pH 8.5, 5 mM MgCl2, 10% DMF) and mixed well using a vortex mixer and pipetting. Diluted Tn5-linker complex was added to the eluate and the mixture was incubated at 55°C for 10 min, then 1.25 μ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 1.0× the volume of AMPure XP beads and eluted into 24 μL of 10 mM Tris-Cl, pH 8.5. Reads were demultiplexed and mapped to the mm10 reference genome with STAR. Cells with detected read counts less than half and greater than 1.8 times of the average count were omitted from the analysis. Fucci-negative cells were also removed. Integration of well-based scRNA-seq data with 10× scRNA-seq data and UMAP dimension were performed using Seurat (Butler et al., 2018). Genes that were expressed in more than five cells and cells expressing at least 200 genes were selected for analysis.

Standard RNA-seq analysis

Request a detailed protocol

Total RNA was prepared using TRIzol reagent (Thermo Fisher Scientific) in accordance with the manufacturer’s protocol. After rRNA was depleted using the NEBNext rRNA Depletion Kit, the RNA sequence library was prepared using the NEBNext Ultra Directional RNA Library Prep Kit (New England Biolabs). Paired-end sequencing was performed with NextSeq500 (Illumina). Sequence reads were quantified for annotated genes using CLC Genomics Workbench (Version 7.5.1; Qiagen). Gene expression values were cut off at a normalization expression threshold value of 3. Differential expression was assessed via empirical analysis with the DGE (edgeR test) tool in CLC Main Workbench, in which the Exact Test of Robinson and Smyth was used (Robinson and Smyth, 2008). An FDR-corrected p value was used for testing statistics for RNA-seq analysis. Previously described lists of TSAs and Aire-dependent TSAs (Sansom et al., 2014) were used for the analysis.

RTOC and RNA-seq analysis

Request a detailed protocol

mCherrylo cells (4 × 104–1 × 105) were sorted from Fucci mice and subsequently reaggregated with trypsin-digested thymic cells (1–2 × 106) from E15.5 wild-type mice. RTOCs were cultured on Nucleopore filters (Whatman) placed in R10 medium containing RPMI1640 (Wako) supplemented with 10% fetal bovine serum, 2 mM L-glutamine (Wako), 1× nonessential amino acids (Sigma-Aldrich), 0.1 pM cholera Toxin Solution (Wako 030-20621), 5 μg/mL insulin solution from bovine pancreas (SIGMA I0516-5ML), 2 nM triiodo-L-thyronine (SIGMA T2877-100MG), 1000 units/mL LIF (nacalai NU0012-1), 0.4 μg/mL hydrocortisone,10 ng/mL EGF (Gibco PMG8041), 1 μg/mL RANKL (Wako), penicillin-streptmycin mixed solution (Nacalai Tesque), and 50 µM 2-mercaptoethanol (Nacalai Tesque) for 5 days. For RNA-seq of RTOC experiments, RamDA-seq was used (Hayashi et al., 2018), which allows RNA-seq analysis of low numbers of cells. Briefly, sorted cells were lysed in TCL buffer (Qiagen). After purification of nucleic acids with Agencourt RNA Clean XP (Beckman Coulter) and subsequent treatment with DNase I, the RT-RamDA mixture containing 2.5× PrimeScript Buffer (TAKARA), 0.6 µM oligo(dT)18 (Thermo), 10 µM 1st NSR primer mix, 100 µg/mL of T4 gene 32 protein, and 3× PrimeScript enzyme mix (TAKARA) were added to the purified nucleic acids for reverse transcription. Samples were added to second-strand synthesis mix containing 2× NEB buffer 2 (NEB), 625 nM dNTP Mixture (NEB), 25 µM 2nd NSR primers, and 375 U/mL of Klenow Fragment (3'–5' exo-) (NEB). After cDNA synthesis and subsequent purification by AMPure XP (Beckman Coulter), sequencing library DNA was prepared using the Tn5 tagmentation-based method. Single-read sequencing was performed using a HiSeq2500 (v4, high out mode). Sequence reads were quantified for annotated genes using CLC Genomics Workbench (Version 7.5.1; Qiagen).

Immunohistochemistry

Request a detailed protocol

The thymus was fixed with 4% paraformaldehyde and frozen in OCT compound. After washing cryosections (5 µm) with PBS, sections were blocked with 10% normal goat serum. Keratin-5 was detected using a combination of a polyclonal rabbit anti-mouse keratin-5 antibody (1:500) and Alexa Fluor 647-donkey-anti-rabbit IgG. Aire was detected using a labeled monoclonal antibody (1:300). Confocal color images were obtained using an LAS X (Leica) microscope.

Immunocytochemistry

Request a detailed protocol

Thymic cell suspensions prepared via Liberase digestion were stained with anti-CD45-PE and anti-TER119-PE. After depletion of labeled CD45+ and TER119+ cells via anti-PE microbeads and a magnetic-activated cell sorting separator, negatively selected cells were stained with antiEpCAM (CD326), anti-CD80, anti-Ly51, and UEA-1. Venus+ CD80hi mTECs were sorted and spun down on glass slides using a cytospin. Slides were then fixed with acetone and stained with anti-Aire antibody and DAPI for nuclear staining. Confocal images were obtained using an LAS X microscope.

Statistical analysis

Request a detailed protocol

Statistically significant differences between mean values were determined using Student’s t-test (***p < 0.001, **p < 0.01, and *p < 0.05). Principle component analysis was performed using the prcomp function in R-project. The sample size was not predetermined by statistical methods, but was based on common practice and previous studies (Akiyama et al., 2016; Akiyama et al., 2014). All replicates are biological replicates. All outliers were included in the data.

Appendix 1

Appendix 1—key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Genetic reagent (Mus musculus)B6.Cg-Aire < tm2Mmat>/RbrcRIKEN BioResource Research CenterBRC No:RBRC03515
Genetic reagent (Mus. musculus)B6;129-Gt(ROSA)26Sor < tm1(Fucci2aR)Jkn>RIKEN BioResource Research CenterBRC No:RBRC06511
Genetic reagent (Mus. musculus)B6(Cg)-Foxn1tm3(cre)Nrm/JJackson LaboratoryIMSR Cat#JAX:018448, RRID:IMSR_JAX:018448
Genetic reagent (Mus. musculus)CAG-Cre transgenic miceProvided by Jun-ichi Miyazaki
AntibodyPurified anti-mouse CD16/32(Rat monoclonal)BioLegendCat#101302, RRID:AB_312801FACS(1:200)
AntibodyAPC/Cyanine7 anti-mouse CD45(Rat monoclonal)BioLegendCat#103116, RRID:AB_312981FACS (1:200)
AntibodyPE Rat anti-mouse CD45(Rat monoclonal)eBioscienceCat#12-0451-82, RRID:AB_465668FACS (1:200)
AntibodyAPC/Cyanine7 anti-mouse TER-119/Erythroid Cells(Rat monoclonal)BioLegendCat#116223, RRID:AB_2137788FACS (1:200)
AntibodyPE anti-mouse TER-119/Erythroid Cells(Rat monoclonal)eBioscienceCat#12-5921-82, RRID:AB_466042FACS (1:200)
AntibodyBrilliant Violet 510 anti-mouse CD326 (Ep-CAM)(Rat monoclonal)BioLegendCat#118231, RRID:AB_2632774FACS (1:400)
AntibodyFITC anti-mouse CD326 Ep-CAM (Rat monoclonal)BioLegendCat#118208, RRID:AB_1134107FACS (1:400)
AntibodyAlexa Fluor 647 anti-mouse Ly-51 (Rat monoclonal)BioLegendCat#108312, RRID:AB_2099613FACS (1:400)
Chemical compound, drugBiotinylated Ulex Europaeus Agglutinin I (UEA I)Vector LaboratoriesCat#B-1065–2FACS (1:800)
Chemical compound, drugStreptavidin PE/Cyanine7 ConjugateeBioscienceCat#25-4317-82FACS (1:800)
Chemical compound, drugStreptavidin APC/Cyanine7 ConjugateBD PharmingenCat#554063 RRID:AB_10054651FACS (1:400)
AntibodyPE anti-mouse CD80 (Armenian hamster monoclonal)eBioscienceCat#12-0801-81, RRID:AB_465751FACS (1:300)
AntibodyPacific Blue anti-mouse CD80 Antibody (Armenian hamster monoclonal)BioLegendCat#104724, RRID:AB_2075999FACS (1:300)
AntibodyAlexa Fluor 647 anti-mouse Aire (Rat monoclonal)eBioscienceCat#51-5934-80IHC (1:100)
AntibodyPurified Rabbit anti-Keratin 5 (rabbit polyclonal)BioLegendCat#905504, RRID:AB_2616956IHC (1:400)
AntibodyAlexa Fluor 647 anti-Rabbit IgG (H + L)(Donkey polyclonal)InvitrogenCat#A-31573, RRID:AB_2536183IHC (1:1000)
Chemical compound, drugLiberase TMRoche DiagnosticsCat#5401127001
Chemical compound, drug7-Aminoactinomycin DCalbiochemCat#129935-1MGCN
Chemical compound, drugSYTOX Blue Nucleic Acid StainInvitrogenCat#S11348
Software, algorithmFlowJo version 10BDFlowJo, RRID:SCR_008520
Software, algorithmCell Ranger v3.0.010× GenomicsCell Ranger, RRID:SCR_017344
Software, algorithmSEURAT version 4.1.0https://github.com/satijalab/seurat/blob/master/vignettes/install.RmdSEURAT, RRID:SCR_007322
Software, algorithmVelocyto version 0.6https://github.com/velocyto-team/velocyto.RVelocyto, RRID:SCR_018167
Software, algorithmpagoda2 version 1.0.9https://github.com/kharchenkolab/pagoda2pagoda2, RRID:SCR_017094
Software, algorithmCell Ranger ATAC version1.1.010× GenomicsCell Ranger ATAC, RRID:SCR_021160
Software, algorithmSignac version 1.5.0https://github.com/timoast/signac/blob/master/vignettes/install.RmdSignac, RRID:SCR_021158
Software, algorithmMonocle3, version 0.2.3https://cole-trapnell-lab.github.io/monocle3/docs/installation/Monocle3, RRID:SCR_018685
Software, algorithmCLC Genomics Workbench Version 7.5.1QIAGENCLC Genomics Workbench, RRID:SCR_011853
Commercial assay or kitChromium Single Cell 3’ Library & Gel Bead Kit v210× GenomicsCat#PN-120237
Commercial assay or kitChromium Single Cell A Chip Kit10× GenomicsCat#PN-120236
Commercial assay or kitChromium i7 Multiplex Kit10× GenomicsCat#PN-120262
Commercial assay or kitChromium Next GEM Single Cell ATAC Library & Gel Bead Kit10× GenomicsCat#PN-1000176
Commercial assay or kitChromium Next GEM Chip H Single Cell Kit10× GenomicsCat#PN-1000162
Commercial assay or kitSingle Index Kit N, Set A10× GenomicsCat#PN-1000212
Commercial assay or kitNEBNext rRNA Depletion KitNew England BiolabsCat#E6310
Commercial assay or kitNEBNext Ultra Directional RNA Library Prep Kit for IlluminaNew England BiolabsCat#E7420
Commercial assay or kitKAPALibraryQuantificationKits Illumina/UniversalNippon GeneticsCat#KK4824
Chemical compound, drugKAPAHiFi DNA PolymeraseNippon GeneticsCat#KK2102
Commercial assay or kitAgilent High Sensitivity DNA KitAgilent TechnologiesCat#5067–4626
Commercial assay or kitMultina DNA-12000SHIMADZUCat#S292-36600-91

Data availability

FASTQ data of RNA-Seq and ATAC-seq are deposited in DDBJ (DRA009125, DRA010209 , DRA012308, DRA012309, DRA012452, and DRA013875).

The following data sets were generated
    1. Akiyama T
    (2019) DDBJ
    ID DRA009125. Single cell analysis of thymic epithelial cells.
    1. Akiyama T
    (2020) DDBJ
    ID DRA010209. RNA-seq analysis of transit-amplifying Aire+ mTECs.
    1. Akiyama T
    (2021) DDBJ
    ID DRA012308. Single cell analysis of transit-amplifying thymic epithelial cells.
    1. Akiyama T
    (2021) DDBJ
    ID DRA012309. RNA-seq analysis of transit amplifying mTEC in reaggregation thymic organ culture.
    1. Akiyama T
    (2021) DDBJ
    ID DRA012452. Single cell ATAC sequence analysis of thymic epithelial cells (Exp.1).
    1. Akiyama T
    (2022) DDBJ
    ID DRA013875. Single cell ATAC sequence analysis of thymic epithelial cells (Exp.2).
The following previously published data sets were used
    1. Amit I
    2. Abramson J
    3. Bornstein C
    4. Nevo S
    5. Giladi A
    6. Kadouri N
    (2018) NCBI Gene Expression Omnibus
    ID GSE103967. Large-scale single cell mapping of the thymic stroma identifies a new thymic epithelial cell lineage.
    1. Dhalla F
    2. Baran-Gale J
    3. Maio S
    4. Chappell L
    5. Hollander G
    6. Ponting CP
    (2019) ArrayExpress
    ID E-MTAB-8105. Single cell RNA-seq of medullary thymic epithelial cells (mTEC).
    1. Kernfeld E
    2. Genga R
    3. Magaletta M
    4. Neherin K
    5. Xu P
    6. Maehr R
    (2018) NCBI Gene Expression Omnibus
    ID GSE107910. Single-cell RNA sequencing resolves cellular heterogeneity throughout embryonic development of the thymus.
    1. Wells KL
    2. Miller CN
    3. Gschwind AR
    4. Phipps JD
    5. Anderson MS
    6. Steinmetz LM
    (2020) NCBI Gene Expression Omnibus
    ID GSE137699. Single cell sequencing defines a branched progenitor population of stable medullary thymic epithelial cells.

References

    1. Lajtha LG
    (1979) Stem cell concepts
    Differentiation; Research in Biological Diversity 14:23–34.
    https://doi.org/10.1111/j.1432-0436.1979.tb01007.x

Decision letter

  1. Shimon Sakaguchi
    Reviewing Editor; Osaka University, Japan
  2. Betty Diamond
    Senior Editor; The Feinstein Institute for Medical Research, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Integrative analysis of scRNAs-seq and scATAC-seq revealed transit-amplifying thymic epithelial cells expressing autoimmune regulator" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Betty Diamond as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. The scATACseq experiment included a relatively low number (n=2) of mice. Understandably this is a complicated experiment, however, the variation between individual animals remains unknown and may influence the interpretation of the results. The authors are advised to show the clustering analysis for the two animals separately in a supplement to confirm that the changes are seen in both animals, not just in one. Similarly, the individual scRNA-seq UMAPs of the three animals should be included in the supplement. The authors should discuss the limitations of small groups in single-cell experiments in Discussion. Technically the integration of scATAC and scRNA seq results may be feasible.

2. Genes expressed in clusters should be provided in the supplement. Without this information, it is difficult to see the feasibility of the cluster annotations and to compare them with other similar studies. The authors should give an unbiased list of all genes expressed in mTECs, in particular those genes that are expressed in TACs. To harmonize the findings in the field, it would be useful for the readers if the authors would compare their cluster-specific genes for the overlap with TEC single-cell clusters from other studies (Bornstein et al., 2018, Dhalla et al., 2020, Baran-Gale et al., 2020, Wells et al., 2020). In addition, Wells et al., 2020 reported mTEC TACs to be mTEClo preceding Aire expression and giving rise to both Aire+ and Ccl21a+ cells. How do the authors reconcile their results (TAC as Aire+) with Wells et al., paper (TACs as Aire-)? How do they interpret the finding that Wells et al., find Ki67 significantly higher in mTEClo than in mTEChi, implying that cell division precedes high the expression of Aire?

4. The candidate cell population for TACs, cluster R1 expresses Aire and the proliferating cell marker Mki67. R1 also expresses Ccl21a. The authors did subclustering of R1 and found that R1A-D express Aire whereas R1E has a higher expression of Ccl21a. The authors note that "Thus, it is possible that TECs expressing cell-cycle-related genes, proposed by scRNA-seq analysis, contain at least two proliferating TECs subsets having different chromatin accessibilities and gene expression profiles." To confirm that both R1A-D and R1E subsets are proliferating TACs, the authors might show the proliferating gene markers in these subsets. Was Mki67 expressed among all R1 subpopulations? Would this argue for the presence of TAC among both Aire+ and Aire- cell populations? Assuming that TACs as a proliferating cell type should express multiple genes associated with cell cycling, the authors focus on Mki67 only, but did R1 TAC express other proliferating cell markers which would support the claim that these are indeed actively dividing cells?

6. The authors focus their study on the CD80high cycling cells (Aire+). Figure 4 show transcription profiles of the isolated cycling CD80+ mTECs. Dots in the "TSA genes" panel (right) don't appear in the "All genes" panel (left). Are those lost dots of TSA genes?

In Figure 4E (left) low expressed genes seem to be skewed towards Venus- mTEChi (in comparison to high expressed genes). A statistical assessment for the comparisons of the TSA, Aire-dep TSA and Aire-indep TSA profiles to the general profile (Figure 4E and 4G), considering expression levels, would confirm the visual assessment.

They should also discuss the reason why m-cherry low cells express a lower level of tissue-specific antigens even though they express Aire? Is Aire expression alone insufficient for TSA expression? Does the transcriptome data provide any mechanistic insight? In Fig4F, Aire expression is similar between Venus+ and Venus-. Is it compatible with Figure 4A showing less Aire-expressing cells in Venus+ than in Venus-?

7. In Figure 6, they showed that RTOC culture of mCherry-low cells produced mCherry-high cells, thereby they claimed that mCherry-low cells differentiated into mCherry-high cells. To substantiate this notion, they should rule out the possibility of selective cell survival of possibly contaminated mCherry-high cells or that transferring adult mTECs into the embryonic cell environment may trigger other signaling pathways that induce the upregulation of the cell cycle and mCherry expression? For example, what about the expression profiles of cell-death/apoptosis-related genes in mCherry low and high cells? What about the cell number of mCherry high cells after RTOC culture using mCherry-high cells alone as a control group? Is the cellularity of survivors comparable to RTOC culture using mCherry-low as shown in figure 6A? Would the same happen if they would transfer these cells to the adult thymus?

8. The authors nicely confirm, by the fine analysis of their scRNA-seq data, that the TAC population contains Aire+ cells and Ccl21+ cells in a visually mutually exclusive manner. However, they don't formally clarify whether the Ccl21-expressing TACs have differences in their chromatin accessibility pattern compared to the Aire-expressing TACs.

It's also worth showing the expression of CD80 in the scRNA-seq UMAP of TACs alone and in the one of all TECs (Figure 2). This would notably allow to determine whether CD80 expression is restricted to Aire-positive TACs or encompasses Aire-negative TACs (Ccl21+).

Also, projection of the R1A-E scRNAseq clusters onto the scATAC-seq UMPA would be enlightening.

9. Trajectory analyses provide a nice confirmation of published results identifying a trajectory from TACs to mTEChi. However, the authors don't discuss whether their data support a potential trajectory from TACs to mTEClo (already suggested/ reported) which seems to be present in Figure 3-Figure sup 2B. Would this mean that some TACs could mature into Ccl21+ mTECs (mTEClo)? and if so, how Aire and Ccl21 are expressed in these TACs? Which TAC sub-cluster do they belong to?

10. In figure S7, the data of fetal thymi showed that proliferating fetal mTECs might have a gene expression profile different from the adult counterpart. One caveat of the interpretation of the data is that the difference between the public data and the authors' data might be attributed to a batch effect because the clusters from two data sets seemed almost completely discrete. They should mention how they processed the data to diminish the risk of such artifacts. Or, it is better to add the data of fetal thymi obtained by the authors themselves, if possible.

Reviewer #1 (Recommendations for the authors):

Regarding figure 4, they should discuss the reason why m-cherry low cells express a lower level of tissue-specific antigens even though they express Aire? Is Aire expression alone insufficient for TSA expression? Does the transcriptome data provide any mechanistic insight?

In figure 6, they showed that RTOC culture of mCherry-low cells produced mCherry-high cells, thereby they claimed that mCherry-low cells differentiated into mCherry-high cells. To substantiate this notion, they should rule out the possibility of selective cell survival of possibly contaminated mCherry-high cells. For example, what about the expression profiles of cell-death/apoptosis-related genes in mCherry low and high cells? What about the cell number of mCherry high cells after RTOC culture using mCherry-high cells alone as a control group? Is the cellularity of survivors comparable to RTOC culture using mCherry-low as shown in figure 6A?

In figure S7, the data of fetal thymi showed that proliferating fetal mTECs might have a gene expression profile different from the adult counterpart. One caveat of the interpretation of the data is that the difference between the public data and the authors' data might be attributed to a batch effect because the clusters from two data sets seemed almost completely discrete. They should mention how they processed the data to diminish the risk of such artifacts. Or, it is better to add the data of fetal thymi obtained by the authors themselves, if possible.

Reviewer #2 (Recommendations for the authors):

– Figure 2A has 2 Ic populations but lacks Ia?

– Figure 1c Y axis is shown as "expression level", this should be rather "accessibility level"?

– The Fucci system incorporates genetically encoded Cherry and Venus probes that highlight G1 and S/G2/M phases of the cell cycle in animal cells. Did the authors control for the transgene inactivation in these mice as in some cases the expression of the transgenes may change?

Reviewer #3 (Recommendations for the authors):

Related to point (1): projection of the R1A-E scRNAseq clusters onto the scATAC-seq UMPA would be enlightening.

https://doi.org/10.7554/eLife.73998.sa1

Author response

Essential revisions:

1. The scATACseq experiment included a relatively low number (n=2) of mice. Understandably this is a complicated experiment, however, the variation between individual animals remains unknown and may influence the interpretation of the results. The authors are advised to show the clustering analysis for the two animals separately in a supplement to confirm that the changes are seen in both animals, not just in one.

Unfortunately, the two mice are not individually labeled, so the data cannot be separated. Therefore, to address the reviewer's concerns, we have performed another scATAC-seq analysis of TEC in mice of the same age (page 4, line 33). As shown in Figure 1—figure supplement 2 and 3, the data were reproducible, confirming clear and similar separation of proliferating TEC clusters.

Similarly, the individual scRNA-seq UMAPs of the three animals should be included in the supplement. The authors should discuss the limitations of small groups in single-cell experiments in Discussion. Technically the integration of scATAC and scRNA seq results may be feasible.

As with the scATAC-seq, UMAP and clustering analysis from the two experiments of scRNA-seq were displayed separately. These biological duplicates were similar, and consistent conclusions were reached. This finding is described on page 5, line 24, and the data are exhibited in Figure 2_supplement_2. Also, a paragraph has been added to the Discussion (page 14, line 19) describing limitations due to the small number of single-cell analyses.

2. Genes expressed in clusters should be provided in the supplement. Without this information, it is difficult to see the feasibility of the cluster annotations and to compare them with other similar studies. The authors should give an unbiased list of all genes expressed in mTECs, in particular those genes that are expressed in TACs.

We added a heatmap figure of the top 5 genes specifically expressed in each cluster in Figure_2_supplement_1. In addition, the list of genes differentially expressed in each cluster was summarized in Supplementary Table 1.

To harmonize the findings in the field, it would be useful for the readers if the authors would compare their cluster-specific genes for the overlap with TEC single-cell clusters from other studies (Bornstein et al., 2018, Dhalla et al., 2020, Baran-Gale et al., 2020, Wells et al., 2020). In addition, Wells et al., 2020 reported mTEC TACs to be mTEClo preceding Aire expression and giving rise to both Aire+ and Ccl21a+ cells. How do the authors reconcile their results (TAC as Aire+) with Wells et al., paper (TACs as Aire-)? How do they interpret the finding that Wells et al., find Ki67 significantly higher in mTEClo than in mTEChi, implying that cell division precedes high the expression of Aire?

Wells and we reported a cell cluster of TECs expressing high levels of cell proliferation markers in droplet scRNA-seq analysis. It is important to link cell subsets defined by scRNAseq with cell populations defined by flow cytometric analysis.

Wells et al., proposed that Ki67 is a marker of TA-TECs for flow cytometric analysis. Although Ki67hi mTECs were assigned as a part of the Aire-negative MHC class IIlo TEC subset, actual data showed that the expression level of MHC class II appears close to the border separating MHCIIlo and MHChi subsets (Figure 2c in Wells et al., eLife 2020; 9: e60188), which was also mentioned by the authors. Moreover, the expression level of the Aire protein in Ki67+ TECs was broad (Figure 2d in Wells et al., eLife 2020; 9: e60188), suggesting that Ki67+ TECs may contain both Aire+ and Aire- mTEC subsets. In contrast, we used Fucci/Aire-GFP dual reporter mice to characterize proliferating mTECs that are mCherrylo and Venus+. CD80 expression of TA-TECs was high, but appears to be close to the border line separating CD80hi and CD80lo subsets. So, expression level of MHC II, CD80, and Aire of proliferating TECs may be slightly lower and broader, compared to mature mTECs. In flow cytometric analysis, definition of cell types based on cell surface markers depends on how the border lines are set. This may be one reason for the apparent inconsistency between their TA-TECs (MHCllloAire-) and our TA-TEC (CD80hiAire+).

Moreover, we noticed that some Venus+ mCherrylo cells (dividing cells in Fucci) are present in the CD80lo fraction, albeit at low frequency (Figure 5—figure supplement 1). This may be consistent with the subcluster analysis in Figure 4, which suggests that the proliferating TEC subset was divided into two; Aire+ subsets (86.2%) and Aire- subsets (13.8%). These data may also support a practical consistency between Wells et al.,’s and our FACS analyses.

It should be noted that not all proliferating cells are transit-amplifying (TA) cells. TA cells are defined as short-lived proliferating cells that link stem cells and mature cell types. Cell fate mapping of possible candidates should be critical for identifying them as “TA”-TEC. Thus, it is important to test the differentiation potential of the proliferating CD80lo fraction into Aire+ and post-Aire mature mTECs in the future.

These points were briefly described in the Discussion (page 13, line 15)

4. The candidate cell population for TACs, cluster R1 expresses Aire and the proliferating cell marker Mki67. R1 also expresses Ccl21a. The authors did subclustering of R1 and found that R1A-D express Aire whereas R1E has a higher expression of Ccl21a. The authors note that "Thus, it is possible that TECs expressing cell-cycle-related genes, proposed by scRNA-seq analysis, contain at least two proliferating TECs subsets having different chromatin accessibilities and gene expression profiles." To confirm that both R1A-D and R1E subsets are proliferating TACs, the authors might show the proliferating gene markers in these subsets. Was Mki67 expressed among all R1 subpopulations? Would this argue for the presence of TAC among both Aire+ and Aire- cell populations? Assuming that TACs as a proliferating cell type should express multiple genes associated with cell cycling, the authors focus on Mki67 only, but did R1 TAC express other proliferating cell markers which would support the claim that these are indeed actively dividing cells?

Thank you for this suggestion. To address these questions, we did sub-cluster the proliferating TEC cluster with a higher resolution (see Materials and methods). Interestingly, the proliferating TEC cluster was divided into 7 subclusters, depending on expression of Aire, Ccl21a, and cell cycle scoring. So, Aire+ clusters and Ccl21a+ clusters were separated, and these two clusters were further separated by expression of cell cycle phase markers (cell cycle scoring) including Mki67. Thus, these data are consistent with the idea that there may be Aire+ and Ccl21a+ TECs in this proliferating TEC cluster. We think these findings are important, and data are shown in new Figure 4 and described on Page 8, line 17.

6. The authors focus their study on the CD80high cycling cells (Aire+). Figure 4 show transcription profiles of the isolated cycling CD80+ mTECs. Dots in the "TSA genes" panel (right) don't appear in the "All genes" panel (left). Are those lost dots of TSA genes?

The axis scales were different between the two panels, which might cause confusion. We corrected the axis. Indeed, dots among the “TSA genes” are present in the “All genes” panels.

In Figure 4E (left) low expressed genes seem to be skewed towards Venus- mTEChi (in comparison to high expressed genes). A statistical assessment for the comparisons of the TSA, Aire-dep TSA and Aire-indep TSA profiles to the general profile (Figure 4E and 4G), considering expression levels, would confirm the visual assessment.

To confirm the data in Figure 4E and G (new Figure 5E and G), GSEA analysis was performed and exhibited in Figure 5_supplement_3. Whereas both Aire dep TSA and Aire-indep TSAs were significantly enriched in Venus- mTEChi, the enrichment (enrichment score) of Aire-dep TSAs was higher than that of Aire-indep TSAs.

They should also discuss the reason why m-cherry low cells express a lower level of tissue-specific antigens even though they express Aire? Is Aire expression alone insufficient for TSA expression? Does the transcriptome data provide any mechanistic insight?

As pointed out by this reviewer, it is possible that AIRE expression is not sufficient for inducing TSAs in mTECs. In another scenario, there may be factors that inhibit the AIRE function during cell division in TA-TECs. This point is described in the Discussion (Page 13, line 24). Unfortunately, we have not gained any useful insights about mechanisms from transcriptomic data. Further comparative analysis (such as epigenetic and chromatin structure analyses) between TSA+ and TSAlo Aire+ TECs will be important to address this issue.

In Fig4F, Aire expression is similar between Venus+ and Venus-. Is it compatible with Figure 4A showing less Aire-expressing cells in Venus+ than in Venus-?

As pointed out by the Reviewer, the fluorescence intensity of Aire GFP in Venus+ mTEChi cells may be lower than that of that of Venus– mTEChi cells in Figure 4A (Figure 5B in revised manuscript). However, as described on Page 9 line 15, a direct comparison between two intensities may be difficult because the fluorescence intensity of Aire GFP in Venus+ mTEChi cells showed a broad peak due to the compensation between GFP and Venus proteins, which have very close fluorescence spectra. In addition to this technical issue, because mTEChi contains both Aire-positive and Aire-negative cells (Post-Aire mTECs), the average expression level of Aire in mTEChi may have become similar to that of Venus+ cells, which contain mainly Aire-positive cells, in the bulk RNA-seq analysis (Figure 5F)

7. In Figure 6, they showed that RTOC culture of mCherry-low cells produced mCherry-high cells, thereby they claimed that mCherry-low cells differentiated into mCherry-high cells. To substantiate this notion, they should rule out the possibility of selective cell survival of possibly contaminated mCherry-high cells or that transferring adult mTECs into the embryonic cell environment may trigger other signaling pathways that induce the upregulation of the cell cycle and mCherry expression? For example, what about the expression profiles of cell-death/apoptosis-related genes in mCherry low and high cells? What about the cell number of mCherry high cells after RTOC culture using mCherry-high cells alone as a control group? Is the cellularity of survivors comparable to RTOC culture using mCherry-low as shown in figure 6A?

This reviewer noted that mCherryhi contamination of sorted cells could preferentially survive in RTOC. We compared expression of pro-apoptotic genes between mCherryhi and mCherrylo. GSEA analysis suggested that pro-apoptotic genes were rather higher in mCherryhi than mCherrylo (Figure_7_figure_supplement_1D). Thus, it is unlikely that mCherryhi could be more resistant to apoptosis, thereby surviving preferentially in RTOC.

Moreover, we have done RTOC using only mCherryhi cells. A purity check of sorted mCherrylo cells indicated that the contamination was below 0.05% (Figure_7_figure_supplement_1C). Therefore, the number of mCherryhi cells mixed with RTOC was changed from 200, 500, and 1000 cells, which would be equivalent to the assumed contamination of approximately 0.2%, 0.5%, and 1 % of mCherryhi cells during cell sorting, respectively. Data suggested that surviving cells are 1, 0, and 5 cells in these experiments. Thus, it is unlikely that detected mCherryhi in RTOC were simply survivors of contaminated mCherryhi. In addition, this experiment rules out the possibility that mCherryhi becomes more proliferative in the embryonic thymus environment. These points are briefly described in Page 11, line 1.

Would the same happen if they would transfer these cells to the adult thymus?

This is a good suggestion because a similar experiment was previously carried out to address the potential of mTEC stem cells (Sekai et al., J. Immunol. Methods 467, 29, 2019). In that experiment, they used 5x 105 cultured cells for the adult thymic transfer. Because we cannot culture TA-TECs to date, to prepare this number of cells from mice, we would need to use 20 thymi per experiment and would then have to sort the TA-TECs. This is quite a difficult experiment in reality. However, our preliminary trial using 1 x 105 TA-TECs was not successful in that we failed to detect transferred cells. Together with the low efficiency of RTOC and the failure of MHC mismatched RTOC, we think some kind of niches are necessary to maintain TA-TECs. Determination of such niches will be an important study in the future.

8. The authors nicely confirm, by the fine analysis of their scRNA-seq data, that the TAC population contains Aire+ cells and Ccl21+ cells in a visually mutually exclusive manner. However, they don't formally clarify whether the Ccl21-expressing TACs have differences in their chromatin accessibility pattern compared to the Aire-expressing TACs.

Integration of TAC-subcluster analysis with scATAC-seq revealed that chromatin accessibility of Aire+ proliferating TECs and Ccl21a+ proliferating cells may not be very different. This was briefly mentioned on Page 8 line 26. However, this point should be carefully concluded because a relatively low number cells of each subcluster was analyzed. So, it is important to compare chromatin structure between these clusters after isolation using a flow cytometer.

It's also worth showing the expression of CD80 in the scRNA-seq UMAP of TACs alone and in the one of all TECs (Figure 2). This would notably allow to determine whether CD80 expression is restricted to Aire-positive TACs or encompasses Aire-negative TACs (Ccl21+).

We demonstrated CD80 expression of the proliferating TEC subcluster in a Feature plot (Figure 4D). Whereas CD80 positive cells are fewer, it is likely that CD80 expressing cells are present among Aire-positive cells, but not among Ccl21+ cells in the proliferating TEC clusters.

Also, projection of the R1A-E scRNAseq clusters onto the scATAC-seq UMPA would be enlightening.

As mentioned above, we projected the proliferating TEC subclusters (R1A-G) onto the scATAC-data. The analysis did not show clear evidence of a difference in chromatin accessibility among these clusters. The result is briefly described in the main text (Page 8, line 26), and the data are presented in Figure 4-figure_supplement_1.

9. Trajectory analyses provide a nice confirmation of published results identifying a trajectory from TACs to mTEChi. However, the authors don't discuss whether their data support a potential trajectory from TACs to mTEClo (already suggested/ reported) which seems to be present in Figure 3-Figure sup 2B. Would this mean that some TACs could mature into Ccl21+ mTECs (mTEClo)? and if so, how Aire and Ccl21 are expressed in these TACs? Which TAC sub-cluster do they belong to?

This point may be related to the question of whether the proliferating TEC subset in scRNA-seq is a single population. We think that velocity analysis of our scRNA-seq data did not show a clear trajectory from the proliferating TEC cluster to Ccl21+ mTECs. However, since the Ccl21a+ cluster is present in the proliferating TEC cluster (R1), it may be possible that these cells differentiate into Ccl21+ mTECs (mTEClo). To date, trajectory analysis of the proliferating Ccl21a cluster has been difficult because of the low cellularity in our scRNA-seq dataset.

10. In figure S7, the data of fetal thymi showed that proliferating fetal mTECs might have a gene expression profile different from the adult counterpart. One caveat of the interpretation of the data is that the difference between the public data and the authors' data might be attributed to a batch effect because the clusters from two data sets seemed almost completely discrete. They should mention how they processed the data to diminish the risk of such artifacts. Or, it is better to add the data of fetal thymi obtained by the authors themselves, if possible.

In Figure S7, a batch effect was removed using the Seurat function, in which selected “anchor” genes commonly expressed in cells are used to remove the batch effect. This function was used to remove the batch effect between our data and those of others, so it turned out to be effective. We also used the Harmony algorithm to remove the batch effect (Figure 8—figure supplement 1C). Again, TEC data of fetal thymus differ from adult TECs. We further integrated fetal thymus data with previously reported scRNA-seq data of TECs (Bornstein et al., Nature 2018; 559:622, Wells et al., eLife 2020;9:e60188, Dhalla et al., eLife. 2020 9:e56221). Data analysis suggested that TECs of fetal thymus seem to be different from adult TEC clusters (Figure 8—figure supplement 2). Overall, it is likely that fetal TECs are considerably different from adult TECs in terms of their gene expression profiles. Unfortunately, we did not perform reproducing experiments for fetal TEC scRNA-seq by ourselves, because it is fairly costly and time-consuming to obtain scRNA-seq data of TECs from various embryonic days.

https://doi.org/10.7554/eLife.73998.sa2

Article and author information

Author details

  1. Takahisa Miyao

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Data curation, Formal analysis, Investigation, Writing – original draft
    Competing interests
    No competing interests declared
  2. Maki Miyauchi

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Data curation, Formal analysis, Investigation, Validation
    Competing interests
    No competing interests declared
  3. S Thomas Kelly

    Laboratory for Cellular Epigenomics, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3904-6690
  4. Tommy W Terooatea

    Laboratory for Cellular Epigenomics, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
  5. Tatsuya Ishikawa

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Data curation, Formal analysis, Investigation
    Competing interests
    No competing interests declared
  6. Eugene Oh

    Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Data curation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8534-0183
  7. Sotaro Hirai

    Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  8. Kenta Horie

    Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  9. Yuki Takakura

    Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Formal analysis, Validation
    Competing interests
    No competing interests declared
  10. Houko Ohki

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  11. Mio Hayama

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Formal analysis, Investigation, Validation
    Competing interests
    No competing interests declared
  12. Yuya Maruyama

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  13. Takao Seki

    Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  14. Hiroto Ishii

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  15. Haruka Yabukami

    Laboratory for Cellular Epigenomics, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  16. Masaki Yoshida

    YCI Laboratory for Immunological Transcriptomics, RIKEN Center for Integrative Medical Sciences, Kanagawa, Japan
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  17. Azusa Inoue

    YCI Laboratory for Metabolic Epigenetics, RIKEN Center for Integrative Medical Sciences, Kanagawa, Japan
    Contribution
    Methodology, Supervision, Validation
    Competing interests
    No competing interests declared
  18. Asako Sakaue-Sawano

    Laboratory for Cell Function Dynamics, RIKEN Center for Brain Science, Saitama, Japan
    Contribution
    Methodology, Supervision
    Competing interests
    No competing interests declared
  19. Atsushi Miyawaki

    Laboratory for Cell Function Dynamics, RIKEN Center for Brain Science, Saitama, Japan
    Contribution
    Methodology, Supervision, Validation
    Competing interests
    No competing interests declared
  20. Masafumi Muratani

    Transborder Medical Research Center, and Department of Genome Biology, Faculty of Medicine, University of Tsukuba, Ibaraki, Japan
    Contribution
    Data curation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0276-8000
  21. Aki Minoda

    Laboratory for Cellular Epigenomics, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Formal analysis, Validation
    Competing interests
    No competing interests declared
  22. Nobuko Akiyama

    Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    Contribution
    Formal analysis, Funding acquisition, Investigation, Supervision, Writing – review and editing
    For correspondence
    nobuko.akiyama@riken.jp
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1085-5996
  23. Taishin Akiyama

    1. Laboratory for Immune Homeostasis, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan
    2. Immunobiology, Graduate School of Medical Life Science, Yokohama City University, Yokohama, Japan
    Contribution
    Formal analysis, Funding acquisition, Investigation, Project administration, Supervision, Validation, Writing – original draft, Writing – review and editing
    For correspondence
    taishin.akiyama@riken.jp
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0341-6154

Funding

Japan Society for the Promotion of Science (Grants-in-Aid for Scientific Research 17H04038)

  • Taishin Akiyama

Japan Society for the Promotion of Science (Grants-in-Aid for Scientific Research 20H03441)

  • Taishin Akiyama

Japan Society for the Promotion of Science (Grants-in-Aid for Scientific Research 17K08622)

  • Nobuko Akiyama

Japan Society for the Promotion of Science (Grants-in-Aid for Scientific Research 20K07332)

  • Nobuko Akiyama

Princess Takamatsu Cancer Research Fund

  • Taishin Akiyama

Uehara Memorial Foundation

  • Taishin Akiyama

NOVARTIS Foundation

  • Taishin Akiyama

Ministry of Education, Culture, Sports, Science and Technology (18H04989)

  • Taishin Akiyama

Ministry of Education, Culture, Sports, Science and Technology (19H04821)

  • Nobuko Akiyama

Japan Science and Technology Agency (CREST JPMJCR2011)

  • Taishin Akiyama

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors declare no competing financial interests. We thank the sequencing staff at the RIKEN Center for Integrative Medical Sciences for assisting with RNA-seq. Computations were performed on the NIG supercomputer at ROIS, National Institute of Genetics and HOKUSAI supercomputer at ISD, RIKEN.

Ethics

All mice were handled in accordance with Guidelines of the Institutional Animal Care and Use Committee of RIKEN, Yokohama Branch . Permission #2018-075.

Senior Editor

  1. Betty Diamond, The Feinstein Institute for Medical Research, United States

Reviewing Editor

  1. Shimon Sakaguchi, Osaka University, Japan

Publication history

  1. Received: September 17, 2021
  2. Preprint posted: October 5, 2021 (view preprint)
  3. Accepted: April 21, 2022
  4. Version of Record published: May 17, 2022 (version 1)
  5. Version of Record updated: June 1, 2022 (version 2)

Copyright

© 2022, Miyao et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,599
    Page views
  • 285
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Takahisa Miyao
  2. Maki Miyauchi
  3. S Thomas Kelly
  4. Tommy W Terooatea
  5. Tatsuya Ishikawa
  6. Eugene Oh
  7. Sotaro Hirai
  8. Kenta Horie
  9. Yuki Takakura
  10. Houko Ohki
  11. Mio Hayama
  12. Yuya Maruyama
  13. Takao Seki
  14. Hiroto Ishii
  15. Haruka Yabukami
  16. Masaki Yoshida
  17. Azusa Inoue
  18. Asako Sakaue-Sawano
  19. Atsushi Miyawaki
  20. Masafumi Muratani
  21. Aki Minoda
  22. Nobuko Akiyama
  23. Taishin Akiyama
(2022)
Integrative analysis of scRNA-seq and scATAC-seq revealed transit-amplifying thymic epithelial cells expressing autoimmune regulator
eLife 11:e73998.
https://doi.org/10.7554/eLife.73998

Further reading

    1. Developmental Biology
    2. Immunology and Inflammation
    Hyun-Woo Jeong, Rodrigo Diéguez-Hurtado ... Ralf H Adams
    Tools and Resources

    The blood-brain barrier (BBB) limits the entry of leukocytes and potentially harmful substances from the circulation into the central nervous system (CNS). While BBB defects are a hallmark of many neurological disorders, the cellular heterogeneity at the neurovascular interface and the mechanisms governing neuroinflammation are not fully understood. Through single cell RNA sequencing of non-neuronal cell populations of the murine cerebral cortex during development, adulthood, ageing and neuroinflammation, we identify reactive endothelial venules (REVs), a compartment of specialised post-capillary endothelial cells (ECs) that are characterized by consistent expression of cell adhesion molecules, preferential leukocyte transmigration, association with perivascular macrophage populations, and endothelial activation initiating CNS immune responses. Our results provide novel insights into the heterogeneity of the cerebral vasculature and a useful resource for the molecular alterations associated with neuroinflammation and ageing.

    1. Immunology and Inflammation
    Shasha Li, Michael D Bern ... Wayne M Yokoyama
    Research Article

    BTB domain And CNC Homolog 2 (Bach2) is a transcription repressor that actively participates in T and B lymphocyte development, but it is unknown if Bach2 is also involved in the development of innate immune cells, such as natural killer (NK) cells. Here, we followed the expression of Bach2 during murine NK cell development, finding that it peaked in immature CD27+CD11b+ cells and decreased upon further maturation. Bach2 showed an organ and tissue-specific expression pattern in NK cells. Bach2 expression positively correlated with the expression of transcription factor TCF1 and negatively correlated with genes encoding NK effector molecules and those involved in the cell cycle. Lack of Bach2 expression caused changes in chromatin accessibility of corresponding genes. In the end, Bach2-deficiency resulted in increased proportions of terminally differentiated NK cells with increased production of granzymes and cytokines. NK cell-mediated control of tumor metastasis was also augmented in the absence of Bach2. Therefore, Bach2 is a key checkpoint protein regulating NK terminal maturation.