Combined multidimensional single-cell protein and RNA profiling dissects the cellular and functional heterogeneity of thymic epithelial cells

Klein, Fabian; Veiga-Villauriz, Clara; Börsch, Anastasiya; Maio, Stefano; Palmer, Sam; Dhalla, Fatima; Handel, Adam E.; Zuklys, Saulius; Calvo-Asensio, Irene; Musette, Lucas; Deadman, Mary E.; White, Andrea J.; Lucas, Beth; Anderson, Graham; Holländer, Georg A.

doi:10.1038/s41467-023-39722-9

Download PDF

Article
Open access
Published: 10 July 2023

Combined multidimensional single-cell protein and RNA profiling dissects the cellular and functional heterogeneity of thymic epithelial cells

Nature Communications volume 14, Article number: 4071 (2023) Cite this article

3613 Accesses
1 Citations
8 Altmetric
Metrics details

Subjects

Abstract

The network of thymic stromal cells provides essential niches with unique molecular cues controlling T cell development and selection. Recent single-cell RNA sequencing studies have uncovered previously unappreciated transcriptional heterogeneity among thymic epithelial cells (TEC). However, there are only very few cell markers that allow a comparable phenotypic identification of TEC. Here, using massively parallel flow cytometry and machine learning, we deconvoluted known TEC phenotypes into novel subpopulations. Using CITEseq, these phenotypes were related to corresponding TEC subtypes defined by the cells’ RNA profiles. This approach allowed the phenotypic identification of perinatal cTEC and their physical localisation within the cortical stromal scaffold. In addition, we demonstrate the dynamic change in the frequency of perinatal cTEC in response to developing thymocytes and reveal their exceptional efficiency in positive selection. Collectively, our study identifies markers that allow for an unprecedented dissection of the thymus stromal complexity, as well as physical isolation of TEC populations and assignment of specific functions to individual TEC subtypes.

Single-cell transcriptional profiling of human thymic stroma uncovers novel cellular heterogeneity in the thymic medulla

Article Open access 17 February 2021

Developmental dynamics of two bipotent thymic epithelial progenitor types

Article Open access 25 May 2022

Integrative scATAC-seq and scRNA-seq analyses map thymic iNKT cell development and identify Cbfβ for its commitment

Article Open access 20 June 2023

Introduction

The thymus is essential for the formation and maintenance of the adaptive immune system. Its stroma provides a unique microenvironment promoting the generation and selection of T lymphocytes tolerant to an individual’s own tissue antigens yet responsive to an unlimited range of pathogens or malignantly transformed cells. Thymic epithelial cells (TEC) constitute the major cellular element of the stromal scaffold^1,2,3. Other cellular components of the stroma are different mesenchymal cell types and endothelial cells^2,4. TEC attract blood-borne lymphoid progenitors, commit them to a T cell fate, provide the molecular cues essential for expansion and differentiation, and shape the T cell antigen receptor (TCR) repertoire via stringent processes of positive and negative selection based on the cells’ antigen specificity^5,6,7.

The TEC compartment is composed of separate cortical (c) and medullary (m) lineages which have typically been defined by the cells’ anatomical location, a limited number of phenotypic markers and several functional characteristics^8,9,10. The surface marker Ly51 and reactivity to UEA1 are used to distinguish between cTEC (Ly51⁺UEA1⁻) and mTEC (Ly51⁻UEA1⁺). Markers such as CD80 and MHCII identify subsets of mTEC such as immature (CD80^loMHCII^lo; mTEC^lo) and mature epithelia (CD80^hiMHCII^hi; mTEC^hi). The latter cells are further differentiated based on the cells’ capacity to express the Autoimmune Regulator (Aire)^11,12. Recent single-cell RNA-sequencing (scRNAseq) uncovered a remarkable TEC heterogeneity which could previously not be appreciated using the few cell surface markers available for flow cytometry^13,14,15. For instance, a scRNAseq analysis of the TEC compartment of 4-week-old mice demonstrated a single cortical TEC type but 4 separate mTEC subtypes, namely immature and a mature mTEC, post-Aire mTEC and tuft-like mTEC¹⁵. Investigations of TEC heterogeneity across the life trajectory (1–52-weeks of age) identified 9 different TEC subtypes whose relative frequencies vary with age¹³. For example, the cTEC compartment is composed of at least two main subtypes, designated perinatal and mature cTEC¹³. Perinatal cTEC represent a major subpopulation early in life (~40% of all TEC in the first week after birth) but their relative frequency rapidly decreases thereafter with only a small fraction of these cells being detected in adult animals. Conversely, mature cTEC increase in frequency and represent the majority of cortical epithelia from 4 weeks of life onwards¹³. Intertypical TEC are characterised by a gene expression profile that includes signatures typical for both cortical and medullary thymic epithelial lineages. They also express genes including Pdpn, Ccl21a, Ly6a, and Plet1 that have previously been associated with mTEC, thought to have a progenitor potential, and localised at the cortico-medullary junction (CMJ)^{16,17,18,19,20}. However, only few of these transcriptionally defined TEC subsets can currently be assigned to any cytometrically characterised TEC subpopulation, (for clarity, we refer to transcriptionally defined TEC clusters as subtypes and cytometrically specified TEC as subpopulations). Moreover, the hitherto absence of suitable and informative cell surface markers to physically isolate most of the TEC subtypes for in vitro analyses and in vivo transfer studies has disallowed to date further functional characterizations of these cells and the physical establishment of direct precursor-progeny relationships at single-cell resolution.

To address this limitation, we sought to screen mouse TEC for the expression of 260 cell surface markers employing massively parallel flow cytometry and the Infinity Flow computational pipeline to infer a co-expression pattern for any of the tested epitopes²¹. This approach identified several novel TEC surface markers that when suitably combined identified perinatal cTEC, intertypical TEC and tuft-like mTEC which had previously only been classified either by their distinct RNA expression profiles or a combination of cell surface and intracellular markers^22,23. The identity of these phenotypically defined TEC subpopulations was verified by simultaneous measurements of mRNA and surface protein expression using CITEseq and, in the case of perinatal cTEC, further characterised functionally, spatially, and developmentally.

Results

Establishment of a cell surface expression atlas across thymic stromal cell subsets

To resolve thymic stroma heterogeneity at a phenotypic level, we sought to identify new cell surface markers that reliably and accurately identify TEC subsets hitherto only defined by the cells’ individual gene expression profiles. For this purpose, we used massively parallel flow cytometry for 260 individual cell surface markers followed by an analysis employing Infinity Flow, a computational machine learning algorithm²¹. Thymic stomal cells were isolated as single cells from 1-, 4-, and 16-week-old mice, physically enriched and subsequently stained for 12 backbone markers that either alone or in combination reliably identified haematopoietic (CD45), different epithelial (EpCAM1, Ly51, UEA1, MHCII, CD40, CD80, CD86, Sca1, AIRE, Podoplanin), endothelial (CD31) and some mesenchymal cells (Sca1, Ly51, Podoplanin; Fig. 1a). As a next step, cells were split into aliquots and stained individually for 260 exploratory markers (Supplementary Table 1). Following sample acquisition, Infinity Flow was used to impute the expression level of each of the exploratory markers tested at single-cell resolution²¹. The resultant predictions of expression are based on non-linear functions of the recorded backbone markers. The observed heterogeneity and co-expression patterns were further analysed and visualised by the single-cell analysis pipeline Seurat²⁴. Hierarchical clustering of the data resulted in 7 clusters for data drawn from 1-week-old mice and in 10 clusters for that of older animals, as illustrated in two dimensions by a Uniform Manifold Approximation and Projection (UMAP) (Supplementary Fig. 1a–c).

**Fig. 1: Infinity Flow analysis reveals TEC heterogeneity.**

At each of the three separate timepoints, the major thymic stromal cell types, epithelia, fibroblasts, pericytes, and endothelial cells, could reliably be identified based on the expression of key markers including EpCAM1 (CD326) identifying TEC, CD140a marking fibroblasts, Ly51 and CD146 singling out pericytes, and CD31 staining endothelial cells. Additional markers identified subpopulations within these cell clusters (Supplementary Fig. 1a–c). Several of the antibody specificities to detect backbone epitopes were also included among the selected 260 exploratory markers (e.g., EpCAM1, CD31, Ly51, and Sca1) which allowed direct comparisons between exploratory and identical backbone markers, thus verifying the utility of the Infinity Flow algorithm. For these markers, we noted highly similar expression profiles, therefore demonstrating the reliability of the computational approach taken (Supplementary Fig. 1d).

