Experimental methodsTissue collection

Inbred male and female C57BL/6J mice were obtained from The Jackson Laboratories or through the National Institute of Aging mouse colony. Foxn1tdTomato mice were generated by crossing Foxn1-cre (Jax 018448) with B6.Cg-Gt(ROSA)26Sortm14(CAG-tdTomato)Hze/J mice (Jax 007914). Foxn1nTnG mice were generated by crossing Foxn1-cre (Jax 018448) with ROSAnT-nG mice (Jax 023537). Foxn1z/z mice were generated as previously described64. As a model of thymic injury, mice were given a sublethal dose of TBI (550 cGy) with a Cs-137 γ-radiation source. Mice were maintained at The Memorial Sloan Kettering Cancer Center, Fred Hutchinson Cancer Center, Walter and Eliza Hall Institute or University of Georgia animal houses. All experiments were performed according to Institutional Animal Care and Use Committee guidelines.

Human thymus tissue was obtained from the archives of the Duke University Department of Pathology as formalin-fixed paraffin-embedded sections. All tissues were used anonymously, only recording patient age, sex and surgical diagnosis. We show images from one 50-year-old female patient. All human tissues were collected according to a protocol approved by the Duke University Institutional Review Board.

Isolation of cells and flow cytometry

The thymus was enzymatically digested following an adapted protocol73,74. In brief, thymi were mechanically dissociated into 1–2-mm pieces. Tissue pieces were incubated with a digestion buffer (either RPMI with 10% FCS, 62.5 μm ml−1 liberase TM, 0.4 mg ml−1 DNase I; or RPMI with 25 mM HEPES, 20 μg μl−1 DNase1 and 1 mg ml−1 collagenase/dispase). Between incubation steps, supernatant containing dissociated cells were transferred to tubes equipped with a 100-μm filter. Cells were pelleted by centrifugation at 400g for 5 min. All steps were performed at 4 °C unless indicated. For sequencing experiments, cell pellets were incubated with anti-mouse CD45 microbeads and CD45+ cells were depleted from cell suspension using magnetic-associated cell sorting (MACS) on LS columns according to the manufacturer’s protocol. Following red blood cell lysis using ACK buffer, the CD45-depleted cell fraction was incubated with an antibody cocktail for 15 min at 4 °C and cells of interest were purified by fluorescent-associated cell sorting (FACS) on a BD Biosciences Aria II using a 100-μm nozzle. Cells were sorted into tubes containing RPMI supplemented with 2% BSA. FACS-purified cells were spun down at 400g for 5 min and resuspended in PBS supplemented with 0.04% BSA for generation of single-cell suspensions.

For flow cytometry and cell sorting, surface antibodies against CD45 (30-F11), CD31 (390 or MEC13.3), TER-119 (TER-119), MHCII IA/IE (M5/114.15.2), EpCAM (G8.8), Ly51 (6C3), PDGFRα (APA5), CD104 (346-11A), L1CAM (555), Ly6D (49-H4), Gp38 (8.1.1), CD26 (H194-112), CD62P (RB40.34), podoplanin (8.1.1), CD62P (RB40.34), CD9 (KMC8) and CD309 (Avas12a) were purchased from BD Biosciences, BioLegend or eBioscience. Ulex europaeus agglutinin 1 (UEA1) was purchased from Vector Laboratories. Antibody against phospho-AKT was purchased from Cell Signaling Technologies; claudin-3 and anti-rabbit secondary antibodies were purchased from Invitrogen (Thermo Fisher); DCLK1 (aa690-720) was purchased from LSBio; GP2 (2F11-C3) was purchased from MBL Life Science; and anti-GFP (Aves GFP-1020) was purchased from AvesLabs. Anti-FOXN1 antibody was a gift from H.-R. Rodewald75. Flow cytometry was performed on a Fortessa X50 or Symphony A6 (BD Biosciences) and cells were sorted on an Aria II (BD Biosciences) using FACSDiva (BD Biosciences). Analysis was performed by FlowJo (Treestar Software). Detail of specific vendors, fluorochromes, catalog numbers, lot numbers dilutions and gating can be found in Supplementary Tables 11, 12.

Cells were isolated as described and depleted of CD45+ cells by MACS depletion (Miltenyi Biotech). CD45− cells were incubated for 5 min with recombinant mouse KGF (100 ng ml−1) when phospho-AKT was assessed. For phospho-AKT staining, cells were fixed and permeabilized in 1.6% paraformaldehyde at 37 °C followed by 90% methanol at 4 °C. After thorough washing to remove all methanol, cells were stained for both intracellular and extracellular antigens simultaneously.

