C57BL/6 mice were maintained in the Max Planck Institute of Immunobiology and Epigenetics. Foxn1-eGFP57, Foxn1-cre58, Rosa26-LSL-EYFP59, Foxn1-s-Fgfr2IIIb60, pLck-cre61, Rosa26-LSL-Cas9-EYFP62 and Foxn1-mCardinal9 transgenic mice have been described previously. The Foxn1-Fgf7 transgene was created by inserting a cDNA fragment corresponding to nucleotides 347–934 in GenBank accession number NM_008008 as a NotI fragment into pAHB14 (ref. 63); in some aged female Foxn1-Fgf7 transgenic mice (FVB/N-tg(Foxn1-Fgf7)1Tbo/Mpie), the two thymic lobes were asymmetric in size and shape; these mice were not included in our analysis. The Foxn1-Fgfr2IIIb transgene was created by inserting a cDNA fragment corresponding to nucleotides 1214–3366 in GenBank accession number NM_201601.2 as a NotI fragment into pAHB14 (ref. 63) and used to generate transgenic mice (FVB/N-tg(Foxn1-Fgfr2)1Tbo/Mpie). The hU6-sgRNAHprt transgene was cloned as a NotI fragment into the Bluescript vector and consists of the human U6 promotor (nucleotides 1–264 in GenBank accession number JN255693) followed by the mouse Hprt target sequence (5′-GATGGGAGGCCATCACATTGG-3′; nucleotides 255–274 in GenBank accession number J00423), the sgRNA backbone (nucleotides 218–139 (reverse complement) in Addgene plasmid 42250) and a short 3′ sequence (TTTTTTGGAA); for injection into fertilized eggs, the construct was linearized with SacI. Transgenic mice were generated on an FVB/N background (FVB/N-tg(hU6-sgRNA-Hprt)1Tbo/Mpie) and subsequently backcrossed to a C57BL/6J background. For timed matings, the day of plug detection was designated as E0.5. Genotyping information is summarized in Supplementary Table 15. Mice were kept in the animal facility of the Max Planck Institute of Immunobiology and Epigenetics under specific pathogen-free conditions. All animal experiments were performed in accordance with the relevant guidelines and regulations, approved by the review committee of the Max Planck Institute of Immunobiology and Epigenetics and Regierungspräsidium Freiburg, Germany (licences 35-9185.81/G-12/85; 35-9185.81/G-16/67).

KGF treatment

At the age of 4 weeks, male mice received nine intraperitoneal injections of KGF (Kepivance, Biovitrum, lot D120961G; 5 mg kg–1 body mass) at days 1, 2, 3, 8, 9, 10, 15, 16 and 17; the mice were killed on day 21.


Embryos for RNA in situ hybridization (ISH) were fixed in 4% paraformaldehyde (PFA) and subsequently embedded in paraffin using standard techniques.


The Applied Biosystems 7500 Fast system was used to detect the signal generated with gene-specific primers combined with 5′-FAM (6-carboxyfluorescein)-labelled hydrolysis probes from Universal Probe Library (Roche). Primer sequences were as follows: Fgf7, 5′-TGGCTGACACCATGACTAGC-3′ and 5′-GGCTACAGGCTGTCGTTTTT-3′ (probe 42); Fgf10, 5′-CGGGACCAAGAATGAAGACT-3′ and 5′-GCAACAACTCCGATTTCCAC-3′ (probe 80); CD31 (Pecam), 5′-CGGTGTTCAGCGAGATCC-3′ and 5′-ACTCGACAGGATGGAAATCAC-3′ (probe 45); Ly51 (Enpep), 5′-TGGACTCCAAAGCTGATCCT-3′ and 5′-TCAGCCCATCTGACTGGAAT-3′ (probe 83). Expression levels were normalized to those of Hprt, using primers 5′-TCCTCCTCAGACCGCTTTT-3′ and 5′-CCTGGTTCATCATCGCTAATC-3′ (probe 95).