The initial expression analysis not only confirmed by flow cytometry the heterogeneity among thymic stromal cell types, but also revealed a dynamic change over time in the relative representation of individual TEC subpopulations (Supplementary Fig. 1a–c). In a second analysis, we focused exclusively on EpCAM1⁺ cells and disclosed in 1-week-old but not older mice three separate cTEC subclusters as defined by the cells’ differential expression of Ly51 and UEA1 thus illustrating a greater heterogeneity of the cTEC population early in life (Fig. 1b–d and Supplementary Fig. 2a–d). In mice 4 weeks of age and older, mTEC with a low surface expression of MHCII (designated mTEC^lo) segregated into 4 separate subclusters based on the differential expression of the surface markers analysed (Fig. 1b–d).

cTEC heterogeneity identified by differential cell surface marker expression

We next queried whether the expression of CD83, CD40, HVEM (CD270), and Ly51 unequivocally classified individual cTEC subpopulations, since their intensity profile differed across cTEC clusters identified in 1-week-old mice (Fig. 2a). The expression of CD40 and HVEM were exclusively restricted to a subcluster designated cTEC I (see below) whereas the two other markers were detected across all cTEC subclusters, but with a stronger signal on the cTEC I (Fig. 2a).

**Fig. 2: Surface expression profile of perinatal cTEC.**

Analysing a previously published scRNAseq dataset of TEC¹³, transcripts for Cd83, Cd40, and Enpep (encoding Ly51) were detected in perinatal cTEC, albeit at various levels (Fig. 2b, c). In contrast, transcripts for Tnfrsf14, the gene encoding HVEM, were detected in only a few TEC but across several clusters, thus failing to unequivocally identify perinatal cTEC. This finding highlighted the limitations of gene expression studies to identify surface markers that matched the cells’ RNA profile. To further assess the relationship of the cTEC I subcluster to TEC subtypes identified by scRNAseq, we generated a score of similarity using SingleR which related the RNA expression profile of individual cells to the computed cell surface expression pattern of cTEC I. This analysis demonstrated the highest similarity score to the pairing of cTEC I with perinatal cTEC (Fig. 2d).

We then aimed to define surface markers that allow the isolation of perinatal cTEC by flow cytometry. We used the presence of markers highly expressed in cluster cTEC I of 1-week-old mice while excluding markers detected on the majority of mature cTEC isolated from 4- to 16-week-old animals (Figs. 1b–d; 2a), as the population of cortical epithelia in older animals only includes perinatal cTEC at a very low frequency. Within the population of cTEC, we identified a subpopulation of cells that concomitantly expressed CD83 and CD40 but were Sca1 negative early in postnatal life (Fig. 2e; Supplementary Fig. 3a). As early as 4 weeks postnatally, TEC with a Sca1⁺ phenotype appeared among CD83⁺CD40⁺ cTEC. These cells were electronically excluded from further analysis as they represent mature cTEC that accumulate with age. The frequency of CD83⁺CD40⁺Sca1⁻ epithelia changed substantially during the life course as the cells’ relative representation progressively increased throughout organogenesis, plateaued in 1-week-old mice (4.6 × 10⁴ ± 1.4 × 10⁴ cells) and subsequently decreased, displaying the lowest representation in 8-week-old animals (2.4 × 10³ ± 2.3 × 10³ cells; Fig. 2e, f; Supplementary Fig. 3b).

Perinatal cTEC displayed in contrast to other cTEC subpopulations, comprising cTEC clusters II and III, higher cell surface levels for HVEM, Ly51 and MHCII which allowed the identification of these cells in combination with a high cell surface expression of CD83 and CD40, in the absence of Sca1. Perinatal cTEC also showed a higher Foxn1 promoter activity as demonstrated in reporter mice where the expression of GFP is under the transcriptional control of the Foxn1 locus (Fig. 2g).

Identification of intertypical and tuft-like TEC

The Infinity Flow analysis of adult mice revealed 4 distinct mTEC^lo clusters. Clusters I and II were observed in 4- and 16-week-old mice and displayed a similar expression profile for most of the 260 exploratory markers (Fig. 1b–d) including a shared expression of Sca1 and CD146 (Fig. 3a, b). To match these two phenotypically defined subpopulations to their corresponding transcriptome-determined TEC subtypes, we probed the RNA profiles of single TEC for the expression of Ly6a/Ly6e (encoding Sca1) and Mcam (encoding CD146). While transcripts for Ly6a/Ly6e were detected especially among intertypical TEC¹³, Mcam-specific RNA was only detected at low levels and in different TEC, but mostly within the intertypical TEC subtype (Fig. 3c). This subtype is characterised by transcriptional features characteristic of both cTEC and mTEC —as phenotypically defined by the conventional surface marker Ly51 and UEA-reactivity— contribute to this unique TEC subtype¹³. Using again the SingleR package, the similarity scores for both 4- and 16-week-old mTEC^lo I and II were calculated to be the highest when matched to the intertypical TEC subtype (Fig. 3d; Supplementary Fig. 4a). The mTEC^lo subpopulation contained cells that co-expressed Sca1 and CD146 and cells with this phenotype increased with postnatal age (Fig. 3e, f; Supplementary Fig. 4b). Although initially only detected among mTEC^lo this subpopulation was increasingly also observed within the cTEC compartment of mice older than 4 weeks of age (Fig. 3e, f).

**Fig. 3: Identification of intertypical TEC within cTEC and mTEC.**

A unique set of surface markers that specifically recognised the mTEC^lo cluster III could not be found. However, the simultaneous expression of CD66a and CD117 in the absence of Sca1 and CD63 positivity identified mTEC^lo cluster IV. This cluster was only detected in adult mice (Fig. 4a; Supplementary Fig. 5a), although the defining 4 cell surface markers could also be detected in the mTEC^lo cluster II of 1-week-old animals (Supplementary Fig. 5a). It is therefore possible that cluster mTEC^lo II of 1-week-old mice represents epithelia that form the separate cluster mTEC^lo IV in older animals. The single-cell transcriptomic analysis only partially matched the phenotypic analysis of mTEC^lo clusters as Ly6a/Ly6e- and Cd63-specific RNA could indeed be detected in the vast majority of TEC (Figs. 3c and 4b) whereas transcripts for Ceacam1 (encoding CD66a) and Kit (encoding CD117) were largely absent in these cells (Fig. 4b). The similarity score revealed the best match between mTEC^lo cluster IV and post-AIRE and tuft-like mTEC (Fig. 4c; Supplementary Fig. 5b).

**Fig. 4: A combination of surface markers to define tuft-like mTEC.**

We next sought to define a phenotypic profile of cells belonging to cluster mTEC^lo IV that would allow their physical isolation by flow cytometry. For this purpose, we identified within the Sca1⁻CD63⁻ mTEC^lo a subpopulation of cells that stained positively for both CD66a and CD117, thus mirroring the features identified for mTEC^lo cluster IV (Fig. 4d). Most of the Sca1⁻CD63⁻CD66a⁺CD117⁺ mTEC^lo cells (~70%) also expressed the serine/threonine-protein kinase Dclk1 (Fig. 4e) which was previously identified as a typical intracellular marker for tuft-like mTEC^{15, 23}. While staining for the tuft-like mTEC marker L1CAM was not successful in our hands, both the Dclk1 negative and positive fractions shared the CD104^lo tuft-like mTEC phenotype as previously described (Supplementary Fig. 5c)¹⁵. Furthermore, the presence of Dclk1 negative cells within tuft-like mTEC is supported by previous scRNAseq results^13,22. We also observed that an absence of Dclk1 in Sca1⁻CD63⁻CD66a⁺CD117⁺ mTEC^lo correlated with a lower surface expression of CD66a and CD117, indicating that these cells are not yet fully differentiated into tuft-like mTEC (Supplementary Fig. 5d). Conversely, the vast majority of Dclk1-positive TEC were detected among Sca1⁻CD63⁻CD66a⁺CD117⁺ mTEC^lo (Supplementary Fig. 5e). In the absence of the transcription factor Pou2f3 Dclk1 expression was absent and Sca1⁻CD63⁻CD66a⁺CD117⁺ mTEC^lo cells were not detected (Fig. 4f, g) thus confirming their identity as tuft-like thymic epithelia^22,23. Using these phenotypic features, we noted the presence of tuft-like mTEC to change over time with the highest frequency and cellularity in 4-week-old mice (Fig. 4h, i). As CD66a⁻CD117⁻ cells (“non-tuft-like”) made up a substantial portion (~40% at 4 weeks) of the relatively few Sca1⁻CD63⁻ mTEC^lo, we sought to further determine their identity via bulk RNAseq. Specifically, we investigated whether these cells shared any transcriptomic features characteristic of tuft-like mTEC. This analysis revealed major transcriptional differences between the non-tuft-like and tuft-like cells. Genes associated with tuft-like cells were enriched in the CD66a⁺CD117⁺ fraction of cells and included Dclk1, Il25, Ceacam1, and Kit (Supplementary Fig. 5f). Similarly, the top genes identified in a previous scRNAseq experiment¹³ defining tuft-like mTEC were expressed at a very high level in these tuft-like cells whilst their transcripts were not present in non-tuft-like cells (Supplementary Fig. 5g). Correlation of the entire non-tuft-like transcriptome to the annotated TEC subsets confirmed this delineation from cells identified as tuft-like mTEC (Supplementary Fig. 5h). Because the CD66a⁻CD117⁻ non-tuft-like mTEC displayed a gene expression profile that did not match any of the known mTEC subpopulations it is likely that these rare cells (~4% of total TEC) represent a mixture of mTEC subpopulations.