Thymic tissue clearing and immunofluorescence

After euthanasia, mice were transcardially perfused with PBS followed by 4% PFA. Thymi were dissected and post-fixed in 4% PFA for 4 h at 4 °C. For confocal imaging, fixed tissue was sectioned at 200 μm using a Leica VT1000 S vibratome. Tissue clearing was performed as previously described76 with some modifications. In brief, tissue was immersed in monomer buffer (4% acrylamide and 0.25% (w/v) azo-initiator (Wako Pure Chemical Industries) in PBS) and incubated at 4 °C overnight. The solution was transferred to a vacuum tube and bubbled with nitrogen gas for 15 min. The gel was set for 3 h at 37 °C with gentle rotation, after which the tissue was transferred to clearing buffer (8% SDS and 50 mM sodium sulfite in PBS) and cleared at 37 °C until turning semi-transparent. To remove SDS, samples were transferred to the following buffers to wash for 1 h each with rotation: (1) 1% SDS, 0.5% Triton-X in PBS; (2) wash buffer (1% BSA and 0.5% Triton-X in PBS) for two washes; and (3) blocking buffer (4% normal serum, 1% BSA and 0.3% Triton-X in PBS) for two washes. Antibodies were diluted in blocking buffer at the dilutions indicated below. The antibodies used were rabbit anti-pan-cytokeratin (Dako, cat. no. Z0622), anti-K5 (BioLegend, cat. no. poly19055), rat anti-mouse K8/18 (Troma-1; Developmental Studies Hybridoma Bank), rabbit anti-K14 (Abcam, cat. no. EPR17350), rat anti-mouse AIRE (WEHI, clone 5H12), rabbit anti-human/mouse DCLK1 (LSBio, cat. no. LS-C100746) and biotinylated UEA1 lectin (Vector Labs, cat. no. B-1065). The secondary antibodies used were Alexa Fluor 647 donkey anti-rabbit IgG (H+L) (Invitrogen, cat. no. A31573), Alexa Fluor 647 goat anti-rat IgG (H+L) (Invitrogen, cat. no. A-21247) and Alexa Fluor 647 streptavidin conjugate (Invitrogen, cat. no. S21374). After staining, samples were washed in PBS with 0.3% Triton-X. For imaging, samples were incubated in EasyIndex optical clearing solution (refractive index, RI = 1.46) (LifeCanvas Technology) at room temperature until turning fully transparent. A table of antibodies and vendors can be found in Supplemeentary Table 13.

Confocal microscopy imaging

Tissue sections were imaged on a Zeiss LSM 880 confocal microscope using a Plan-Apochromat ×25/0.8 multi-immersion objective at a voxel size of 0.22 μm in XY and 2 μm in Z.

Light-sheet microscopy imaging

Whole thymic lobes were scanned using a Zeiss Z.1 Light-sheet microscope. The detection objective was an EC Plan-Neofluar ×5/0.16. Stacks were acquired at a resolution of 0.915 μm in XY and approximately 4.9 μm in Z. Dual-side images were fused using the maximum intensity option.

Image presentation

All images shown are processed using Imaris v.9.7.1 (Bitplane). Regions of interest in tissue sections are presented as 10-μm Z-projections.

Volume calculation of thymic regions

The total volumes of entire right thymic lobes were calculated using Imaris by generating the lobe surface from the tdTomato channel. Medullary regions were defined by high GFP+ cell density and Krt14+ and HD-TEC regions in aged thymus were defined by compacted GFP+ TECs. To identify medullary and HD-TEC regions in the images, we developed a pipeline in ImageJ (v.2.3.0/1.53f)77. For medulla, GFP and keratin 14 channels were combined. For HD-TEC regions, only the GFP+ channel was used. The images were filtered using two-dimensional (2D) median and three-dimensional (3D) Gaussian filtering and then binarized using a 2D min–max filter with thresholds set according to fluorescence intensity. The resulting image was used to extract the medullary or HD-TEC surface. The cortical volume was calculated as total thymic lobe volume minus the medullary and HD-TEC volume.

Segmentation of TEC nuclei

