Developmental trajectory and evolutionary origin of thymic mimetic cells

Animals

C57BL/6 mice were maintained at the Max Planck Institute of Immunobiology and Epigenetics. Foxn1^−/− (ref. ²⁷), Foxn1:cre⁵², Rosa26-LSL-EYFP⁵³, Foxn1:mCardinal⁵⁴, Ascl1^fl/fl (ref. ⁵⁵), Foxn1:BlFoxn4, Foxn1:CmFoxn4, Foxn1:CmFoxn1³⁸ and Foxn1:Fgf7²⁰ transgenic mice have been described previously. The Foxn1:Bmp4 transgene was created by T. Schlake and B.K. by inserting a cDNA fragment corresponding to nucleotides 497–1729 in GenBank accession number NM_007554.3 as a NotI fragment into pAHB14³⁴. The Δ3ex2 Foxn1 deletion mutant (internal designation Chi6) transgene was generated by deletion of nucleotides 504–728 of mouse Foxn1 cDNA (Genbank accession number NM_008238.2) and insertion as a NotI fragment into pAHB14³⁴. To generate transgenic mice, constructs were linearized and injected into FVB pronuclei according to standard protocols. The Foxn1^Δ3ex2 mice were bred to a Foxn1-deficient background. Genotyping information is summarized in Supplementary Table 2. Mice carrying the original nu mutation (CByJ.Cg-Foxn1nu/J) were purchased from Charles River and used for snRNA-seq experiments. Mice were analysed at the age of 4–6 weeks, unless otherwise stated.

The zebrafish line carrying an internal deletion of the foxn1 gene has been described⁵⁶. Adult zebrafish (3 months of age) were used for experiments.

Mice and zebrafish were kept in the animal facility of the Max Planck Institute of Immunobiology and Epigenetics under specific pathogen-free conditions (mice: 14 h light, 10 h dark; temperature 22 ± 2 °C; relative humidity 55 ± 10%; zebrafish: 13 h light, 11 h dark, water temperature 28 °C). 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 the Regierungspräsidium Freiburg, Germany (mice: licenses 35–9185.81/G-12/85; 35–9185.81/G-16/67; zebrafish: license 35–9185.81/G-14/41). All strains are made available from the corresponding author upon request, subject to standard material transfer agreements.

Ammocoete larvae of L. planeri (body length, 8–10 cm) were caught by a licensed fisherman from the wild in the Freiburg region (Riedgraben, March-Neuershausen) under permission by the local governmental authority (Landratsamt Breisgau-Hochschwarzwald, license 420.1.13-2024-034414). Juvenile cat shark (S. canicula) specimens were kindly supplied by Markéta Kauka (Max Planck Institute for Evolutionary Biology, Plön, Germany); juvenile bamboo sharks (Chiloscyllium punctatum) were purchased from a local pet shop. Upon arrival at the laboratory, lampreys and sharks were euthanized using 0.02% tricaine methanesulfonate. For RNA isolation, tissues were removed under a dissection microscope and dissolved in TRI reagent (T9424, Sigma-Aldrich); for histological studies, dissected tissues were fixed in 4% neutral formalin.

Sample sizes were based on our experience and accepted practice in the respective fields, balancing statistical robustness, resource availability and animal welfare. No statistical methods were used to predetermine sample size. Provided the transgenic status and age matched the experimental requirements, mice were randomly assigned to experimental groups, irrespective of sex. Blinding was not possible because the thymus phenotype—the transgenic status of the respective mouse—is evident from flow cytometry, imaging analysis or genotyping information.

RNA in situ hybridization

Tissues were fixed with modified Davidson’s fluid and dehydrated with Li-ethanol and ethanol in a stepwise manner and finally embedded in paraffin; RNA ISH on paraffin sections was performed using DIG-labelled probes⁵⁷. Double ISH was carried out as follows. DIG- and fluorescein-labelled RNA antisense probes were simultaneously hybridized to RNA in tissue sections. The DIG-labelled probe was detected first, either with an alkaline-phosphatase-conjugated anti-DIG antibody (1:2,000 dilution in maleic acid buffer (MAB); 100 mM maleic acid, pH 7.5, 150 mM NaCl, 2 mM Levamisol, 1% blocking reagent (Roche), 0.1% Tween-20) for chromogenic detection, or with an anti-DIG-POD antibody (1:300 dilution in MAB) for fluorescent detection. The presence of the DIG-hapten was revealed by staining with BM Purple (Roche) or by Cy3 fluorescence, using the Tyramide Signal Amplification Plus system (AkoyaBioscience). The sections were washed several times in PBS; the fluorescein-labelled probe was detected by a peroxidase-conjugated anti-Fluorescein-POD (1:300 dilution in MAB) and revealed by Cy5 fluorescence.