RNA ISH on paraffin sections was performed using DIG-labelled probes as described64. Sequence coordinates in GenBank accession numbers were as follows: Foxn1, nucleotides 2181–3584 in XM_006532266.3; Fgf7, nucleotides 153–877 in NM_008008.3; Fgf10, nucleotides 859–1570 in NM_008002.3; Fgfr1, nucleotides 761–1614 in NM_001079909.2; Fgfr2, nucleotides 328–800 in EF143340; Fgfr2_exon3b, nucleotides 1819–1964 in NM_201601.2; Fgfr2_exon3c, nucleotides 2169–2306 in NM_010207.2; Hspb1, nucleotides 224–632 in NM_013560.2; Trpm5, nucleotides 499–962 in NM_020277.2.


Thymi were fixed in 4% PFA, washed in PBS, incubated in 20% sucrose overnight and embedded in OCT. Sections of 8–10 µm were dried overnight at room temperature and before staining were moistened in PBS followed by a 30-min blocking step (PBS with 0.5% BSA, 0.2% Tween and anti-mouse IgG (1:50)). Antibody staining was performed at room temperature in staining buffer (PBS with 0.5% BSA, 0.2% Tween and 3% serum). Sections were stained for 2 h with primary antibodies (Supplementary Table 6) and then for 45 min with secondary antibodies and streptavidin. Sections were washed with PBS between incubations. After staining, sections were mounted in Fluoromount G.

Image analysis

Images were acquired on Zeiss microscopes (Axioplan 2 or Imager Z1 with ApoTome attachment) equipped with AxioCam MRc 5 cameras.

Flow cytometry

To generate single-cell suspensions for analytical and preparative flow cytometry of TECs, the procedures described in refs. 18,65 were followed. Relevant staining reagents are listed in Supplementary Table 16. The enzymatic cocktail required to liberate TECs destroys the extracellular domains of the CD4 and CD8 surface markers (but not that of the CD45 molecule); hence, when analysis of thymocyte subsets was desired, thymocyte suspensions were prepared in parallel by mechanical liberation, achieved by gently pressing thymic lobes through 40-µm sieves. To isolate thymic mesenchymal and endothelial cells, the cell suspension of total thymocytes was depleted of CD45+ cells; the EpCAMCD45 cell population was stained for Ly51 and CD31 to purify EpCAMCD31Ly51+ mesenchymal and EpCAMCD31+Ly51 endothelial cells. Cell sorting and analytical flow cytometry were carried out using MoFlow and Fortessa instruments, respectively (both from Dako Cytomation-Beckman Coulter); flow cytometry experiments were carried out using FACSDiva and FlowJo software. The fraction of Foxn1-expressing cells was determined by eGFP fluorescence emanating from the Foxn1-eGFP transgene57, which faithfully recapitulates acute levels of Foxn1 expression18,34. The thymopoietic index was calculated by dividing the total number of thymocytes by the number of TECs.

Single-cell RNA amplification and library preparation

scRNA-seq was performed using the CEL-Seq2 method45 with several modifications44. A fivefold volume reduction was achieved using a nanolitre-scale pipetting robot (Mosquito HTS, TTP Labtech)66. EpCAM+CD45 TECs were sorted into 384-well plates containing 240 nl of primer mix67 and 1.2 μl of mineral oil (Sigma-Aldrich) replacing Vapor-lock. Sorted plates were centrifuged at 2,200g for 3 min at 4 °C, snap-frozen in liquid nitrogen and stored at −80 °C before further processing. To convert RNA into cDNA, 160 nl of reverse-transcription reaction mix including 0.4 µM template-switch oligonucleotide (5′-AAGCAGTGGTATCAACGCAGAGTGAATrGrGrG-3′; adapted from ref. 68) and 2.2 μl of second-strand reaction mix were added.

For the purpose of simultaneous transcriptome and barcode analysis, the volume of each well was equally split in half; that is, 1.1 µl per well was transferred to a new 384-well plate. The original plate (including the mineral oil) was used for barcode analysis, while the copy of the plate was used for analysis of individual cell transcriptomes.

