Research ArticleT CELLS

Developmental bifurcation of human T follicular regulatory cells

See allHide authors and affiliations

Science Immunology  28 May 2021:
Vol. 6, Issue 59, eabd8411
DOI: 10.1126/sciimmunol.abd8411

Reconstructing TFR cell development

Antibody responses, including against self-antigens, are tightly controlled within germinal centers (GCs) by balancing the activity of T follicular regulatory (TFR) and helper (TFH) cells. Kumar et al. reconstructed the developmental trajectory of TFR cells using flow cytometry indexed single-cell transcriptomics of follicular T cells isolated from human peripheral blood, lymph nodes, and tonsils. After distinguishing TFR cells from other follicular T cell populations based on gene signature, they used transcriptomics to support a model in which mature GC TFR cells arise from regulatory T cells, with circulating TFR cells representing a separate developmental pathway. These results suggest that mature TFR cells primarily arise independently from their less mature counterpart in peripheral blood, providing further insight into how human TFR cells develop.


Germinal centers (GCs) are anatomic structures where B cells undergo affinity maturation, leading to production of high-affinity antibodies. The balance between T follicular helper (TFH) and regulatory (TFR) cells is critical for adequate control of GC responses. The study of human TFH and TFR cell development has been hampered because of the lack of in vitro assays reproducing in vivo biology, along with difficult access to healthy human lymphoid tissues. We used a single-cell transcriptomics approach to study the maturation of TFH and TFR cells isolated from human blood, iliac lymph nodes (LNs), and tonsils. As independent tissues have distinct proportions of follicular T cells in different maturation states, we leveraged the heterogeneity to reconstruct the maturation trajectory for human TFH and TFR cells. We found that the dominant maturation of TFR cells follows a bifurcated trajectory from precursor Treg cells, with one arm of the bifurcation leading to blood TFR cells and the other leading to the most mature GC TFR cells. Overall, our data provide a comprehensive resource for the transcriptomics of different follicular T cell populations and their dynamic relationship across different tissues.


Understanding of the collaborations between T and B cells leading to germinal center (GC) reactions was revolutionized by the discovery of follicular helper T (TFH) cells (13). TFH cells have unique access to the B cell follicles and are crucial for isotype switching and affinity maturation (4). In addition, identification of a dedicated population of Foxp3+ cells, termed T follicular regulatory (TFR) and shown to be endowed with suppressive function, provided further insight into the regulation of the GC response (58). The mechanisms leading to the development and maturation of TFH and TFR cells have been studied in great detail in mice (915). However, there are clear differences between the development of murine and human follicular cells. For example, murine TFH cells require interleukin-6 (IL-6) for their differentiation, whereas human TFH cells do not rely on IL-6; instead, they require IL-23, IL-12, and transforming growth factor–β (TGFβ) (16). Conversely, TGFβ is inhibitory for TFH function in mice (17). Furthermore, the study of human TFH and TFR maturation poses challenges because of the inadequacy of in vitro experimental systems and poor access to lymphoid tissue with different stages of normal follicular cell maturation.

In vitro experiments have been important to clarify the polarization of different effector T cell subsets, such as T helper 1 (TH1), TH2, and regulatory T (Treg) cells. However, the nature of follicular T cell maturation requires an orchestrated interaction with different cell types: First, antigen presentation by dendritic cells initiates the follicular differentiation program; and later, interactions with B cells within a different niche (the T-B border) lead to a mature follicular phenotype (4, 9, 1820). The reliance on distinct microenvironments for T follicular cell maturation has precluded the development of optimal in vitro assays. In humans, most studies have focused on circulating TFH and TFR cells at different time points after immunization (21). However, these studies are insufficient to examine the full maturation process that takes place in lymphoid tissues because mature follicular cells expressing BCL6 cannot be found in circulation. In addition, the poor access to lymphoid tissue in humans and lack of known cell surface receptors for sorting a near-pure population of TFR cells, identified as CXCR5+Foxp3+, from lymphoid tissue create a major hurdle for more detailed analysis of these cells (22).

We previously used samples from B cell–deficient individuals and paired samples of tonsil and blood to establish that most human circulating TFR cells have an immature phenotype (23). Here, we use tissues with increasing frequency of the most mature follicular T cells to capture the different maturation stages of human TFR cells and used single-cell transcriptomics [single-cell RNA sequencing (scRNA-seq)] to dissect the heterogeneity of follicular T cells across those tissues. We tested the hypothesis that, as TFR cells differentiate from precursor thymic Treg cells, they follow a bifurcated trajectory toward circulating TFR cells and the most mature tissue TFR cells. Our results demonstrate that circulating TFR cells do not appear to be precursors of mature TFR cells but rather represent a second arm of the developmental bifurcation arising from Treg cells that seed lymphoid tissues. Our results clarify the interpretation of apparently contradictory results regarding the position of TFR cells in human GCs, possibly arising from studying TFR cells at distinct maturation stages (22, 24). In addition to providing further insight into human TFR cell development, the datasets generated here using TFH and TFR cells collected from different tissues and maturation stages may serve as a valuable resource for future studies seeking to further understand human TFH and TFR biology.


Distinct anatomic locations have different frequency of mature follicular T cells

In this study, we sought to investigate the maturation and tissue adaptation of human TFH and TFR cells using samples from distinct human lymphoid tissues. Although time course experiments would be optimal for these studies, accessing lymphoid tissues longitudinally remains challenging in humans. Here, we examined follicular T cells collected from pediatric tonsils, iliac lymph nodes (LNs) from kidney transplant recipients (before transplantation), and peripheral blood from healthy donors to capture follicular T cells at different maturation stages (fig. S1, A and C). Whereas TFR cells in peripheral blood have immature characteristics as previously described (23), tonsils are enriched for GCs and have a high frequency of the most mature follicular cells defined as PD-1+ICOS+ (Fig. 1, A and B, and fig. S1B). Conversely, human iliac LNs are heterogeneous in terms of GC frequency, most frequently presenting none or very few (<3) GCs (Fig. 1, C and D). Overall, the frequency of mature follicular cells appeared to mirror the expected GC content of the tissue. We therefore compared these tissues with different GC content to obtain insights about the development of follicular T cells, taking advantage of the power of single-cell transcriptomics to reconstruct differentiation trajectories (25).

