Age-related epithelial defects limit thymic function and regeneration

The thymus is essential for establishing adaptive immunity yet undergoes age-related involution that leads to compromised immune responsiveness. The thymus is also extremely sensitive to acute insult and although capable of regeneration, this capacity declines with age for unknown reasons. We applied single-cell and spatial transcriptomics, lineage-tracing and advanced imaging to define age-related changes in nonhematopoietic stromal cells and discovered the emergence of two atypical thymic epithelial cell (TEC) states. These age-associated TECs (aaTECs) formed high-density peri-medullary epithelial clusters that were devoid of thymocytes; an accretion of nonproductive thymic tissue that worsened with age, exhibited features of epithelial-to-mesenchymal transition and was associated with downregulation of FOXN1. Interaction analysis revealed that the emergence of aaTECs drew tonic signals from other functional TEC populations at baseline acting as a sink for TEC growth factors. Following acute injury, aaTECs expanded substantially, further perturbing trophic regeneration pathways and correlating with defective repair of the involuted thymus. These findings therefore define a unique feature of thymic involution linked to immune aging and could have implications for developing immune-boosting therapies in older individuals.

The thymus is essential for establishing adaptive immunity yet undergoes age-related involution that leads to compromised immune responsiveness.The thymus is also extremely sensitive to acute insult and although capable of regeneration, this capacity declines with age for unknown reasons.We applied single-cell and spatial transcriptomics, lineage-tracing and advanced imaging to define age-related changes in nonhematopoietic stromal cells and discovered the emergence of two atypical thymic epithelial cell (TEC) states.These age-associated TECs (aaTECs) formed high-density peri-medullary epithelial clusters that were devoid of thymocytes; an accretion of nonproductive thymic tissue that worsened with age, exhibited features of epithelial-to-mesenchymal transition and was associated with downregulation of FOXN1.Interaction analysis revealed that the emergence of aaTECs drew tonic signals from other functional TEC populations at baseline acting as a sink for TEC growth factors.Following acute injury, aaTECs expanded substantially, further perturbing trophic regeneration pathways and correlating with defective repair of the involuted thymus.These findings therefore define a unique feature of thymic involution linked to immune aging and could have implications for developing immune-boosting therapies in older individuals.
Thymic T cell differentiation requires the close interaction between thymocytes and the supporting stromal microenvironment, which is composed of highly specialized TECs, endothelial cells (ECs), mesenchymal cells like fibroblasts (FBs), dendritic cells, innate lymphoid cells and macrophages 1 ; however, thymic function is not static over lifespan with a well-described decline in function that accelerates upon puberty and is characterized by tissue atrophy, disrupted stromal architecture, reduced export of new naive T cells and, ultimately, diminished responsiveness to new antigens 2 .In addition to its chronic functional decline with age, the thymus is also extremely sensitive to acute damage induced by routine insults such as stress and infection, but also more severe damage such as that caused by common cancer therapies such as cytoreductive chemo-or radiation therapy 2,3 .Although the thymus harbors tremendous capacity for endogenous repair, this regenerative capacity also declines with age 2,3 .This deficiency manifests most prominently after the thymic damage caused by the conditioning regimens required for hematopoietic cell transplantation (HCT), which leads to prolonged T cell lymphopenia, an important contributor to transplant-related morbidity and mortality due to infections and malignant relapse 2,3 .In fact, both pre-and post-transplant thymic function can be positive prognostic indicators of HCT outcomes 4 .Recent advances in single-cell technology have provided new insights into the heterogeneity of TECs in young and aged mice and in humans [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] ; how TECs orchestrate T cell differentiation and how dysfunction in these processes is linked to autoimmunity and immunodeficiency; however, perhaps as a consequence of this complexity, the mechanisms underlying thymic involution and regeneration after damage remain poorly understood 2,3 .
Using single-cell and spatial transcriptomics, lineage-tracing and advanced imaging, we report here the age-associated emergence of unique TEC states linked with tissue degeneration.aaTECs formed atypical high-density epithelial clusters that were devoid of thymocytes; an accretion of nonfunctional thymic tissue that worsened with age and exhibited features of partial epithelial-to-mesenchymal Article https://doi.org/10.1038/s41590-024-01915-9 all ECs showed enrichment for lymphatic markers (Supplementary Fig. 3c), consistent with the scarcity of lymphatic ECs observed in sections or flow cytometric analysis of thymi from adult mice 28 .Expression of Plvap and Cldn5 delineated capEC 29 , whereas vECs demonstrated the highest expression of Vwf and Vcam1 (Fig. 1c,d, Extended Data Fig. 2c and Supplementary Table 1), genes that correlate with the vessel diameter.Venous ECs expressed high levels of P-selectin (Selp) (Fig. 1c,d, Extended Data Fig. 2c and Supplementary Table 1), indicative of thymic portal ECs (TPECs), which have been linked to homing of hematopoietic progenitors in the thymus 30 .Notably, Bmp4, which is produced by thymic ECs and is involved in regeneration after acute insult 31 , was expressed highest by vECs (Fig. 1c,d, Extended Data Fig. 2c and Supplementary Table 1).
Comparison of the populations from young (2-mo) and aged (18-mo) mice revealed considerable overlay, with no major new cell populations emerging with age in the endothelial or fibroblast compartments, although transcriptional differences were observed within existing cell populations suggesting an age-associated change in cell state rather than cell type (Fig. 1e and Supplementary Table 2).In contrast, within the TEC compartment we observed two distinct age-associated epithelial cell types, referred henceforth as aaTEC1 and aaTEC2 (Fig. 1e).These were apparent only in samples from 18-mo mice and mapped to only one previously published dataset 12 , likely due to the age and cell selection protocols used (Extended Data Fig. 3b).Assessing quantitative changes in cell subsets with age, we found that the two aaTEC populations exhibited the greatest changes with age, but we also observed a considerable increase in early prog (Fig. 1f).Using our sequencing dataset as well as previous publications to design flow cytometry panels to examine these stromal populations 6 , we found little change in the frequencies of EC subsets or their number with age, although there were increased numbers of medFBs, PCs and vSMCs (Extended Data Fig. 4a-c).This age-related increase in mesenchymal-lineage cells may reflect that observed in other tissues where fibrosis is a key driver of functional decline 32 .Within the epithelial compartment, we observed decreases in all mTECs, including TEC tuft cells, but in contrast we observed a marked increase of an atypical TEC population with age that expressed EpCAM but was negative for Ly51 and UEA1 (features of cTECs or mTECs, respectively) (Fig. 1g), consistent with the loss of mTEC-or cTEC-specific markers within the aaTEC1 subset.This population of Ly51/UEA1 double negative (DN)-TECs gradually emerges and expands across lifespan (Fig. 1h and Extended Data Fig. 4d).Differential expression analysis between each aaTEC subset versus all other TECs identified claudin-3 (Cldn3) and podoplanin (Pdpn) as potentially distinctive, in conjunction with established markers (Fig. 1i and Extended Data Fig. 4e).Flow cytometric analysis confirmed their utility, resolving populations of both aaTEC1 (Epcam + MHCII + Ly51 − UEA1 − Cldn3 + ) and aaTEC2 (Epc am − MHCII + PDGFRa − Pdpn + ) subsets that increased with age (Fig. 1j,k and Extended Data Fig. 4f).Overall, these data reveal that the most prominent shift among stromal cell types during thymic involution was the emergence of two unique epithelial cell types that lack typical TEC features.transition (EMT).The accumulation of aaTECs in the involuted thymus was exacerbated by acute injury and was associated with diminished regenerative capacity compared to young mice.We found evidence that aaTECs drew tonic signals away from other TECs, acting as a 'sink' for epithelial regeneration cues such as FGF and BMP signaling.These structural and functional changes to the thymic epithelium could be linked to molecular changes in the fibroblast compartment and specifically, their age-related upregulation of programs associated with inflammaging.These data define a key feature of the involuted thymus that limits organ function and restricts regenerative capacity after acute injury; findings that may have important implications for the development of therapeutic strategies to improve thymic function in aged individuals.