Sequence coordinates in GenBank accession numbers for the probes were as follows:

• M. musculus: Foxn1, nucleotides 2181–3584 in XM_006532266.3; Foxj1, nucleotides 1406–1917 in NM_008240.3; Hnf4a, nucleotides 1469–1920 in NM_008261.3; Foxi1, nucleotides 1067–1576 in NM_023907.4; Pax9, nucleotides 1260–1818 in NM_011041.3.

• S. canicula: RAG1, nucleotides 3400–3900 in XM_038808256; FOXN1, nucleotides 101–600 in XM_038813127; FOXN4, nucleotides 2001–2500 in XM_038815571; CD3E, nucleotides 1–594 in KY434199; FOXJ1, nucleotides 1501–2000 in XM_038820881; FOXI1, nucleotides 1231–1730 in XM_038822317; MYOG, nucleotides 620–1120 in XM_038820992; POU2F3, nucleotides 1401–1900 in XM_038779605; TTR, nucleotides 180–680 in XM_038808583.

• L. planeri: CDA1, nucleotides 8–496 in MG495252; TTR, nucleotides 148–593 in XM_061573254; MyHC1, nucleotides 5013–5609 in AB126173.

Image analysis

Images were acquired on Zeiss microscopes (Axioplan 2) equipped with an Mrc 5 camera; in some figure panels, Cy5 signals were converted to false (yellow) colour for better visualization.

Flow cytometry and cell sorting

Single-cell suspensions of TECs for preparative flow cytometry were obtained as described^31,58. Thymic epithelial cells have the surface phenotype EPCAM⁺CD45⁻; thus, cell surface staining was performed using anti-EPCAM (G8.8), conjugated with APC (1:1,000, BioLegend) or anti-EPCAM (G8.8), conjugated with biotin (1:1,000, BioLegend), in combination with streptavidin, conjugated with eFluor 450 (1:1,000, eBioscience), and anti-CD45 (30-F11), conjugated with PE Cy7 (1:2,000, BioLegend) at 4 °C in PBS supplemented with 0.5% BSA and 0.02% NaN₃. In order to differentiate cells with past and acute Foxn1 expression, triple-transgenic Foxn1:cre; Rosa26-LSL-eYFP; Foxn1:mCardinal mice were used for cell sorting. Cells with past expression of Foxn1 were sorted as EPCAM⁺YFP⁺ mCardinal⁻ cells, whereas cells with acute Foxn1 expression were sorted as EPCAM⁺YFP⁺ mCardinal⁺ cells. EPCAM⁺CD45⁻ cells (after negative enrichment using anti-CD45 magnetic-activated cell sorting (MACS) beads and anti–Ter-119 MACS beads, Miltenyi Biotec) were sorted directly into TRI reagent (T9424, Sigma-Aldrich). Cell sorting was carried out using the MoFlow instrument (Dako Cytomation-Beckman Coulter) controlled with the Summit (5.5) software. Analytical flow cytometry was performed for TECs as follows: anti-EPCAM (G8.8), conjugated with APC (1:1,000, BioLegend); anti-Ly51 (alias BP-1; 6C3), conjugated with PE (1:1,600, eBioscience); UEA1, conjugated with FITC (1:1,000, Vector Labs) or UEA1, conjugated with biotin (1:600, Vector Labs), in combination with streptavidin, conjugated with eFluor 450 (1:1,000, eBioscience). When analysis of haematopoietic fractions was desired, thymocyte suspensions were prepared in parallel by mechanical liberation, best achieved by gently pressing thymic lobes through 40 μm sieves. Cell surface staining (anti-CD45 (30-F11), conjugated with PE/Cy7 (1:2,000, BioLegend); anti-CD4 (GK1.5), conjugated with FITC (1:1,000, BioLegend); anti-CD8a (53-6.7), conjugated with APC (1:800, eBioscience); anti-TCRβ (H57-597), conjugated with PE (1:400, eBioscience); anti-CD19 (eBio1D3), conjugated with PerCP/Cy5.5 (1:500, eBioscience) or PE/Cy7 (1:1,000, eBioscience); anti-B220 (alias CD45R; RA3-6B2), conjugated with biotin (1:200, eBioscience); anti-IgM (II/4.1), conjugated with PE (1:300, eBioscience), anti-CD93 (alias C1qRp; AA4.1), conjugated with APC (1:300, eBioscience); streptavidin conjugated with eFluor 450 or FITC (1:1,000, eBioscience)) was performed at 4 °C in PBS supplemented with 0.5% BSA and 0.02% NaN₃. Flow cytometry experiments were evaluated using FACSDiva (8.0.2) and FlowJo (9.3.1) software. The relevant gating strategies^20,37,38 are shown in Supplementary Fig. 4.