Tuft-like mTEC originate from mTEC that have once expressed the tissue restricted antigen Csnb¹⁵. We utilised Csnb^Cre::Rosa26^LSL-YFP reporter mice¹⁵ that allow in vivo fate mapping within the mTEC linage to test whether Sca1⁻CD63⁻CD66a⁺CD117⁺ mTEC^lo originate from a Csnb-expressing precursor. In keeping with the previous study, we identified 70–80% of the Sca1⁻CD63⁻CD66a⁺CD117⁺ mTEC^lo mTEC to be YFP labelled (Supplementary Fig. 5i), suggesting that they represent bona fide tuft-like mTEC. The mTEC^lo compartment is composed of medullary epithelia that either have not yet expressed the transcriptional facilitator AIRE or, alternatively, belong to a group of cells that have differentiated from AIRE-positive, mature mTEC. Tuft-like mTEC have previously been shown to derive, at least in part, from AIRE-positive precursors²² and be enriched for expression of the surface glycoprotein Tspan8, an AIRE-enhanced tissue-restricted antigen^14,25. We therefore tested whether Sca1⁻CD63⁻CD66a⁺CD117⁺ mTEC^lo (i.e. tuft-like) mTEC are positive for the expression of Tspan8. As many as 60% of tuft-like mTEC expressed Tspan8, further validating their identity and identifying these cells to contain post-AIRE mTEC (Supplementary Fig. 5j). The combined expression of YFP and Tspan8 was only detected in a small fraction of Sca1⁺CD146⁺ mTEC^lo intertypical TEC (Supplementary Fig. 5i,j), suggesting they belong to a TEC developmental stage that does not yet promiscuously express tissue-specific antigens.

CITEseq validates novel TEC markers

We next concurrently determined the abundance of surface proteins and mRNA expression in individual thymic stromal cells using Cellular Indexing of Transcriptomes and Epitopes by sequencing (CITEseq)²⁶. This approach was also taken to confirm unequivocally that the combined use of the newly described cell surface markers identifies specific and previously defined TEC subtypes. For this purpose, CD45⁻Ter119⁻ cells were isolated from thymi of 1- and 16-week-old mice and stained with oligonucleotide-coupled antibodies each directed against EpCAM1, MHCII, UEA1, CD80, CD86, CD40, CD83, HVEM, CD73, Sca1, CD63, CD117, CD200, CD54, CD49a, CD274, CD9, Ly6G/Ly6C, Podoplanin, or CD31, respectively. Labelled cells were further subjected to single-cell sequencing along with estimating the abundance of antibody-derived tags (ADTs)²⁶.

Cell clustering resulted in the identification of 12 clusters with data computed either from single-cell gene expression profiling or, alternatively, from ADT, and displayed by means of t-distributed stochastic neighbour embedding (tSNE) (Supplementary Fig. 6a–c). The cell-type annotation performed was based on the gene expression profiles derived from the Immunological Genome Project (ImmGen)²⁷ and confirmed the identity of individual clusters as endothelial cells, fibroblasts, stromal cells, and epithelial cells, respectively. Hence, the chosen combination of selected surface markers was sufficient to identify individual stromal cell types (Supplementary Fig. 6d, e).

Analysing the captured gene expression profiles of only TEC identified 8 separate clusters (A-H), whereas examining the ADT data recognised 9 (1–9) clusters (Fig. 5a, b). Comparing the two approaches demonstrated a nearly pairwise relationship (Fig. 5c). The few exceptions observed concerned on one hand the clusters D and E which represented a mixture of clusters 4 and 5 and on the other hand cluster G which split into clusters 7 and 8 as a result of the clusters’ differential expression of Ly6C/Ly6G (Fig. 5d). The limited number of antibodies used in the CITEseq analysis identified three cTEC cluster (defined as UEA1 non-reactive cells:1 corresponding to cluster A [1/A], 2/B, and 3/C), three mTEC^lo (MHCII^loCD80^lo: 4/D, 5/E, and 9/H) and three mTEC^hi subpopulations (MHCII^hiCD80^hi: 6/F, 7/G, and 8) (Fig. 5d–f and Supplementary Fig. 6f).

**Fig. 5: CITEseq validates new TEC markers.**

We next assigned CITEseq-defined TEC cluster identities to those subtypes we have previously classified¹³. Clusters 1/A and, to a lesser extent, 2/B corresponded to perinatal cTEC. Cluster 2/B was further related to mature cTEC and cluster 3/C to mature cTEC and intertypical TEC (Fig. 5d, f). Clusters 4/D and 5/E related to intertypical TEC while cluster 6/F was most similar to both mature mTEC and proliferating TEC. Clusters 7/G and 8/G displayed a high similarity to mature mTEC (Supplementary Fig. 6f) whereas cluster 9/H linked to tuft-like mTEC. Notably, the transcriptional signatures characteristic of post-Aire mTEC and neural (n) TEC could not be detected. Furthermore, we analysed the data for the expression of gene transcripts characterising the recently described mimetic mTEC subsets (designated tuft, microfold, enterocyte/hepatocyte, neuroendocrine, ciliated, ionocyte, keratinocyte, and muscle cells)²⁸. Transcripts identifying “tuft cells” were detected in cluster 9/H (Supplementary Fig. 7) whereas mRNA signatures characterizing microfold-, enterocyte/hepatocyte-, and keratinocyte-like TEC were detected in clusters G/7-8 within the mature mTEC population. Genes associated with neuroendocrine, ciliated, ionocyte, and muscle cells among specified TEC subsets could not be robustly identified.

Replicating the changes of specific TEC subtypes with age, clusters 1/A and 2/B were more abundant in 1-week-old mice whereas the frequencies of clusters 3/C, 4/D, 5/E, 7/G, and 9/H were increased in 16-week-old animals (Fig. 5d; Supplementary Fig. 6g).

TEC in cluster 1/A displayed the highest CD83, CD40, and HVEM protein expression among CITEseq defined clusters, thus confirming the cells’ identity as perinatal cTEC. The expression of these markers was reduced in cluster 2/B and completely absent in cluster 3/C, suggesting the former to represent a developmentally intermediate cell state between perinatal and mature/intertypical-like cTEC (Fig. 5d, g). Furthermore, cTEC maturation was paralleled by a decrease in Foxn1 transcription and an increase in CD73, CD49a, and Sca1 protein expression (Supplementary Fig. 6h).

The differential CD117, CD63, and Sca1 protein expression (as measured by ADT) identified cluster 9/H as tuft-like mTEC (CD117⁺CD63⁻Sca1⁻; see above and Fig. 5d, h) and thus confirmed the flow cytometric definition and gating strategy used to identify these cells as both accurate and practical (Fig. 4g–i). This conclusion was further corroborated by the detection of Dclk1 and Ceacam1 transcripts in cluster 9/H (Fig. 5h) and a high similarity score with the tuft-like mTEC subtype (Fig. 5d, f)¹³.

ADT-based detection of Sca1 protein expression matched to cells with a transcriptional signature of intertypical TEC within the cTEC (cluster 3/C) and mTEC^lo subpopulations (clusters 4 + 5/D +E; Fig. 5d, f, h). Hence, intertypical TEC could unequivocally be identified by Sca1 expression alone. As the transcriptional signature identifying intertypical TEC was spread across three CITEseq-defined clusters (3/C, 4/D and 5/E; Fig. 5d, f), the detection of CD146 expression appeared to deconvolute TEC heterogeneity further since fractions of Sca1⁺ cTEC and Sca1⁺ mTEC^lo stained positively for CD146⁺ (Fig. 3e).

The ADT-based documentation of surface markers identified individual TEC subtypes. However, the corresponding gene expression profiles were on their own insufficient to recognize these cells, not least because of the occasional discrepancy between surface protein and RNA expression (Supplementary Fig. 6h, i). CITEseq could therefore validate the utility of the selected, novel surface markers and the gating strategies chosen. Together they identified four TEC subpopulations that correspond to a specific transcriptionally defined cluster and two subpopulations that represent a mixture of two related clusters, namely UEA1⁻CD83⁺CD40⁺Sca1⁻ perinatal cTEC (cluster 1/A), UEA1⁻CD83⁻CD40⁻Sca1⁻ mature cTEC (2/B), UEA1⁻CD83⁻CD40⁻Sca1⁺ intertypical-like cTEC (3/C), UEA1⁺MHCII^loCD80^loSca1⁺ intertypical mTEC (4 + 5/D + E), UEA1⁺MHCII^hiCD80^hi mature mTEC (6 + 7 + 8/F + G), and UEA1⁺MHCII^loCD80^loSca1⁻CD63⁻CD66a⁺CD117⁺ tuft-like mTEC (9/H) (Supplementary Fig. 8).