Emergence of atypical epithelial cell populations with age
Consistent with its well-described decline in function over lifespan 2,3,21,22 , we found that total thymic cellularity declined in female mice across 18 months of age (Extended Data Fig. 1a) coincident with morphological changes, such as a relative decrease in the cortical-to-medullary ratio (Extended Data Fig. 1b,c).Quantification of the major structural cell lineages (TECs, ECs and FBs) by flow cytometry revealed little alteration in ECs or FBs, but a diminished TEC compartment that mirrors the overall loss of thymic cellularity, with a more severe loss in medullary TECs (mTECs) compared to cortical TECs (cTECs) despite the observed cortical thinning (Extended Data Fig. 1b-d).
Primary stromal cell lineages were defined based on transcription of lineage-specific genes: TECs with Epcam, H2-Aa; ECs with Pecam1, Cdh5; FBs with Pdgfra, alongside less abundant stromal cell types, including mesothelial cells (MECs) (Upk3b and Nkain4); vascular smooth muscle cells (vSMCs) (Acta2 and Myl9); pericytes (PCs) (Myl9 and Acta2) (Fig. 1a and Extended Data Fig. 2a); and extremely rare nonmyelinating Schwann cells (nmSCs) (Gfap, Ngfr (p75) and S100b) 25 .To more precisely define the heterogeneity of the major CD45 − structural compartments (epithelial, endothelial and fibroblast), we then subsampled and reanalyzed each major stromal cell population separately and used public subset signatures (Supplementary Table 1) provided collectively within ThymoSight to precisely define compartment heterogeneity.Using this approach, we annotated steady-state clusters within each cell lineage with respect to publicly available datasets (Fig. 1c,d and Extended Data Fig. 2b-d).Unsupervised clustering analysis distinguished the Pdgfra-expressing mesenchyme into three main groups, two of which were consistent with mouse capsular or human interlobular FB signatures (capsFB; including but not limited to expression of Dpp4, Smpd3 and Pi16) and mouse medullary or human perilobular FBs (medFB; Ptn and Postn) 9,23,26 (Fig. 1c,d, Extended Data Fig. 2b and Supplementary Table 1).We also identified an intermediary subset of FB (intFB) marked by Inmt and Gpx3, which did not map to the public gene signatures (Fig. 1c,d and Extended Data Fig. 2b).