Isolation of cell nuclei from thymus tissue

Thymus tissues were recovered under a dissection microscope; thymocytes were mechanically liberated by applying gentle pressure on the tissue on a 40-μm sieve and extensive washing to deplete as many haematopoietic cells as possible. For nude mice, thymic rudiments from two individuals were pooled to constitute one sample. The resulting tissue remnants were transferred into a 1.5 ml Eppendorf tube containing 150 µl of freshly prepared ice-cold lysis buffer (10 mM Tris-HCl pH 8, 10 mM NaCl, 3 mM MgCl₂, 0.1% non-denaturing detergent Igepal CA-630, 0.1% Tween-20, 1% BSA), supplemented with 1 U μl⁻¹ Protector RNase inhibitors (Roche, 3335402001). Tissues were homogenized using a RNAse-free disposable pestle for 1.5 ml tubes (Fisher Scientific, 12141364) with a circular motion until the majority of clumps were homogenized. Nuclei quality was assessed under a phase-contrast microscope. Subsequently, 500 µl of ice-cold wash buffer (10 mM Tris-HCl pH 8, 10 mM NaCl, 3 mM MgCl₂, 0.1% Tween-20, 1% BSA) with 0.2 U μl⁻¹ of Protector RNase inhibitors was added, mixed by inversion, and spun briefly to collect liquid from the cap of the tube. The nucleus suspension was then filtered through a 70 µm Flowmi Cell Strainer (H13680-0070, Bel-Art) and centrifuged at 300g for 5 min at 4 °C. The supernatant was removed, and the nuclear pellet was resuspended in 200 μl of 0.5% BSA in PBS, supplemented with 0.2 U μl⁻¹ of Protector RNase inhibitors and filtered again through a 40-µm Flowmi Cell Strainer (H13680-0040, Bel-Art). Nuclei were quantified by Trypan blue staining using the Countess 3 automated cell counter (Thermo Fisher Scientific). Only Trypan blue-positive nuclei were counted (to exclude fat droplets). Each nucleus suspension was normalized to a concentration of 800 nuclei per µl using the resuspension buffer.

PCR with reverse transcription

RNA was isolated from total thymus tissue using TRI reagent (T9424, Sigma-Aldrich). Thymus tissue was isolated from specimens of unspecified sex from juvenile brown-banded bamboo shark (C. punctatum)¹³ and foxn1^+/+ and foxn1^−/− adult zebrafish⁵⁶ on the ikzf1:EGFP background⁴³. For each RNA extraction from zebrafish thymi, three organs were pooled and a total of three such pools were processed. RNA was quantified using the Qubit RNA HS Assay Kit (ThermoFisherScientific, Q32852) and the Qubit 4 Fluorometer (ThermoFisherScientific, Q33226). RNA quality was checked by determining the 18S/28S rRNA ratio using the Fragment Analyzer RNA Kit (ThermoScientific, DNF-471-0500) and the 5200 Fragment Analyzer System (ThermoScientific, M5310AA). cDNA was prepared using random hexamers and the SMARTScribe Reverse Transcriptase (Clontech).

The primers used for amplification from bamboo shark cDNA and the respective amplicon sizes are listed in Supplementary Table 3; the primers used for amplification from zebrafish cDNA and the respective amplicon sizes are listed in Supplementary Table 4. In the PCR reaction, primers were used at a concentration of 10 pmol μl⁻¹. For comparative analyses, the amounts of input cDNA were calibrated by comparison to zebrafish ef1a expression levels. Amplicons, independently generated for three thymic samples each per genotype, were pooled prior to gel electrophoresis. For gel source data, see Supplementary Fig. 8.

Bulk RNA-seq of TECs