Fig. 1 TFH and TFR cell phenotype in different human tissues correlates with GC activity.

The content of adult peripheral blood (n = 22), human tonsils (n = 32), and iliac LNs (n = 15) of CXCR5+Foxp3 TFH and CXCR5+Foxp3+ TFR cells were examined. (A) Representative flow cytometry plots of PD-1 and ICOS expression in TFH and TFR cells from different tissues. TFH, TFR, and Treg gates are indicated in the first plot. The quantified results are displayed in fig. S1. (B) Expression of PD-1 and CXCR5 by Foxp3+ cells from the three tissues. (C) Distribution of the 15 LNs analyzed on the basis of GC content. (D) Representative histological sections of GCs with different GC content (GCs are identified by black arrowheads).

Iliac LNs have TFH and TFR cells at different maturation stages

To capture TFH and TFR cells at different stages of maturation, we sorted TFH (320 cells), TFR (512 cells), and Treg (320 cells) cells from two human peripheral blood samples, two iliac LNs, and two tonsils from independent donors for single-cell transcriptomic analysis (Fig. 2A). We followed a “permissive” sorting strategy to prevent exclusion of cells within intermediate stages of differentiation. The sorting of human TFR cells based on surface markers is challenging for two reasons: (i) Many CD25+ cells within the tonsil are Foxp3 (23), and (ii) the most mature TFR cells down-regulate CD25 (2628). We decided to use CD25 to enrich for TFR cells despite the limitation that this sorting strategy would exclude the most mature Foxp3+CD25 TFR cells, since identifying this very low frequency population in our scRNA-seq analysis was not likely to be feasible. Therefore, we hypothesized that studying TFR cell maturation remained possible even in the absence of the most mature CD25 stage. We sorted TFR, Treg, and TFH cells (sorting strategy depicted in fig. S2) and used the Smart-seq2 method (29) to perform indexed scRNA-seq, a plate-based method that generates a dataset where individual cell transcriptomes are indexed to their fluorescence-activated cell sorting (FACS) profile. This approach allowed us to confirm the flow cytometry phenotype of individual cells and investigate their corresponding transcriptional signature.

Fig. 2 scRNA-seq datasets of TFH, TFR, and Treg cells from blood, LN, and tonsil.

(A) Experimental design for scRNA-seq data generation, with indication of the number of cells that passed quality control (two donors per tissue). (B) UMAP visualization of cells from all three tissues that passed quality control. The cells are represented in two main clusters. (C) Differential gene expression between the two clusters defined in (B). Genes that are significantly differentially expressed in clusters 0 and 1 are represented in blue and red, respectively. Genes were considered significant with an adjusted P < 0.05. (D) UMAP projection of the cells defined in (B), plotted based on the tissue of origin. Almost all cells isolated from peripheral blood belong to cluster 0 (yellow), whereas most of the tonsil cells belong to cluster 1, with only a few cells in cluster 0. Cells isolated from LNs are distributed in both clusters. (E) UMAP projection of the cells defined in (B) labeled with the cell type enriched for by sorting. (F) UMAP projection of the cells as in (E), plotted based on their tissue of origin as in (D). These data show that almost all cells from tonsil in cluster 1 were sorted as TFH and TFR, whereas tonsil Tregs are excluded from this cluster. Furthermore, cells from LN in cluster 1 also fell within TFH and TFR sorting gates.

Cells collected from all patients and tissues were analyzed together for computational analysis. We obtained scRNA-seq data from 721 cells passing quality control criteria for downstream analysis (Fig. 2A and fig. S3, A to C). A Uniform Manifold Approximation and Projection (UMAP) of filtered cells shows that those cells are separated into two main clusters (Fig. 2B). Differential gene expression between those two main clusters (Fig. 2C) allowed their assignment as mature versus immature and circulating cells (4). Mature cells (cluster 1) display CXCR5, PDCD1 [programmed cell death protein 1 (PD-1)], ICOS, and SH2D1A SLAM (signalling lymphocytic activation molecule)–associated protein (SAP) as up-regulated genes, whereas immature circulating cells (cluster 0) showed higher expression of CCR7, SELL (CD62L), SELPLG, and S1PR1 (Fig. 2C).

We then overlaid the information of tissue of origin of these cells (Fig. 2D) and the identity of these cells based on the sorting information (Fig. 2, E and F) to gain further insight into the characteristics of each cluster. We found that cells sorted from blood were mostly excluded from cluster 1. Furthermore, Treg cells from all three tissues were also excluded from cluster 1. Cluster 1 contained virtually all cells sorted as TFH and TFR from the tonsil and a proportion of TFR and TFH cells from LNs. These results suggest that, as anticipated, LNs contain a mixture of TFH and TFR cells at different stages of maturation, whereas the GC-rich tonsils have a profile of mature follicular T cells (TFH and TFR).

Maturation stages of TFH cells across tissues

Taking advantage of the segregation of cells based on their maturation states (as described in Fig. 2F), we first used our data to examine the maturation of TFH cells because the maturation of TFH cells has been better characterized than TFR. We aimed to first recreate the expected maturation trajectory of TFH cells as validation for our approach before investigating TFR maturation. For this, we analyzed only the transcriptome of cells sorted as TFH cells from the different tissues. Although TFH cells were now analyzed in the absence of other cell subsets, they still group in two major clusters (Fig. 3A). Cluster 1 represents cells expressing higher levels of genes associated with greater TFH cell maturity, such as CXCL13, PDCD1, and TOX2 (Fig. 3B). In addition, we also observed the segregation of circulating and tonsil TFH cells into clusters 0 and 1, respectively, with LN TFH cells represented in both clusters (Fig. 3C and fig. S4A).

Fig. 3 Maturation stages of TFH cells across different tissues.