aaTECs form high-density clusters of epithelial cells
To study epithelial structures in the involuted thymus, we created a reporter strain, Foxn1 Cre × Rosa26 nTnG (Foxn1 nTnG ), where all cells express a nuclear-localized tdTomato reporter, except those that have activated Foxn1 (a master transcription factor in TEC 33 ), which instead, express a nuclear-localized green fluorescent protein (GFP) (Extended Data Fig. 5a).The nuclear localization of the Foxn1 nTnG reporter enabled    cellular resolution of TEC by light-sheet imaging of whole cleared thymic lobes from young and old mice.As expected, we observed a relatively low density of GFP + cells in the thymic cortex compared to the subcapsular and medullary regions (Fig. 2a).At the whole-tissue level, the medulla formed a highly complex and interconnected structure in the young thymus that degenerated into isolated islets upon involution (Fig. 2b and Supplementary Video 1a,b).Another notable feature of the involuted thymus, entirely absent in the young, was the emergence of zones of very high-density GFP + TEC (HD-TEC) clusters (Fig. 2a,b).These HD-TEC clusters formed band-like structures that were associated with the medulla of the involuted thymus (Fig. 2b, Extended Data Fig. 5b,c and Supplementary Video 1b).Although the volume of cortex and medulla, and the number of cTECs and mTECs (calculated from whole-tissue imaging) declined with age, HD-TEC clusters emerged to comprise a substantial volume and number of TECs (Fig. 2c).Spatial transcriptomics using Visium comparing thymi from 2-mo to 18-mo mice demonstrated aaTEC1 and aaTEC2 signatures (generated from the 20 genes most differentially expressed within each subset   EpCAM + tdTom -GFP + TEC Lineage g, Vein plots describing the continuous transition of 18-mo early prog , mTEC1, mTEC prol and aaTEC subsets to their predicted descendants (represented by diagonal flows) and the dynamic relative frequencies (vein width on the y axis) of these TEC subsets in the thymus over the binned pseudotime.h, Expression of thymocyte markers Thy1 and Lck overlaid on the 18-mo spatial transcriptomics dataset.Outline represents thymocyte-poor area overlaid onto heatmap showing aaTEC1 or aaTEC2 signatures.i, Two representative images in 12-18-mo male and female Foxn1 nTnG mice showing tdTomato and GFP expression with HD-TEC areas highlighted, with few or no tdTomato + cells.Scale bar, 50 μm.j, Human tissue sections from a 50-year-old woman.Shown are consecutive sections with H&E, cytokeratin or CD1a staining.k, aaTEC1 and aaTEC2 gene signatures (top 20 marker genes from our mouse data converted to human orthologs; Supplementary Fig. 3 and Supplementary Table 3) were overlaid on human thymic epithelial cells (n TEC = 40,144) from single-cell sequencing datasets generated and published elsewhere 9,19,20 .Summary data represents mean ± s.e.m. and each dot represents an individual biological replicate.Statistics were generated (b-d) using a two-tailed Mann-Whitney test.

Article
https://doi.org/10.1038/s41590-024-01915-9compared to all other TECs) formed clusters that were only detected in the involuted thymus and the medullary distribution (as assessed by hematoxylin and eosin (H&E) staining) of these signatures resembled that of HD-TEC clusters (Fig. 2d and Supplementary Table 3).To determine whether the HD-TEC structures were composed of aaTECs, we surveyed the expression of keratin subunits distinguishing aaTEC from other TEC populations (Fig. 2e,f).We found that most GFP + TEC in these regions were K5 + , some were K8 + and none were K14 + (Fig. 2e,f).Notably, the prospective aaTEC1 marker claudin-3 was capable of clearly demarcating HD-TEC regions in the aged thymus (Fig. 2e,f).Taken together, these data strongly suggest that the DN-TECs identified through flow cytometry, HD-TECs identified by imaging approaches and aaTECs identified through single-cell RNA sequencing (scRNA-seq) are the same cell populations.

aaTECs are derived from FOXN1-expressing epithelial cells
Given the paucity of expression of TEC genes by aaTECs, yet their clear association with other TEC populations (Fig. 1a), we sought to establish the lineage relationship between aaTECs and TECs.Flow cytometric analysis of Foxn1 nTnG mice confirmed the emergence of two atypical GFP + TEC-derived populations with age: one population that was EpCAM + but lacked canonical cTEC or mTEC markers, mirroring the DN-TEC/ aaTEC1 population described above, and another that was EpCAM − , consistent with the aaTEC2 phenotype (Fig. 3a).Further analysis of these populations in Foxn1 nTnG mice confirmed that claudin-3 expression was increased in EpCAM + UEA1 − cells in 18-mo mice, while podoplanin was capable of marking EpCAM − aaTEC2 (Fig. 3b-d).These findings were supported by single-cell sequencing analysis of 20-mo nonhematopoietic thymic stroma from Foxn1 Cre × R26-fl-Stop-fl tdTomato (Foxn1 tdTom ) reporter mice, in which tdTomato is expressed in all cells with a history of Foxn1 expression.Integration of scRNA-seq data of 1,093 cells from 20-mo Foxn1 tdTom TECs with the broader dataset of 7,412 (18-mo and 20-mo) wild-type (WT) cells revealed that transcription of the tdTomato reporter was detected in all TEC clusters, including both aaTEC1 and aaTEC2 (Fig. 3e).Together, these orthogonal lineage-tracing approaches confirm the thymic epithelial origin of both aaTEC1 and aaTEC2 subsets.We next assessed the lineage relationships among TEC subsets by performing unbiased RNA velocity analysis.Consistent with previous reports, this approach demonstrated a clear lineage trajectory in young mice stemming from the mTEC prol and continuing into differentiated mTEC2 and mimetic cell lineages 5,7 (Fig. 3f,g, Extended Data Fig. 2d and Supplementary Video 2a,b).This relationship was preserved in the involuted thymus, with the additional progression of mTEC1 and early prog cells into aaTEC1, whereas aaTEC2 were derived from early prog and aaTEC1 (Fig. 3f,g, Extended Data Fig. 2d and Supplementary Video 2a,b) 8 .Although these data suggest that aaTECs are terminally differentiated subsets, given their expression of TEC precursor markers such as claudin-3 and Plet1 (Extended Data Fig. 4e), an alternative hypothesis is that aaTECs could instead represent a stalled progenitor cell differentiation stage.

aaTEC niches do not support T cell development
Closer inspection of our Visium sequencing data highlighted a marked lack of thymocyte transcripts in aaTEC zones (Fig. 3h).Consistent with this observation, HD-TEC clusters in the involuted thymus of Foxn1 nTnG mice excluded other tdTomato + cells (Fig. 3i, Extended Data Fig. 5b and Supplementary Video 3), which given its ubiquitous expression in all non-TECs, will largely mark thymocytes.These data suggest that aaTEC microenvironments were thymocyte 'deserts' that do not support thymocyte differentiation and are deprived of thymic crosstalk factors.Notably, staining of sections of aged human thymus also reveals the presence of epithelial-rich thymocyte-poor regions (Fig. 3j), consistent with previously published findings 34 .Finally, we found that aaTEC gene signatures were apparent in previously generated human TEC scRNA-seq datasets, but the utility of prospective markers such as claudin-3 and podoplanin for the identification of aaTECs on human tissue is still to be determined 9,19,20 (Fig. 3k, Extended Data Fig. 5d-f and Supplementary Table 3).

aaTEC form microenvironmental 'scars' associated with EMT
We next sought to understand the molecular alterations that may underlie thymic involution and how these relate to aaTECs.Despite the profound dysregulation that occurs during thymic involution, we found only relatively minor changes in the expression of key epithelial and thymocyte growth factors, largely restricted to fibroblasts (Extended Data Fig. 6a).cTECs are a crucial population involved with early thymocyte development, as well as regulating thymic size via expression of Foxn1 and its downstream targets such as Dll4 (refs.1,3,21,35).We did not observe any change in expression with age of Foxn1, Dll4 or other genes involved in cTEC identity and function, including Psmb11, Prss16, Enpep (which encodes Ly51) or Ly75 (which encodes CD205); although this is likely a technical limitation due to the low number of cTECs captured through single-cell sequencing of the whole CD45 − compartment.To identify concerted age-dependent transcriptional changes within each stromal population, we used gene set enrichment analysis (GSEA) 36,37 coupled with network enrichment analysis using Cytoscape 38 to integrate the GSEA results into networks sharing common gene sets.Pathways with significant enrichment across at least one population with age were broadly grouped into eight biological categories corresponding with the hallmarks of aging (Fig. 4a, Extended Data Fig. 6b and Supplementary Tables 4 and 5) [39][40][41] .Consistent with reports demonstrating a link between mitochondrial function, aging and senescence 42,43 , we found a broad decrease in the transcription of genes within pathways associated with mitochondrial function and metabolism across most cell populations and in fibroblasts upregulation of genes associated with immune function (including antigen processing and presentation).Together, these findings support the role of these cells and pathways in the induction of inflammaging or the senescence-associated secretory phenotype (SASP) 44,45 .Notably, changes in pathways associated with proteostasis were largely restricted to epithelial cell populations (Fig. 4a, Extended Data Fig. 6b and Supplementary Tables 4 and 5).  4 and 5) and Cytoscape network analysis was used to integrate enriched pathways (false discovery rate (FDR) ≤ 0.05) sharing a core set of genes.Dotplot of top five pathways within each category (Supplementary Table 5).b, GSEA pathway enrichment within aaTEC1 or aaTEC2 subsets (generated by comparing aaTEC1 and aaTEC2 to all other TECs; Supplementary Table 6).c, Heatmap of 8,795 genes within cTEC, mTEC1, aaTEC1, aaTEC2 and medFB subsets ranked by cadherin-1 (encoded by Cdh1) expression.d, Scatter-plot of Cdh1 and Vim with cTEC, mTEC1, aaTEC1, aaTEC2 and medFB subsets.e, Scatter-plot of Cdh1 and Vim transcription overlaid with expression of epithelial and mesenchymal genes and EMT known regulators.f, Expression of key epithelial genes and thymopoietic factors by various 18-mo TEC subsets, including aaTECs.g, Heatmap of AIRE-(left) and FEZF2-dependent/independent (right) genes reported previously 49 .Heatmap shows scaled normalized gene expression.h, CellChat interaction overview summarizing number of interactions between grouped populations in 2-mo and 18-mo mice.i, CellChat interaction analysis between stromal cell populations with early prog , mTEC1, aaTEC1 or aaTEC2 as cellular receivers (see also Extended Data Fig. 8b).Matrix represents all significantly enriched pathways targeting either early prog , mTEC1, aaTEC1 or aaTEC2 (color-coded by the receiver population) and split by the type of CellChat signaling (secreted, cell-cell and ECM).j, Levels of PTN and MK in thymus at 2-mo (n = 5), 12-mo (n = 5) or 18-mo (n = 5) female C57BL/6 mice.Summary data represents mean ± s.e.m. and each dot represents an individual biological replicate.Statistics for j were generated using the Kruskal-Wallis test with Dunn's correction.ECM, extracellular matrix.Because the approach above involved direct comparison of cell populations from young and aged mice, aaTECs (which are not present in significant numbers in young mice) could not be assayed.Therefore, we inspected enriched pathways of aaTECs versus all other TECs within the aged setting.GSEA revealed loss of antigen presentation within aaTEC1s (Fig. 4b and Supplementary Table 6), consistent with the progressive loss of thymic epithelial function upon differentiation into aaTECs.For aaTEC2s, one of the most highly enriched pathways was the hallmark EMT gene signature (Fig. 4b and Supplementary Table 6).Analysis of gene expression across 8,795 EMT-related genes in descending order according to expression of E-cadherin (Cdh1), as previously described 46 , further suggested that aaTEC2 lay in a liminal zone between epithelial and mesenchymal identity (Fig. 4c).Scatter-plots based on E-cadherin (Cdh1) and vimentin (Vim) transcriptional expression (as prototypical epithelial and mesenchymal markers, respectively) overlaid with expression of archetypal epithelial or mesenchymal genes, further suggested that aaTEC2 lost epithelial features and partially gained mesenchymal traits (Fig. 4d,e).These observations are consistent with a partial EMT (pEMT) state that parallels a senescence-associated pEMT tubular epithelial population identified in kidney fibrosis 47 , although we do not exclude the possibility that some aaTEC2s undergo a full EMT 48 .
The unusual morphology and microenvironment formed by aaTECs prompted us to further assess their function relative to typical TECs.Consistent with scRNA-seq analyses, imaging and flow cytometry confirmed that aaTECs did not represent prominent, known mimetic cells, lacking expression of markers of tuft cells, M cells or corneocytes (Extended Data Fig. 7a-c).We found that aaTECs expressed low levels of key TEC mediators of T cell differentiation, including Dll4, Psmb11, Kitl, Ccl25, Cxcl12 and Il7 (Fig. 4f).Neither aaTEC populations expressed Fezf2, Aire, key NF-κB target genes or the receptors Tnfrsf11a (which encodes RANK), Ltbr or Cd40, which are associated with driving mTEC differentiation (Fig. 4f and Extended Data Fig. 7d,e).Accordingly, we found little or no expression of AIRE-dependent or FEZF2-dependent tissue-restricted antigen (TRA) expression 49 in aaTEC1s or aaTEC2s, being largely restricted to mTEC2 and mTEC prol populations (Fig. 4g and Extended Data Fig. 8a).
Using CellChat 50 to identify cell-cell interactions between stromal populations, we found a clear shift with age where signals were redistributed with the emergence of the two aaTEC populations (Fig. 4h).Focusing on specific signaling pathways influencing aaTEC1s and aaTEC2s, as well as their putative upstream mTEC1s and early prog precursors (as well as their sources), we found enrichment for growth factors such as FGF, BMP and EGF, as well as factors associated with promoting EMT such as PTN (pleiotrophin), MDK (midkine) and ANGPTL (angiopoietin-like) coming from across the stromal compartment (Fig. 4i and Extended Data Fig. 8b) [51][52][53][54][55][56] .Consistent with this observation, the amounts of PTN and MDK increased in the thymus with age (Fig. 4j) and, along with many other TECs, aaTECs expressed abundant EMT-related integrins and Sdc4, which can act as a receptor for PTN, MDK and ANGPTL4 (ref.57) (Extended Data Fig. 8c).Given the link between EMT, fibrosis and senescence, and the well-described loss of TRA expression with age 5 , these data suggest that aaTECs form unique nonfunctional microenvironments in the involuted thymus associated with senescence and EMT 32,58,59 .Although not fibrotic themselves, these aberrant high-density aaTEC 'scars' mirror the mesenchymal scarring found in other tissues with age and may be responsible for diminishing overall thymic function 32,41 .