Perinatal cTEC present an enhanced potential for positive selection

We next sought to localize perinatal cTEC within the thymus stromal architecture. Because Ly51 abundance was higher on perinatal cTEC in comparison to other cortical epithelial populations (Fig. 2g), we used this differential to localize perinatal cTEC on thymus tissue sections (Fig. 6a). Quantification of the Ly51 signal intensity in immunohistology detected these cells in close proximity to the medulla with gradual increase of the Ly51 signal but invariable cytokeratin 8 (K8) staining across the cortex from subcapsular region to the inner and eventually deep cortex (Fig. 6a, b).

**Fig. 6: Perinatal cTEC display an increased capacity for positive selection.**

The newly assigned surface markers enabled the isolation and functional testing of individual cTEC populations ex vivo. We therefore investigated the capacity of perinatal (CD83⁺CD40⁺Sca1⁻) and non-perinatal (CD83⁻CD40⁻) cTEC to effect positive thymocyte selection. We co-cultured these cells in thymic epithelial cell cultures (TECx) together with CD69⁻CD4⁺CD8⁺ (i.e., pre-selection double-positive) thymocytes for two days (for gating see Supplementary Fig. 9). Thymocytes were then analysed for phenotypic features associated with positive thymic selection, i.e., the upregulation of TCR and CD5 (Fig. 6c). The number of total thymocytes and those with a CD5^hiTCRβ^hi phenotype were significantly increased in TECx composed of perinatal cTEC when compared to aggregates composed of other cortical epithelia (Fig. 6d). Taken together, these results identified perinatal cTEC to be juxtaposed to the medulla and particularly efficient in positively selecting thymocytes.

The number of perinatal cTEC significantly decreases with age¹³. We therefore explored whether this variation was paralleled by a change in the efficiency to impose positive thymocyte selection. We thus monitored and compared thymocyte maturation in 4- and 16-week-old thymi and classified their sequential maturational stages according to TCRβ and CD69 expression (i.e., stage 0: TCRβ⁻CD69⁻ → stage 1: TCRβ⁺CD69^+/− → stage 2: TCRβ⁺CD69⁺ → stage 3: TCRβ⁺CD69⁻). The frequency of pre-selection thymocytes (i.e., stage 0) was increased, whereas the relative abundance of cells with a post-selection phenotype (stage 2 and 3) was significantly reduced in older animals (Fig. 6e, f). These in vivo results indicated a compromised capacity of older mice to positively select thymocytes, which correlated with a decrease in the availability of perinatal cTEC (Fig. 2e, f).

Crosstalk with thymocytes induces maturation of perinatal cTEC

We finally investigated whether thymic crosstalk^8,29 could explain the inverse correlation between the decreased frequency of perinatal cTEC with age and the expansion of thymocytes after birth. We therefore first determined the frequency of perinatal cTEC in Rag2-deficient (Rag2^−/−) mice, which have a hypoplastic thymus secondary to a thymocyte developmental arrest at the DN3a stage. We found a high fraction of perinatal cTEC in these mice that was not influenced by age (Fig. 7a). Hence, thymocytes at developmental stages up to the beta-checkpoint did not influence the age-related changes in perinatal cTEC frequencies.

**Fig. 7: Crosstalk with thymocytes facilitates cTEC maturation.**

To probe whether thymocytes at later developmental stages, especially unselected CD4⁺CD8⁺ (double positive, DP) thymocytes, controlled the frequency of perinatal TEC, Rag2^−/− mice were injected with antibodies directed against CD3ε. This treatment results in a substantial increase in pre-selection DP thymocytes^30,31,32. Four weeks after antibody or control injections, the thymus of actively treated Rag2^−/− mice contained an abundance of DP thymocytes which correlated with numerical and phenotypic changes in the cTEC compartment (Fig. 7b–e). The latter were marked by a reduction in perinatal cTEC, parallel to an increase in non-perinatal cTEC, specifically Sca1⁺ cells (Fig. 7d, e) which corresponded to intertypical TEC according to our CITEseq data (Fig. 5d, f, h). Taken together, these results identified the abundance of and/or signalling by pre-selection DP thymocytes as the mechanism by which the frequency of perinatal cTEC was controlled.

Discussion

Single-cell transcriptomic analyses have uncovered an unexpected heterogeneity within many cell populations of a seemingly identical phenotype. Cells of the thymic stromal compartment constitute no exception to this observation^{4, 13,14,15,22}. The apparent lack of suitable cell surface markers identifying unequivocally TEC subpopulations identical to individual TEC subtypes precludes the isolation of live TEC and their ex vivo functional analysis. Here we report that this limitation has been substantially overcome. We describe novel cell surface markers that identify the comparable subtype of scRNAseq-defined perinatal, mature and intertypical cTEC, and mature, intertypical and tuft-like mTEC.

Predicting the cell surface phenotypes from corresponding scRNAseq profiles is challenging as technical limitations detecting low transcript copy numbers and the acknowledged disparity between transcript detection and protein expression render this attempt difficult. For example, a comparison of 7 scRNAseq methods uncovered that high-throughput methods, including the widely used 10x Chromium, have lower sensitivities in comparison to the low-throughput methods Smart-seq2 and CEL-Seq2 when capturing rare transcripts^33,34. We therefore opted for an alternative method and stained TEC for the expression of hundreds of cell surface markers. This screening approach of massive parallel flow cytometry combined with Infinity Flow analysis discovered surface markers previously not inferred to be expressed by TEC. CITEseq which combines the detection of promising candidate markers and single-cell transcriptomic profiles finally established the accuracy of the cell surface markers chosen to identify TEC subtypes.

By applying the newly identified surface markers we show that perinatal cTEC are enriched within the TEC scaffold towards the cortico-medullary junction (CMJ) and that these cells are particularly efficient in positively selecting maturing thymocytes. As expected from their role in shaping the TCR repertoire, we find perinatal cTEC typically juxtaposed to thymocytes with an activated phenotype (i.e., CD69⁺). The age-dependent decline in the frequency of perinatal cTEC is noted both when using flow cytometry and scRNAseq to classify these cells. The actual pace by which this regression is observed differs, however, between these two methods, demonstrating a seemingly faster kinetic for perinatal cTEC identified by their RNA expression profile¹³. This may be explained by differences in the half-lives of specific transcripts and their corresponding proteins. Nonetheless, a post-natal decrease of these cells to an almost complete absence early in adulthood is expected to compromise the robustness of thymopoiesis and possibly the efficiency by which thymocytes are positively selected. The observed decrease in perinatal cTEC correlates with other compositional changes within the epithelial scaffold and may constitute an intrinsic driver for thymus senescence. This understanding is consistent with scRNAseq data that shows the quiescence of a population of medullary precursor cells and correlates this alteration with an impaired maintenance of the medullary TEC compartment¹³. In parallel, these age-related changes link to less efficient T cell selection, a decreased self-antigen representation, a decreased TCR repertoire diversity, and a reduced frequency of thymus-resident naïve T cells (this report and ref. ¹³).

Intertypical TEC express the glycosyl phosphatidylinositol-anchored cell surface protein Sca1 independently whether they are positive for the cortical marker Ly51⁺ or reactive with the lectin UEA-1, a general feature of medullary TEC. Here, we now show that intertypical TEC can be further split into two sizable subpopulations based on their expression of the cell adhesion molecule CD146. Because oligonucleotide labelled anti-CD146 antibodies were not available for the CITEseq analysis and transcripts for this marker are typically lowly expressed among intertypical TEC, a distinction between CD146 negative and positive cells was not possible when analysing scRNAseq data. Intertypical TEC may contain progenitors with a developmental bias towards the mTEC lineage¹³. Indeed, a recent report provides further support of this contention since the gene expression profile of intertypical TEC largely overlaps with a heterogeneous progenitor population which has been claimed to act as mTEC biased postnatal TEC progenitor²⁰. Our profiling of Tspan8 expression and lineage tracing furthermore suggest that the majority of Sca1⁺CD146⁺ intertypical TEC relate to immature mTEC that have not yet fully differentiated to express collectively a broad range of tissue restricted antigens.