(A) UMAP projection of all TFH cells from the three tissues. The cells fall into two major clusters. (B) Differential gene expression between the two clusters defined in (A), highlighting some genes involved in TFH differentiation and maturation. (C) UMAP projection as defined in (A), with TFH cells separated based on the tissue from where they were sorted. Tonsil TFH cells are almost entirely contained in cluster 1, whereas TFH cells from the LN are split between the two clusters. (D) Arranging the TFH cells in a continuum based on their expression profiles shows a pseudo-space relationship of TFH cells across the tissues, with LV0 representing an increased activation/maturation state of these cells. (E) Down-regulation and (F) up-regulation of gene expression of selected genes involved in TFH maturation, along the pseudo-space axis (LV0). Cells have the same color scheme as in (D). The dotted line represents the calculated event switch of gene expression as estimated in pseudo-space. Loess curve (in blue) represents the pattern of gene expression change across the cells. (G) A timeline of the estimated event switch of genes up- and down-regulated along with the TFH maturation across the pseudo-space (upward/purple: up-regulated; downward/yellow: down-regulated; length (k) indicates “activation strength,” representing how fast the gene is up- or down-regulated along the trajectory).

To model the transcriptome changes in the maturation of TFH cells, we obtained a pseudo-time relationship of all TFH cells from the different tissues by calculating latent variables (LVs) using the Bayesian-Gaussian Process Latent Variable Modeling (BGPLVM) (30), adapted from a study of Treg tissue adaptation (31). Automated relevance determination (ARD) determined LV0 as the most relevant and LV4 as the second most relevant LVs (fig. S4B). We observed that LV0 (fig. S4B) can be used to describe the maturation states of TFH cells. Using this method, we were able to align all TFH cells along a maturation/activation pseudo-time (or pseudo-space as cells are from different tissues) (Fig. 3D). We found that TFH cells are ordered from the most immature, where circulating TFH cells predominate, toward intermediate maturation states with LN and some tonsil TFH cells until reaching the most mature states where tonsil TFH cells are mostly represented (Fig. 3D).

This approach also allowed us to investigate gene expression changes in relation to the calculated pseudo-space axis LV0. We found that genes that are well known to be highly expressed in circulating cells, e.g., CCR7 or SELL, are progressively down-regulated with respect to LV0 (Fig. 3E). Conversely, other genes are progressively up-regulated, such as ICOS, BCL6, and CXCL13 (Fig. 3F). Furthermore, it is possible to infer the timing of the relative switch of gene up-/down-regulation by determining the inflection point in gene expression relative to LV0. For example, CXCL13 is up-regulated at a later maturation stage (the calculated inflection point is marked with a dotted line) after the down-regulation of GPR183 (EBI2). The sequential switch of gene expression, across the pseudo-space, for 15 genes implicated in TFH maturation is represented (Fig. 3G), along with information about additional genes and their switch in expression across pseudo-space (fig. S4C). Together, these data show that it is possible to use a pseudo-time algorithm to infer the developmental trajectory of follicular T cells, and therefore, we sought to use a similar approach for TFR cells.

Identification of tissue TFR cells

Although it was possible to directly evaluate the maturation of TFH cells across different tissues, TFR cells could not be assessed directly in the same manner because our permissive sorting strategy was likely to include non-TFR cells. In human peripheral blood, most CD4+CD25+CXCR5+ cells are Foxp3+, whereas in tonsils, there is a large proportion of CXCR5+CD25+Foxp3 cells that do not represent bona fide TFR cells (fig. S5) (23). However, one of the advantages of scRNA-seq is that the transcriptional information can be used to reassign cells sorted as tonsil CXCR5+CD25+ cells as either bona fide TFR or “contaminating T follicular cells” based on gene signature.

To identify bona fide TFR cells for downstream analysis of TFR maturation, we used a comparative approach involving the same method applied in two complementary ways: First, we analyzed cells sorted as TFH and TFR together to highlight the distinct regulatory phenotype and identify cells enriched for this regulatory gene signature and then, inversely, by analyzing cells sorted as Treg and TFR together to highlight the distinct follicular phenotype and identify cells enriched for a follicular gene signature. The combination of the two approaches is expected to allow the independent identification of the same group of cells—bona fide TFR cells expressing both a follicular and regulatory gene signature.

When we analyzed cells sorted as TFH and TFR from the different tissues (blood, LN, and tonsil) using unbiased clustering, we identified five clusters (Fig. 4A). Cells sorted as TFR from tissues (LNs and tonsils) were mainly split into clusters 0 and 3, whereas most tissue TFH cells were contained in cluster 4. Sorted blood TFH and TFR cells were in two adjacent clusters 1 and 2, respectively (Fig. 4A and fig. S6A). We anticipated that cells sorted as TFR from the blood would be homogeneous (almost all are Foxp3+), whereas we expected less than half of the cells from tissues to be TFR cells. Differential gene expression across all clusters identified cluster 3 as the bona fide TFR cells enriched for a transcriptional signature including IL1R2, BCL6, PRDM1, and FOXP3 (Fig. 4B). The contaminating follicular T cells that were captured when sorting tissue CD25+ cells fell in cluster 0 and displayed a transcriptional signature similar to the recently described IL-10-producing T follicular (IL10 TF) cells (32), including high expression of IL10 and IL21. These IL10 TF cells were also previously characterized as CD25+, further supporting their shared identity (32). Cluster 4 expressed high levels of BCL6, CXCL13, and PDCD1, consistent with their classification as tissue TFH cells. Last, the transcriptional profile of clusters 1 and 2 shows a phenotype consistent with T cells found in circulation, including the expression of CCR7, SELL, and S1PR1 genes. Higher expression of IL2RA and FOXP3 in cluster 2 supports their assignment as blood TFR cells. We noticed that cluster 3, which constitutes bona fide TFR cells, also included a few sorted TFH cells. Therefore, we considered the tissue TFR population as the cells sorted as TFR (from LN and tonsil) after removing the cells assigned as IL10 TF (which belong to cluster 0). We then reexamined the profile of both TFH and TFR cell subsets from each tissue and found gene expression characteristics in agreement with the known profile of these cells (Fig. 4C and fig. S6B).

Fig. 4 Identification of bona fide tissue TFR cells.