TECs derived from mouse thymus were sorted directly into TRI reagent (T9424, Sigma-Aldrich), while for zebrafish total thymus tissue was used. RNA isolation was performed according to standard protocols. RNA was quantified using the Qubit RNA HS Assay Kit (ThermoFisherScientific, Q32852) and the Qubit 4 Fluorometer (ThermoFisherScientific, Q33226). RNA quality was checked by determining the 18S/28S rRNA ratio using the Fragment Analyzer RNA Kit (ThermoScientific, DNF-471-0500) and the 5200 Fragment Analyzer System (ThermoScientific, M5310AA). For each zebrafish library, RNA from three animals was pooled and three such pools were sequenced. Libraries were prepared using the Ultra RNA Library Prep Kit (NEB). Samples were run on Illumina HiSeq 2500, HiSeq 3000 or NovaSeq 6000 instruments and sequenced to a depth of 10 × 10⁶ to >150 × 10⁶ reads per sample. Transcriptomes were analysed on the Galaxy platform⁵⁹ using Trim Galore! version 0.4.3.1 (developed by Felix Krueger at the Babraham Institute), HISAT2 version 2.1.0⁶⁰ and featureCounts version 1.6.1.0⁶¹.

snRNA-seq of thymic tissue

The Chromium GEM-X Single Cell 3′ v4 protocol (CG000731, Rev B) was followed starting from step 1.1 according to the manufacturer’s guidelines. The Chromium GEM-X Single Cell 3′ Kit v4 (PN-1000686) and the Chromium GEM-X Single Cell 3′ Chip Kit v4 (PN-1000690) were used to process each sample. A total of 36.6 µl of normalized nucleus suspension (see ‘Isolation of cell nuclei from thymus tissue’) was added to the master mix to target a recovery of 20,000 nuclei. The ready GEM-X chip was loaded onto a 10x Genomics Chromium Xo instrument, and final libraries were completed as per 10x Genomics guidelines. Libraries were quantified using the Qubit High Sensitivity DNA Assay (Invitrogen, Q32851), and their molar concentrations were calculated on the basis of their size distribution using the Fragment Analyzer NGS 1–6,000 bp High Sensitivity DNA kit. Libraries were pooled, cleaned of adapter dimers and denatured according to Illumina guidelines. Libraries were sequenced in paired-end mode with a read length of 10 × 100 × 100 × 10 bp (i5 index, R1, R2, i7 index) on a NovaSeq 6000 instrument (Illumina), aiming for a minimum sequencing depth of 20,000 read pairs per nucleus.

Analysis of snRNA-seq data

Counts were generated with 10x Genomics Cell Ranger (8.0.1) and background noise was reduced with CellBender (0.3.0, –fpr 0.01)⁶². Droplets were excluded if they exhibited low total UMI counts (<500), high proportions of mitochondrial counts (>5%), or abnormal numbers of detected genes (beyond 1 median absolute deviation (MAD) from the median) or complexity (ratio of detected genes over total count; beyond 3 MADs from the median). Doublets were removed with scDblFinder⁶³. After initial dimensionality reduction and clustering of the data, each cell was scored for mimetic and background marker gene sets with AUCell²¹. Clusters predominantly comprising cells scoring highly for signatures compatible with expected contaminants (thymocytes and adipocytes for control and nude samples, respectively) were excluded from further analysis. The reduced dataset was clustered again. Parameters were selected to maximize the mean silhouette width. Subclustering was performed on clusters selected manually based on their mean silhouette width, again performing parameter sweeps and selecting parameters based on the total within cluster sum of squares and mean silhouette widths of subclusters. The final clusters were annotated by investigation of dominant signature scores and marker genes.

Assessment of mimetic signatures in bulk RNA-seq data

Bulk RNA-seq data from all samples were processed jointly. Genes were included if their counts exceeded 5 in at least 50% of samples. Differential expression of genes between samples was assessed with voom⁶⁴ and limma with treat(., lfc = log2(1.2), robust = TRUE)^65,66 to generate the t-statistic. Relative enrichment of gene signatures was determined with the competitive gene set test camera²⁵ with cameraPR, employing the t-statistic as output by treat as the ranking statistic. The output of camera is a two-sided and Benjamini–Hochberg adjusted P value. Results are reported as log₁₀(adj. P), conditionally multiplied by −1 for signatures with an upwards direction of change. Although we were mainly interested in TEC signatures^2,20, we also collected background signatures of unrelated cell types from PanglaoDB⁶⁷, Tabula Muris⁶⁸ (bladder, spleen, kidney, liver, marrow, muscle, lung, non-myeloid brain, from FACS) and MSigDB^69,70 (M8: mouse cell-type signature gene sets, except for Tabula Muris Senis signatures). From the collection of signatures, those that comprised up to 10% less or 10% more genes than the smallest and largest TEC or mimetic signature, respectively, were retained as background signatures (see Supplementary Fig. 1). To investigate mouse-derived signatures in RNA-seq samples from D. rerio, we performed orthology transfer using data from the Ensembl database (release 112)^71,72, allowing many-to-many relationships.

Identification of mimetic cells in scRNA-seq data