We further specify the cell surface phenotype for tuft-like mTEC (L1CAM⁺CD104⁺¹⁵). These cells share transcriptional (e.g. expression of IL25, Trmp5, Dclk1, and IL17RB) and morphological characteristics with gut epithelial tuft cells^13,15,22, play a function in central T cell tolerance induction²² and control both the homeostasis of type 2 innate lymphoid cells and the generation of type 2 natural killer T cells^15,22,23. The lineage tracing of these cells further shows that the majority of tuft-like mTEC derive from AIRE expressing TEC, a finding in keeping with previously reported observations^15,22. However, whether all tuft-like mTEC differentiate from mature mTEC remains an open question because the labelling method needed to draw this conclusion (i.e., inducible, Aire-dependent tracing) may not be completely effective²². Interestingly, the Csnb lineage tracing identified an increased frequency of labelled tuft-like cells in comparison to mature mTEC where labelling is initiated. While we have no unequivocal explanation for this increase, we nevertheless conclude that the majority of tuft-like mTEC (at least 60%) are the progeny of mature medullary epithelia. There is, however, room to speculate that all tuft-like mTEC may be derived in this way, since any contributions from another Csnb-negative (i.e., non-labelled) precursor would dilute the frequency of labelled tuft-like cells, a result that we did not observe.

To probe the utility of the new set of cell surface markers to phenotype altered TEC scaffolds, we next analysed the composition of the thymic epithelia in FOXN1^Δ505/WT mice which express a dominant negative mutation of FOXN1 and consequently show substantial defects in TEC differentiation³⁵. A previous scRNAseq-based analysis of these animals revealed a relative enrichment of perinatal and mature cTEC against a reduction of tuft-like mTEC whereas the frequency of intertypical TEC remained unchanged. The flow cytometric analysis of the TEC scaffold in FOXN1^Δ505/WT mice identifies the same variations and thus maps accurately to the transcriptional analysis of these cells, thus demonstrating that the phenotypic and gene expression-based analyses draw comparable conclusions (Supplementary Fig. 10).

With the approach taken, five of the previously defined nine TEC subtypes¹³ can now be unequivocally identified using cell surface markers. The still small number of discriminatory cell surface markers so far identified likely accounts for this minor limitation. Implementing more markers in the screening process may identify additional cell surface markers that will identify the remaining TEC subtypes for which we have not yet identified an unambiguous cell surface marker profile. Alternatively, the use of intracellular markers may be informative in identifying the remaining TEC subtypes, namely proliferating TEC, post-AIRE mTEC, and nTEC. However, an obvious drawback for this approach will be that TEC identified in this fashion will be non-viable and can therefore not be used for functional studies in vitro or after transfer in vivo.

It is important to establish precursor-progeny relationships for specific TEC subtypes now that we can identify and purify specific subpopulations. For example, the potential can be tested whether CD146⁺ intertypical TEC give rise to mature mTEC, competent to effect negative selection. Another effort could be directed in dissecting the molecular requirements of tuft-like mTEC controlling the development of type 2 lymphoid cells employing ex vivo functional assays. Finally, the screening workflow described here will also be valuable in identifying novel biomarkers apt to monitor changes in cell subpopulations and their functions resulting from spontaneous or engineered changes in gene function.

Taken together, we have identified novel surface markers that enable the isolation and functional assessment of different TEC subpopulations that correspond to previously identified subtypes so far only defined by their transcriptome. This is accomplished by combining a high throughput screening workflow with a computational expression projection followed by unsupervised clustering.

Methods

Mice

Animals were maintained under specific pathogen–free conditions at the University of Oxford Biomedical Science facilities under local and United Kingdom Home Office regulations and permissions. They were housed in cages grouped to a maximum of five animals per cage on wood litter. Environmental enrichment such as tunnel, gnawing and nesting material was provided. The animal house was maintained under artificial lighting (12 h) between 7:00 am and 7:00 pm, in a controlled ambient temperature of 22 °C ± 2 °C, and relative humidity between 30% and 70%. Experiments performed were approved by the United Kingdom Home Office regulations and co-housed age- and gender-matched wild-type C57BL/6 mice were used in all experiments as a reference for genetically modified animals.

Csnb^Cre mice¹⁵ were crossed to the Rosa26^YFP mouse line³⁶ to induce lineage tracing in the mature mTEC compartment.

Mice heterozygous for a Foxn1 allele with a single nucleotide loss at position 1470 (designated FOXN1^Δ505/WT) were generated at the Genome Engineering Facility of the MRC Weatherall Institute of Molecular Medicine, University of Oxford as previously described³⁵.

Rag2^−/− mice were bred and maintained in the mouse facility of the Department of Biomedicine at the University of Basel in accordance with permissions and regulations of the Cantonal Veterinary Office of Basel-Stadt.

For timed pregnancies 7- to 14-week-old mice were mated over-night and separated early next morning. For pregnant females the mating was considered E0.5 that morning.

TEC isolation

Isolated thymi were cleaned from adipose tissue, separated into the two lobes, and subsequently subjected to three rounds of enzymatic digestion with Liberase (2.5 mg/ml, Roche, Cat no: 5401127001) and DNaseI (10 mg/ml, Roche, Cat no: 10104159001) diluted in PBS (Gibco, Cat no: 70011044) at 37 °C. After filtration through a 100-μm cell strainer and resuspension in FACS buffer (PBS supplemented with 2% FBS), cell number was determined using a CASY cell counter (Innovatis). For most analyses CD45⁺ hematopoietic cells were depleted by incubation with anti-CD45 beads (Miltenyi) as per manufacturer’s recommendations and subsequently subjected to the AutoMACS separator (Miltenyi) “depleteS” program.

Flow cytometry

Cells were counted and stained in FACS buffer containing antibodies of interest for 30 min at 4 °C in the dark. For the identification of dead cells an additional staining with propidium iodide (PI, Sigma, Cat no: P4864) or Zombie red (Biolegend, Cat no: 423110) was used. For intra-cellular staining, cells were fixed and permeabilised after cell-surface staining using the Cytofix/Cytoperm (BD Biosciences, Cat no: 554714) or the Transcription Factor Staining Buffer Set (Invitrogen, Cat no: 00-5523-00) according to the manufacturer’s protocol. Cells were analysed and sorted on a BD FACSAria III instrument (BD Biosciences). Cells were sorted into FACS buffer. Cell purities of at least 95% were confirmed by post-sort analysis. The following antibodies were used: CD4-APCCy7 (1:400, GK1.5, Biolegend), CD4-FITC (1:1000, GK1.5, Biolegend), CD5-PerCPCy5.5 (1:400, 53-7.3, Biolegend), CD8α-AF700 (1:800, 53-6.7, Biolegend), CD40-FITC (1:200, 3/23, Biolegend), CD45-AF700 (1:400, 30-F11, Biolegend), CD63-PE (1:400, NVG-2, Biolegend), CD66a-APC (1:400, Mab-CC1, Biolegend), CD66a-FITC (1:400, Mab-CC1, Biolegend), CD69-PECy5 (1:800, H1.2F3, Biolegend), CD69-PECy7 (1:200, H1.2F3, Biolegend), CD73-BV421 (1:400, TY/11.8, Biolegend), CD80-PECy5 (1:2000, 16-10A1, Biolegend), CD80-BV605 (1:400, 16-10A1, Biolegend), CD83-Bio (1:200, Michel-19, Biolegend), CD83-PE (1:400, Michel-19, Biolegend), CD104-FITC (1:400, 346-11A, Biolegend), CD117-BV421 (1:200, 2B8, BD Biosciences), CD146-APC (1:800, ME-9F1, Biolegend), EpCMA1-PerCPCy5.5 (1:800, G8.8, Biolegend), Dclk1 (1:1000, ab31704, Abcam), HVEM-APC (1:400, HMHV-1B18, Biolegend), Ly51-PECy7 (1:400, 6C3, Biolegend), MHCII-APC/Fire750 (1:2000, M5, Biolegend), UEA1-Cy5 (1:400, Vector Laboratories, in-house labelling), Sca1-BV510 (1:800, D7, Biolegend), TCRβ-FITC (1:400, H57-597, Biolegend), TCRβ-PE (1:1000, H57-597, Biolegend), Tspan8-APC (1:400, FAB6524A, R&D Systems). Biotinylated antibodies were detected using Streptavidin-BV785 (1:500, Biolegend) and unlabelled Dclk1 using anti-rabbit IgG AF647 (1:4000, Invitrogen). Data was analysed using FlowJo (version 10).

Massively parallel flow cytometry