(A) UMAP visualization of cells sorted as TFH and TFR cells (based on cell surface markers) from all tissues (blood, LN, and tonsil), falling into five clusters. (B) Heatmap of all DEGs across the clusters identified in (A). The gene expression patterns confirm the identity of cluster 3 cells as the bona fide tissue TFR population, whereas cluster 0 defines cells previously described as CXCR5+CD25+ IL10 TF cells. Cluster 4 shows a gene expression pattern corresponding to lymphoid tissue TFH cells. Clusters 1 and 2 correspond to blood TFH and TFR cells, respectively. (C) Violin plots of the representative genes in TFH and bona fide TFR cells, from blood, LN, and tonsil, after removing the IL10 TF cells contaminating the sorted tissue TFR cells. (D) UMAP projection of sorted TFR and Treg cells from all tissues in five clusters. (E) Heatmap of all DEGs across the clusters represented in (D). The gene expression pattern of cluster 0 cells was ascribed to blood TFR cells; cluster 1 to IL10 TF cells; cluster 2 to blood Treg cells; cluster 3 to bona fide tissue TFR cells; and, last, cluster 4 to tissue Treg cells. (F) Violin plots of gene expression in bona fide TFR and Treg cells from the different tissues after removing the contaminating IL10 TF cells. (B and E) Genes were considered significant with an adjusted P < 0.05.

Using the same approach, we next analyzed cells sorted as Treg and TFR to highlight the follicular phenotype of TFR cells and independently identify bona fide TFR cells based on transcriptional signature (Fig. 4D). Similar as before, this approach revealed cluster 1 as IL10 TF cells, with the associated markers IL10 and IL21, whereas cluster 3 represented bona fide TFR cells (Fig. 4E and fig. S7A). Again, removing the IL10 TF cells, it was possible to confirm the follicular signature of TFR cells distinct from Treg cells (Fig. 4F and fig. S7B). We found an overlap in the identity of bona fide TFR cells from the above two approaches; however, a few cells were still classified differently in both analyses. We therefore identified the IL10 TF cells that were common in both analyses and segregated them from the sorted tissue TFR cells. For all downstream analysis, we considered these sorted tissue TFR cells, without IL10 TF cells, as the bona fide tissue TFR cells.

Since we were now able to assign the identity of mature follicular T cells sorted from lymphoid tissue as TFH, TFR, and IL10 TF, we further analyzed the distinct features of these populations relative to each other. First, we used the index sort data to locate the assigned TFH, TFR, and IL10 TF cells in the flow cytometry datasets to investigate their phenotype in relation to the sorting strategy. We found that using conventional staining with CXCR5, ICOS, CD127, and CD25, TFR and IL10 TF cells were not segregated (fig. S8), indicating that an alternative gating strategy is needed to separate these two populations. To investigate this, we performed an unbiased comparison of cell type–specific transcriptomes to identify the key genes expressed by each cell subset (Fig. 5, A to D). Cluster 0 corresponded to bona fide TFR cells with the characteristic expression of FOXP3 and IL1R2. TFH cells (cluster 1) displayed higher expression of BCL6, CD40LG, and CXCL13. Last, IL10 TF cells (cluster 2) shared with TFR cells high expression of IL2RA and displayed high expression of IL10, IL21, and LAG3 (Fig. 5, B to D). Consistent with the phenotype described for IL10 TF cells, the major producers of IL-10, we observed high expression of CTLA4 on these cells (32).

Fig. 5 Comparison of the three major mature CXCR5+ populations in lymphoid tissues.

(A) UMAP visualization of cells assigned as TFH, TFR, and IL10 TF with a mature profile (within cluster 1 of Fig. 2B). (B) Heatmap of a selection of significantly DEGs among the three subpopulations from (A). Genes were considered significant with an adjusted P < 0.05; the selection represents the top DEGs. (C) Violin plots of genes showing a mature phenotype of all cell subsets, albeit with varied expression. (D) Violin plots of genes of interest, namely genes differentially expressed by TFH, TFR, and IL10 TF cells. (E) Histograms representative of the phenotype of TFH (Foxp3CD25CXCR5+ICOS+), TFR (Foxp3+CXCR5+ICOS+), and IL10 TF (FOXP3CD25+CXCR5+PD-1+) cells and (F) respective quantification analyzed by flow cytometry. The gating strategy and fluorescence minus one data are shown in fig. S7. One-way analysis of variance (ANOVA) (nonparametric) with Dunn’s multiple comparison test, *P < 0.05 and **P < 0.01. Data from four tonsils.

We also investigated the expression of other markers, namely, CADM1 and CXCR6, previously described as able to discriminate our populations of interest (32) and genes identified as differential from our data. However, we did not find CADM1, to be differentially expressed between IL10 TF and TFR cells, transcriptionally or at protein level. Although transcripts for CXCR6 were very lowly expressed at the single-cell level in our dataset, we found that a large proportion of TFR cells were CXCR6+ by flow cytometry (Fig. 5, E and F, and fig. S9, A to C). Last, as IKZF3 is highly expressed by all IL10 TF cells (Fig. 5D), we observed that TFR and IL10 TF cells display an opposite preference for the expression of IKZF3 and IKZF2 (Helios): Although virtually all IL10 TF cells express IKZF3, TFR cells show a preference for IKZF2 expression. We confirmed this specificity at the protein level, where we also found that cells positive for Helios (IKZF2) tend not to express IKZF3 (Fig. 5, E and F, and fig. S9D). We conclude that the addition of CXCR6 staining can improve the purity of TFR cells sorted from tonsils (fig. S9E).

For independent validation of our results, we investigated the tissue TFR transcriptional profile in a published scRNA-seq analysis of TFH cells sorted from tonsils (33). In this dataset, cells were enriched for CXCR5+PD1+ T cells, and scRNA-seq was performed using the 10x Chromium platform in contrast to the Smart-seq2 method we used. We analyzed this dataset to identify TFR cells, although these cells were expected to be a minority in this nonenriched pool of cells. UMAP projections of these data (fig. S10A), along with expression profiles of FOXP3 and IL2RA, identified one cluster (cluster 3) likely to contain the CD25+ IL10 TF cells and the TFR cells (fig. S10B). We then reanalyzed the cells from this cluster, leading to the identification of three subclusters (fig. S10C) including one cluster (cluster 0) showing a transcriptional signature of TFR cells with high expression of IL1R2 and FOXP3 and a second cluster (cluster 1) showing a profile characteristic of CD25+ IL10 TF cells (with IL10, IL21, and RUNX2) similar to our datasets (fig. S10D). Overall, our approach allowed us to distinguish CXCR5+Foxp3+ TFR cells from CXCR5+CD25+Foxp3 IL10 TF cells based on gene expression.

Maturation of TFR cells along a bifurcated trajectory