Signature gene sets for mimetic cells were derived from single-cell differential gene expression results from Supplementary Table 2 in Michelson et al.². The gene lists were filtered (adjusted P < 0.01, |log₂(FC)| >1, expressed in at least 10% of cells) and used as input to AUCell²¹ to score signatures in each cell. Based on the resulting histograms of areas under the curve (AUCs), thresholds were defined manually for each mimetic signature (see Extended Data Fig. 1a). A cell that met thresholds for multiple signatures was disambiguated by ordering all eligible signatures and choosing the highest-ranked one, first by the cell’s rank within the signature AUCs (that is, a signature is ranked higher if the cell’s position on the histogram is further to the right), and second (in case of ties) by the z-score standardized AUC. The signatures ‘Tuft1’ and ‘Tuft2’, as well as ‘Skin_basal’ and ‘Skin_keratinized’ led to the identification of an overlapping set of cells which were therefore collapsed into unified Tuft and Skin populations (Extended Data Fig. 1b). Using the final signatures, mimetic cells were identified in the uniform manifold approximation and projection (UMAP) graphs calculated from single-cell datasets²⁰.

Shifts in population composition were evaluated with scCODA²⁴. The model was run multiple times, using each of the major populations (early and postnatal progenitors, mTEC, cTEC, Aire-stage and unassigned) as the reference. Detected changes were minimal for the ‘unassigned’ population and it was therefore used as the final reference. A population was deemed ‘overall credibly changed’ if more than half the models detected a credible change.

Lineage tracing of mimetic cell lineages

CRISPR–Cas9 barcodes for single cells²⁰ were binarized (barcode present or absent) and a logistic regression model was fitted with fixed covariates ‘signature’ and ‘sample’. A P value for the ‘signature’ term was determined by likelihood ratio test with a reduced model. Proportions of Foxn1-expressing cells between populations were compared with a binomial test, implemented in the findMarkers function from the scran package⁷³.

Determination of fate probabilities in scRNA-seq data

scRNA-seq data from barcoded and unbarcoded samples from three time points (embryo, newborn, 4 week) were integrated with fastMNN⁷⁴. The integrated dataset was re-clustered and clusters with expression of Pth/Chga (ectopic parathyroid) and Cd3/Cd4/Cd8a (thymocytes) were excluded. Diffusion pseudotime was calculated, starting from the cell with the highest AUC for the early progenitor signature in the embryonic sample⁷⁵. The pseudotime and integrated representation of the data were used for generation of a PseudotimeKernel and ConnectivityKernel with CellRank^29,30, respectively. Both kernels were combined and analysed with a GPPCA estimator⁷⁶ to identify between 10 and 16 macrostates (gppca.fit(n_states = [10, 16])). The optimal number of macrostates based on the minChi criterion⁷⁶ was 12. Macrostates were manually labelled according to their composition of annotated cell types and defined as initial (early progenitors) or terminal (Aire-stage and mimetic populations) states. Intermediate states were not retained. Fate probabilities were calculated with the remaining macrostates and the overlap of cells exceeding fate probability thresholds in pairs of fates was quantified by the Jaccard index. Because the automatically determined macrostates did not reflect all mimetic populations, we additionally ran CellRank after manually identifying terminal cells, picking representatives with the highest signature AUC for each population of interest. Fate probabilities and Jaccard indices were also calculated for these data.

Data handling and statistics

Box plots encapsulate the first to third quartile, a line indicates the median. Whiskers extend to the furthest point with a distance of up to 1.5 times the interquartile range from the boxes. For the data analyses presented in Figs. 2a,b and 3e and Extended Data Figs. 6 and 8a, differences in means between groups were compared by ANOVA. Separate ANOVAs were conducted for each dependent variable. The resulting P values were corrected for multiple testing with the conservative Bonferroni method. For tests with adjusted P < 0.05, we rejected the null hypothesis of equal means and performed pairwise two-sided t-tests between all (Figs. 2a,b and 3e and Extended Data Fig. 6a) groups, or between selected groups (Extended Data Fig. 8, F₂ samples only). P values from pairwise t-tests were corrected for multiple testing with the conservative Bonferroni method. For Extended Data Figs. 2d and 4d, groups were compared via a two-sided t-test without (Extended Data Fig. 2d) or with (Extended Data Fig. 4d) Welch correction (GraphPad Prism 9.5.1). For LOESS (Fig. 3a), we tested span parameters between 0.4 and 0.9 in increments of 0.01 with tenfold cross-validation. The minimal mean squared error was achieved with a span of 0.53. The numbers of biological replicates are indicated in the figure panels or figure legends.

Reporting summary

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