aaTECs limit thymic repair following acute injury in aged mice
The thymus is extremely sensitive to injury but has substantial capacity for repair.The ability of the thymus to regenerate is thought to decline with age yet the mechanisms of this deficit are poorly understood 2,3,60 .We found that 1-2-mo mice subjected to sublethal total body irradiation (TBI) had fully recovered thymic cellularity by day 28 after TBI (Fig. 5a).By contrast, aged mice exhibited a significant delay in the restoration of thymic cellularity and did not approach pre-damage levels until approximately day 42 (Fig. 5a).Plotting thymic size after damage relative to baseline cellularity showed that restoration of thymic cellularity during the regenerative phase was impaired in aged mice (Fig. 5b).Histomorphological analysis correlated with these findings, showing similar depletion of thymocytes, accumulation of cellular debris and granulation tissue formation observed 1 day after injury in both cohorts, but evidence of sustained fibrosis, dystrophic calcification and occasional dyskeratotic epithelial cells only in 18-mo mice, suggesting a relative dysfunction in the regenerative process with age (Extended Data Fig. 9a).To assess the stromal compartments that orchestrate thymic regeneration, we performed flow cytometric analysis of endothelial, mesenchymal and epithelial stromal cell populations before and after regeneration.Although there were few major differences in the early response of young versus aged mice in most populations (Extended Data Fig. 9b), one notable exception was the re-emergence of EpCAM + MHCII + Ly51 − UEA1 − DN aaTEC1s and EpC AM − MHCII + PDGFRα − PDPN + aaTEC2s, which are depleted but rapidly restored as a prominent feature of the regenerated thymus of aged mice (Fig. 5c,d and Extended Data Fig. 9b).Similarly, whole-organ imaging analysis of aged Foxn1 nTnG mice 28 days after TBI revealed prominent high-density aaTEC regions after damage (Fig. 5e and Supplementary Video 4).Quantification revealed that, despite the thymus remaining smaller than pre-damage, aaTECs comprised an increased volume and number, indicating a substantial preferential increase in aaTECs coincident with the impaired thymic regeneration of aged mice (Fig. 5e,f).
To explore comprehensively the mechanisms of the damage response and impaired regeneration in aged thymic stroma, we analyzed 58,309 CD45 − nonhematopoietic cells from the thymus of 2-mo and 18-mo by scRNA-seq at days 1, 4 and 7 after TBI and annotated them based on steady-state (day 0) subset signatures before integration (Fig. 5g-i, Extended Data Fig. 10a, Supplementary Fig. 1 and Supplementary Table 3).Age-associated TEC populations remained a major feature of the aged thymus after damage (Fig. 5j).RNA velocity analysis implied that, in 2-mo mice, the re-emergence of differentiated conventional mTEC populations (mTEC2, mimetic cells) stemmed largely from the mTEC prol populations (Fig. 5k).In contrast, after damage in 18-mo mice there was a skewing in the inferred direction and rather than mTEC prol cells driving differentiation toward conventional differentiated mTECs, mTEC1 and early prog cells drove differentiation toward aaTEC1 and aaTEC2 (Fig. 5k), consistent with findings by others 5 .Therefore, the emergence of aaTECs may perturb normal TEC differentiation in the early stages of regeneration.