TEC nuclei were identified in confocal images of thymus sections using the spot detection function in Imaris. Total TEC spots were filtered and then TEC subsets were segmented according to the shortest distance to the indicated surface or the section edge. Medullary or HD-TECs were defined by the shortest distance to indicated surface ≤0 μm, subcapsular TECs were defined as shortest distance to section edge ≥−25 μm and the remaining spots were defined as cTECs.

Mean cell density was calculated by dividing the number of specific TEC subsets to the volume of different thymic regions.

Quantification of TECs

The number of nuclei in the various TEC subsets in the right lobe was calculated by multiplying the mean cell densities ascertained by confocal analysis of the slices from the left lobe by the volumes determined by light-sheet imaging of the right lobe.

Tissue preparation for sequencing

The scRNA-seq of FACS-sorted cell suspensions was performed on Chromium instrument (10x Genomics) following the user-guide manual (CG00052 Rev E) and using the Single Cell 3′ Reagent kit (v2). The viability of cells before loading onto the encapsulation chip was 73–98%, as confirmed with 0.2% (w/v) Trypan blue stain. Each sample, containing approximately 8,000 cells, was encapsulated in microfluidic droplets at a final dilution of 66–70 cells per µl (a multiplet rate ~3.9%). Following a reverse transcription step, the emulsion droplets were broken, barcoded cDNA purified with DynaBeads and amplified by 12 cycles of PCR: 98 °C for 180 s, 12× (98 °C for 15 s, 67 °C for 20 s, 72 °C for 60 s) and 72 °C for 60 s. The 50 ng PCR-amplified barcoded cDNA was fragmented with the reagents provided in the kit, purified with SPRI beads and the resulting DNA library was ligated to the sequencing adaptor followed by indexing PCR: 98 °C for 45 s; 12 × 98 °C for 20 s, 54 °C for 30 s, 72 °C for 20 s and 72 °C for 60 s. The final DNA library was double-size purified (0.6–0.8×) with SPRI beads and sequenced on an Illumina NovaSeq platform. Sequencing for Foxn1lacz and Foxn1tdTom was performed on an Illumina NextSeq.

Visium spatial gene expression slides were permeabilized at 37 °C for 12–18 min and polyadenylated. Messenger RNA was captured by primers bound to the slides. Reverse transcription, second-strand synthesis, cDNA amplification and library preparation proceeded using the Visium Spatial Gene Expression Slide and Reagent kit (10x Genomics, PN 1000184) according to the manufacturer’s protocol. After evaluation by real-time PCR, cDNA amplification included 11–12 cycles; sequencing libraries were prepared with eight cycles of PCR. Indexed libraries were pooled equimolar and sequenced on a NovaSeq 6000 in a PE28/120 run using the NovaSeq 6000 S1 Reagent kit (200 cycles; Illumina).

Library preparation and sequencing

After preparing our single-cell suspension solution, we utilized the library preparation and next-generation sequencing services offered by the University of Georgia’s Genomics and Bioinformatics Core to generate our scRNA-seq library. Ten thousand thymic stromal cells were loaded onto a 10x Genomics Chromium 3′ Single Cell Gene Expression Solution v3 microfluidics chip (10x Genomics) to generate an Illumina sequencer-ready library. Sequencing was then performed on an Illumina NextSeq 500/550, using four flow lanes that resulted in four BCL files that were shared with us using Illumina’s Basespace online platform.

Computational analysisMapping of single-cell and spatial transcriptome libraries

The scRNA-seq FASTQ files were processed with Cell Ranger (v.7.0.1) and Visium libraries were processed with Space Ranger (v.1.3.1) from 10x Genomics. All samples were mapped to the mouse mm10-2020-A genome assembly, except for the Foxn1tdTom dataset that was mapped to a custom mouse mm10-2020-A, including the sequences for the tdTomato gene and WPRE element (custom genome FASTA and index files for the tdTomato-WPRE sequence were downloaded from GSE125464).

Single-cell RNA-seq and spatial transcriptomics quality control and initial analysis

The Cell Ranger and Space Ranger-generated filtered_feature_bc_matrix.h5 files were processed following the guidelines on the shunPykeR GitHub repository78, an assembled pipeline of publicly available single-cell analysis packages put in coherent order, which allow data analysis in a reproducible manner and seamless usage of Python and R code. Genes that were not expressed in any cell, and also ribosomal and hemoglobin genes, were removed from downstream analysis. Each cell was then normalized to a total library size of 10,000 reads and gene counts were log-transformed using a pseudo-count of 1. Principal-component analysis (PCA) was applied to reduce noise before data clustering. To select the optimal number of principal components to retain for each dataset, the knee point (eigenvalues smaller radius of curvature) was used as a guide. Leiden clustering79 was used to identify clusters within the PCA-reduced data.