Identification of the bona fide TFR cells based on their transcriptional profile enabled us to next examine TFR cell development and maturation. Previous work demonstrated that murine TFR cells develop preferentially from thymic-derived Treg cells (8, 26, 34). We also previously found that human TFR cells appear to develop from thymic Treg cells after birth following exposure to immunizing antigens (23). Although we found that blood and tissue Treg cells (from tonsil and LN) fall within the same cluster (Fig. 2F), we anticipated some heterogeneity within Treg cell populations, as has been suggested (35). We confirmed that tonsil and blood Treg cell populations have clear differences in their transcriptome and surface marker phenotype, including CTLA-4, CD62L, and CCR7. The differential gene expression analysis of tissue and blood Treg cells showed up-regulation of FOS, JUN, IFITM1, and NFKBIA in tissue Treg cells. (fig. S11). As TFR differentiation occurs in lymphoid tissues (22), we studied the relationship between tissue Treg cells and circulating and lymphoid tissue TFR cells (Fig. 6, A and B). We first compared the expression of genes known to define each cell type and found the observed gene expression pattern to be in line with known signatures (Fig. 6C). We also looked at differentially expressed genes (DEGs) between blood and tissue TFR cells (Fig. 6, D and E) that showed an up-regulation of circulating markers in blood TFR cells, whereas maturation markers were preferentially up-regulated in tissue TFR cells.

Fig. 6 Transcriptional relationship between tissue and circulating TFR cells and Treg cells.

(A) UMAP projection of tissue Treg and circulating and tissue TFR cells. (B) The same UMAP coordinates of tissue Tregs and circulating and tissue TFR cells, plotted on the basis of the tissue of origin. Note that some TFR cells from LNs are considered as having circulating/immature characteristics based on the transcriptional assignment as described in Fig. 4. (C) Violin plots of genes characteristic of the three cell populations. (D) Differential gene expression between blood and tissue TFR cells, with a representation of genes that are significantly up-regulated in mature tissue TFR cells. (E) Genes significantly up-regulated in blood TFR cells. Genes were considered significant with an adjusted P < 0.05.

We next aimed to examine the hypothesis that tissue and circulating TFR cells represent the end points of two divergent developmental trajectories arising from precursor thymic-derived conventional Treg cells in lymphoid tissue (Fig. 7A). Therefore, we used two independent strategies to calculate the ordering of cells along a pseudo-time, to define their putative lineage trajectory. We first calculated LVs (as described for TFH cells) that could define the expected maturation trajectory of these cells in silico (fig. S12). We found that the two most significant LVs were able to describe the tissue and cell-type transcriptional differences among the three populations of cells: LV0 showed a pattern of “gradual” tissue segregation among the cells sorted from either blood or tissue, whereas LV1 highlighted the cell-type differences between Treg and TFR cells (Fig. 7B). As an independent strategy, we generated diffusion maps to calculate an independent pseudo-time ordering (Fig. 7C) (36), along with Slingshot for trajectory inference of these cells, which allows for multiple branching trajectory inference (37). To identify possible differentiation trajectories, we tested the different cell subsets (i.e., tissue Treg, blood TFR, and tissue TFR) as the putative points of origin for trajectory inference with Slingshot. This led to identification of two possible trajectories: Assuming tissue Treg cells as the point of origin led to the inference of a bifurcated trajectory (Fig. 7D), whereas the other two points of origin led to the same linear trajectory but with opposite direction (Fig. 7E).

Fig. 7 A bifurcated maturation trajectory of TFR cells from tissue Treg cells.

(A) Diagram representing the tested hypothesis, with TFR maturation along a bifurcated trajectory from tissue Treg cells toward immature blood TFR cells in one direction and toward fully mature TFR cells in secondary lymphoid organs in the other direction. (B) All three cell populations were placed in a calculated pseudo-space where the LV0 axis segregates blood and tissue cells, whereas the LV1 axis describes variation between Treg cells and the two TFR populations. (C) Diffusion maps showing segregation among the three cell populations. (D) Predicted bifurcating trajectory, calculated using Slingshot when tissue Treg cells were considered as the point of origin, overlaid on the diffusion map. (E) Predicted linear trajectory calculated with Slingshot when tissue TFR or circulating TFR cells were considered as the point of origin (the trajectory is the same in both cases, with opposite direction). (F) Vector field representing the RNA velocity calculated using Velocyto overlaid upon the diffusion maps. Arrows (vectors) on some tissue Treg cells point toward their future state as either blood or tissue TFR cells. (G and H) Spliced (s) and unspliced (u) normalized expression value of key genes, used to calculate RNA velocity. (G) Select genes up-regulated in the direction of tissue TFR cells and (H) select genes up-regulated toward circulating TFR cells. (I) Flow cytometry contour plots showing the absence of CD25+ cells among the most mature (CXCR5++PD-1++) Foxp3+ TFR cells in human tonsils (n = 4).

Next, to identify the most likely trajectory based on the biological state of the cells, we used RNA velocity calculations (38, 39). RNA velocity uses the expression of spliced and unspliced RNA from single cells to predict the future states of individual cells, which facilitate the prediction of the trajectory direction at individual points (i.e., the cells) without the need to input a point of origin. By calculating and projecting the RNA velocity vector field onto the diffusion map, we visualized the direction of the predicted future state of each cell and assessed how these data conform to the three inferred Slingshot trajectories (Fig. 7F). We observed that RNA velocity supports the calculated trajectory with a point of origin among tissue Treg cells, bifurcating toward blood and tissue TFR cells (Fig. 7, D and F).