For barcode analysis, all transcripts in each well were amplified by template-switch PCR (5′-AAGCAGTGGTATCAACGCAGAGT-3′ (template-switch oligonucleotide primer) and 5′-CAGAGTTCTACAGTCCGA-3′ (short P5 primer)), followed by amplification of the scar region of Hprt transcripts (5′-GCCGGTAATACGACTCACTATAGGGAGTTCTACAGTCCGACGATC-6bp UMI-6bp cell barcode-CCAGTTAAAGTTGAGAGATCATCTCC-3′ (Hprt barcoding primer), 5′-GCCTTGGCACCCGAGAATTCCATAGCGATGATGAACCAGGTTATGACC-3′ (Hprt short P7 primer)) using the PrimeSTAR GXL system (Takara Bio). Libraries were completed by addition of full-length adaptors by PCR (RP1, RPI1-48 TruSeq Small RNA Illumina adaptor sequences). Libraries from 96 cells were pooled before clean-up. In total, 32 libraries (E16.5 time point), 32 libraries (P0 time point), 28, 40 and 80 libraries (three P28 time points), 32 libraries (1-year time point), 44 libraries (Foxn1-Fgf7 transgenic mice, P28 time point) and 60 libraries (Foxn1-Fgf7 transgenic mice, 1-year time point) were sequenced on the MiSeq sequencing system (2 × 300 bp) at a depth of ≥4,000 reads per cell.

Transcriptomes were generated as described previously66; for the P28 wild-type dataset from a non-barcoded mouse, the whole sample volume was used for transcriptome generation. cDNAs from 96 cells were pooled before clean-up and in vitro transcription, generating four libraries from one 384-well plate. For all purification steps, 0.8 μl of AMPure/RNAClean XP beads (Beckman Coulter) was used per 1 μl of sample, including in library clean-up.

In total, 32 libraries (E16.5 time point), 32 libraries (P0 time point), 28, 40 and 80 libraries (three P28 time points), 28 libraries (non-barcoded dataset of P28 time point), 32 libraries (1-year time point), 44 libraries (Foxn1-Fgf7 transgenic mice at P28 time point) and 60 libraries (Foxn1-Fgf7 transgenic mice at 1-year time point) (each library was generated by pooling 96 cells) were sequenced on the Illumina HiSeq 2500 or NovaSeq 6000 sequencing system (paired-end multiplexing run, high-output mode) at a depth of ~170,000 reads per cell.

Quantification of transcript abundance

Paired-end reads were aligned to the transcriptome using bwa (version 0.6.2-r126) with default parameters69. The transcriptome contained all gene models based on the mouse ENCODE VM9 release downloaded from the UCSC genome browser, comprising 57,207 isoforms, with 57,114 isoforms mapping to fully annotated chromosomes (1 to 19, X, Y, mitochondria). All isoforms of the same gene were merged to a single gene locus. Furthermore, gene loci overlapping by >75% were merged to larger gene groups. This procedure resulted in 34,111 gene groups. The right mate of each read pair was mapped to the ensemble of all gene loci and to the set of 92 ERCC spike-ins in the sense direction70. Reads mapping to multiple loci were discarded. The left read contained the barcode information: the first six bases corresponded to the cell-specific barcode followed by six bases representing the unique molecular identifier (UMI). The remainder of the left read contained a poly(T) stretch. The left read was not used for quantification. For each cell barcode, the number of UMIs per transcript was counted and aggregated across all transcripts derived from the same gene locus. On the basis of binomial statistics, the number of observed UMIs was converted into transcript counts71.

scRNA-seq data analysis

Clustering analysis and visualization of all datasets were performed with the VarID algorithm46. Cells with a total number of transcripts of <1,000 (1-year wild-type dataset), <1,500 (wild-type and Foxn1-Fgf7 transgenic P28 datasets) and <3,000 (wild-type E16.5, wild-type P0 and Foxn1-Fgf7 transgenic 1-year datasets) were discarded, and the count data of the remaining cells were normalized by downscaling. Notably, before normalization, cells yielding transcriptomes containing >2% Kcnq1ot1 transcripts, a previously identified marker of low-quality cells, were removed from the analysis72. Moreover, transcripts correlating to Kcnq1ot1 with a Pearson’s correlation coefficient of >0.65 were also removed. Furthermore, mitochondrial genes as well as Hprt were excluded. The following genes and correlating gene groups were removed from the analysis using the CGenes parameter (all datasets): Jun, Fos and predicted genes with Gm identifiers. For the wild-type P28 datasets, Malat1, Xist and Neat1 were excluded using the FGenes parameter. A pruned k nearest neighbour (kNN) matrix was inferred using the pruneKnn function of VarID with default parameters except alpha (set to 1) and no_cores (set to 10). The UMAP representation was used for cell cluster visualization73. Transition probabilities between clusters were computed using the transitionProbs function of VarID with the P value set to 0.001. Differentially expressed genes between two subgroups of cells were identified similarly to in a previously published method74. First, negative binomial distributions reflecting the gene expression variability within each subgroup were inferred on the basis of the background model for the expected transcript count variabilitycomputed by the RaceID3 algorithm44. Using these distributions, a P value for the observed difference in transcript counts between the two subgroups was calculated and corrected for multiple testing with the Benjamini–Hochberg method.