Cells were isolated and CD45 depletion plus backbone staining were performed as described. The surface backbone panel included antibodies directed against CD45 (1:400, 30-F11, AF700, Biolegend), EpCAM1 (1:800, G8.8, PerCPCy5.5, Biolegend), Ly51 (1:400, 6C3, PECy7, Biolegend), MHCII (1:2000, M5, BV510, Biolegend), CD40 (1:200, 3/23, PECy5, Biolegend), CD80 (1:400, 16-10A1, BV605, Biolegend)), CD86 (1:800, GL-1, BV650, Biolegend), Sca1 (1:1000, D7, BV785, Biolegend), Podoplanin (1:200, 8.1.1, APC, Biolegend), CD31 (1:1000, 390, FITC, Biolegend), the Ulex europaeus agglutinin I (UEA1) lectin labelled with biotin (1:1000, Vector Laboratories), followed by secondary streptavidin-BV421 (1:1000, Biolegend) staining and Zombie red (1:1000, Biolegend) staining. Subsequently, the stained cells were distributed across the three 96-well plates provided with the LEGENDScreen kit (Biolegend, Cat no: 700009), each well containing a unique PE-labelled exploratory antibody as well as isotype controls and blanks. PE-labelled antibodies targeting GP2, Tspan8, CD177 and F3 were used as additional exploratory surface antibodies. Due to the low cell numbers obtained after CD45 depletion only ¼ of the recommended quantity of exploratory antibodies was used. Plates were incubated at 4 °C for 30 min in the dark. Thereafter, fixation was performed using the Cytofix buffer (BD Biosciences, Cat no: 554714) for 1 h at 4 °C in the dark. As an additional backbone marker, cells were stained intracellularly for anti-AIRE (1:400, 5H12, AF750, Invitrogen) in Cytoperm buffer (BD Biosciences, Cat no: 554714) and one well stained with anti-FOXN1³⁷ (1:1600, 2/41, PE, kind gift from Hans-Reimer Rodewald) as an additional exploratory marker, over-night at 4 °C in the dark. The next day cells were resuspended in 100 μl FACS buffer before analysis.

Infinity Flow and single-cell clustering and expression analysis

For the Infinity Flow computational analysis of the LEGENDScreen datasets, the acquired fcs files were gated on CD45 negative cells or specifically on EpCAM1⁺ TEC using the FlowJo software. The newly exported fcs files were then used as the dataset for the Infinity Flow pipeline as recently published²¹. The augmented data matrices generated during this process were then further analysed using the Seurat package for hierarchical clustering of the cells and differential expression analysis²⁴, mostly following the workflow presented in https://satijalab.org/seurat/articles/pbmc3k_tutorial.html. We used default parameters, except for the data normalization method, normalization.method = “CLR”. Values below zero were set to zero to allow for log normalization. Markers were filtered by hand to exclude T-cell related and focus on stromal cell related genes (Supplementary Table 1). The top 20 PCA dimensions were used for clustering and UMAP projections and a clustering resolution of 0.275 was used.

We compared the Infinity Flow data matrices with the scRNAseq dataset of reference ¹³ by identifying the most closely related genes for each Infinity Flow protein. However, some of the antibodies bind to protein complexes, and here we chose the most abundant transcript related to such a complex—for example, we chose H2-Ab1 RNA transcripts for MHC class II protein detection. Furthermore, UEA1 detection was identified via Fut1 RNA expression, since FUT1 synthesises the glycan target of UEA1; for the complete assignment see Supplementary Table 2. We compared the Infinity Flow fluorescence values with the scRNAseq normalised log counts. Clusters from each dataset were then compared using the SingleR package in R³⁸, with the Wilcox ranked sum test (using the SingleR option de.method = “wilcox”).

Histological analyses

Frozen thymus tissue sections (7 µm) were fixed in acetone and stained using antibodies specific for CD69 (1:100, H1.2F3, Biolegend), Ly51 (1:200, 6C3, Biolegend), K8 (1:500, TROMA-1, NICHD supported Hybridoma Bank), K14 (1:500, Poly19053, Biolegend). Images were acquired using a Leica DMi8 microscope.

Thymic epithelial cell cultures (TECx)

Perinatal cTEC (CD45⁻EpCAM1⁺MHCII⁺Ly51⁺CD83⁺CD40⁺) and non-perinatal cTEC (CD45⁻EpCAM1⁺MHCII⁺Ly51⁺CD83⁻CD40⁻) were sorted from the thymi of 2-week-old C57BL/6 mice and put in co-cultures with CD69⁻ DP thymocytes sorted from the same thymi, respectively. Cells were transferred in a 1:1 TEC to DP ratio into 1.5 mL tubes containing 1 mL Iscove’s modified Dulbecco’s medium (IMDM, Gibco, Cat no: 12440053) supplemented with 10% FBS, 100 units/mL penicillin and 100 μg/mL streptomycin (Sigma, Cat no: P4333-100ML) and 1× GlutaMAX supplement (Gibco, Cat no: 35050061). Co-cultures were maintained at 37 °C in a humidified atmosphere containing 10% CO₂ for 48 h and then analysed by FACS. As a control DP cells were also cultured without the addition of TEC.

Anti-CD3 injections

6- to 7-week-old Rag2^−/− animals were injected intraperitoneally with 50ug of anti-CD3ε (clone KT3, in-house produced) or HBSS. Four weeks post injection thymi were analysed for the appearance of DP thymocytes and for changes within their cTEC compartment.

Cellular indexing of transcriptomes and epitopes by sequencing (CITEseq)

Cells were isolated from six thymi of 1-week- and three thymi of 16-week-old C57BL/6 mice and depleted of CD45⁺ cells by AutoMACS. Subsequently cells were stained for CD45 (1:400, 30-F11, AF700, Biolegend), EpCAM1 (1:800, G8.8, PerCPCy5.5, Biolegend), Ly51 (1:400, 6C3, PECy7, Biolegend), Ter119-FITC (1:400, TER119, Biolegend) and with PI. In addition cells were stained with antibodies coupled to oligonucleotides directed against CD9 (MZ3, Totalseq-A, Biolegend), CD40 (3/23, Totalseq-A, Biolegend), CD49a (HMα1, Totalseq-A, Biolegend), CD54 (YN1/1.7.4, Totalseq-A, Biolegend), CD63 NVG-2, Totalseq-A, Biolegend), CD73 TY/11.8, Totalseq-A, Biolegend), CD83 (Michel-19, Totalseq-A, Biolegend), CD117 (2B8, Totalseq-A, Biolegend), CD146 (human with cross reactivity to mouse) (P1H12, Totalseq-A, Biolegend), CD200 (OX-90, Totalseq-A, Biolegend), CD274 MIH6, Totalseq-A, Biolegend), HVEM (HMHV-1B18, Totalseq-A, Biolegend), Ly6D (49-H4, Totalseq-A, Biolegend), Ly6C/Ly6G (Gr1) (RB6-8C5, Totalseq-A, Biolegend), MadCAM1 (MECA-367, Totalseq-A, Biolegend), Podoplanin (8.1.1, Totalseq-A, Biolegend), CD80 (16-10A1, Totalseq-A, Biolegend), CD86 (GL-1, Totalseq-A, Biolegend), MHCII (M5, Totalseq-A, Biolegend), Sca1 (D7, Totalseq-A, Biolegend), CD31 (390, Totalseq-A, Biolegend), EpCAM1 (G8.8, Totalseq-A, Biolegend), CD36 (HM36, Totalseq-A, Biolegend), CD133 (15-2C11, Totalseq-A, Biolegend), CD157 BP-3, Totalseq-A, Biolegend), CD300LG (ZAQ5, Totalseq-A, Biolegend), and the Ulex europaeus agglutinin I (UEA1) lectin labelled with biotin, followed by secondary staining with streptavidin-PE (Totalseq-A, Biolegend) coupled to an oligonucleotide. CD45⁻Ter119⁻EpCAM1⁺ and CD45⁻Ter119⁻EpCAM1⁻ cells were sorted in a 70–30% ratio into a 1.5 mL tube containing FACS buffer for the 1-week-old and 16-week-old samples, respectively. For both timepoints an estimate of 28,000 total cells were loaded on two wells of a 10x Genomics Chromium Single Cell Controller. After single-cell capture cDNA and library preparation were performed according to the manufacturer’s instructions using a Single-Cell 3’ v3 Reagent Kit (10x Genomics) with the changes as described in²⁶ to capture cDNA and produce libraries from antibody-derived oligos (ADT). Sequencing was performed on one lane of the Illumina NovaSeq 6000 system with a mix of 90% cDNA library and 10% ADT library resulting in 151nt-long paired-end reads.

The dataset was analysed by the Bioinformatics Core Facility, Department of Biomedicine, University of Basel. cDNA reads were aligned to ‘mm10’ genome using Ensembl 102 gene models with the STARsolo tool (v2.7.9a) with default parameter values except the following parameters: soloUMIlen=12, soloBarcodeReadLength=0, clipAdapterType=CellRanger4, outFilterType=BySJout, outFilterMultimapNmax=10, outSAMmultNmax=1, soloType=CB_UMI_Simple, outFilterScoreMin=30, soloCBmatchWLtype=1MM_multi_Nbase_pseudocounts, soloUMIfiltering=MultiGeneUMI_CR, soloUMIdedup=1MM_CR, soloCellFilter=None. ADT libraries were also processed using the STARsolo tool with default parameters except soloCBmatchWLtype=1MM_multi_Nbase_pseudocounts, soloUMIfiltering= MultiGeneUMI_CR, soloUMIdedup=1MM_CR, soloCellFilter=None, clipAdapterType=False, soloType=CB_UMI_Simple, soloBarcodeReadLength=0, soloUMIlen=12, clip3pNbases = 136.