CD45− TBI series

The quality of the single cells was computationally assessed based on total counts, number of genes and mitochondrial and ribosomal fraction per cell, with low total counts, low number of genes (≤1,000) and high mitochondrial content (≥0.2) as negative indicators of cell quality (Supplementary Fig. 1). Cells characterized by more than one negative indicator were considered as low-quality cells. Although cells were negatively sorted before sequencing for the CD45 marker, a small number of CD45+ cells (expressing Ptprc), and also a few parathyroid cells (expressing Gcm2), were detected within our dataset (Supplementary Fig. 1). To remove bad-quality cells and contaminants in an unbiased way, we assessed the metrics at the cluster level rather than on individual cells. Leiden clusters with a low-quality profile and/or a high number of contaminating cells were removed. Finally, cells marked as doublets by scrublet80 were also filtered out. Overall, a total of 12,497 cells, representing 13.3% of our data, were excluded from further analysis (Supplementary Fig. 1).

After removal of low-quality and doublet cells, PCA (ncomps=45) and unsupervised clustering analysis was applied to the steady-state CD45− slice of the data using top highly variable genes (nhvgs = 3,500) and using Leiden (resolution = 0.3). Batch effect correction was performed using harmony81 with default parameters and using sample (Supplementary Fig. 1) as the batch key. Major cell lineages (epithelium, endothelium and fibroblast) were annotated based on canonical markers (Extended Data Fig. 2). Each major lineage was then sliced and reanalyzed in a similar fashion (epithelium, nhvgs = 3,500, ncomps = 50, harmony_key = ‘sample’, resolution = 1.4; endothelium, nhvgs = 3,500, ncomps = 30, harmony_key = ‘sample’, resolution = 0.1; fibroblast, nhvgs = 3,500, ncomps = 50, harmony_key = ‘sample’, resolution = 0.5) to interrogate these linages heterogeneity to a higher degree (Fig. 1a,c). Similarly for the steady-state CD45− data, the TBI CD45− slice of the data (nhvgs = 3,500, ncomps = 65, harmony_key = ‘sample’, resolution = 0.7) and their subsequent epithelial (nhvgs = 3,500, ncomps = 35, harmony_key = ‘sample’, resolution = 1.3), endothelial (nhvgs = 3,500, ncomps = 35, harmony_key = ‘sample’, resolution = 0.2) and fibroblast (nhvg s= 3,500, ncomps = 30, harmony_key = ‘sample’, resolution = 0.4) lineages were reanalyzed and annotated separately (Extended Data Fig. 10); however, when highly variable genes were calculated in the TBI setting, the day 1 TBI part of the data was excluded from the calculation due to the presence of a high number of inflammatory response genes. Finally, the steady-state and TBI slice annotations were transferred on the complete dataset (nhvgs = 3,500 (no day 1), ncomps = 65, harmony_key = ‘sample’) shown in Fig. 5i.

Foxn1tdTom data

The quality of the single cells was computationally assessed as described for the CD45− TBI series (nhvgs = 3,500, ncomps = 25, harmony_key = ‘sample’). A total of 4,062 cells, representing 23.0% of the data were excluded from further analysis. The epithelial cell lineage was sliced and reanalyzed further (nhvgs = 3,500, ncomps = 45, harmony_key = ‘sample’) to allow identification of smaller epithelial cell populations present in the epithelial compartment of the CD45− TBI series.

Foxn1LacZ data

The quality of the single cells was computationally assessed as described for the CD45− TBI series (nhvgs=3,500, ncomps = 35, harmony_key = ‘sample’). A total of 5,594 cells, representing 40.0% of the data, were excluded from further analysis. The epithelial cell lineage was sliced and reanalyzed further (nhvgs = 3,500, ncomps = 45) to allow identification of smaller epithelial cell populations present in the epithelial compartment of the CD45− TBI series.

Differential expression analysis

Differential expression analysis for comparisons of interest was performed using the sc.tl.rank_gene_groups() function from scanpy82 with the Wilcoxon rank-sum method83. In all cases, differentially expressed genes (DEGs) were considered statistically significant if the FDR-adjusted P value was ≤0.05. No fold change threshold was applied.

