
95% of researchers rate our articles as excellent or good
Learn more about the work of our research integrity team to safeguard the quality of each article we publish.
Find out more
ORIGINAL RESEARCH article
Front. Immunol. , 03 March 2025
Sec. T Cell Biology
Volume 16 - 2025 | https://doi.org/10.3389/fimmu.2025.1553535
This article is part of the Research Topic Thymus Research and Development: A New Look to the Past, Current Knowledge, and Future Perspectives View all 9 articles
Thymus-committed regulatory T cells (Tregs) are essential for immune homeostasis. Recent findings stress their heterogeneity, suggesting possible alternate routes for thymic Treg development with unique features in humans, namely the clear evidence of Treg commitment at the double-positive (DP) stage and the presence of a significant population of CD8 single-positive (SP) FOXP3pos Tregs. Here, we present a dedicated analysis strategy to a spectral flow cytometry-based study of thymus from children and aged adults (≥ 74-years-old), to further elucidate Treg development and heterogeneity in the human thymus. We applied an unsupervised analysis pipeline to data generated from 6 high-dimensional panels, taking advantage of a common backbone of 11 markers, and we were able to map thymocytes along T cell maturation stages. Generating UMAP and FlowSOM cluster coordinates from the backbone, we projected all other markers onto these, characterizing clusters with the information of all markers. Focusing this analysis on events inside a putative total Treg gate, we could portray rarer subsets of human thymic Tregs and investigate their trajectories using pseudotime analysis. We uncover clusters within human DP thymocytes uniquely expressing FOXP3 or CD25, a DP-branching trajectory towards a CD103posCD8SP Tregs endpoint, and define trajectories towards CD4SP Tregs, including towards a cluster of CXCR3posCD4SP Tregs, that may consist of thymic resident or recirculating Tregs, and do not expand in the elderly. Our flow cytometry approach separates Treg populations with likely distinct functions and facilitates the design of future studies to unravel the complexity of human regulatory T cells.
Regulatory T cells (Tregs) expressing FOXP3 are essential players in T cell homeostasis, regulating autoimmune and other inflammatory responses (1, 2). Whilst Tregs may differentiate from conventional (Tconv) counterparts, the main functions in controlling autoimmunity are accomplished by thymic-generated Tregs (3, 4).
The thymic Treg heterogeneity has become a focus of recent research, driven by the development of novel animal models and advancements in single-cell transcriptomic methodologies (5–9). The uncovering of alternate pathways in Treg agonist selection, characterized by a different order in the acquisition of the Treg hallmarks CD25, the α-chain of the high affinity IL-2 receptor, and FOXP3, the lineage-associated transcription factor (5, 9), led to an update in the two-step TCR affinity-based concepts of their development (10). Current models consider that milder agonist signals drive the generation of a FOXP3-first precursor on the first step, and stronger signals generate a CD25-first precursor. The second step is mediated by γc-cytokine signals, namely IL-2 and IL-15, but a role for IL-4 and Lymphotoxin has recently been described in the regulation of FOXP3-first precursor maturation (5, 11–13). Tregs generated via each of these routes express different TCR repertoires and thus participate in the regulation of different responses (5, 12), underscoring the relevance of elucidating the signals driving their differentiation. Adding to the heterogeneity of Treg in the thymus, recirculating Tregs have also been described (14). Although FOXP3 reporter mouse models have been instrumental to shed light on many of these processes, translating these findings to human biology remains a significant challenge.
Two important considerations need to be taken into account in mouse vs. human Treg development (5, 12, 15). First, FOXP3 expression is found at much higher frequencies in human double-positive (DP) thymocytes than in mice (6, 15–17), hence human Tregs may mainly commit at DP (18, 19) rather than at the single positive (SP) stage. Second, CD8posFOXP3pos Tregs, whose relevance is gaining attention (20), are found at much higher frequencies in the human thymus compared to mice (15, 18, 21). The intricacies in human thymic Treg development and heterogeneity are far from being elucidated (22).
The advent of spectral flow cytometry enabled the use of more complex panels with overlapping spectra, resulting in an increase in the number of markers analyzed simultaneously (23, 24). When combined with emerging technologies such as CITE-Seq, single-cell transcriptomics, and spatial transcriptomics, these advances drive to a deeper understanding of thymic biology (6–8, 25, 26). However, many of these methodologies are resource-intensive, and spectral flow cytometry is likely to remain the primary approach for a substantial proportion of studies due to its accessibility and versatility. To handle the high-dimensional data generated, unsupervised clustering algorithms are being increasingly used to identify potential subpopulations based on the simultaneous expression of multiple markers. Furthermore, trajectory inference methods are providing valuable insights into cellular differentiation pathways (27, 28).
Here, we applied a flow cytometry-based approach to investigate human thymic Treg heterogeneity, using a common backbone strategy and adapting existing R-based algorithms. We integrated data from 6 different spectral flow cytometry panels to analyze thymic Tregs expressing FOXP3 and/or CD25, allowing us to define subpopulations and investigate lineage trajectories. In addition, we compared the relative distributions of the identified subsets in the thymus from children in parallel with older individuals. Altogether, we describe a strategy for integrating flow cytometry data from multiple panels combined with trajectory inference, enabling us to help decipher Treg development in the human thymus.
Human thymic specimens, from newborns to 15-year-old children, were obtained from routine thymectomy performed during pediatric corrective cardiac surgery at Hospital de Santa Cruz - Unidade Local de Saúde de Lisboa Ocidental, following parental written informed consent. Thymic specimens from adults older than 74-years-old were obtained from patients undergoing cardiac surgery for coronary artery disease or severe aortic stenosis at Hospital de Santa Maria – Unidade Local de Saúde de Santa Maria, after patient written informed consent. Samples used in the manuscript are listed in Supplementary Table S1. The thymic samples would be otherwise discarded. The study was approved by the Ethical Boards of the Centro Académico de Medicina de Lisboa, and of the Hospital de Santa Cruz, Portugal.
Thymic samples were maintained at 4°C and processed within less than 24h upon collection. Total thymocytes were recovered through tissue dispersion using a 70μM filter (BD Biosciences) and a syringe plunger and subsequently isolated on a Ficoll-Paque™ Plus (Cytiva) density gradient, as previously described (18). Thymic tissue with blood was discarded to avoid cross-contamination with peripheral lymphocytes.
Thymocytes isolated from 8-month-old infants (n=3) were simultaneously stained for 6 different panels (panels A to F), sharing a common backbone (Figure 1A; Supplementary Table S2). A backbone was defined from markers sharing the same fluorochrome for the same antigen (clone). Thymocytes from older patients (≥ 74-year-old, n=3) were processed only for panel F, as the obtained cell numbers were insufficient to do more than one panel. Surface staining was performed on fresh thymocytes for 25 minutes at room temperature, and were fixed, permeabilized, and stained using the eBioscience™ Foxp3/Transcription Factor Staining Buffer Set (eBioscience), as per manufacturer’s instructions, for intracellular staining. Samples were acquired on a 3-laser Cytek™ Aurora (Cytek Biosciences, Fremont, California) spectral flow cytometer and analyzed as outlined in Figure 1A.
Figure 1. Backbone strategy for panel integration in flow cytometry analysis of human thymic Tregs. (A) The scheme depicts on the left the experimental design, flow cytometry panel markers, and gating strategy to select total putative Tregs, and on the right the unsupervised analysis pipeline for file processing and panel integration. (B) UMAP of human thymic Tregs showing scaled expression of each backbone marker and a density plot (bottom right). (C) UMAP of human thymic Tregs (n=311.224), colored by clusters (left) with the heatmap of scaled median expression of backbone markers from each cluster with gray bars showing the numbers of cells on each cluster along the respective proportion within total Treg thymocytes (right) (see also Supplementary Figures S1-S5). Created in BioRender. Biorender, I. (2025) https://BioRender.com/j67h095.
All FCS files were unmixed and manually compensated in SpectroFlo software (v3.2.1; Cytek Biosciences, Fremont, California) or FlowJo (v10.8.1; BD Biosciences, Franklin Lakes, New Jersey), followed by the gating strategy outlined in Supplementary Figure S1 using FlowJo (v10.8.1). Thymocytes were gated based on forward and side scatter, doublets were excluded by plotting the area and height of forward and side scatter, live cells were gated by live/dead staining, and non-thymocytes were excluded by lineage markers (CD11c, CD14, CD16, CD19 and CD123). TCRγδ T cells were also excluded, included in the lineage channel in panels A-E, or on a separate channel in panel F. For Treg cells, a further gate based on FOXP3 and CD25 expression was applied. Gated samples were exported as FCS files and imported into R (v4.3.3; R Foundation for Statistical Computing, Vienna, Austria) as flowSets (flowCore v2.14.2). Arcsinh transformation was applied to all samples (FlowVS v1.34.0). After transformation, automatic quality control based on flow rate, signal acquisition and dynamic range was performed (flowAI v1.32.0). To correct for technical intersample variance, normalization was applied using fdaNorm (FlowStats v4.14.1). Upon pre-processing, samples were exported as FCS files (flowCore v2.14.2).
Pre-processed FCS files (samples from panels A to F) were uploaded to R and downsampled to the minimum number of events across all samples (322.259 for total thymocytes, 17.442 for Treg and 248.748 for Tregs in the comparison of young and older human thymi). Downsampling was performed randomly and the resulting flowSets were exported as FCS files. Files were reimported as flowSets, one for each panel, and subsetted to backbone markers only. FlowSets were merged, and FCS files for each sample only with backbone markers were exported. Batch correction was performed for technical panel variance. We compared normalization by warp and gauss with batch removal through Harmony, as implemented by the FlowCT (v1.0.0) package. Harmony was selected as the best overall method.
For UMAP (uwot v0.1.16), backbone marker expression levels (batch corrected with Harmony) were used (Supplementary Figures S2, S3). The number of neighbors and minimum distance were tested and defined. Visualization was achieved by plotting the marker expression values scaled from 0 to 1 on the UMAP (CATALYST v1.26.1 and ggplot2 v3.5.1). Cells were clustered using FlowSOM (CATALYST v1.26.1) and visualized in a cluster tree (clustree v0.5.1). The number of clusters was decided based on tree stabilization, and phenotypically similar clusters were manually merged to avoid over-clustering (Supplementary Figure S4). The single-cell experiment was exported as an R Data Format (RDS). Clusters of cells lacking marker expression, which separated early in the tree, were excluded. The diffusion map was implemented as in the destiny package (v3.16.0).
Trajectory inference was performed with Slingshot (v2.10.0) on a diffusion map representation of the flow cytometry data, which captures continuous transitions between cell states. We used the clusters previously established using FlowSOM and defined cluster 14 as the starting cluster based on its high expression of CD1a, CD4 and CD8 and lower CD3 expression, and then Slingshot was used to fit smooth trajectories through the high-dimensional space. Final clusters (9, 13, 16, 1 and 5) were defined based on their high expression of CD45RA and/or their location at the extremities of the diffusion map. Finally, each cell is assigned a pseudotime value, indicating its relative position along the inferred trajectory. Cells at earlier stages have lower pseudotime values, while more differentiated cells have higher values. Visualization was achieved by plotting the inferred pseudotime values and scaled marker expression on the diffusion map (ggplot2). The single-cell experiment was exported as an RDS.
The visualization of the markers not included in the backbone was performed by importing the FCS files after event down-sampling from a particular panel containing all markers’ information and the RDS file (before cluster removal and diffusion map generation), followed by subsetting to the particular panel. UMAP coordinates and cluster codes were extracted and provided to the complete dataset. Events from previously removed clusters were again excluded. Diffusion map coordinates, lineages, and curves were extracted from the RDS and applied to the dataset. Marker visualization was conducted as previously described, both on the UMAP and the diffusion map.
Fresh thymocytes, isolated from 7-month-old to 1-year-old infants (n=3), were used to assess cytokine production after stimulation with phorbol myristate acetate (PMA; 50ng/mL, Sigma-Aldrich) plus ionomycin (500ng/mL; Sigma-Aldrich), in the presence of brefeldin A (10µg/mL; Sigma-Aldrich) for 4 hours, followed by surface staining, fixation, permeabilization, and intracellular staining, as previously described (29).
All software and algorithms used can be found in Supplementary Table S3, together with the relevant references.
Spectral flow cytometry allows the expansion of the number of markers in a single-tube. When combined with the potential of the new analysis algorithms it offers powerful tools for detailed, single-cell analysis, yet to be fully explored. We took profit of existing algorithms in the R platform and designed an approach to help decipher thymic Treg heterogeneity in the human thymus. In order to identify and characterize relevant stages from intermediate developmental stages to mature Treg subpopulations, we used a 3-laser spectral flow cytometer and an unsupervised analysis approach. We developed a panel-integration strategy, based on a common backbone approach and an R-studio pipeline (Figure 1A and methods section), adapting existing algorithms. Three thymuses from 8-month-old children undergoing cardiac reconstructive surgery were analyzed (Figure 1A). We selected events inside a total thymocyte gate defined by the expression of FOXP3 and/or CD25 to include all putative Treg-lineage cells (Figure 1A; Supplementary Figure S1). The common backbone of 11 markers, using the same antibody clones and fluorochromes, was designed to contain, besides the Treg-lineage, markers able to discriminate major T cell developmental stages (CD1a, CD3, CD45RO, CD45RA, CD69, CD31, CD4, and CD8), and CXCR3, a differentiation marker potentially associated with recirculating cells (Figure 1A; Supplementary Figure S2). Each panel (A-F) included other non-backbone markers potentially informative on thymic Treg subsets, allowing us to gain information on a total of 44 markers (Figure 1A). The analysis pipeline is also outlined in Figure 1A. The initial “raw”-FCS files were transformed, cleaned, and normalized using available R-packages (detailed in Methods section) before downsampling, upon which two different FCS files were created: one (green in Figure 1A) including only data from the backbone markers, used to generate UMAP coordinates and FlowSOM clustering, and a second one (purple in Figure 1A), including all marker information from each tube that can then be projected in the UMAP coordinates and the identified clusters.
The UMAP generated with information from backbone markers allowed us to evaluate the putative total Tregs across the main developmental stages in the human thymus (Figure 1B). UMAP regions became apparent, such as cells expressing just one of the Treg-lineage markers (FOXP3posCD25neg-low, CD25posFOXP3neg-low), and DP, CD4SP and CD8SP thymocytes. Also observed was a peninsula of CXCR3 expressing CD4SP Tregs (Figure 1B). The density plot allowed us to visualize the relative paucity of immature cells, with lower CD3 expression and high expression of CD1a, CD4 and CD8α, as well as of CD8SP Tregs (Figure 1B). These were not observable in the UMAP performed in parallel using total thymocytes (Supplementary Figure S2), highlighting the advantage of the focused putative Treg analysis. Importantly, the expression of backbone markers from each panel was virtually equal in all panels (Supplementary Figure S3), indicating the reproducibility of signal from the different panels and illustrating the successful panel-integration.
Next, we used the information of the backbone-marker-only FCS files to apply FlowSOM clustering (see Methods). The clustree visualization led us to eliminate two small clusters, completely separated in the tree from the onset, and merge other minor clusters (Supplementary Figure S4). All 16 clusters were present in each of the three individual samples with only slight frequency variations (Supplementary Figure S5). The 16 identified clusters were then projected onto the UMAP with the scaled expression of the backbone markers in each cluster depicted in the heatmap in Figure 1C. In addition to the classical Treg clusters, expressing high levels of CD25 and of FOXP3 (clusters 12, 13, 15, 16) several emerged that deserved exploration. Clusters 4 to 9 lacked FOXP3 expression or expressed it at very low levels, and thus was difficult to exclude the presence of Tconv. These included CD4SP (clusters 4-5), CD8SP (clusters 8-9), and DP thymocytes (clusters 6-7), the latter displaying the highest CD1a expression, in accordance with an immature state. Conversely, clusters 11-14 express very low to no CD25 and are either CD8SP (clusters 12-13) or DP (clusters 11 and 14). From the latter, cluster 14 stands out as expressing the lowest levels of CD3, highest of CD1a and low CD69, denoting the possibility of representing the most immature cells in the gate. In contrast, cluster 11 appears to be gaining expression of CD25. The CD4SP cluster 1 appeared separated in the UMAP, displaying marked CXCR3 expression along with CD45RO, low levels of CD45RA or CD31, but interestingly, some CD1a and CD8α, suggesting the inclusion of immature cells (Figures 1B, C). Overall, the unsupervised analysis of the backbone markers’ expression unveils important heterogeneity within putative total Tregs in the human thymus, raising questions concerning relationships between subpopulations and developmental trajectories.
We have then searched to use the information from all other markers to shed light on these aspects, taking advantage of the panel-integration strategy. Our methodology allowed us to display the expression of all other markers from the six panels projected onto the UMAP (Figure 2A; Supplementary Figure S6) and then associate their expression with the FlowSOM-derived clusters (Figure 2B). The DP cluster 14 expressed the highest levels of the proliferative marker Ki-67, the Wnt signaling-associated transcription factor TCF1, and the cortical-associated chemokine receptor CCR9, together with the lowest expressions of HLA-ABC, CD27, CD28 or CCR7, attesting the more immature nature of these cells (Figure 2; Supplementary Figure S6). Interestingly, this phenotype is partially shared by the DP cluster 6, which did not express FOXP3 (Figure 2; Supplementary Figure S6). However, this CD25-only cluster featured much lower levels of Ki-67, and some CD69 expression alongside the brightest HLA-DR (Figure 2; Supplementary Figure S6), possibly betraying a stronger TCR-signal strength perceived.
Figure 2. Characterization of total Tregs in the human thymus. (A) Scaled expression of additional markers not included in the backbone markers projected on the UMAP coordinates generated with backbone markers. (B). Heatmap of scaled median expression of the backbone and all other additional markers in each of the 16 clusters defined using FlowSOM (see also Supplementary Figure S6).
The other low-null FOXP3 clusters, CD4SP (clusters 4-5) and CD8SP (clusters 8-9) thymocytes, differed mostly in the symmetrical expression of CD45RO and CD45RA, suggesting that these represent consecutive thymocyte developmental states in each of the SP populations, with RUNX3 being more expressed in the latter (Figures 1, 2), as expected due to its association with the CD8 lineage (30). Importantly, they all expressed low-null CTLA-4, an important suppressor marker (31), and intermediate levels of the α-chain of the IL-7 receptor, CD127 (Figure 2), further suggesting that these clusters may include Tconv lineage cells. On the other hand, the clusters expressing FOXP3 and low-null levels of CD25 (clusters 10-13) were all CD103-expressing (Figure 2). These included DP and CD8SP thymocytes, supporting a role for CD103 in CD8 Treg development from the DP stage (15, 32).
FOXP3posCD25pos CD4SP (clusters 15-16) and CD8SP (clusters 12-13) populations represented a significant fraction of events, as expected (Figure 1C; Supplementary Figure S5), differing in the gradual changes in expression levels of different markers and on the higher levels of PD-1 and RUNX3 within the CD8SP clusters (Figure 2; Supplementary Figure S6). These subtle differences contrast with the dramatic immunophenotype found in the CXCR3-expressing cluster 1, displaying the highest levels of the Treg lineage markers CD25, FOXP3, and CTLA-4 in parallel with CD127, activation-associated markers such as CD95 and ICOS, as well as the adhesion molecule ICAM-1 (Figure 2), in the absence of significant expression of other chemokine receptors associated with T cell differentiation (Supplementary Figure S6). Of note, this cluster also expressed T-bet (Figure 2), linked to CXCR3 expression, but there was no evidence of IFN-γ production in a complementary stimulation assay (Supplementary Figure S7). Intriguingly, some of these cells also expressed, albeit at low levels, CD8 (α and β) and CD1a, suggesting the inclusion of some developing thymocytes in this cluster and thus that not all these cells are CD4SP (Figure 1, 2). Of remark, although the manual analysis recapitulates all these findings, as illustrated in Supplementary Figure S8, our panel-integration approach allows a deep Treg profiling in the human thymus, raising relevant questions about their developmental relationships.
Our results further support alternate paths of Treg development in the human thymus (5–7, 12), leading us to take advantage of our approach to explore available trajectory algorithms. We employed trajectory inference with pseudotime projection (see Methods) onto diffusion maps defined according to coordinates generated with the backbone marker information (Figure 3A). This allows the projection of all markers from panels A to F (Supplementary Figure S9), as well as the 16 identified clusters (Figure 3B). The high expression of CD1a and low expression of CD3 (Figure 3A), corresponding to cluster 14 (Figure 3B), presents as a tip of the diffusion map and was used as the point of origin for the pseudotime trajectory inference analysis (Figure 3C). The CD45RA-expressing clusters and/or clusters located at extremities (clusters 9, 13, 16 and 5) of the diffusion map were considered endpoints, plus an additional one corresponding to the CXCR3-expressing cluster 1, also at an extremity (Figure 3). Overall, we visualized five trajectories (Figure 3C) that could be further characterized by projecting markers from panels A to F (Figure 4).
Figure 3. Trajectory inference for total Tregs in the human thymus. (A) Projections of the scaled expression of the backbone markers in the diffusion map generated with backbone data from the merged panels; the distribution of events is also shown (bottom right). (B) Diffusion map projection generated with backbone marker expression, colored by FlowSOM clusters. (C) Pseudotime trajectory inference, showing the five identified trajectories (lines) (see also Supplementary Figure S9).
Figure 4. Analysis of inferred Treg trajectories in the human thymus. Projection of each inferred Treg trajectory labeled at the column top in the diffusion map generated with backbone data, showing events in each lineage colored by pseudotime values (first row) and the projection of the scaled expression of selected markers, including backbone and non-backbone markers in the subsequent rows.
Three trajectories corresponded to CD4SP endpoints, on the right side of the diffusion map and collecting a large part of events (Figure 3A, density plot, and Figure 3C), and to two CD8SP endpoints, on the left side of the diffusion map (Figures 3A, C).
Focusing on the CD8SP clusters, it becomes apparent that the determinant markers in these two divergent trajectories were the expression of FOXP3 and CD103 versus CD25 expression (Figure 4). In addition, CD8SP CD25-only lacked CTLA-4 and featured a less mature phenotype, with sustained expressions of CD1a and CD28, and higher levels of CD38 at the endpoint (Figure 4 and Supplementary Figure S9). This visualization and analysis further interrogate the Treg-lineage commitment of CD25-only CD8SP.
Concerning CD4SP endpoints, three different trajectories were identified, including a CD25-only (Figures 3C, 4). However, contrasting to CD8SP, a clear marker distinction was not found, which suggests a possible overlap of truly CD25-only with CD25pos CD4SP that later acquired FOXP3 expression. Finally, the CXCR3pos CD4SP trajectory showed more concordance with the one expressing FOXP3 than the CD25-only trajectory (Figure 4).
Thus, combining trajectory inference with our panel-integration approach further distinguishes different Treg populations that likely follow distinct developmental pathways in the human thymus.
Next, we investigated whether the subset distribution was altered with age, providing further clues for understanding human thymic Treg biology. For this, we analyzed thymic tissue from three representative age groups, namely collected from individuals aged 74 years or older (n=3) undergoing cardiac surgery, older children and male adolescents (5-15 years old, n=4), and children below 1-year-old (n=5) (Figure 5). As expected, the total number of cells recovered in the older thymi was low, imposing the selection of only one of the panels in Figure 1 (Panel F) to ensure the acquisition of a high number of putative Tregs (Supplementary Figure S10). We applied the same unsupervised analysis pipeline, minus the panel-integration process as only panel F was used. The UMAP generated with total thymocytes (Supplementary Figure S11A), and with total Tregs (Supplementary Figure S11B) allowed the visualization of markers and the identification of areas corresponding to the subsets and trajectories from the analysis using the 6 panels. FlowSOM clustering led us to identify 18 clusters (Supplementary Figure S12) that were projected onto the UMAP (Figure 5A). The clusters that likely correspond to the 5 trajectory endpoints from Figures 3 and 4 are manually annotated in Figure 5A. The analysis of the cluster distribution projected onto the UMAPs of the concatenated putative Tregs from individuals according to the three age classes showed relative stability during childhood and adolescence with marked alterations with aging (Figure 5B). This was even clearer in the individual analysis of cluster distribution provided in Figure 5C. Aging was associated with an apparent expansion of CD25-only cells with a unique phenotype (clusters 3, 4, and 5), rare in the younger samples. These cells were CD25low, lacking FOXP3 and CTLA-4, expressing CD45RO in the full absence of CD45RA, and high levels of CD127, CD69 (clusters 3-4) and HLA-DR (cluster 4), a phenotype compatible with locally activated or recirculating cells. Additionally, the most immature CD1a expressing clusters (clusters 1; 6-8; 12; 17-18) appeared reduced, even when considering the relative expansion of clusters 3-5. Also striking was the almost absence of CD8SP Tregs expressing FOXP3 with low levels of CD25 in the older thymi (clusters 9 and 11). The lack of CD103 in panel F precluded a detailed analysis of these CD8SP Tregs. Finally, the CXCR3pos CD4SP Tregs were divided in this analysis into two clusters (clusters 1 and 2), mostly differing in the expression of CD8α and CD8β, as well as of CD1a (Figure 5A). Interestingly, cluster 1, expressing these markers, was virtually absent in the older thymus samples, and the cluster 2 was not expanded. These findings favor a recent thymic origin for at least some of these CXCR3pos CD4SP cells and argue against the recirculating hypothesis.
Figure 5. Comparison of total Tregs from young and older human thymus. (A) UMAP was generated using total Tregs from the thymus of ≤ 1 year old children (n=5), older children of 5-15 years of age (n=4), and older adults ≥ 74 years old (n=3) analyzed with panel F (Figure 1), followed by FlowSOM clustering (n=241.574); the 18 considered clusters are projected onto the UMAP (left) with the heatmap displaying the scaled median expression and the gray bars depicting the number of events on each cluster with their respective proportion (right); manual annotation highlights the clusters corresponding to the endpoints identified in trajectory inference analysis. (B) Cluster projection on the UMAP of Tregs according to age groups of thymus donors. (C) Stacked bar graphs depict the distribution of cluster frequencies in each thymus with the age group identified at the top. (see also Supplementary Figures S11, S12).
In conclusion, all the identified subsets were conserved in the human thymus from infancy to old age. The study of older subjects represents an important strategy to better understand Treg biology, their heterogeneity and developmental paths in the human thymus.
The heterogeneity of the thymic Tregs and its biological implications have been increasingly recognized (4, 5, 7–9, 12, 28, 33), raising important questions concerning mechanisms involved and the need to uncover markers for their distinction. We designed an experimental approach to profile Tregs in the human thymus leveraging the potential of spectral flow cytometry and the advances in algorithms for unsupervised analysis. Focusing our analysis of putative total Tregs, including all FOXP3 and/or CD25 expressing thymocytes, we uncover Treg subsets associated with distinct developmental paths. Moreover, we show that these subsets were conserved during aging.
The spectral flow cytometry could be applied flexibly and using fewer resources than CyTOF and sequencing-based methods (23, 24, 34). This is particularly relevant for subpopulation sorting for mechanistic and fate-mapping studies, ultimately required for functional validation. Our backbone marker approach allowed leveraging information from 44 potentially co-expressed markers distributed into six distinct panels to help decipher Treg heterogeneity in the human thymus using a 3-laser flow cytometer. The common backbone markers, allowing the classification of thymocytes along known developmental stages and the total Treg identification, were used to generate a UMAP where the expression of the other markers was projected. Importantly, we have not attempted an “imputation” of values from markers not included in the same panel or tube (35) to characterize populations or subsets according to non-backbone marker expression. Instead, we generated backbone-based UMAP and FlowSOM cluster coordinates to project and infer cluster characterization using all markers across 6 panels, a strategy similar to one previously used with CyTOF data (34). While this methodology does not allow the attribution of an expression value for a marker not present in the panel (36, 37), it characterizes additional markers from the different panels using the UMAP coordinates generated with the information of backbone markers common to all panels. It also allowed us to investigate the expression of non-backbone markers across FlowSOM clusters and a pseudotime trajectory again defined by the backbone marker information. Notably, our strategy can likely be applied to more panels with more markers and in multiple biological contexts.
The analysis performed within the putative total Treg significantly increased the resolution, granting the identification of relevant clusters present in thymic tissue from donors of different ages. The inclusion of thymocytes only expressing CD25 or FOXP3 was decided based on our previous data on human Treg development (15, 38) and to fit with the FOXP3-first and CD25-first precursors described in mice (5, 12). However, we were aware of possible contamination by conventional thymocytes, as both markers can be upregulated following activation (39). Our study does not attempt to provide definitive answers into the relationships between cells classified in different clusters, but the marker information and trajectory inference analysis provide significant contributions. We understand the limitations in trajectory inference analysis, and possible biases introduced when choosing start and endpoints. However, our choices were founded on solid knowledge on thymic development and our analysis was intended to generate hypothesis to be further probed. Thymocytes expressing only one of these Treg-lineage markers were found along the T cell development, namely in DP clusters expressing higher levels of CD1a and in mature clusters defined by high levels of CD45RA in CD4SP and CD8SP. It is, therefore, tempting to consider that the CD1a expressing DP thymocytes correspond to the two alternate pathways previously suggested in line with an earlier Treg commitment in the human thymus (5, 6, 12, 13, 18, 40). However, the current algorithms did not retrieve pseudotime trajectories initiating in the two populations, a limitation expected to be solved soon given the rapid development of applications and the expanding number of markers in the panels. It is also important to contextualize our data to the recent studies using CITE-seq and single-cell transcriptomics of human thymocytes, which either do not detect or consider a DP Treg-committed population (25, 33), detect but do not discuss its implications (7) or detect and discuss, leaving open Treg commitment at more than one-time point across pseudotime (6).
Of particular interest in our analysis are CD8SP Treg clusters, as, possibly due to their rarity, these cells are not considered in trajectory inference studies with human thymocytes (6, 28, 33), and are found annotated only in the most recent human thymus atlas (8). We investigated two different CD8SP endpoints, clearly diverging halfway in the trajectory. One of the branches expressed CD25 and negligible levels of FOXP3 and CTLA-4 throughout the trajectory and includes a fraction of the CD25-only DP cluster (cluster 6) that failed to map to other trajectories. This branch is likely to include non-Treg cells. In contrast, the other CD8SP trajectory expressed FOXP3 and CD103. CD103, an integrin that dimerizes to integrin β7 and binds to E-cadherin (41), is already expressed at the DP stage, apparently marking CD8 Treg-lineage choice, as previously suggested (18). Thus, CD103 represents an relevant marker to sort and functionally investigate CD8 Tregs in the human thymus.
Interestingly, although we were able to identify mature (CD45RApos) CD4SP Tregs expressing either CD25 or FOXP3, it was not possible to identify a marker or marker combination clearly separating these trajectories, possibly reflecting greater fluidity in CD4 Treg commitment and the variable times of induction found for CD4SP Tregs (6).
Finally, the trajectory leading to CXCR3pos CD4SP presents significant overlap to both other lineages, however, it appears to represent a different fate, with a distinctive phenotype characterized by the highest expression of the bona fide Treg markers, such as FOXP3, CD25 and CTLA-4, as well as the highest expression of CXCR3, T-bet and ICAM-1 and lowest levels of TCF-1. This phenotype is in concordance with the phenotype of a population of cells described as recirculating in studies from both mice and humans (7, 8, 14, 33, 42). A recent study using a recently described novel trajectory inference methodology named tviblindi and different datasets, also identified a population with the same phenotype at the end of a trajectory (28). Nevertheless, it must be noted that these cells could also be thymic-resident and that the expression of CD45RO would also agree with an intermediate developmental stage. In line with this possibility is the observation of CD8α, CD8β, and CD1a expression in this cluster, albeit at low levels. Although CD8 reexpression has been associated with activated CD4 T cells in tissues (43), the expression of CD1a is intriguing and could be compatible with the coexistence within this population of cells at different developmental stages. Interestingly, in our comparative analysis of samples from children, adolescents, and older individuals, the CXCR3-expressing population subdivides into two clusters, differing in the expression of CD8 and CD1a. The cluster expressing CD8 and CD1a was almost absent in the samples from older individuals, as are clusters corresponding to most immature developmental stages, suggesting that it may represent immaturity rather than an activated state. In contrast to the lack of increase of CXCR3pos CD4SP Tregs in the older thymi, we do observe a marked expansion in clusters featuring an activated or possibly recirculating phenotype, namely expressing CD25, CD127 and negligible levels of FOXP3. These clusters may represent Tconvs and were rare in the thymic tissue of children and adolescents. We were careful during tissue processing to avoid any blood-derived contamination, but a possible relation of these cells with the adipose tissue observed in the thymic samples from older individuals cannot be excluded and remains to be investigated. The function of the CD25-only populations in the human thymus should be investigated in future studies. Our results support a function far beyond being precursors of bona fide FOXP3 Tregs. Additionally, we showed that the expression of CD25 not accompanied by FOXP3 was associated with distinct phenotypes in CD8SP and CD4SP. It would be important to define if the CD25-only populations have regulatory properties or if they play a role in the maturation of conventional thymocytes and may impose any functional imprint.
In conclusion, we investigated human thymic Treg heterogeneity using a spectral flow cytometry approach, that can be applied to other population-focused studies. The simplicity of the approach takes advantage of the many channels offered by spectral flow cytometry and of existing unsupervised analysis algorithms to visualize the expression of markers, projected into UMAPs and FlowSOM clusters defined by a backbone. Based on the acquisition of large numbers of thymocytes, we subcluster putative total Tregs and identify Treg subpopulations, likely to represent relevant alternate paths of development for CD4 and CD8 Tregs, in addition to possible thymic resident locally activated or recirculating populations. While the hypotheses generated by our analysis need fate-mapping studies with human thymocytes for validation, our flow cytometry data may provide handles for the design of these experiments.
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
The studies involving humans were approved by the Ethical Boards of Centro Académico de Medicina de Lisboa and Hospital Santa Cruz. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.
BM: Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing, Visualization. MP-P: Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing, Writing – original draft. NM: Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing. EJ: Investigation, Writing – review & editing. EC: Investigation, Writing – review & editing. TV: Investigation, Methodology, Writing – review & editing, Resources. MA: Investigation, Writing – review & editing, Methodology, Resources. RA: Investigation, Writing – review & editing, Methodology, Resources. AA: Investigation, Methodology, Writing – review & editing, Conceptualization, Formal analysis, Funding acquisition, Supervision, Writing – original draft, Validation. AS: Investigation, Writing – review & editing, Conceptualization, Funding acquisition, Supervision, Writing – original draft, Project administration, Validation.
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by national funding through FCT -Fundação para a Ciência e a Tecnologia, I.P., under the project PTDC/MED-IMU/0938/2020, DOI:10.54499/PTDC/MED-IMU/0938/2020 (Sousa, A.E., 2020). NM received fellowships funded by Janssen-Cilag Farmaceutica. MP-P (2023.02835.BD) and EC (2022.153366.BD) received PhD scholarships from FCT.
We would like to thank all the study participants and their legal representatives for their collaboration; the teams and directors of cardiothoracic surgical services from Hospital de Santa Cruz/ULS Lisboa Ocidental and Hospital de Santa Maria/ULS de Santa Maria for vital cooperation to obtain the thymus samples; the scientific advice of Diana Santos, Helena Nunes-Cabaço and Tomas Kalina.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.
The author(s) declare that no Generative AI was used in the creation of this manuscript.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1553535/full#supplementary-material
1. Sakaguchi S, Mikami N, Wing JB, Tanaka A, Ichiyama K, Ohkura N. Regulatory T cells and human disease. Annual Reviews. (2020) 38:541–66. doi: 10.1146/annurev-immunol-042718-041717
2. Dikiy S, Rudensky AY. Principles of regulatory T cell function. Immunity. (2023) 56:240–55. doi: 10.1016/j.immuni.2023.01.004
3. Sakaguchi S, Yamaguchi T, Nomura T, Ono M. Regulatory T cells and immune tolerance. Cell. (2008) 133:775–87. doi: 10.1016/j.cell.2008.05.009
4. Sjaastad LE, Owen DL, Tracy SI, Farrar MA. Phenotypic and functional diversity in regulatory T cells. Front Cell Dev Biol. (2021) 9:715901. doi: 10.3389/fcell.2021.715901
5. Owen DL, Mahmud SA, Sjaastad LE, Williams JB, Spanier JA, Simeonov DR, et al. Thymic regulatory T cells arise via two distinct developmental programs. Nat Immunol. (2019) 20:195–205. doi: 10.1038/s41590-018-0289-6
6. Heimli M, Flåm ST, Hjorthaug HS, Trinh D, Frisk M, Dumont K-A, et al. Multimodal human thymic profiling reveals trajectories and cellular milieu for T agonist selection. Front Immunol. (2023) 13:1092028. doi: 10.3389/fimmu.2022.1092028
7. Park J-E, Botting RA, Domínguez Conde C, Popescu D-M, Lavaert M, Kunz DJ, et al. A cell atlas of human thymic development defines T cell repertoire formation. Science. (2020) 367:eaay3224. doi: 10.1126/science.aay3224
8. Yayon N, Kedlian VR, Boehme L, Suo C, Wachter BT, Beuschel RT, et al. A spatial human thymus cell atlas mapped to a continuous tissue axis. Nature. (2024) 635:708–18. doi: 10.1038/s41586-024-07944-6
9. Herppich S, Toker A, Pietzsch B, Kitagawa Y, Ohkura N, Miyao T, et al. Dynamic imprinting of the treg cell-specific epigenetic signature in developing thymic regulatory T cells. Front Immunol. (2019) 10:2382. doi: 10.3389/fimmu.2019.02382
10. Lio C-WJ, Hsieh C-S. A two-step process for thymic regulatory T cell development. Immunity. (2008) 28:100–11. doi: 10.1016/j.immuni.2007.11.021
11. Owen DL, La Rue RS, Munro SA, Farrar MA. Tracking regulatory T cell development in the thymus using single-cell RNA sequencing/TCR sequencing. J Immunol. (2022) 209:1300–13. doi: 10.4049/jimmunol.2200089
12. Santamaria JC, Borelli A, Irla M. Regulatory T cell heterogeneity in the thymus: impact on their functional activities. Front Immunol. (2021) 12:643153. doi: 10.3389/fimmu.2021.643153
13. Borelli A, Santamaria JC, Zamit C, Apert C, Chevallier J, Pierre P, et al. Lymphotoxin limits Foxp3+ regulatory T cell development from Foxp3lo precursors via IL-4 signaling. Nat Commun. (2024) 15:6976. doi: 10.1038/s41467-024-51164-5
14. Thiault N, Darrigues J, Adoue V, Gros M, Binet B, Perals C, et al. Peripheral regulatory T lymphocytes recirculating to the thymus suppress the development of their precursors. Nat Immunol. (2015) 16:628–34. doi: 10.1038/ni.3150
15. Caramalho Í, Nunes-Cabaço H, Foxall RB, Sousa AE. Regulatory T-cell development in the human thymus. Front Immunol. (2015) 6:395. doi: 10.3389/fimmu.2015.00395
16. Cupedo T, Nagasawa M, Weijer K, Blom B, Spits H. Development and activation of regulatory T cells in the human fetus. Eur J Immunol. (2005) 35:383–90. doi: 10.1002/eji.200425763
17. Darrasse-Jèze G, Marodon G, Salomon BL, Catala M, Klatzmann D. Ontogeny of CD4+CD25+ regulatory/suppressor T cells in human fetuses. Blood. (2005) 105:4715–21. doi: 10.1182/blood-2004-10-4051
18. Nunes-Cabaço H, Caramalho Í, Sepúlveda N, Sousa AE. Differentiation of human thymic regulatory T cells at the double positive stage. Eur J Immunol. (2011) 41:3604–14. doi: 10.1002/eji.201141614
19. Vanhanen R, Leskinen K, Mattila IP, Saavalainen P, Arstila TP. Epigenetic and transcriptional analysis supports human regulatory T cell commitment at the CD4+CD8+ thymocyte stage. Cell Immunol. (2020) 347:104026. doi: 10.1016/j.cellimm.2019.104026
20. Liston A, Aloulou M. A fresh look at a neglected regulatory lineage: CD8+Foxp3+ Regulatory T cells. Immunol Lett. (2022) 247:22–6. doi: 10.1016/j.imlet.2022.05.004
21. Cosmi L, Liotta F, Lazzeri E, Francalanci M, Angeli R, Mazzinghi B, et al. Human CD8+CD25+ thymocytes share phenotypic and functional features with CD4+CD25+ regulatory thymocytes. Blood. (2003) 102:4107–14. doi: 10.1182/blood-2003-04-1320
22. Pala F, Notarangelo LD, Bosticardo M. Rediscovering the human thymus through cutting-edge technologies. J Exp Med. (2024) 221:e20230892. doi: 10.1084/jem.20230892
23. Nolan JP, Condello D. Spectral flow cytometry. CP Cytometry. (2013) 63:1.27.1–1.27.13. doi: 10.1002/0471142956.cy0127s63
24. Robinson JP. Flow cytometry: past and future. BioTechniques. (2022) 72:159–69. doi: 10.2144/btn-2022-0005
25. Chopp LB, Gopalan V, Ciucci T, Ruchinskas A, Rae Z, Lagarde M, et al. An integrated epigenomic and transcriptomic map of mouse and human αβ T cell development. Immunity. (2020) 53:1182–1201.e8. doi: 10.1016/j.immuni.2020.10.024
26. Cordes M, Canté-Barrett K, Van Den Akker EB, Moretti FA, Kiełbasa SM, Vloemans S, et al. Single-cell immune profiling reveals novel thymus-seeding populations, T cell commitment, and multi-lineage development in the human thymus. Science Immunology. (2022) 7(77):eade0182. doi: 10.1101/2022.02.18.481026
27. Saelens W, Cannoodt R, Todorov H, Saeys Y. A comparison of single-cell trajectory inference methods. Nat Biotechnol. (2019) 37:547–54. doi: 10.1038/s41587-019-0071-9
28. Stuchly J, Novak D, Brdickova N, Hadlova P, Iksi A, Kuzilkova D, et al. Deconstructing complexity: A computational topology approach to trajectory inference in the human thymus with tviblindi. eLife (2024) 13:RP95861. doi: 10.7554/eLife.95861.1
29. Nunes-Cabaço H, Ramalho-dos-Santos A, Pires AR, Martins LR, Barata JT, Sousa AE. Human CD4 T cells from thymus and cord blood are convertible into CD8 T cells by IL-4. Front Immunol. (2022) 13:834033. doi: 10.3389/fimmu.2022.834033
30. Taniuchi I, Osato M, Egawa T, Sunshine MJ, Bae S-C, Komori T, et al. Differential requirements for runx proteins in CD4 repression and epigenetic silencing during T lymphocyte development. Cell. (2002) 111:621–33. doi: 10.1016/S0092-8674(02)01111-X
31. Wing K, Onishi Y, Prieto-Martin P, Yamaguchi T, Miyara M, Fehervari Z, et al. CTLA-4 control over foxp3+ Regulatory T cell function. Science. (2008) 322:271–5. doi: 10.1126/science.1160062
32. Kutleša S, Wessels JT, Speiser A, Steiert I, Müller CA, Klein G. E-cadherin-mediated interactions of thymic epithelial cells with CD103+ thymocytes lead to enhanced thymocyte cell proliferation. J Cell Sci. (2002) 115:4505–15. doi: 10.1242/jcs.00142
33. Morgana F, Opstelten R, Slot MC, Scott AM, Van Lier RAW, Blom B, et al. Single-cell transcriptomics reveals discrete steps in regulatory T cell development in the human thymus. J Immunol. (2022) 208:384–95. doi: 10.4049/jimmunol.2100506
34. Bendall SC, Simonds EF, Qiu P, Amir ED, Krutzik PO, Finck R, et al. Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum. Science. (2011) 332:687–96. doi: 10.1126/science.1198704
35. Mocking TR, Duetz C, Van Kuijk BJ, Westers TM, Cloos J, Bachas C. Merging and imputation of flow cytometry data: A critical assessment. Cytometry Pt A. (2023) 103:818–29. doi: 10.1002/cyto.a.24774
36. Pedreira CE, Costa ES, Barrena S, Lecrevisse Q, Almeida J, Van Dongen JJM, et al. Generation of flow cytometry data files with a potentially infinite number of dimensions. Cytometry Pt A. (2008) 73A:834–46. doi: 10.1002/cyto.a.20608
37. Abdelaal T, Höllt T, Van Unen V, Lelieveldt BPF, Koning F, Reinders MJT, et al. CyTOFmerge: integrating mass cytometry data across multiple panels. Bioinformatics. (2019) 35:4063–71. doi: 10.1093/bioinformatics/btz180
38. Caramalho I, Nunes-Silva V, Pires AR, Mota C, Pinto AI, Nunes-Cabaço H, et al. Human regulatory T-cell development is dictated by Interleukin-2 and -15 expressed in a non-overlapping pattern in the thymus. J Autoimmun. (2015) 56:98–110. doi: 10.1016/j.jaut.2014.11.002
39. Kmieciak M, Gowda M, Graham L, Godder K, Bear HD, Marincola FM, et al. Human T cells express CD25 and Foxp3 upon activation and exhibit effector/memory phenotypes without any regulatory/suppressor function. J Transl Med. (2009) 7:89. doi: 10.1186/1479-5876-7-89
40. Nunes-Cabaço H, Ribot JC, Caramalho Í, Serra-Caetano A, Silva-Santos B, Sousa AE. Foxp3 induction in human and murine thymus precedes the CD4+ CD8+ stage but requires early T-cell receptor expression. Immunol Cell Biol. (2010) 88:523–8. doi: 10.1038/icb.2010.4
41. Hardenberg J-HB, Braun A, Schön MP. A yin and yang in epithelial immunology: the roles of the αE(CD103)β7 integrin in T cells. J Invest Dermatol. (2018) 138:23–31. doi: 10.1016/j.jid.2017.05.026
42. Kirberg J, Bosco N, Deloulme J-C, Ceredig R, Agenès F. Peripheral T lymphocytes recirculating back into the thymus can mediate thymocyte positive selection. J Immunol. (2008) 181:1207–14. doi: 10.4049/jimmunol.181.2.1207
Keywords: thymus, human immunology, regulatory T cells, T cell development, spectral flow cytometry
Citation: Moleirinho B, Paulo-Pedro M, Martins NC, Jelagat E, Conti E, Velho TR, Abecasis M, Anjos R, Almeida ARM and Sousa AE (2025) A backbone-based flow cytometry approach to decipher regulatory T cell trajectories in the human thymus. Front. Immunol. 16:1553535. doi: 10.3389/fimmu.2025.1553535
Received: 30 December 2024; Accepted: 11 February 2025;
Published: 03 March 2025.
Edited by:
Valentin P. Shichkin, OmniFarma LLC, UkraineReviewed by:
Robert Cody Sharp, University of Florida, United StatesCopyright © 2025 Moleirinho, Paulo-Pedro, Martins, Jelagat, Conti, Velho, Abecasis, Anjos, Almeida and Sousa. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Afonso R. M. Almeida, YWZvbnNvLmFsbWVpZGFAZ2ltbS5wdA==; YWZvbnNvLmFsbWVpZGFAbWVkaWNpbmEudWxpc2JvYS5wdA==
†These authors have contributed equally to this work and share first authorship
‡These authors share last authorship
Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.
Research integrity at Frontiers
Learn more about the work of our research integrity team to safeguard the quality of each article we publish.