Further analysis steps were performed using R (v4.1.2). Note that cell filtering was done based only on the analysis of the gene expression, not ADT abundance. Cells were considered as high-quality cells if they had at least 2000 UMI counts, which is the threshold derived from the distribution of UMI counts across cells, forming a data set of 9953 cells.

Multiple Bioconductor (v3.14) packages including DropletUtils (v1.14.2), scDblFinder (v1.8.0), scran (v1.22.1), scater (v1.22.0), scuttle (1.4.0) and batchelor (v1.10.0) were applied for the further analysis of the dataset mostly following the steps of the workflow presented at https://bioconductor.org/books/release/OSCA/. Normalised³⁹ log-count values for the gene expression were used to construct a shared nearest-neighbour graph⁴⁰, which nodes, i.e., cells, were clustered by ‘cluster_louvain’ method from the R igraph package⁴¹. Counts reflecting the ADT abundance in cells were also log-normalised and clustered in a similar manner.

The data set was subjected to the cell-type annotation using the Bioconductor package SingleR (v1.8.1) and samples from the Immunological Genome Project (ImmGen) provided by the Bioconducter package celldex (v1.4.0) as the reference. Clusters of cells mostly assigned to ‘Epithelial cells’ (5834 cells) were filtered (Supplementary Fig. 6d). Note that one of the clusters (cluster A, Supplementary Fig. 6a) was excluded at this step, because it was mostly composed of cells with elevated percentage of reads mapping to mitochondrial and ribosomal genes and lower number of counts.

The gene expression of filtered cells was re-analysed by removing the batch effect formed by the combination of the sample of origin and the number of counts per cell (cells with >12,000 counts and cells with 12,000 counts) and re-clustered (Fig. 5a, b). Cells were also subjected to the cell-type annotation using scRNAseq transcriptional profiles of single TEC as the reference data set¹³ (Fig. 5d, f). The scoreMarkers function of the scran package was applied to find marker genes of clusters 1–3. The standardised log-fold change across all pairwise comparisons ‘mean.logFC.cohen’>1 was used as the significance threshold defining the set of marker genes.

A t-SNE dimensionality reduction was used for visualizing single cells on two dimensions. T-SNE coordinates were calculated using the runTSNE function from the scater package and default parameters. For the visualization of cells based on the gene expression, coordinates of principal components and 2000 most variable genes with excluded mitochondrial and ribosomal genes were used as the input. For the visualization of cells based on the ADT abundance, coordinates of principal components and all ADTs were used as the input.

Bulk RNA sequencing

Triplicates of Sca1⁻CD63⁻CD66a⁺CD117⁺ tuft-like and Sca1⁻CD63⁻CD66a⁻CD117⁻ non-tuft-like mTEC cells were sorted from thymi of 6-week-old C57BL/6 mice into trizol (Invitrogen, 15596026). Subsequently samples were submitted to ultra-low input bulk RNAseq. Reads were trimmed using Trimmomatic (version 0.36) to remove adapter sequences and aligned to the mouse genome (mm10) using STAR (version 2.7.3a)^42,43. HTSeq (version 0.12.4) was used to assign reads to genes with the option “intersection-nonempty”⁴⁴. Differentially expressed genes were identified using edgeR (FDR < 0.05)⁴⁵. Spearman correlation coefficients were calculated between Sca1⁻CD63⁻CD66a⁻CD117⁻ non-tuft-like mTEC bulk RNAseq samples and scRNAseq data from reference ¹³ for all differentially expressed genes with log₂ fold change ≥1.

Statistical analysis

GraphPad Prism (version 9) was used to perform all statistical analyses, except for the bulk RNAseq and CITEseq datasets. The statistical tests used are described in the figure legends, and exact p-values are shown within each figure. Non-significant differences are not specified.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

The CITEseq and the bulk RNAseq datasets have been deposited in the Gene Expression Omnibus database under accession numbers GSE215418 and GSE226128, respectively. Source data are provided with this paper.

References

Anderson, G. & Takahama, Y. Thymic epithelial cells: working class heroes for T cell development and repertoire selection. Trends Immunol. 33, 256–263 (2012).
Article CAS PubMed Google Scholar
James, K. D., Jenkinson, W. E. & Anderson, G. Non-epithelial stromal cells in thymus development and function. Front. Immunol. 12, 634367 (2021).
Article CAS PubMed PubMed Central Google Scholar
Anderson, G., Jenkinson, E. J., Moore, N. C. & Owen, J. J. MHC class II-positive epithelium and mesenchyme cells are both required for T-cell development in the thymus. Nature 362, 70–73 (1993).
Article CAS PubMed ADS Google Scholar
Handel, A. E. et al. Developmental dynamics of the neural crest-mesenchymal axis in creating the thymic microenvironment. Sci. Adv. 8, eabm9844 (2022).
Article CAS PubMed PubMed Central Google Scholar
Klein, L., Kyewski, B., Allen, P. M. & Hogquist, K. A. Positive and negative selection of the T cell repertoire: what thymocytes see (and don’t see). Nat. Rev. Immunol. 14, 377–391 (2014).
Article CAS PubMed PubMed Central Google Scholar
Kadouri, N., Nevo, S., Goldfarb, Y. & Abramson, J. Thymic epithelial cell heterogeneity: TEC by TEC. Nat. Rev. Immunol. 20, 239–253 (2020).
Article CAS PubMed Google Scholar
Zlotoff, D. A. & Bhandoola, A. Hematopoietic progenitor migration to the adult thymus. Ann. N. Y Acad. Sci. 1217, 122–138 (2011).
Article CAS PubMed PubMed Central ADS Google Scholar
Abramson, J. & Anderson, G. Thymic epithelial cells. Annu Rev. Immunol. 35, 85–118 (2017).
Article CAS PubMed Google Scholar
Laufer, T. M., DeKoning, J., Markowitz, J. S., Lo, D. & Glimcher, L. H. Unopposed positive selection and autoreactivity in mice expressing class II MHC only on thymic cortex. Nature 383, 81–85 (1996).
Article CAS PubMed ADS Google Scholar
Derbinski, J., Schulte, A., Kyewski, B. & Klein, L. Promiscuous gene expression in medullary thymic epithelial cells mirrors the peripheral self. Nat. Immunol. 2, 1032–1039 (2001).
Article CAS PubMed Google Scholar
Herzig, Y. et al. Transcriptional programs that control expression of the autoimmune regulator gene Aire. Nat. Immunol. 18, 161–172 (2017).
Article CAS PubMed Google Scholar
Gray, D. H. et al. Developmental kinetics, turnover, and stimulatory capacity of thymic epithelial cells. Blood 108, 3777–3785 (2006).
Article CAS PubMed Google Scholar
Baran-Gale, J. et al. Ageing compromises mouse thymus function and remodels epithelial cell differentiation. Elife 9, e56221 (2020).
Article CAS PubMed PubMed Central Google Scholar
Dhalla, F. et al. Biologically indeterminate yet ordered promiscuous gene expression in single medullary thymic epithelial cells. EMBO J. 39, e101828 (2020).
Article CAS PubMed Google Scholar
Bornstein, C. et al. Single-cell mapping of the thymic stroma identifies IL-25-producing tuft epithelial cells. Nature 559, 622–626 (2018).
Article CAS PubMed ADS Google Scholar
Mayer, C. E. et al. Dynamic spatio-temporal contribution of single beta5t+ cortical epithelial precursors to the thymus medulla. Eur. J. Immunol. 46, 846–856 (2016).
Article CAS PubMed PubMed Central Google Scholar
Michel, C. et al. Revisiting the road map of medullary thymic epithelial cell differentiation. J. Immunol. 199, 3488–3503 (2017).
Article CAS PubMed Google Scholar
Onder, L. et al. Alternative NF-kappaB signaling regulates mTEC differentiation from podoplanin-expressing precursors in the cortico-medullary junction. Eur. J. Immunol. 45, 2218–2231 (2015).
Article CAS PubMed Google Scholar
Ulyanchenko, S. et al. Identification of a bipotent epithelial progenitor population in the adult thymus. Cell Rep. 14, 2819–2832 (2016).
Article CAS PubMed PubMed Central Google Scholar
Nusser, A. et al. Developmental dynamics of two bipotent thymic epithelial progenitor types. Nature 606, 165–171 (2022).
Article CAS PubMed PubMed Central ADS Google Scholar
Becht, E. et al. High-throughput single-cell quantification of hundreds of proteins using conventional flow cytometry and machine learning. Sci. Adv. 7, eabg0505 (2021).
Article CAS PubMed PubMed Central ADS Google Scholar
Miller, C. N. et al. Thymic tuft cells promote an IL-4-enriched medulla and shape thymocyte development. Nature 559, 627–631 (2018).
Article CAS PubMed PubMed Central ADS Google Scholar
Lucas, B. et al. Diversity in medullary thymic epithelial cells controls the activity and availability of iNKT cells. Nat. Commun. 11, 2198 (2020).
Article CAS PubMed PubMed Central ADS Google Scholar
Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e3529 (2021).
Article CAS PubMed PubMed Central Google Scholar
Rattay, K., Meyer, H. V., Herrmann, C., Brors, B. & Kyewski, B. Evolutionary conserved gene co-expression drives generation of self-antigen diversity in medullary thymic epithelial cells. J. Autoimmun. 67, 65–75 (2016).
Article CAS PubMed Google Scholar
Stoeckius, M. et al. Simultaneous epitope and transcriptome measurement in single cells. Nat. Methods 14, 865–868 (2017).
Article CAS PubMed PubMed Central Google Scholar
Heng, T. S., Painter, M. W. & Immunological Genome Project, C. The Immunological Genome Project: networks of gene expression in immune cells. Nat. Immunol. 9, 1091–1094 (2008).
Article CAS PubMed Google Scholar
Michelson, D. A., Hase, K., Kaisho, T., Benoist, C. & Mathis, D. Thymic epithelial cells co-opt lineage-defining transcription factors to eliminate autoreactive T cells. Cell 185, 2542–2558.e2518 (2022).
Article CAS PubMed PubMed Central Google Scholar
Hollander, G. A. et al. Developmental control point in induction of thymic cortex regulated by a subpopulation of prothymocytes. Nature 373, 350–353 (1995).
Article CAS PubMed ADS Google Scholar
Jacobs, H. et al. CD3 components at the surface of pro-T cells can mediate pre-T cell development in vivo. Eur. J. Immunol. 24, 934–939 (1994).
Article CAS PubMed Google Scholar
Shinkai, Y. & Alt, F. W. CD3 epsilon-mediated signals rescue the development of CD4+CD8+ thymocytes in RAG-2-/- mice in the absence of TCR beta chain expression. Int Immunol. 6, 995–1001 (1994).
Article CAS PubMed Google Scholar
Levelt, C. N., Wang, B., Ehrfeld, A., Terhorst, C. & Eichmann, K. Regulation of T cell receptor (TCR)-beta locus allelic exclusion and initiation of TCR-alpha locus rearrangement in immature thymocytes by signaling through the CD3 complex. Eur. J. Immunol. 25, 1257–1261 (1995).
Article CAS PubMed Google Scholar
Zheng, G. X. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).
Article CAS PubMed PubMed Central ADS Google Scholar
Ding, J. et al. Systematic comparison of single-cell and single-nucleus RNA-sequencing methods. Nat. Biotechnol. 38, 737–746 (2020).
Article CAS PubMed PubMed Central Google Scholar
Rota, I. A. et al. FOXN1 forms higher-order nuclear condensates displaced by mutations causing immunodeficiency. Sci. Adv. 7, eabj9247 (2021).
Article CAS PubMed PubMed Central ADS Google Scholar
Srinivas, S. et al. Cre reporter strains produced by targeted insertion of EYFP and ECFP into the ROSA26 locus. BMC Dev. Biol. 1, 4 (2001).
Article CAS PubMed PubMed Central Google Scholar
Rode, I. et al. Foxn1 protein expression in the developing, aging, and regenerating thymus. J. Immunol. 195, 5678–5687 (2015).
Article CAS PubMed Google Scholar
Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat. Immunol. 20, 163–172 (2019).
Article CAS PubMed PubMed Central Google Scholar
Lun, A. T., Bach, K. & Marioni, J. C. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Genome Biol. 17, 75 (2016).
Article PubMed Google Scholar
Xu, C. & Su, Z. Identification of cell types from single-cell transcriptomes using a novel clustering method. Bioinformatics 31, 1974–1980 (2015).
Article CAS PubMed PubMed Central Google Scholar
Blondel, V. D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. Journal of Statistical Mechanics: Theory and Experiment 2008, P10008 (2008).
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).
Article CAS PubMed PubMed Central Google Scholar
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Article CAS PubMed Google Scholar
Anders, S., Pyl, P. T. & Huber, W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169 (2015).
Article CAS PubMed Google Scholar
Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).
Article CAS PubMed Google Scholar