High-resolution lineage tracing

The lineage tracing method developed here was based on CRISPR–Cas9-mediated scarring in exon 3 of the Hprt gene (Fig. 2a). Repair of Cas9-induced double-strand breaks results in a number of different sequence outcomes, which we refer to here as barcodes because they indelibly mark all cell progeny. Most barcodes carry deletions (Fig. 2b), preventing secondary modifications of their sequence; however, in rare instances, we recorded cells with two barcodes of similar sequence, presumably originating from ongoing modifications of the target sequence or from relaxed X-chromosome inactivation. In bulk and single-cell analyses of male mice, barcodes can be unambiguously read out on both the DNA and transcript level from the single X chromosome. In female mice, which carry two X chromosomes, DNA analysis at the level of cell populations yields ambiguous results; however, for single-cell analysis, the phenomenon of dosage compensation enables unambiguous barcode analysis at the RNA level. For the bulk analyses described here, we therefore used only male mice, whereas single-cell transcriptome analyses were conducted in both male and female mice. In the hU6-sgRNAHprt; Foxn1-cre; Rosa26-LSL-Cas9-EYFP triple-transgenic mice used here, TECs are marked in early embryogenesis, as soon as Foxn1 expression begins; our previous analysis indicated that, at the onset of Foxn1 expression, the thymic rudiment harbours ~4,000 epithelial cells9, thus placing an upper limit on the number of barcodes that can be observed in TECs at later developmental stages. However, the observed number of barcodes was three- to fourfold lower, presumably because the outcome of the repair process is not random; hence, some barcodes, although independently generated, are identical in sequence. The different frequencies of barcode generation must be taken into account when reconstructing lineage relationships; rare barcode sequences are more informative than frequently generated barcodes. On the basis of previous experiments using our Foxn1-cre transgenic line and the accessibility of the Rosa26 locus in TECs34, we assume that the overwhelming majority of TECs (>95%) are barcoded early in embryogenesis. However, we cannot exclude the possibility that some cells in the thymic rudiment are recruited into the Foxn1-positive lineage at later stages of development. If so, this would however only lessen the ability to identify clonal relationships across all time points.

Single-cell barcoding analysis

Paired-end fastq files were used for the identification of scar sequences in single cells. The first six bases of the left read contained the UMI information, followed by six bases representing the cell barcode. The remainder of the left read contained the Hprt scar sequence. The right mate of the paired-end reads also contained the overlapping Hprt scar sequence; that is, both reads contained the full scar sequences. Primers (forward, 5′-GCTCGAGATGTCATGAAGG-3′; reverse, 5′-GGGGGGCTATAAGTTCTT-3′) were used to extract the targeted region of the Hprt gene containing the edited sequence. Because both the left and right mates of the paired-end reads contain the full scar sequences, only sequences that appeared in both paired-end reads were used for further analysis. Moreover, only cells yielding ≥200 reads were included in the analysis; cells were excluded from further analysis if more than one scar sequence was detected in male cells (the threshold for the second sequence was set at ≥10% of the major sequence) and if more than two scar sequences were detected in female cells (the threshold for the third sequence was set at ≥10% of the second sequence).

Barcodes shared by cTECs and mTECs