aaTECs co-opt growth signals at baseline and during repair
After acute damage, we and others have demonstrated the importance of epithelial growth factors such as BMP4 and keratinocyte growth factor (KGF) for endogenous thymic repair 31,61 .Indeed, we found a broad upregulation of endogenous epithelial regenerative factors after damage such as Fgf7, Fgf10, Fgf21 and Bmp4, as well as thymopoietic factors such as Flt3l and Kitl (Fig. 6a and Supplementary Tables 7-10).In contrast, stroma from aged mice did not upregulate these regenerative programs to the same extent or, in many instances, at all (Fig. 6a and Supplementary Tables 7-10).cTECs are crucial as targets for epithelial regenerative cues via expression of proteins such as FOXN1 and its downstream targets such as Dll4, Ccl25 and Cxcl12 (ref.31).Consistent with this, Foxn1 expression increased in cTECs after damage, which was maintained in aged mice, which was also reflected in its downstream chemokine targets Ccl25 and Cxcl12 (Fig. 6a,b and Supplementary Tables 7-10).In contrast, increased expression of the key regeneration mediator Dll4 and Psmb11 (a gene encoding for β5t, a proteasome subunit critical for CD8 + T cell selection), were abrogated or decreased,  9b).An aging index was generated by calculating the ratio of aged to young AUC for each indicated population (n = 10 per cell type).e, Whole-tissue imaging of 12-mo male Foxn1 nTnG mice at baseline or 28 days after TBI (550 cGy).f, Total volume and volume of cortex, medulla and aaTEC regions, as well as the number of cTECs, mTECs and aaTECs of the thymic right lobe.g-j, scRNA-seq was performed on CD45 − cells isolated from 2-mo or 18-mo thymus at baseline (day 0) and days 1, 4 and 7 after TBI.UMAP of 81,241 CD45 − cells annotated by age cohort (g), day after TBI (h) or structural cell subset mapped from Fig. 1c and Extended Data Fig. 10a (i).j Associated frequency analysis of all TEC subsets after TBI within each age cohort.k, RNA velocity analysis on all TEC subsets at days 0, 1, 4 and 7 after TBI in 2-mo or 18-mo mice.respectively, in aged mice (Fig. 6b and Supplementary Tables 7-10).Receptors for epithelial growth factors were broadly expressed across TEC subsets, including both aaTEC1s and aaTEC2s (Fig. 6c); however, using CellChat 50 to infer interactions based on expression of receptorligand pairs, as well as downstream signaling effectors, demonstrated a clear skewing of interactions with age toward aaTECs for two prominent growth cues for TECs: FGF and BMP signaling (Fig. 6d).These data suggest that aaTECs draw these pro-growth factors away from conventional TECs with age.This was supported by active signaling through FGF receptors in response to KGF in cTECs, mTECs and aaTECs, but not fibroblasts, which should not respond to KGF (Fig. 6e).

Reduction in FOXN1 activity favors aaTEC differentiation
FOXN1 is the master regulator of thymic epithelial lineage, crucial for TEC differentiation and function.It is also important for thymic regeneration after damage, acting as a key target downstream of these regeneration pathways 2,3 .FOXN1 protein could be readily detected in cTEC and mTEC populations but not aaTEC1s or aaTEC2s (Fig. 7a).
Consistent with this, there was little transcription of Foxn1 by aaTEC subsets compared to putative precursor populations (cTEC, early prog and mTEC1; Figs.4f and 7b,c).Using Dynamo 62 to predict the differentiation potential (negative stem cell potential) of cells within the 18-mo TEC dataset, we found that Foxn1 expression correlated only with differentiation away (and not toward) aaTECs (Fig. 7d,e).
To reveal quantitative insights into the regulatory capacity of Foxn1, we first computed the RNA Jacobian using Foxn1 as both regulator and target gene (Fig. 7f).The response heatmap showed that self-activation of Foxn1 followed an almost linear trend with few intermediate plateaus, suggesting endogenous Foxn1 self-induction and differential regulation of Foxn1 targets based on expression level, consistent with previously published reports of the self-regulation of FOXN1 (ref.63).Next, we applied in silico perturbation analysis of Foxn1 to assess the impact of this major regulator on cell fate outcomes.We first simulated the impact of Foxn1 deactivation and found that it diverted differentiation away from mTEC prol and mTEC2 subsets and toward the aaTECs (Fig. 7g).In contrast, computationally activating Foxn1 expression reinforced the transition of cTEC, early prog and mTEC1     subsets toward the normal differentiation trajectory and away from the aaTEC compartment (Fig. 7g).To test these in silico perturbation predictions, we turned to the Foxn1 Z/Z mouse model expressing a hypomorphic allele that causes early loss of Foxn1 in TECs and premature thymic involution 64 .scRNA-seq of TECs from 6-mo Foxn1 Z/Z mice and controls was integrated and compared to our 2-mo and 18-mo datasets.
The TEC subset frequency in 6-mo Foxn1 WT/WT mice was comparable to 2-mo WT mice for the majority of populations, with the emergence of small but detectable populations of aaTEC1 and aaTEC2 subsets (Fig. 7h and Extended Data Fig. 10b); however, in 6-mo Foxn1 Z/Z mice we observed a large expansion of aaTEC1s and aaTEC2s collectively composing approximately 50% of TEC; an expansion even greater than 18-mo WT mice (Fig. 7h and Extended Data Fig. 10b).Collectively, these data suggest that loss of FOXN1 expression in TECs favors differentiation into an aaTEC fate.

Discussion
Here, we define age-associated changes to the thymic microenvironment in the involuting thymus that impairs function in two ways.First, atypical aaTECs form high-density epithelial clusters, devoid of thymocytes.The accretion of aaTEC regions directly contributes to the loss of functional thymic tissue with age and, given that these regions expand after damage, exacerbates injury-induced atrophy.Second, we found evidence that the emergence of aaTECs perturbs the network of growth factors supporting stromal cell function and thymocyte differentiation, likely constituting an additional impediment to thymic function.Notably, similar features of epithelial-rich, thymocyte-devoid regions can also be observed in the human thymus 34 and our data offer strategies for further interrogation of aaTECs, including the validation of claudin-3 expression by aaTEC1s and podoplanin by aaTEC2s.Despite the relatively early emergence of aaTECs, their genesis seems to be linked to hallmarks of aging.Genetic approaches demonstrated that aaTEC populations were derived from Foxn1 + precursors, yet both had lost expression of canonical markers of cTECs and mTECs.RNA velocity analysis suggest that both aaTEC1s and aaTEC2s derive from mTEC1s, consistent with its progenitor-like phenotype 5,7,65 , as well as a mTEC1-like cell that shares a signature with a recently described progenitor 8 ; however, more sophisticated lineage-tracing approaches will be required to more directly test the precursor-progeny relationships.Nevertheless, our data strongly suggest that loss of FOXN1 is a key driver for the emergence of aaTECs.FOXN1 is crucial for many aspects of TEC biology including their generation, maintenance and regeneration 33 .FOXN1 expression declines with age, a process that has been implicated in contributing to thymic involution 33,64 .Age-associated