Generation of public gene signatures to characterize our steady-state subsets

We used the sc.tl.score_genes() function from scanpy81 (that calculates averaged scores based on cluster specific genes; scores are subtracted with a randomly sampled reference gene set) to generate gene signatures based on markers provided in the literature23,26,27 (Extended Data Fig. 2 and Supplementary Table 1) to assist annotation of our steady-state epithelial, endothelial and fibroblast subsets in the CD45− TBI series data (Supplementary Fig. 3).

Mapping of our steady-state subsets onto the TBI, Foxn1 tdTom and Foxn1 LacZ data

We used scanpy’s sc.tl.score_genes() function with the top 20 DEGs from the steady-state defined subsets (Wilcoxon, FDR ≤ 0.05, sorted in descending order by Wilcoxon z-score; Supplementary Table 3) to generate unique cell type subset signatures, which we mapped to the respective lineage subsets in the TBI (days 1, 4 and 7; Extended Data Fig. 10), Foxn1tdTom and Foxn1LacZ data.

Public datasets reanalysis

Re-analysis of single-cell transcriptome datasets from public nonhematopoietic mouse and human thymic samples (CD45− populations) were processed as described in the ‘Single-cell RNA-seq and spatial transcriptomics analysis’ section of the Methods. The Data Availability section provides a complete list of the raw count data files used as the entry point for each dataset reanalysis.

ThymoSight