DNA was isolated from sorted EpCAM+CD45Ly51+UEA-1 (cTEC) and EpCAM+CD45Ly51UEA-1+ (mTEC) populations from each mouse, and the region of exon 3 of the Hprt gene was amplified using the following primers: 5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCTTTCATAGAGACAAGGAATGTGTCC-3′ (forward, P5-DD302)and 5′-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTAGTTGATTATGTAGCAtAGTTTGACAAG-3′ (reverse, P7-DD305). Libraries were sequenced at a depth of ~250,000 reads per sample on the MiSeq sequencing system (2 × 300 bp).

Next, a table containing the counts of all barcodes across cTECs and mTECs for each mouse (n = 33 mice) was constructed, and the frequency distribution of barcodes was determined (Fig. 2d). To quantify enrichment of shared barcodes between the cTEC and mTEC samples for a given mouse, we first extracted the set of all barcodes Bi for sample i that were observed no more than twice in all animal samples; then, we determined the number of barcodes within set Bi co-occurring in another sample j and divided this by the number of barcodes in Bi to compute the co-occurrence probability of rare barcodes, termed Pij:

$${P}_{{ij}}=frac{sum _{kin {B}_{i}}{chi }_{{B}_{j}}left(kright)}{sum _{kin {B}_{i}}{chi }_{{B}_{i}}left(kright)}$$

Here χT(x) denotes the indicator function; that is, χT(x) = 1 if xTχT(x) = 0 otherwise.

If i and j denote mTEC and cTEC samples from the same mouse, we would like to test whether Pij is relatively increased compared with cases where i and j denote samples from different mice.

An increased co-occurrence probability in corresponding samples compared with all other mice is indicative of a common barcode repertoire and is interpreted to mean that cTEC and mTEC populations in a mouse arise from common progenitors marked by particular barcodes; conversely, a similar ratio across all animals suggests the random occurrence of rare barcodes and argues against a common origin. The observation that some of the rare barcodes are not shared by cTEC and mTEC populations from the same mouse can be explained by the low frequency of these rare barcodes, resulting in sampling dropouts, but may also have a biological explanation, for example, if either mTEC or cTEC progeny derived from the same progenitor have died out. In this context, it is worth noting that, without prefiltering based on barcode frequency, the fraction of barcodes shared by mTEC and cTEC populations was >50%.

To quantify a single enrichment value Em for a given animal m, the ratios calculated for a corresponding pair of cTEC and mTEC samples, mcTEC and mmTEC, from the same mouse were divided (with cTECs or mTECs as a reference) by the average of the ratios for pairings involving either the respective cTEC or mTEC sample and a sample from any other animal, excluding the samples from animal m. The enrichment value Em is then calculated as the maximum of these ratios with mTECs or cTECs from the same animal as a reference:

$$begin{array}{ccc}{E}_{m} & = & {rm{max }}(frac{{P}_{{m}_{{rm{cTEC}}}{m}_{{rm{mTEC}}}}}{frac{1}{N-2}sum _{jne {m}_{{rm{cTEC}}}{vee jne m}_{{rm{mTEC}}}}{P}_{{m}_{{rm{cTEC}}}j}},\ & & frac{{P}_{{m}_{{rm{mTEC}}}{m}_{{rm{cTEC}}}}}{frac{1}{N-2}sum _{jne {m}_{{rm{cTEC}}}{vee jne m}_{{rm{mTEC}}}}{P}_{{m}_{{rm{mTEC}}}j}})end{array}$$

Here N is the total number of samples and indicates logical conjunction. In the summation of the denominators, we exclude pairings involving the mTEC or cTEC samples for animal m.

The enrichment values of shared rare barcodes for each mouse were then plotted for each individual mouse grouped by age bin in Fig. 2f. Note that the analysis of bulk populations is most robust with respect to sampling when the numbers of cTECs and mTECs are approximately equal, as is the case for the P14 time point9. In situations where only a small number of barcodes are recovered, a diminished degree of co-occurrence is probably the result of sampling dropouts and the associated reduced statistical power.

Lineage analysis