Article
https://doi.org/10.1038/s41590-024-01915-9 TECs do not express FOXN1 nor its downstream targets.Moreover, in mice expressing a FOXN1 hypomorph 64 we found that accelerated emergence of aaTECs accompanies early thymic atrophy; however, although our velocity analysis suggests that aaTECs are the downstream products of differentiation from mTEC1s and early prog and may represent end-stage epithelial cells arising and accumulating in the involuted thymus, an alternative hypothesis could be that aaTECs are progenitor cells that have become blocked.In fact, expression of markers known to be expressed by TEC precursors such as claudin-3 and Plet1 could support this alternative interpretation and further studies will be needed to test these hypotheses [66][67][68] .
The molecular drivers of aaTECs likely also involve additional factors.For instance, we found considerable changes in the stromal microenvironment beyond TECs, including significant changes in fibroblasts that mirror many of the main hallmarks of aging 39,41 , reflected by a loss in mitochondrial, metabolic and proteostasis programs, and increase in pathways associated with inflammaging or the SASP 39 .This latter finding likely also reflects the broader role of inflammaging and SASP in driving EMT in aged tissues 58,59 .This link is especially notable given that one of the main features of the aging thymus is the replacement of functional tissue with fat 69,70 , which can directly drive loss of thymic function 71 .Furthermore, there is evidence that the emergence of fat in the human thymus may be triggered by EMT 72 .These data also support recent observations that age-associated changes in tissues are organ-and cell-lineage specific (for example senescence-associated inflammaging impacting on muscle regeneration) 32,39,40 .Our findings are therefore consistent with the concept that aaTECs represent a thymus-specific manifestation of these programs, and that age-related changes in other cells of the thymic microenvironment, in particular fibroblasts, could also contribute to aaTEC emergence; however, our data also suggest that the emergence of aaTECs impairs thymic function beyond the replacement of functional tissue, by competing with conventional TECs for growth signals.This effect seemed especially important after acute damage where expansion of aaTECs correlated with impairment to thymic regeneration, a significant clinical problem given the relationship between T cell reconstitution and clinical outcomes following HCT 4 .Notably, there is also evidence of epithelial-rich areas devoid of thymocytes in human thymus found only in aged and dysregulated tissue (such as myasthenia gravis) 34 ; however, without a direct way to deplete these cells (or block their emergence to begin with) further studies will be needed to elucidate the specific contributions of aaTECs in thymic aging and responses after damage.
These observations suggest that the defective response and recovery to acute damage with age could be due to: (1) failure of aged stromal subsets to orchestrate regenerative programs; (2) upregulation of a partial EMT program leading to expansion of aaTEC; and (3) the co-opting of regenerative factors by aaTEC from typical differentiated mTEC subsets.These features seem to be at least partially driven by changes within the fibroblast compartment, and particularly their upregulation of genes associated with inflammaging/SASP.In summary, these studies highlight unique stromal cell responses to age-and stress-related thymic atrophy.Furthermore, the discovery of aaTECs, along with the functional changes in fibroblasts with age consistent with inflammaging/SASP, provide therapeutic targets for improving T cell immunity more broadly.Age-associated TECs therefore constitute a nexus of stromal cell dysfunction in thymic involution and impaired regeneration.https://doi.org/10.1038/s41590-024-01915-9 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 Foxn1 lacz and Foxn1 tdTom 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 analysis
Mapping 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 Foxn1 tdTom 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.h5files were processed following the guidelines on the shunPykeR GitHub repository 78 , 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 clustering 79 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 scrublet 80 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).
Foxn1 LacZ data.The quality of the single cells was computationally assessed as described for the CD45 − TBI series (n hvgs =3,500, n comps = 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 (n hvgs = 3,500, n comps = 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 scanpy 82 with the Wilcoxon rank-sum method 83 .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 scanpy 81 (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 literature 23,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), Foxn1 tdTom and Foxn1 LacZ 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-7,10-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.
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.
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) 87uilt 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) 50was 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), mTEC prol /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 tool 88 .Low-quality reads and adaptor contaminants were removed https://doi.org/10.1038/s41590-024-01915-9using Trimmomatic 89 (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 STAR 90 aligner using default parameters and-runThreadN set to 32 to increase execution speed.The STAR-aligned files were then used as input to the feature-Counts 91 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 edgeR 92 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 log 10 -transformed raw bulk counts (per biological replicate) and log 10 -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.Extended Data Fig. 5 | Validation of aaTEC identification and imaging.a, Generation of Foxn1 nTnG mice .ROSA nT-nG (nT/nG) mice 73 were intercrossed with Foxn1 Cre mice 74 .Representative flow cytometric plots of TEC from 11 weeks old WT, nT/nG and Foxn1 nTnG show specific detection of GFP in nearly all TEC only in the latter strain.Quantification of the relative proportions of TEC expressing the reporters are shown in the bar graph on the right (n = 2 to 3 from 2 experiments).b, Representative confocal images of thymic sections from 12-mo Foxn1 nTnG mice with high-density TECs located in peri-medullary region.c, Representative confocal images of thymic sections from 12-mo Foxn1 nTnG mice stained with antipan-keratin, with high-density TEC regions highlighted.d, UMAPs integrating all human CD45 − non-hematopoietic cells from published datasets 9,19,20 (ThymoSight) annotated by age and dataset (n = 115,536).e, UMAPs integrating all human thymic epithelial cells from public datasets 9,19,20 (ThymoSight) annotated by age and dataset (n = 40,144).f, Signatures of our mouse epithelial cell subsets (Supplementary Table 3), including aaTEC, overlaid onto the integrated human TEC data derived from 9,19,20 .Extended Data Fig. 10 | Mapping of pre-existing thymic stromal sequencing datasets after acute damage.a, Broad structural cell subsets were annotated based on expression of canonical markers such as Pdgfra, Epcam, H2-aa, Pecam, and Cdh5 and steady-state subset signatures(top 20 marker genes).Leiden clustering and cell subset signatures derived from Supplementary Table 3 across endothelial, fibroblast and epithelial cells isolated at days 1, 4, and 7 after TBI.b, Stacked barplots comparing TEC subsets frequency in the 6-mo Foxn1 Z/Z mice and controls to our own TEC subsets in 2-mo and 18-mo wild-type mice at steady state.♀: female mice, ♂: male mice.≤ ≥ α μ μ