Download references

Acknowledgements

We thank Emilie Cosway and Sonia Parnell for their technical assistance and Lilly von Muenchow for critical reading of the manuscript. Further we would like to acknowledge Grozdan Cvijetic for his help to set up the Infinity Flow computational analysis pipeline. We also thank Jakub Abramson for sharing the Csnb^Cre mice. This research was supported in part by the generosity of an anonymous donor to Stanford University establishing the DiGeorge Syndrome Research Fund (G.A.H.), Swiss National Science Foundation grant IZLJZ3_171050; 310030_184672 (G.A.H.), Medical Research Council grant MR/S036407/1 (G.A.H.), Wellcome Trust grant 105045/Z/14/Z and 211944/Z/18/Z (G.A.H.); Swiss National Science Foundation Early Postdoc.Mobility Fellowship (P2BSP3_188183) and Postdoc.Mobility Fellowship (P500PB_206823) to F.K.; NIHR Clinical Lectureship to F.D.; work in the G.A. lab is supported by an MRC Programme Grant (MR/T029765/1).

Author information

Authors and Affiliations

Department of Paediatrics and Institute of Developmental and Regenerative Medicine, University of Oxford, Oxford, UK
Fabian Klein, Clara Veiga-Villauriz, Stefano Maio, Fatima Dhalla, Adam E. Handel, Mary E. Deadman & Georg A. Holländer
Department of Biomedicine, University of Basel, Basel, Switzerland
Anastasiya Börsch
Mathematical Institute, University of Oxford, Oxford, UK
Sam Palmer
Nuffield Department of Clinical Neurosciences, University of Oxford, Oxford, UK
Adam E. Handel
Paediatric Immunology, Department of Biomedicine, University of Basel and University Children’s Hospital Basel, Basel, Switzerland
Saulius Zuklys, Irene Calvo-Asensio, Lucas Musette & Georg A. Holländer
Institute for Immunology and Immunotherapy, Medical School, University of Birmingham, Birmingham, UK
Andrea J. White, Beth Lucas & Graham Anderson
Department of Biosystems Science and Engineering, ETH Zurich, Basel, Switzerland
Georg A. Holländer

Authors

Fabian Klein
View author publications
You can also search for this author in PubMed Google Scholar
Clara Veiga-Villauriz
View author publications
You can also search for this author in PubMed Google Scholar
Anastasiya Börsch
View author publications
You can also search for this author in PubMed Google Scholar
Stefano Maio
View author publications
You can also search for this author in PubMed Google Scholar
Sam Palmer
View author publications
You can also search for this author in PubMed Google Scholar
Fatima Dhalla
View author publications
You can also search for this author in PubMed Google Scholar
Adam E. Handel
View author publications
You can also search for this author in PubMed Google Scholar
Saulius Zuklys
View author publications
You can also search for this author in PubMed Google Scholar
Irene Calvo-Asensio
View author publications
You can also search for this author in PubMed Google Scholar
Lucas Musette
View author publications
You can also search for this author in PubMed Google Scholar
Mary E. Deadman
View author publications
You can also search for this author in PubMed Google Scholar
Andrea J. White
View author publications
You can also search for this author in PubMed Google Scholar
Beth Lucas
View author publications
You can also search for this author in PubMed Google Scholar
Graham Anderson
View author publications
You can also search for this author in PubMed Google Scholar
Georg A. Holländer
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

F.K., C.V.-V., S.M., F.D., S.Z., L.M., I.C.-A., M.E.D., A.J.W., and B.L. performed experiments; F.K., C.V.-V., S.P., A.B., F.D., A.E.H., S.Z., and L.M. analysed data; F.D. provided CsnbCre mice; G.A. provided the thymi from Pou2f3-/- mice; F.K. and G.A.H. conceived the project and wrote the manuscript.

Corresponding author

Correspondence to Georg A. Holländer.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary Information

Peer Review File

Reporting Summary

Source data

Source Data

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Klein, F., Veiga-Villauriz, C., Börsch, A. et al. Combined multidimensional single-cell protein and RNA profiling dissects the cellular and functional heterogeneity of thymic epithelial cells. Nat Commun 14, 4071 (2023). https://doi.org/10.1038/s41467-023-39722-9

Download citation

Received: 23 October 2022
Accepted: 21 June 2023
Published: 10 July 2023
DOI: https://doi.org/10.1038/s41467-023-39722-9

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.