Last, we examined the RNA velocity of individual genes in the cell subsets represented on the diffusion maps by examining the levels of unspliced and spliced transcripts (Fig. 7, G and H). We found evidence of up-regulation of IL1R2 in tissue TFR cells (abundant unspliced mRNA compared with spliced mRNA), whereas it is absent in blood TFR cells. On the contrary, we observed a down-regulation of IL2RA (encoding CD25) in tissue TFR cells (where most transcripts are already spliced), whereas blood TFR cells maintain a high content of unspliced IL2RA transcripts (Fig. 7G), consistent with the loss of CD25 by the most mature CXCR5++PD-1++ tissue TFR cells (Fig. 7I). The coordinated reduction of IL2RA expression with up-regulation of IL1R2 in maturing tissue TFR cells is especially relevant, as this has been independently confirmed by flow cytometry by another group (27). We also observed a consistent pattern of TFR cells up-regulating other known follicular markers, such as ICOS, BCL6, and TRAF3, in line with their known maturation state in tissues that has previously been mostly described in mice (10, 40). On the other arm of the bifurcation, genes such as LEF1, FOXP1, and CCR7 were all up-regulated in blood TFR cells, in line with their immature and circulating characteristics (Fig. 7H) (13, 22). Examples of RNA velocity of other cell surface markers and transcription factors are also represented (fig. S13A). In addition, we show a selection of surface markers that are likely to change in fully mature tissue TFR cells (fig. S13B). These results support a divergent maturation trajectory of blood and tissue TFR cells from precursor Treg cells.


Our results show that, in vivo, various tissues with different GC content can be used to capture the distinct differentiation states of human follicular T cells. The frequency of mature follicular T cells (presumably GC resident), displaying a PD1+ICOS+ phenotype, was increased in tissues with greater frequency of GCs such as the tonsil. On the other extreme, circulating T follicular cells—and especially TFR cells—have been reported to have a more immature profile (22, 23). It is likely that a different proportion of GCs can be found in LNs from other locations, particularly from the mesenterium, which is very immunologically active. However, we anticipate that the differentiation of TFH and TFR cells follows a similar sequence of events in those different locations as results in mice described a generally conserved TFH and TFR differentiation in different tissues (4).

We took advantage of scRNA-seq to model the maturation of TFH and TFR cells from blood, iliac LNs, and tonsils, thereby capturing different maturation states of these cells. We found that T follicular cells (both TFH and TFR) cluster based on their maturation states, with iliac LNs containing T follicular cells representing intermediate maturation states split between the most immature (blood) and the most mature (presumably representing GC-resident tonsil cells). We show that the well-studied maturation of TFH cells can be modeled using pseudo-time, which allows for the identification of sequential gene expression changes as TFH cells move along their maturation trajectory from the follicle to the GC. Our study provides a view of TFH gene expression dynamics across human tissues and reveals that the gene expression changes associated with human TFH maturation largely overlap with those previously reported to occur in mice. Pseudo-time analysis was also used in mice to study TFH and TH1 differentiation, further supporting a similarity between the differentiation of murine and human cells (4, 9, 41). Our datasets will be available to researchers interested in identifying previously unknown transcriptional information and their modulation along the TFH differentiation trajectory in humans, providing a valuable resource for the community.

A key challenge that we sought to address in our study was the identification of bona fide TFR cells among the mixed population of sorted CXCR5+CD25+ cells. We identified the transcriptome of bona fide TFR cells using two complementary comparisons: first, comparing all cells sorted as CXCR5+ (TFH and enriched-for-TFR) to identify a putative cluster with regulatory characteristics (the bona fide tissue TFR) and, second, by comparing all cells sorted as CD25+ (Treg and enriched-for-TFR) to identify a putative cluster with follicular characteristics (again, presumably the bona fide tissue TFR cells). This dual approach led us to identify the bona fide tissue TFR cells and also to ascribe the remaining CXCR5+CD25+ follicular T cells as the recently described IL10 TF cells (32). In addition, we could confirm the validity of our approach using an independently published dataset of tonsil TFH cells (33).

We and others have shown that TFR cells originate predominantly from thymic Treg cells in mice (8, 26, 34). This led us to investigate the maturation trajectories from tissue Treg cells toward blood and tissue TFR cells. RNA velocity calculations have the advantage of providing a trajectory direction for each cell (represented as a vector), thus estimating the future cell state independently of a point of origin or termination. This approach, combined with trajectory inference with Slingshot, allowed us to identify the bifurcated maturation trajectory from tissue Treg cells toward blood TFR cells or tissue TFR cells. We cannot exclude the hypothesis that some circulating TFR cells can be recruited into the tissue TFR pool. Conceptually, it is appealing to speculate that circulating TFR cells constitute a pool of cells that can access the tissue as a humoral response unfolds. Our analysis, however, did not identify a trajectory originating from blood TFR toward tissue TFR when the blood TFR subset was “forced” as the point of origin. As a consequence, such recruitment of circulating TFR cells into the tissue, if present, is obscured by the dominant differentiation of tissue TFR cells from local Treg cells. Nevertheless, we would like to propose that circulating TFR cells are not a by-product of tissue TFR maturation but may represent a pool of circulating TFR cells from multiple past immunizations, which can be recruited at a low level into ongoing GC reactions, bringing a diverse self-reactive repertoire (34). Similarly, it has been described that the development of TFH cells may also lead to blood TFH that can contribute to subsequent humoral responses elicited by the same antigen (42), although our previous data has shown that in B cell–deficient individuals, the number of circulating TFR cells is maintained, whereas circulating TFH cells are significantly reduced (23).

TFR cells progressively lose CD25 expression as they reach their most mature state (22, 2628), and loss of CD25 expression coincides with the up-regulation of IL1R2 (24, 27). Although analyzing the most mature CD25 TFR cells was infeasible, and only CD25+ TFR cells were captured in our analysis, we were able to confirm that the tissue TFR cells we examined already captured a signature of increased IL1R2 expression that was particularly notable when we directly assessed the TFR maturation trajectories. Our findings regarding the spliced and unspliced IL1R2 and IL2RA (CD25) mRNA in TFR cells are remarkably in line with the expected TFR terminal maturation, including the changes in surface markers that occur, measured using flow cytometry (27). We could observe that tissue TFR cells have spliced IL2RA transcripts but lack unspliced IL2RA, suggesting that they were actively reducing the production of new CD25 molecules, unlike blood TFR cells that maintain high mRNA levels of unspliced IL2RA. Conversely, IL1R2 was starting to be up-regulated in the most mature TFR cells, where almost all transcripts were unspliced, representing recent gene up-regulation. Thus, although we could not sort the most mature CD25IL1R2+ TFR cells, our transcriptional data successfully models the dynamics of the emergence of this phenotype.

Our results support that human TFR cells arise from a bifurcated maturation trajectory. We further provide an extensive resource on the transcriptional activity of human TFH, TFR, and Treg cells from different lymphoid tissues.