Fig. 1 |
Fig. 1 | Emergence of atypical epithelial populations with age.a, Uniform Manifold Approximation and Projection (UMAP) of 22,932 CD45 − thymic cells from 2-mo and 18-mo female C57BL/6 mice, annotated by cell type subset and outlined by cell compartment (epithelial; fibroblast; endothelial; MEC; vSMC/PC; nmSC).b, ThymoSight integration of public data for murine nonhematopoietic thymic stromal cells, including our own dataset (n = 297,988) annotated by publication source and outlined by cell type and compartment.c, Violin plots highlighting key genes marking individual subsets within individual structural compartments (fibroblast, endothelium and epithelium).d,e, UMAPs of individual structural compartments color-coded by cell type subset (d) and age cohort (e).n EC = 1,661; n FB = 13,240; n TEC = 6,175.f, Scaled change in frequency for each individual structural cell subset with age.g, Gating strategy and quantities for cell populations within the epithelial lineage (based on previous work 6 ) in 2-mo (n = 10) and 18-mo (n = 10) mice.First, based on a CD45 − EpCAM + parent gate, tuft cells were identified by expression of L1CAM, then all other TECs were assessed for expression of conventional TEC markers UEA1 and Ly51.Within the UEA1 hi Ly51 lo mTEC population CD104 + MHCII lo cells were identified as mTEC1.Cells that were deemed as non-mTEC1 were then fractionated based on MHCII and Ly6D.h, Concatenated flow cytometry plots and graphs highlighting the frequency of Ly51 − UEA1 − (DN-TECs) across lifespan (gated on CD45 − EpCAM + MHCII + cells).i, Violin plots of aaTEC1 and aaTEC2 novel markers.j,k, Flow cytometry plots (j) and quantities (k) for aaTEC1 and aaTEC2 populations in 2-mo (n = 15), 12-mo (n = 10) or 18-mo (n = 13) male and female C57BL/6 mice.Summary data represent mean ± s.e.m.; each dot represents an individual biological replicate.Statistics were generated using a two-tailed Mann-Whitney test comparing within individual subsets (g) or Kruskal-Wallis (k) test with Dunn's correction.

Fig. 2 |
Fig. 2 | Age-associated TECs form distinct high-density peri-medullary 'scars'.a, Representative images of thymus in 2-mo and 24-mo male Foxn1 nTnG mice with medulla marked by dotted line.HD-TEC regions are apparent in aged but not young mice.b, Images of whole-tissue light-sheet imaging showing central medulla (magenta) in 2-mo and 24-mo male Foxn1 nTnG mice with HD-TEC (cyan) only apparent in aged mice.Medulla surface was defined on the basis of the frequency of GFP + cells and tdTomato expression (with high tdTom correlating with medullary regions) and confirmed by confocal imaging of KRT14.c, Quantification of total thymus volume, volume of cortex, medulla and HD-TEC regions and the number of mTECs, cTECs and HD-TECs calculated from wholetissue and confocal imaging in 2-mo (n = 3) and 18-mo (n = 3) mice.Summary data represent the mean ± s.e.m.; statistics were derived from independent biological replicates (individual animals) using a two-tailed unpaired t-test.d, Visium spatial sequencing performed on 2-mo or 18-mo C57BL/6 thymus.Displayed are H&E sections, cortex and medulla identified by Leiden clustering and heatmaps of aaTEC1 and aaTEC2 signatures (top 20 differentially expressed genes for each subset versus other TECs; Extended Data Fig. 3 and Supplementary Table3) overlaid.e, Expression of pan-keratin and keratin subunits 8, 5, 14 and claudin-3 on HD-TEC regions of thymus from 12-18-mo male Foxn1 nTnG mice.Scale bar, 50 μm.f, Transcriptional expression of keratin subunits and claudin-3 in the epithelial scRNA-seq dataset.n TEC = 6,175.

Fig. 4 |
Fig.4| Age-associated TEC regions are non-functional and associated with EMT.a, GSEA pathway analysis was performed for each subset based on differentially expressed genes in 18-mo versus 2-mo mice (Supplementary Tables4 and 5) and Cytoscape network analysis was used to integrate enriched pathways (false discovery rate (FDR) ≤ 0.05) sharing a core set of genes.Dotplot of top five pathways within each category (Supplementary Table5).b, GSEA pathway enrichment within aaTEC1 or aaTEC2 subsets (generated by comparing aaTEC1 and aaTEC2 to all other TECs; Supplementary Table6).c, Heatmap of 8,795 genes within cTEC, mTEC1, aaTEC1, aaTEC2 and medFB subsets ranked by cadherin-1 (encoded by Cdh1) expression.d, Scatter-plot of Cdh1 and Vim with cTEC, mTEC1, aaTEC1, aaTEC2 and medFB subsets.e, Scatter-plot of Cdh1 and Vim transcription overlaid with expression of epithelial and mesenchymal genes and EMT known regulators.f, Expression of key epithelial genes and thymopoietic Article https://doi.org/10.1038/s41590-024-01915-9

Fig. 6 |
Fig. 6 | Muted transcriptional response to irradiation in the aged thymic stroma linked to aaTEC disruption of trophic signals.a, Differential expression after damage within stromal cell subsets of key thymopoietic and epithelial growth factors in 2-mo or 18-mo female C57BL/6 mice at days 4 and 7 after TBI (compared to day 1).b, Differential expression of cTECs associated genes at days 4 and 7 after TBI in 2-mo and 18-mo mice.c, Expression of receptors for epithelial growth factors across stromal cell subsets.d, Chord diagraminteraction analysis of FGF and BMP signaling pathways in 2-mo or 18-mo mice at baseline and at day 7 after TBI.e, CD45 − nonhematopoietic stromal cells were isolated from 18-mo female mice and incubated for 5 min with KGF (100 ng ml −1 ) when cells were stained for phosphorylated AKT by flow cytometry (n = 6 thymuses isolated from individual mice).Summary data represent mean ± s.e.m. and each dot represents an individual biological replicate.Statistics were generated using a two-tailed paired t-test.

Fig. 7 |
Fig. 7 | FOXN1 loss of function accelerates aaTEC emergence with age.a, FOXN1 expression across TEC subsets by flow cytometry.b-d, UMAP of 18-mo TEC data at steady state, color-coded for TEC subset (b), Foxn1 expression (c) and differentiation potential (d) with projected RNA velocity data.e, Scatter-plot of differentiation potential (d) versus Foxn1 expression levels color-coded by subset annotation.f, Response heatmap showing the RNA Jacobian element for Foxn1 self-induction (df Foxn1 /dX Foxn1 ) versus the Foxn1 expression (Foxn1 (M s )).g, In silico perturbation analysis of Foxn1 in the 18-mo epithelium at steady state and accompanied cell fate diversions.Velocity arrows in UMAPs show cell fate directionality in the unperturbed dataset (left) and after in silico suppression (middle) or induction (right) of Foxn1.h, scRNA-seq was performed on CD45 − cells from 6-mo Foxn1 Z/Z mice and controls.UMAPs of 3,594 cells color-coded by sample and mapped to our TEC subsets.