To quantify the enrichment of shared barcodes between different populations within the TEC compartment at single-cell resolution, we first identified the cell clusters representing early progenitors, postnatal progenitors, and mature cTECs and mTECs in the scRNA-seq analysis of the individual mice. For a given mouse, we determined the barcode repertoire and counted the number of cells for each barcode; in the rare instances where a pair of barcodes was observed, we used the combination of barcodes for quantification. Cells with more than two barcodes were discarded, as this situation is the result of cell doublets and/or sequencing errors. For each barcode (or pair of quantified barcodes), we next determined its frequency for each TEC population in a given sample. These frequencies were compared with background barcode frequencies derived from the barcode distribution quantified from the bulk DNA sequencing data (n = 33; Fig. 2d), averaging across all samples. We considered sampling dropouts as a reference background model for technical variability and, therefore, assessed significant over-representation of barcodes on the basis of the estimated probability mass of a binomial distribution with a probability parameter informed by the barcode frequency derived from the bulk sequencing experiments. It is well known that UMI-based abundance derived from scRNA-seq data can be modelled by a negative binomial distribution without explicitly modelling zero inflation and that variability in genes (or barcodes) with low expression is described well by the binomial noise component71.

If nb,i is the number of times barcode b was observed in bulk DNA sequencing sample i, the background frequency of barcode fb was calculated as

$${f}_{b}=frac{1}{N}sum _{iin S}frac{{n}_{b,i}}{sum _{kin B}{n}_{k,i}}$$

Here B denotes the set of all barcodes, S denotes the set of all samples and N is the size of S.

Assuming binomial sampling statistics, we then calculated the P value Pb,i for the observed number of cells cb,i with barcode b among all barcode-carrying cells in sample i obtained by scRNA-seq as the right tail probability of the binomial distribution with background probability fb and the total number of barcode-carrying cells as parameters:

$${P}_{b,i}=mathop{sum }limits_{j={c}_{b,i}}^{{C}_{i}}B,(j|,{f}_{b},,,{C}_{i},{rm{w}}{rm{i}}{rm{t}}{rm{h}},{C}_{i})=sum _{kin B}{c}_{k,i}$$

These P values were further corrected for multiple testing by the Benjamini–Hochberg method, considering all P values for all detected barcodes in a given dataset. For cells with two barcodes, we calculated P values accordingly after multiplying the background frequencies of the co-occurring barcodes. Barcodes were considered informative if their P value indicated a significant deviation from the expected frequency (Pb,i < 0.001).

Gene set analysis

Population-specific gene sets were derived by performing differential gene expression analysis of clusters representing early progenitors (c5), postnatal progenitors (c1 and c6), cTECs (c4) and mTECs (c12, c2, c7, c9, c10 and c18) versus all other clusters from data for 4-week-old (P28) mice using the diffexpnb function of the RaceID3 package44. Genes with adjusted P < 0.01 and log2(fold change) > 1 were included in gene sets, from which genes with Gm and Rik identifiers were excluded. Overlapping genes between early progenitor and cTEC gene sets were excluded from the early progenitor gene set. Genes included in the gene sets of each of the four populations are listed in Supplementary Tables 14.

Although the transcriptional profiles of the progenitor clusters differed from those of the mature TEC subsets by the expression of heat shock protein genes, these genes were not included in the final lists, as they did not distinguish between the early and postnatal progenitors. The final gene sets were analysed for enriched biological processes using the Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 Analysis Wizard annotation tool75,76. In the representation of gene set ratios, clusters expressing T cell progenitor-related genes (representing thymic nurse cells) and parathyroid-associated genes (representing ectopic parathyroid tissue) were excluded.

scRNA-seq data comparisons

The present data were compared with publicly available scRNA-seq data for TECs isolated from mice of different ages1. To do this, the raw count matrices and metadata describing the nine subtypes of TECs were obtained through the Bioconductor data package MouseThymusAgeing (https://bioconductor.org; https://doi.org/10.18129/B9.bioc.MouseThymusAgeing). Data normalization, dimensionality reduction and visualization with UMAP were then performed using the default parameters of the scRNA-seq data analysis CRAN package Seurat version 3 (ref. 77).

Statistical analysis and reproducibility

Two-tailed t tests were used to determine the significance levels of differences between the means of two independent samples, considering equal or unequal variance as determined by the F test. For multiple tests, the conservative Bonferroni correction was applied. For all analyses, several biological replicas were studied; numbers of replicas are indicated in the figures and/or figure legends. No statistical methods were used to predetermine sample sizes; blinding and randomization were not used. 

Reporting summary

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

Source link