Study design

We aimed to investigate the development and maturation of human TFR and TFH cells in secondary lymphoid tissues. We used a combination of different tissues to have access to follicular cells at different stages of maturation. We then performed single-cell transcriptomics to assign each cell to its proper cellular identity and pseudo-time and RNA velocity algorithms to reconstruct the developmental pathway of the different populations based on the individual transcriptomes. We also performed flow cytometry to validate our transcriptomics and identify surface markers to enrich for human TFR.

Human samples

Fresh peripheral blood samples were collected from adult healthy volunteers. Tonsils were collected from healthy children submitted to tonsillectomy because of tonsil hypertrophy. Children with any clinical condition under any drug treatment or submitted to tonsillectomy because of chronic tonsillitis were excluded. Iliac LNs were collected from kidney transplant recipients before transplantation (fig. S1C). Patients with immune-mediated disease or under immune-modulatory treatment were excluded. These studies were approved by the Lisbon Academic Medical Center Ethics Committee (ref. no. 505/14). Informed consent was obtained from all adult volunteers, parents, or legal guardians.

Cell isolation

Peripheral blood mononuclear cells (PBMCs) were isolated from blood samples by Ficoll-gradient medium (Histopaque-1077, Sigma-Aldrich) using SepMate tubes (STEMCELL Technologies). Lymphocytes from tonsils and LNs were also isolated by Ficoll-gradient medium after mechanical disruption. Before cell sorting, samples were enriched for CD4+ T cells using the MojoSort Human CD4 T Cell Isolation Kit (BioLegend). Cells were stained with anti-CD127 (eBioRDr5, eBioscience), anti-CD25 (BC96, eBioscience), anti-CD4 (OKT4, BioLegend), anti-CD45RA (HI100, eBioscience), anti-CD45RO (UCHL1, BioLegend), anti-CXCR5 (J252D4, BioLegend), anti-Foxp3 (PCH101, eBioscience), anti-ICOS (C398.4A, BioLegend), and anti–PD-1 (EH12.2H7, BioLegend). Cell sorting was performed in FACS Aria IIu and Aria III instruments (BD Biosciences), and index sorting was performed at the Flow Cytometry Facility of Instituto Gulbenkian de Ciência.

Flow cytometry

PBMCs and/or lymphocytes from tonsils, prepared as above, were stained with the following antibodies: anti-CD4 (OKT4, BioLegend), anti-CD185/CXCR5 (J252D4, BioLegend), anti-CD278/ICOS (C398.4A, BioLegend), anti-CD25 (BC96, eBioscience), anti–CD279/PD-1 (EH12.2H7, BioLegend), anti-Foxp3 (PCH101, eBioscience), anti-CD45RO (UCHL1, BioLegend), anti-CCR7 (150503, R&D Systems), anti-CD62L (DREG-56, BioLegend), anti- human leukocyte antigen (HLA)–DR (L243, BioLegend), anti–CD152/CTLA-4 (L3D10, BioLegend), anti-CD186/CXCR6 (K041E5, BioLegend), anti-CD223/LAG3 (T47-530, BD Biosciences), anti–CADM-1 (3E1, MBL), anti-CD195/CCR5 (RU0, BD Biosciences), anti-IKZF2/Helios (22F6, BioLegend), and IKZF3/Aiolos (14C4C97, BioLegend). Intracellular Foxp3 and CTLA-4 staining was performed using the Foxp3 Fix/Perm Kit (eBioscience) according to the manufacturer’s instructions. For cell viability staining, the Live/Dead Fixable Aqua Dead Stain Kit (Life Technologies) was used. Samples were acquired on a BD LSRFortessa X-20 cytometer and further analyzed on FlowJo software (TreeStar).

RNA processing for single-cell sequencing

To obtain mRNA libraries from single cells, Smart-seq2 protocol was performed with minor modifications (29). In summary, single cells were sorted into 4 μl of lysis buffer [0.2% (v/v) Triton X-100 (Sigma-Aldrich)], deoxynucleotide triphosphate, and External RNA Controls Consortium (ERCC) RNA Spike-In (1:128 M dilution) in 96-well plates, spun down, and immediately frozen at −80°C. SMARTScribe Reverse Transcriptase (50 U; Clontech) were used for reverse transcription, and 25 cycles were used for cDNA amplification. Sequencing libraries were performed with the Nextera XT DNA Sample Preparation Kit (Illumina) using one-fourth of the volumes indicated by Picelli et al. (29), and they were sequenced on an Illumina HiSeq 2000 sequencer v4 chemistry (paired-end 75–base pair reads).

RNA expression quantification and quality control

Gene expression quantification was carried out using Salmon version 0.11.3 (43), with the default parameters against a modified version of human reference genome GRCh38 containing the ERCC sequences. Quality control steps were done in the following way: Percentage of ERCC expression was calculated for each cell and separated from the true gene expression matrix. Cells used for analysis passed the following filtering criteria: number of genes between 1000 and 3500 with more than 40% of expected mapped reads and ERCC sequences expression of less than 30%. Furthermore, we observed a bias of high ERCC expression in four plates showing 50% less number of genes expressed compared with other plates. Therefore, cells from these four plates were also removed. Last, the cells with less than 5.75 (log10) of total sequenced reads were removed.

Downstream analysis

All the downstream analysis was performed using the R package Seurat v.3.1.1 (44). SCTransform (45) command was used for normalization, finding variable features and scaling data, along with regressing out the sequencing depth and percentage of mitochondrial expression for each cell. For all the analysis, the top 10 principal components (PCs) were used to calculate the Shared Nearest Neighbor (SNN) graph and for UMAP calculations. These PCs were deduced on the basis of significant PCs found using JackStraw analysis through command “JackstrawPlot” along with command “ElbowPlot.” For 10x Chromium public datasets, 20 and 10 PCs were used, respectively. For UMAP, “RunUMAP” was used. Depending on the analysis, different resolutions were used for cluster assessment using the command “FindClusters”: clustering of all cells, resolution 0.3; clustering of TFH cells, resolution 0.8; clustering of TFH and TFR cells, resolution 1.3; and clustering of TFR and Treg cells, resolution 1.2. For public datasets, we used resolutions of 0.2 and 0.8 for global and TFR assessment, respectively. For the analysis of mature TFR, TFH, and IL10 TF from tissue, default normalization using functions NormalizeData and highly Variable genes using FindVariableGenes with top 2000 variable genes and selection method “vst” were used. Data were scaled for all genes. Top 10 PCs were used to calculate SNN grah and for UMAP calculations. FindClusters function for clustering of cells with default resolution (0.8) was used. Differential gene expression was calculated using “FindAllMarkers” with default parameters and were considered significantly differentially expressed if it passed the Bonferroni-adjusted P value of less than 0.05.