Extended Data Fig. 1 |
General stromal features of thymic aging.a, Thymic cellularity of female C57BL/6 mice at 2, 6, 9, 12 or 19+ months of age.b-c, Representative images of hematoxylin and eosin (H&E) stained mouse thymi from 2, 6, 8-9, 13, and 19+mo female C57BL/6 mice (b) used to calculate ratio of cortical (dark) to medullary (light) region (c).In (c), each dot represents a biological replicate.d, Flow cytometric analysis of enzymatically digested thymus and absolute cell numbers for major cell types (TECs; cTECs and mTECs; ECs and FBs) in 2, 6, 9, 12, and 18+ mo female C57BL/6 mice.Summary data represents mean ± SEM; each dot represents an individual biological replicate; statistics were generated for a, c, and d using one-way ANOVA with the Dunnett correction for multiple comparisons.Extended Data Fig. 2 | Mapping of pre-existing thymic stromal sequencing datasets.a, Broad structural cell subsets were annotated based on expression of canonical markers such as Pdgfra, Epcam, H2-aa, Pecam, and Cdh5.b, Leiden clustering of our fibroblast population (n FB =13,240) and signatures for murine capsular-medullary and human perilobular-interlobular fibroblasts based on previously published datasets 9,26 .c, Leiden clustering of our endothelial population (n EC =1,661) and signatures for arterial, capillary, venular, and lymphatic endothelial cells based on previously published datasets 27 .d, Leiden clustering of our thymic epithelial population (n TEC =6,175) and signatures of previously published literature and overlaid on our sequencing dataset 5-7,13 .Extended Data Fig. 3 | Thymosight: Integration of thymic sequencing datasets.a-b, UMAPs of (a) all mouse non-hematopoietic thymic stroma cells (ThymoSight integration of public data 5-18 and ours; n = 297,988) annotated by age and sequenced population, and (b) all mouse thymic epithelial cells (n = 205,625; subset of CD45 − ThymoSight data) annotated by publication source, sequenced population, and TEC subset.Extended Data Fig. 4 | Quantification of non-epithelial stromal cells subsets and aaTEC gating.a-c, Concatenated flow cytometry plots and quantities for cell populations within the fibroblast (a), endothelial (b), and "other" (c) cell lineages in 2-mo (n = 10) and 18-mo (n = 10) mice.c, Concatenated flow cytometry plots and quantities for pericytes (PC), vascular smooth muscle cells (vSMC) and mesothelial cells (MEC) (n = 10/age).d, Frequency and numbers of DN-TEC across lifespan: 2-mo (n = 14), 6-mo (n = 5), 9mo (n = 15), 12-mo (n = 5), and 18+mo (n = 18).e, Violin plots with extensive list of aaTEC1 and aaTEC2 markers.f, Gating strategy for aaTECs.aaTEC1 were first gated on CD45 − TER119 − then PDGFRα - CD31 − cells.EpCAM + MHCII + cells were gated as the whole TEC compartment, then mTECs and cTECs were excluded by taking the UEA1 − Ly51 − double negative fraction and gating on CLDN3.aaTEC2 were also first gated on CD45 − TER119 − then PDGFRα -CD31 − cells.EpCAM − MHCII + cells were then gated and PDPN + PDGFRβ - were classed at aaTEC2.Summary data represents mean ± SEM; each dot represents an individual biological replicate.Statistics for a-c were generated using two-tailed Mann-Whitney tests comparing within individual populations and for d using the Kruskal-Wallis test with Dunns correction.

Fig. 6 |
Age related changes in gene expression in thymic stromal cells.a, Differential expression of key epithelial genes and thymopoietic factors with age.b, As in Fig.4a, GSEA pathway analysis was performed for each subset based on differentially expressed genes within each population between 2-mo and 18-mo mice (Supplementary Table

4
-5) and Cytoscape network analysis was used to integrate enriched pathways (FDR≤0.05)sharing a core set of genes.Dotplot of top 5 pathways within each category.Individual pathways are listed.Extended Data Fig.7| Comparison of aaTECs with conventional TECs and mimetic cells.a, 3D reconstruction and representative images of highdensity TEC region from 12-mo Foxn1 nTnG mice stained with DCLK1 or UEA1 to highlight tuft cells and M-like cells, respectively.b, Flow cytometry plots showing proportion of selected mimetic cells (tuft, corneocyte and M-cells) in 2-mo (n = 6) and 18-mo (n = 8) Foxn1 nTnG mice.Mimetic cells were first gated on EpCAM + GFP + cells, then mTECs (UEA1 hi Ly51 lo ) were assessed for the mimetic cell markers DCLK1 (tuft cells), GP2 (microfold cells), and Ly6D (corneocytes).Bar graph shows quantification of mimetic cell numbers.c, Flow cytometry plots showing mimetic cell frequency in 18-mo Foxn1 nTnG mice (n = 8) gated on EpCAM + GFP + UEA1 − Ly51 − DN-TECs.Bar graph shows quantification of mimetic cells comparing mTECs (as in Extended Data Fig.7b) and DN-TECs.d, Aire and Foxn1 expression in TEC subsets.e, Representative confocal images of thymic sections from 12-mo Foxn1 nTnG mice stained with anti-AIRE, with the medulla or high-density TECs highlighted.Scale bar: 50μm.Summary data represents mean ± SEM; each dot represents an individual biological replicate.Statistics for b-c were generated using two-tailed Mann-Whitney tests comparing within individual mimetic cell subsets.Extended Data Fig.8| aaTEC function and interactome with age.a, As in Fig.4g, heatmap of AIRE-(left) and FEZF2-(right) dependent/independent genes from72 .Heatmap shows scaled normalized gene expression.Individual genes are listed.b, As in Fig.4i, CellChat chord diagrams showing outgoing signals from all stromal cell populations towards early prog , mTEC1, aaTEC1 or aaTEC2 cellular receivers.Chord diagrams are color-coded by the sender population and split by the type of CellChat signaling (secreted, cell-cell and ECM).Specific outgoing signals per sender are listed in color-matching boxes on the side of each plot.c, Violin plot of receptors for putative EMT factors Midkine and Pleiotrophin.Extended Data Fig. 9 | Acute thymic recovery following ionizing radiation.a, Morphologic alterations in the context of acute thymic involution after TBI and thymic reconstitution in 2-mo and 18-mo mice.Top and bottom rows represent low-and high-power images from each timepoint, respectively.Annotations: cortex (*), medulla (**), adipocytes (#), areas of dystrophic calcification ( †), areas of dense fibrosis (arrowhead).b, Kinetics of recovery for the epithelial, endothelial and fibroblast defined subsets on day 0 (n = 10, 2-mo and 10, 18-mo), 1 (n = 10, 2-mo and 10, 18-mo), 4 (n = 10, 2-mo and 10, 18-mo) and 7 (n = 10, 2-mo and 10, 18-mo) after TBI in 2-mo and 18-mo mice.Total cellularity for each subset.Statistics compare across ages for each timepoint.Asterisks denote when recovery in either age cohort on a specific timepoint was significantly different to the other cohort.Summary data represents mean ± SEM.Statistics for b were generated using a two-way ANOVA with Šídák correction.