ThymoSight is an R Shiny app that we have developed to allow interactive exploration of all mouse and human publicly available single-cell datasets of the nonhematopoietic thymic stroma. Mouse datasets included are from refs. 5,6,7,10,11,12,13,14,15,16,17,18 and our own data are from this manuscript. Human datasets included are from refs. 9,19,20. ThymoSight also provides dataset metadata fields (if available/applicable) such as tissue, age, stage, sorted cell population, gender, genotype, treatment, linked publication, mapped annotation based on our own subset signatures and original annotation. The app.R code that launches the app together with the Python notebooks used to create consistent annotation fields, reanalyze and integrate the public datasets with ours have been submitted on GitHub (https://github.com/FredHutch/thymosight). The server hosting the interactive app can be accessed at www.thymosight.org.

Single-cell RNA-seq meta-analysis Integration of CD45− steady-state data with the Visium data

We used scanorama (v.1.7.2)84 to integrate our scRNA-seq datasets with our spatial transcriptomic data. Integration was performed between age-matched data at steady state and with default parameters using scanpy’s example tutorial (https://scanpy.readthedocs.io/en/stable/tutorials/spatial/integration-scanorama.html#integrating-spatial-data-with-scrna-seq-using-scanorama).

RNA velocity analysis

Velocyto (v.0.17.17)85 was used to generate loom files, which we subsequently merged with our already-annotated single-cell object. We performed RNA velocity analysis within the thymic epithelium compartment of our data using scVelo (v.0.2.4)86 in stochastic mode.

Pseudotime analysis

We used Dynamo’s (v.1.4.0)62 VectorField() function with the given parameters (basis = ‘umap’, M = 1,000, MaxIter = 170, pot_curl_div = True) to calculate the vector field and to estimate the negative of the single-cell potential (ddhodge potential; Dynamo’s version of pseudotime) of the thymic epithelia in 18 mo mice at steady state.

Cell fate prediction analysis

We used Dynamo’s topography (basis = ‘umap’) and fate (interpolation_num = 100, direction = ‘forward’, inverse_transform = false, average = false) functions to create fate prediction animations for our 18-mo epithelial dataset at steady state, setting each of our epithelial subsets as the progenitor population each time. To visualize the fate transition animation results in a static format we leveraged CellRank’s (v.2.0.3.dev10+g4ae88b9)87 built plot_single_flow() module using the already-calculated Dynamo’s ddhodge potential (binned) to create vein plots resembling fate transition and relative frequency of the epithelial subsets.

Pathway enrichment analysis

Pathway enrichment analysis was performed with GSEA (v.4.3.2)37 according to the gene list and rank metric provided. The GSEA preranked module was used to predict pathway enrichment in threshold-free comparisons: (1) 18-mo versus 2-mo subsets at steady-state and (2) aaTEC1s and aaTEC2s versus other TECs. We created rankings for all DEGs using the Wilcoxon z-score in descending order. Predicted pathways with an FDR ≤ 0.05 were considered significantly enriched.

Network analysis

Network analysis of the significantly enriched GSEA pathways from comparisons of interest was performed using Cytoscape (v.3.10.0)38. We used the EnrichmentMap() function to organize enriched pathways (FDR ≤ 0.05) with a high overlap of genes (default cutoff similarity of 0.375) in the same network allowing for a simplified and intuitive visualization of the distinct processes that are significantly represented in each subset at steady state. This facilitated interpretation of the enriched pathways from the plethora of comparisons and allowed categorization of all resulting pathways into networks based on the overlap of the genes contributing to the pathway’s enrichment (Supplementary Table 5). Manual inspection of the resulting networks allowed allocation of network-related annotations. Individual pathways that were not part of an existing network were manually annotated to the existing categories based on their biological function or grouped under ‘singlet’.

Cell–cell interaction analysis

CellChat (v.1.4.0)50 was used with default parameters to predict cell–cell interactions between all CD45− subsets within the 2-mo and 18-mo cohorts at steady state and at days 1, 4 and 7 after damage against the complete CellChat database. Cell subsets with fewer than 20 cells were excluded from the interactome analysis. For comparisons between the individual TBI time points and age cohorts, individual CellChat objects were integrated using the mergeCellChat() function. For the circus plots shown in Fig. 4h, some of the cell type subsets were grouped together: ECs (aEC, capEC and vEC), FBs (capsFB, intFB, medFB, nmSC and Fat), mTECprol/mTEC2s and mimetic (basal, tuft, neuro, goblet and M cell).

CD45− bulk RNA-seq preprocessing and downstream data analysis Quality control, alignment and gene count quantification

Quality control of the raw read files (FASTQ) was performed using the FastQC tool88. Low-quality reads and adaptor contaminants were removed using Trimmomatic89 (default parameters for paired-end reads) and post-trimmed reads were reassessed with FastQC to verify adaptor removal and potential bias introduced by trimming. The quality control-approved reads were aligned to the GRCm38.p5 (mm10) mouse genome assembly (GENCODE; M12 release) with STAR90 aligner using default parameters and–runThreadN set to 32 to increase execution speed. The STAR-aligned files were then used as input to the featureCounts91 tool (default parameters) to quantify gene expression levels and construct the count matrix.

Low gene count removal and library size normalization

The raw count matrix was converted to a DGEList object in R using the readDGE() function from the edgeR92 package. Lowly expressed genes were removed using the filterByExpr() function for the groups of interest before comparison with the scRNA-seq datasets.

Bulk RNA-seq versus scRNA-seq

scRNA-seq sample reproducibility was verified using bulk RNA-seq data for the CD45− sorted populations. Comparison between bulk and scRNA-seq CD45− transcriptional profiles was performed by computing Pearson’s correlation between log10-transformed raw bulk counts (per biological replicate) and log10-transformed averaged raw single-cell counts (per technical replicate) for the relevant datasets across the TBI timeframe (Supplementary Fig. 2).

Statistics and reproducibility

All statistics were calculated, and display graphs were generated, in GraphPad Prism.

Specific statistical tests used have been highlighted in the figure legends but briefly, statistical analysis between two groups were performed on biological replicates (individual mice) with a two-tailed Mann–Whitney or two-tailed t-test and, where appropriate, a two-tailed paired t-test. Statistical comparison between three or more groups was performed on biological replicates (individual mice) with a Kruskall–Wallis test with Dunn’s correction, one-way analysis of variance with Dunnett’s correction or two-way analysis of variance with Šídák correction.

The imaging studies in Fig. 2a,b were performed independently three times (n = 1–3 mice per experiment). For the images in Fig. 2e, PanK was performed independently three times (n = 1–3 mice per experiment), Krt5 was performed independently four times (n = 1–3 mice per experiment), Krt8 was performed once (n = 2 mice), Krt14 was performed independently seven times (n = 1–3 mice per experiment) and claudin-3 was performed independently four times (n = 1–4 mice per experiment); with one section imaged for each mouse. The studies described in Fig. 3i were performed independently ten times (n = 1–3 mice per experiment) with one section imaged per animal. Figure 3b was performed independently three times (n = 4 per group). In Extended Data Fig. 7b, staining for FOXN1 was performed independently twice (n = 7 per experiment). In Extended Data Fig. 9a, DCLK and UEA1 were performed once (n = 2 mice per experiment) with one section imaged per animal.

Reporting summary

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