Pseudo-space analysis

LVs were calculated using BGPLVM (30), implemented in python package GPy (46) as “GPy.models.BayesianGPLVM” as previously used (31, 41). BGPLVM was used on normalized counts and only considered highly variable genes. The method was run with six LVs and based on their ARD. The LV with the most significant ARD served as the pseudo-space axis. Gene expression switch across the pseudo-space axis was calculated using the R package “Switchde” (47), and genes that showed an absolute correlation of more than 0.25 with the significant LV were used for the analysis. Genes with a q value of <0.05 and with a t0 within LV range were kept for further interpretation.

Trajectory analysis

Diffusion Maps were created using R package Destiny (48) on the normalized counts using all genes. Diffusion component (DC) 1 and 2 were used for further analysis. The R package Slingshot (37) was used for trajectory analysis with the command slingshot, and option clusterLabels were defined as the sorted population and reducedDim as DC1 and DC2. Different start.clus were used to identify the possible lineage trajectories.

RNA velocity analysis

To get the tissue TFR cells for velocity estimates, we excluded the IL10 TF cells and used sorted tissue TFR cells, along with sorted blood TFR and tissue Treg cells. These cells were aligned against the reference genome using STAR (49) against human genome GRCh38, and a bam file was created per cell. Velocyto command line with Smart-seq2 parameter was used to create the loom file for these cells (38). Further RNA velocity estimations were done using the R package “velocyto.R” using the guided tutorial for Smart-seq2 samples (38). For filtration, for each subpopulation, we used 5 as the average expression for spliced reads, 1 for unspliced reads, and 0.5 for spanning reads. Quantile fit was performed on top/bottom 5% of cells. Command “gene.relative.velocity.estimates” was used with k = 15 using spliced and unspliced reads. Command “show.velocity.on.embedding.cor” was used to visualize the RNA velocity on the calculated diffusion maps.


Fig. S1. Quantification of TFH, TFR, and Treg cells from the three tissues used in the study.

Fig. S2. Characteristics of the cell populations from the three tissues used in the study.

Fig. S3. Heterogeneity of the single-cell transcriptomes from the different donors.

Fig. S4. Maturation stages of TFH cells.

Fig. S5. Proportion of Foxp3+ cells within different gates of tonsil.

Fig. S6. Assignment of bona fide TFR cells with comparison of TFH and TFR populations.

Fig. S7. Assignment of bona fide TFR cells with comparison of Treg and TFR populations.

Fig. S8. Flow cytometry profile of cells assigned as TFH, TFR, and IL10 TF.

Fig. S9. Protein validation by flow cytometry.

Fig. S10. Validation of the assignment of bona fide TFR cells with a public dataset from tonsil TFH cells.

Fig. S11. Treg heterogeneity between blood and tonsil.

Fig. S12. The top LVs.

Fig. S13. Calculated RNA velocity of a selection of genes.

Table S1. Differential scRNA-seq gene expression table for mature and immature clusters of all cells (text file).

Table S2. Differential scRNA-seq gene expression table for mature and immature clusters of TFH cells (text file).

Table S3. Differential scRNA-seq gene expression table for identifying bonafide TFR in sorted TFH and TFR cells (text file).

Table S4. Differential scRNA-seq gene expression table for identifying bonafide TFR in sorted TFR and Treg cells (text file).

Table S5. Differential scRNA-seq gene expression table for tissue TFH, TFR, and IL10 TF cells (text file).

Table S6. Differential scRNA-seq gene expression table for blood and tissue TFR cells. (text file).

Table S7. Differential scRNA-seq gene expression table for blood and tissue Treg cells. (text file).

Table S8. Raw data file (Excel spreadsheet).


Acknowledgments: We would like to thank J. Tosello, N. Morais, J. Carneiro, and K. Serre for the helpful discussions and A. Jinat and L. Mamanova for scRNA-seq experimental support. Funding: This work was funded by ENLIGHT-TEN project, European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement 675395 to S.K.; UIDB/04046/2020 and UIDP/04046/2020 grants from FCT, Portugal to M.G.-C.; IGR-Curie 1428 Clinical Investigation Center, Labex DCBIOL (ANR-10-IDEX-0001-02 PSL* and ANR-11-LABX0043), SIRIC (INCa-DGOS-Inserm_12554, projects 2011 to 2017: INCa-DGOS-Inserm_4654) to E.P.; and PTDC/IMI-IMU/7038/2014, EJPRD/0003/2019, and LISBOA-01-0145-FEDER-007391, projeto cofinanciado pelo FEDER através POR Lisboa 2020–Programa Operacional Regional de Lisboa, do PORTUGAL 2020, e pela Fundação para a Ciência e a Tecnologia to L.G. Author contributions: S.K., V.R.F., L.G., S.A.T., M.G.-C., and E.P. conceptualized the project. T.G. and S.K. planned the data analysis. S.K. analyzed and visualized the data. V.R.F. and R.J.M. designed the experiments. V.R.F., A.P.B., and A.Á.-D. collected samples. F.R., D.E., A.Á.D., and R.J.M. performed experiments, M.M. performed the FACS sorting. S.K., F.R., and L.G. wrote the manuscript. E.S. and M.G.-C. provided resources. L.G. supervised the study and acquired funding. Competing interests: R.J.M. has been employed by AstraZeneca since October 2018 and is a stockholder; this had no bearing on this work. In the past 3 years, S.A.T. has consulted for Genentech and Roche and is a member of the Scientific Advisory Boards at Biogen, Foresite Labs, and GlaxoSmithKline. The authors declare that they have no other competing interests. Data and materials availability: RNA-seq data for this work have been deposited in the Array Express database at EMBL-EBI ( under accession number E-MTAB-10364.

Stay Connected to Science Immunology

Navigate This Article