Identity and Diversity of Human Peripheral Th and T Regulatory Cells Defined by Single-Cell Mass Cytometry

Human CD3⁺CD4⁺ Th and T regulatory (Treg) cells are involved in the effector or modulatory function of the immune response (1–5). Th1 cells enhance cytotoxic immune responses against intracellular pathogens and are involved in autoimmunity (1, 3). Th2 cells enhance IgE class-switching and immune responses against helminths and are involved in allergy (1, 3). Th17 cells enhance inflammatory immune responses against some bacteria and fungi and are involved in chronic inflammation (2, 3). T follicular helper (Tfh) cells support B cell maturation and humoral immunity and, like other Th subsets, may be involved in primary or acquired immunodeficiencies (5). In addition to the Th responses, peripherally induced T regulatory type 1 (Tr1) cells and thymic-derived FOXP3-expressing Treg cells, are essential to suppress undesired immune responses and maintain immune homeostasis (4, 6–8). Immune dysregulation may result in alterations in the frequency, function, or location of any one of these CD3⁺CD4⁺ T cell subsets, which determine the severity of different pathologies (4, 9, 10).

Different pathogens can induce diverse T cell responses (3, 11). With high encounter of pathogens over time, humans share a large repertoire of Th1, Th2, Th17, and Tfh cell populations within the memory compartment, as well as Tr1 and Treg cell populations. Intense investigation of Th cell lineages based on their unique cytokine profiles has helped to define combinations of surface and chemokine receptors, as well as transcription factors, specific for each population (12, 13). Likewise, Treg cells were initially identified by their immune suppressive function, followed by intensive investigation for specific surface markers and transcription factors. Table I (14–27) details each Th, Treg, and Tr1 cell population’s defining markers, function, and the additional surface receptors’ and transcription factors’ expression attributed to each main population.

The collective characterization of alterations in different subsets in immune-mediated diseases remains difficult because of the large phenotypic diversity (11, 12) across all CD3⁺CD4⁺ T cell populations, as well as the lack of a comprehensive analysis in healthy conditions. In addition, the concept of lineage definition and stability has been largely challenged by evidence of T cells with common molecular characteristics across different subsets (28), especially in pathologic conditions. This can occur by transition or reprogramming between distinct CD3⁺CD4⁺ T cell subsets. In addition, the possible plasticity between these different Th and Treg cell populations adds a further layer of complexity to the understanding of the immune landscape. Therefore, elucidating the interrelationship between CD3⁺CD4⁺ Th and Treg cell populations is critical to provide a deeper understanding of the immune system, which will pave the way for a comprehensive analysis in immune-mediated diseases (4, 11, 28, 29).

Cytometry by time-of-flight (CyTOF) overcomes the limitations of fluorophore tags by detecting metal-tagged Abs attached to single cells, which allows simultaneous multiparameter interrogation (30). The advent of this technology has opened opportunities to study functional and phenotypic single-cell diversity in depth (31–33), including specific CD3⁺CD4⁺ T cell subsets (34, 35); it also has significant promise for the advancement of clinical diagnosis and prognosis (36, 37).

In this study, we used CyTOF to simultaneously analyze CD3⁺CD4⁺ Th1 (CD183 [CXCR3]⁺ CD194 [CCR4]⁻), Th2 (CCR4⁺CXCR3⁻), Th17 (CD196 [CCR6]⁺ CD161⁺), Tfh (CD185 [CXCR5]⁺), Treg (CD25^hiCD127^lo/−FOXP3⁺), and Tr1 (CD223 [LAG3]⁺ CD49b⁺) cells from the peripheral blood of healthy donors. By evaluating differential expression of transcription factors and chemokine, activation, and inhibitory receptors, we further characterized each baseline population and identified new subsets within. Using data-driven techniques to challenge the distinction of defined human Th, Treg, and Tr1 cell populations, we found shared phenotypes among these subsets, suggesting a continuum of Th and Treg cell populations. Altogether, we present a precise and novel approach to show the multiplicity of CD3⁺CD4⁺ Th, Treg, and Tr1 cells in healthy human blood. The purpose of this study is to establish a comprehensive visualization of the dynamic composition of CD3⁺CD4⁺ T cell subsets at the same time in healthy donors, independently of their functional characterization, which could be used as a baseline to study patients with immune dysregulation.

Sterile X-VIVO 15 with gentamicin and l-glutamine (Lonza), supplemented with 5% (w/v) sterile human serum from male AB plasma (Sigma-Aldrich), was used as cell culture medium (CM). Cell staining medium (CSM) consisted of PBS with 0.05% (w/v) BSA (MACS BSA Solution; Miltenyi Biotec) and 0.02% (w/v) sodium azide. Unlabeled carrier protein–free Abs were purchased from BD Biosciences (San Jose, CA), BioLegend (San Diego, CA), eBioscience (San Diego, CA), and R&D Systems (Minneapolis, MN). Our in-house Ab conjugated with lanthanide isotopes was made using Maxpar X8 Ab labeling kits (Fluidigm), according to the manufacturer’s instructions, except for diluting the final Ab conjugate with PBS-based Ab Stabilizer (Boca Scientific). Isotope-labeled Abs used for mass cytometric staining are listed below. The Abs were purchased from Fluidigm unless otherwise noted. The following Abs were used for surface staining: 141Pr-CCR6 (G034E3), 143Nd-CD45RA (HI100), 145Nd-CD4 (RPA-T4), 148Nd-LAG3 (polyclonal; R&D Systems), 149Sm-CCR4 (205410), 155Gd-CD62L (DREG-56; BioLegend), 161Dy-CD49b (AK-7; BD Biosciences), 163Dy-CXCR3 (G025H7), 164Dy-CD161 (HP-3G10), 165Ho-TIGIT (MBSA43; eBioscience), 166Er-ICOS (DX29; BD Biosciences), 167Er-CD226 (11A8; BioLegend), 168Er-CD8α (SK1), 169Tm-CD25 (2A3), 170Er-CTLA-4 (14D3), 171Yb-CXCR5 (51505), 173Yb-CD3 (UCHT1; BioLegend), 175Lu–PD-1 (EH12.2H7), and 176Yb-CD127 (A019D5). The following Abs were used for intracellular staining: 147Sm-pSTAT5 (47), 151Eu-RORC2 (AFKJS-9; eBioscience), 152Sm-GATA3 (TWAJ; eBioscience), 153Eu-pSTAT1 (58D6), 158Gd-pSTAT3 (4/P-STAT3), 160Gd–T-bet (4B10), and 162Dy-FOXP3 (PCH101). The selection of these 26 markers was based on data in the literature reporting chemokine receptors, transcription factors, activation markers, and the proportion of Th cell and Treg cell subsets identified by flow cytometry in human peripheral blood, as described in Table I.

Blood samples from healthy adults of both genders (median age 56 y, range 30–69 y for 9 of 11 donors tested) were purchased from the Stanford Blood Center or obtained from volunteers within the laboratory, after informed consent, in accordance with the Stanford University Institutional Review Board–approved protocol (IRB-34131).

PBMCs were isolated by Ficoll-Hypaque (GE Healthcare) density gradient and immunostained fresh or processed further for long-term storage. CD4⁺ T cells were isolated by negative selection using magnetic beads (EasySep Human CD4+ T Cell Enrichment Kit, STEMCELL Technologies). Isolated CD4⁺ T cell and PBMC samples were cryopreserved in stocks of 10 × 10⁶ cells per milliliter of FBS with 10% (v/v) DMSO for ≥4 mo. Cells were thawed in warm CM and rested overnight in CM at 37°C, 5% CO₂ before staining.

Human PBMCs and isolated CD4⁺ T cells were stained for mass cytometry, as described by Bodenmiller et al. (38). Samples were analyzed individually, without barcoding or stimulation. Each sample was resuspended to 2 × 10⁶ cells per 200 μl of PBS in FACS tubes and incubated for 1–2 min at room temperature (RT) with 1000:1 parts Cisplatin-¹⁹⁸Pt (Fluidigm). The reaction was blocked with 1:5 parts CM, followed immediately by fixation. Every wash with CSM was performed at 800 × g, RT. We normalized cell volume (100 μl) before adding surface Ab mixture. After permeabilization (15 min on ice, with a final concentration ≥ 90% [v/v] methanol), we washed the cells twice with CSM and normalized cell volume again before adding intracellular Ab mixture. Only 20-min RT incubations of Cell-ID Intercalator-Ir (Fluidigm) were performed.

Cells were passed through a Falcon Tube Strainer Cap (Fisher) and diluted to 0.5 × 10⁶ cells/ml EQ Four Element Calibration Beads (Fluidigm) prediluted with 1:9 parts deionized water. Cells were injected at 45 μl/min into a CyTOF 2 mass cytometer (Fluidigm) outfitted with a Super Sampler (Victorian Airship and Scientific Apparatus) and operated with software v6.0.626. The signal intensity was normalized to EQ calibration beads using the same software, and beads were automatically removed.

Data collected in .fcs file format were normalized for intra- and interfile signal drift using CyTOF software v6.0626 and then analyzed on Cytobank (39). Our gating scheme excluded calibration beads by gating on 140 Ce⁻ events and excluded cell aggregates using DNA signal (193Ir versus 191Ir) and event length, as described (40). Residual non-CD3⁺CD4⁺ T cells were removed by gating CD8⁻CD3⁺CD4⁺ T cell events in purified CD4⁺ T cell and PBMC samples, either fresh or frozen. To identify a given Th cell or Treg cell subset, we use a combination of one to three markers without excluding any of the markers used to define the other subsets.

Every heat map represents differential marker expression between cell populations by normalizing the mean marker intensity of a given cell population to CD3⁺CD4⁺ T cells (arcsinh ratio). This differential expression was then transformed to z-score during visualization in R (ggplots package) for all target populations, including the original CD3⁺CD4⁺ T cell distribution (data not shown).

We applied viSNE (41) in Cytobank to each fresh or frozen sample cohort, using the default settings for number of iterations, perplexity, and theta. We used equal downsampling of total CD3⁺CD4⁺ T cells or proportional downsampling of Th, Treg, and Tr1 cell populations. In each case, 60,000–80,000 events were analyzed from frozen or fresh samples, across all donors. The parameters used in our viSNE analysis were selected based on the optimal separation of the baseline populations on the viSNE map. The 12 markers sufficient to distinguish the Th and Treg cell subsets were CD161, CCR6, CCR4, CXCR3, CXCR5, CD25, CD127, CD45RA, FOXP3, GATA3, T-bet, and RORC2. Additionally, LAG3 and CD49b were necessary to distinguish the Tr1 cell subset.

The distinction and overlap of the different subsets within each baseline population was calculated based on the viSNE analysis (Fig. 2B). Gates were manually drawn based on the position of each subset on the tSNE1 and tSNE2 axes with polygon gate on Cytobank; an illustrative representation of these gates was created with Adobe Photoshop CS6 by coloring the viSNE map blue (ungated cells) and gray (gated cells) and by highlighting the edges of each subpopulation of the gated cells in black. We then calculated the percentage of events from each subset that fell into another subset gate, to develop an overlap matrix. The scaled Venn diagram of the overlapping subsets (Fig. 3B) was drawn such that the frequency of each subset of CD3⁺CD4⁺ T cells (at a minimum of 1%) is represented by the circle diameter.

To maintain consistency for the events used for viSNE and FlowSOM (42), we concatenated the generated viSNE files using Cytobank’s concatenation tool (FCSConcat), before importing these files into the R (v3.3.0) environment. We digitally labeled each donor and Th, Treg, and Tr1 cell populations to track the origin of each cell in later visualizations by creating a designated channel for each parameter.

Each file was converted into a flow object using flowCore from Bioconductor (http://bioconductor.org/packages/release/bioc/html/flowCore.html) and then analyzed with the FlowSOM algorithm (42). Considering the expression of the same 12 or 14 markers used in the previous viSNE analysis, 100 clusters were created using a self-organizing map based–algorithm and then visualized in a minimum spanning tree (MST). The MST shows similar clusters adjacent to each other; however, because most of the edges have similar weights, no unique MST configuration was found for our dataset. In addition, the distance between clusters does not provide any information.

We then applied principal component analysis (PCA) (stats package and ggbiplot package for visualization in R) to the resulting 100 clusters (nodes). The PCA plot displays a simplified representation of these nodes, grouping each Th and Treg cell population inside ellipses that represent the normal distribution of ≥68% of the nodes within each population, capturing all data points within 1 SD from the mean. We have also calculated the correlation of the parameters with each component (determined by the vectors coming from the center of the PCA plot). The length of each vector represents the contribution of each variable to each principal component. FlowSOM and PCA were run in triplicate to validate these findings. The structure of our model remained consistent, and it showed reproducible overlap across all iterations of PCA.

For consistency and efficiency in our data workflow, we used Shiny, a web application framework for R, to integrate multiple R packages for our mass cytometry data analysis in a user-friendly and automated manner.

We compared the frequency of each baseline population of CD3⁺CD4⁺ T cells between fresh and frozen samples using the Welch unpaired t test available in GraphPad Prism v7. Statistical significance was determined with α = 0.05. Each population was analyzed individually, without assuming a consistent SD.

Similarity among Th1, Th2, Th17, Tfh, and Treg cell populations was determined by the number of events from each population within shared nodes after clustering. Only nodes (N) with ≥120 events (≥0.15%) of the target population (X) were considered statistically relevant, when calculating similarity between X and the four other populations (a–d). The proportion of a_i–d_i was weighted by the frequency of X (x_i) for each node, and the sum of these weighted proportions was used to calculate the percentage similarity of X to a–d. Tr1 cells were calculated in the same way, except Tr1 cell similarity was calculated for five other populations (Th1, Th2, Th17, Tfh, and Treg cells) as opposed to the previous four. The following three equations summarize our calculation for similarity:

Our initial CyTOF analysis focused on validating the chemokine and surface receptor frequencies for circulating CD3⁺CD4⁺ Th1, Th2, Th17, Tfh, Treg, and memory Tr1 cell subsets using the previously identified markers (2–8) (Table I). The gating strategy to identify frequencies of CXCR3⁺, CCR4⁺, CCR6⁺, CD161⁺, and FOXP3⁺ cells within the CD3⁺CD4⁺ T cell population is shown in Supplemental Fig. 1A and 1B. Frequencies were within range of those in previous reports (43–49) (Supplemental Table I) and were comparable between fresh and frozen samples, with the exception of CD161 and FOXP3, which were better captured on fresh samples.

Next, we used one to three markers to identify Th1 (CXCR3⁺CCR4⁻), Th2 (CCR4⁺CXCR3⁻), Th17 (CCR6⁺CD161⁺), Tfh (CXCR5⁺), Treg (CD25^hiCD127^lo/−FOXP3⁺) (Supplemental Fig. 1C), and memory Tr1 (LAG3⁺CD49b⁺) cell populations (Supplemental Fig. 1D). The frequencies of these cell populations are shown in Supplemental Table I and are within range of what was previously reported using flow cytometry (24, 26, 50). The frequencies of Th1, Th17, and Treg cell populations were found, by CyTOF, to be significantly different between fresh and frozen samples as a result of the lower expression of CXCR3, CD161, and FOXP3, respectively, in frozen samples. Th2, Th17, and Tfh cells were predominantly memory cells, as defined by the lack of expression of CD45RA (86.9 ± 6.67% of Th2, 88.3 ± 3.93% of Th17, and 76.9 ± 4.75% of Tfh cells were CD45RA⁻, n = 8), whereas Th1 cells and Treg cells were more equally distributed between naive and memory compartments (61.9 ± 9.34% of Th1 cells and 57.7 ± 10.2% of Treg cells were CD45RA⁻, n = 8). These data show that CyTOF analysis recapitulated the frequencies of CD3⁺CD4⁺ T cell subsets observed in total PBMCs by standard flow cytometry (24, 26, 43–48, 50). These data also highlight that Th1, Th17, and Treg cell populations are better identified in fresh samples.

We next measured the expression of the master transcription factors T-bet, GATA3, RORC2, and FOXP3 in unstimulated Th1, Th2, Th17, Tfh, and Treg cell populations (referred to hereafter as baseline populations) using frozen isolated CD4⁺ T cell samples (Fig. 1A). We found the highest expression of T-bet, GATA3, and FOXP3 in Th1, Th2, and Treg cell populations, respectively. No specific RORC2 expression was detected in the nonactivated Th17 population. This finding may be due to the fact that RORC2 expression is activation dependent and more transient compared with other transcription factors. In addition to these transcription factors, we included p-STAT1, p-STAT3, and p-STAT5 in our analysis of frozen CD4⁺ T cell samples to monitor active signaling pathways. No basal p-STAT1 or p-STAT3 expression was observed among the baseline populations. p-STAT5 expression was elevated within Treg cells (25), but it was not significantly different from that of other baseline populations (data not shown).

We next included, in our analysis of frozen isolated CD4⁺ T cell samples, surface markers relevant to further characterize the baseline populations, such as the coinhibitory receptors PD-1, TIGIT, and CTLA-4, the costimulatory receptor ICOS, and the T cell activation and adhesion markers CD226, LAG3, CD49b, and CD62L. The differential marker expression of each baseline population for each individual (n = 8) is shown in Fig. 1B. Th1 cells had elevated expression of CD279 (PD-1) (7/8 individuals) and CD226 (5/8 individuals). Th2 cells had elevated expression of CD226 (7/8 individuals) and LAG3 (6/8 individuals). Th17 cells had elevated expression of CD278 (ICOS), CD226, LAG3, and CD49b (8/8 individuals for all). Tfh cells had elevated expression of PD-1 (8/8 individuals), TIGIT (8/8 individuals), ICOS (7/8 individuals), and CD62L (6/8 individuals). Treg cells had elevated expression of TIGIT and CTLA-4 (8/8 individuals for both). Lastly, there was elevated expression of PD-1, ICOS, CD226, and CD62L (3/3 individuals for all) within the memory Tr1 cell population. Of note, CXCR3 was expressed in Tfh cells (7/8 individuals) and Tr1 cells (3/3 individuals), and CCR4 was expressed in Treg cells (6/8 individuals) and Tr1 cells (3/3 individuals). In addition, Th17 cells expressed CXCR3 and CCR4 in all donors tested, but the intensity of expression varied among donors. Furthermore, CD161 was expressed on memory Tr1 cells (3/3 individuals), but these cells did not coexpress CCR6 like the Th17 cell population. Interestingly, Th1, Th2, and Th17 cells displayed more activation markers, whereas Tfh cells, Treg cells, and Tr1 cells were more positive for coinhibitory and costimulatory receptors, likely indicative of their effector and regulatory functions, respectively. We found no substantial differences in marker expression among baseline populations between fresh and frozen samples.

We further characterized these baseline populations to determine whether the observed proportional changes are due to the presence of small populations of positive cells among the whole or to overall differences in expression within populations.

We next analyzed CD3⁺CD4⁺ T cells using viSNE (41). We initially considered the differential expression of 23 markers, including surface markers, chemokine receptors, transcription factors, and signaling and activation molecules. Then, by sequentially excluding these groups of markers from the analysis, one at a time, we determined that only 12 markers are sufficient to separate the five baseline populations. The CD3⁺CD4⁺ T cell viSNE landscape, based on the expression of these 12 markers, plus two markers that distinguish Tr1 cells (LAG3 and CD49B), is shown in Supplemental Fig. 2.

Ultimately, using viSNE, we analyzed CD3⁺CD4⁺ T cells (frozen samples) based on their unique combined expression of CD161, CCR6, CCR4, CXCR3, CD25, CD127, CXCR5, and CD45RA plus the transcription factors FOXP3, T-bet, GATA3, and RORC2. Fig. 2A shows an example of the workflow used to further characterize the different baseline population subsets by viSNE. We asked whether baseline populations, as traditionally gated, were a single population of cells by showing where they lay on the viSNE map (Fig. 2B). We found a total of 17 subsets, indicating that each baseline population occupies more than one area on the viSNE map; therefore, each population encompasses more than one subset. The frequency of the 17 subsets was quite consistent among donors (Supplemental Fig. 3A). The unselected CD3⁺CD4⁺ T cells were mostly naive T cells, as indicated by the CD45RA expression level (Fig. 2A).

Fig. 3A shows the differential expression of 20 markers for each of the 17 subsets across all individuals (n = 8). Indeed, within each baseline population, we identified multiple subsets with differential expression of transcription factors, costimulatory molecules, and inhibitory or activating receptors. For example, within the Th2 population, we identified subsets with differential expression of GATA3 and FOXP3, which were clearly noticeable at the single-cell level (Supplemental Fig. 3B). A subset of Treg cells expressing GATA3 under baseline conditions and after activation has been described in mice (51, 52). The coexpression of transcription factors GATA3 and FOXP3 demonstrates the interrelationship between Treg and Th2 cell populations. In contrast to what was reported previously (53–56), we did not detect overlapping Th1 or Th17 features in the Treg cell population. The distinction and overlap of the different subsets were calculated based on the viSNE analysis described in Fig. 2B and in 2Materials and Methods (Data analysis). Fig. 3B shows that the space occupied on the viSNE map by Th1 subset 2 completely encompassed Tfh cell subset 3, a Th1-Tfh cell subset previously described as Tfh1 (15). Th2 subset 1 and Treg cell subset 1 were superimposable on the viSNE map, indicating that they are identical. Similarly, Th17 subset 2 encompassed Tfh cell subset 2 and Th1 cell subset 1. The latter Th17-Th1 cell subset demonstrates the high degree of plasticity of Th17 cells. This finding is consistent with previous reports showing the ability of Th17 cells to acquire Th1 cell characteristics (57, 58).

Altogether, these Th and Treg cell subset phenotypes and their overlap highlight specific interrelationships between the baseline populations when expanded phenotyping is used simultaneously. Thus, using additional defining markers, and after eliminating those overlapping subsets, we found a total of 15 unique CD3⁺CD4⁺ Th and Treg cell subsets. A summary of defining, shared, and differentially expressed markers for each subset identified is shown in Table II. We list subsets with bona fide baseline population phenotypes, as well as the other subsets sharing multiple Th and Treg cell phenotypes. In conclusion, by viSNE analysis, which simultaneously engages all known markers, we established a foundation on which to understand the relationship among the baseline populations and their heterogeneity.

Given the overlap observed with extended phenotyping, we used an unsupervised clustering analysis, which is a more quantitative method, alternative to viSNE, not based on manual gating, to confirm our findings and identify consistent interrelationships among the baseline populations. Using the FlowSOM algorithm, we clustered CD3⁺CD4⁺ T cells (frozen samples) based on phenotypic similarity among baseline populations (see 2Materials and Methods; Data analysis). Supplemental Fig. 4A displays 100 nodes visualized in an MST, each representing a cluster of cells that share expression of the same 12 parameters used in our viSNE analysis. Analysis of these nodes by PCA (Supplemental Fig. 4B) shows a consistent phenotypic overlap among baseline populations. For example, overlap of the Th2 subset with the Treg cell subset (Fig. 3B) is reproduced in the PCA analysis (Supplemental Fig. 4B). In summary, this unsupervised analysis confirms the specific overlap among baseline populations previously observed by viSNE.

We then quantified similarity among the Th1, Th2, Th17, Tfh, and Treg cell populations after clustering to identify baseline population interrelationships. The percentage similarity between each baseline population represents the proportion of cells that shares the same phenotype. Because each node represents a single phenotype, we could count the number of phenotypically similar cells to calculate similarity (see 2Materials and Methods; Statistical analysis). Supplemental Fig. 4C displays the similarity between each pair of baseline populations. We found that Th1, Th2, and Tfh cell populations shared high similarity with Th17 cells. Also, the Th2 population has the highest similarity with the Treg cell population.

The similarity for each baseline population overall corresponds with the phenotypic organization of baseline populations displayed by PCA. Even small patterns of similarity are still apparent by PCA, such as the Treg cell population being least similar to Th17, or Th1 being least similar to Th2. The largest distinction among baseline populations within the PCA was the grouping of Th2-Treg cell populations and Th1-Th17-Tfh cell populations on opposite sides (Supplemental Fig. 4B). We found that the variables with the greatest contribution to this distinction were CCR4, CXCR3, and CD161. In conclusion, these findings reflect the population interrelationships in human peripheral blood.

To characterize Tr1 cells and their relationship with the other baseline populations, we proceeded with viSNE and unsupervised clustering analysis, including LAG3 and CD49b (which identify Tr1 cells), using fresh PBMC gated CD3⁺CD4⁺ T cells (Fig. 4). Tr1 cells were detected as a distinct population in three healthy donors, with a mean frequency of 3.05 ± 1.94% of CD3⁺CD4⁺CD45RA⁻ T cells (Fig. 4A). In addition to the Tr1 cell subset, this analysis confirms the detection of the 15 subsets, previously identified in frozen samples, within the baseline populations shown in Figs. 2B and 3A with a similar expression profile. Moreover, we identified two additional populations using this 14-parameter analysis (12 markers plus LAG3 and CD49b) on fresh CD3⁺CD4⁺ T cells: a CXCR5⁺ Th2 subset with low GATA3 expression and a CXCR5⁺CD127^low Th17 subset (Fig. 4A), for a total of 18 unique subsets.

The unsupervised clustering analysis using 14 markers showed similarity between Th and Treg cell populations (Fig. 4B) in fresh and frozen samples that was comparable to that previously observed with 12-marker analysis in frozen samples (Fig. 4B versus Supplemental Fig. 4C). The similarity between each cell population remained consistent among fresh and frozen samples, because the pattern of overlap between Th and Treg cell populations remained the same. These data suggest that the baseline populations and their similarity are reproducible among donors and conditions.

Fig. 4C displays an MST of 100 nodes for the unsupervised clustering and the calculated similarity of Tr1 to the Th and Treg cell populations. Tr1 cells were most similar to Th1 (34.46%) and Tfh (30.90%) cells, followed by Th2 (18.11%) and Th17 (15.53%) cells. In contrast, Tr1 cells had no substantial similarity (0.61%) to Treg cells. In conclusion, based on CD49b and LAG3 expression, we could depict Tr1 cells by viSNE and by unsupervised clustering analyses in fresh and frozen samples.

CD3⁺CD4⁺ Th1, Th2, Th17, Tfh, Treg, and Tr1 cell populations have been identified as major players in determining the type of immune response, in health or immune-mediated diseases (1–5). The identification of these subsets is based on specific surface receptors, signaling pathways, transcription factors, and functional properties. Using CyTOF and high-dimensional single-cell analyses, we achieve a systems-level characterization by visualizing the CD3⁺CD4⁺ T cell subsets at the same time and by grouping them in phenotypically distinct populations using unsupervised clustering. Expanded phenotyping using known surface receptors and master transcription factors for Th1, Th2, Th17, Tfh, Treg (defined in this study as baseline populations), and Tr1 cells was performed on purified or gated CD3⁺CD4⁺ T cells from fresh or frozen samples obtained from healthy donors. Our data show that we could identify distinct baseline populations, as previously described by flow cytometry (24, 26, 43–48, 50). However, we also identified heterogeneity within each baseline population, with a total of 18 subsets identified by viSNE. Among these 18 CD3⁺CD4⁺ T cell subsets, we found overlap of a Treg cell subset with a Th2 cell subset, as well as overlap of a Th17 cell subset with Th1 and Tfh cell subsets, reducing the number of the distinct subsets to 16 which were found independently of donor variability. The clustering analysis further showed novel reproducible similarities among certain subsets. Overall, our in-depth characterization of CD3⁺CD4⁺ T cell subsets reveals unprecedented diversity, but also interrelationships among different subsets.

CyTOF analysis confirms that T-bet and FOXP3 are constitutively expressed in unstimulated Th1 cells and Treg cells, respectively, as previously described (8, 16). Additionally, we found elevated GATA3 expression within circulating Th2 cells, which has not been reported. ViSNE analysis shows low, but detectable, RORC2 expression in four of five Tfh cell subsets identified, including the CD161⁺CCR6⁺ Tfh cell subset. RORC2 expression has been reported in germinal center CXCR5⁺ICOS⁺Bcl-6⁺ Tfh cells, associated with the early development of human Tfh and Th17 cells (17). Our results suggest that Tfh cells may retain the ability to express RORC2 outside the germinal center, and they align with the description of Th17-like (CXCR3⁻CCR6⁺) Tfh cells in peripheral blood (5, 59).

Our extended characterization of each baseline population using 23 markers (Fig. 3) confirmed the elevated expression of PD-1 and CD226 in Th1 cells (60–62), TIGIT, CD62L, ICOS, and PD-1 in Tfh cells (21–23, 63), TIGIT and CD152 (CTLA-4) in Treg cells (64, 65), and CD226, ICOS, and PD-1 in Tr1 cells (26, 66). Of note, we found that unstimulated Treg and memory Tr1 cell populations differ in their expression of TIGIT and PD-1, with the former being strongly positive for TIGIT and the latter positive for PD-1. In this analysis, the Treg cell population did not express ICOS. However, the viSNE data show that a CD45RA⁻ Treg cell subset does in fact express ICOS, as previously described for memory FOXP3 Treg cells (67, 68). These data allude to the unique divergence among the Th, Treg, and Tr1 cell populations in blood.

Our viSNE analysis confirms, in a single snapshot of the CD3⁺CD4⁺ T cell immune landscape, the distinction of the different baseline populations, but, at the same time, highlights the heterogeneity within each baseline population, including the existence of new previously undescribed subsets. Of the 16 identified Th, Treg, and Tr1 cell subsets using viSNE, we also described overlapping phenotypes among Th cell and Treg cell subsets. We found a CCR6⁺ Th2 subset, which matches previously described CCR6⁺CCR4⁺ Th17 cells (69). Similarly, within Th17 cells we identified a Th1-like Th17 subset that matches the previously described CXCR3⁺T-bet⁺ Th17 cell population (70) and a Th17-like Tfh cell subset that matches the previously described CCR6⁺ Tfh cell population (59). Among Tfh cells, we saw a Th1-like Tfh cell subset that matches the previously described CXCR3⁺ Tfh or Tfh1 cell population (15). Lastly, among Treg cells, we found Th2-like Treg and CXCR3⁺ Treg cell subsets, which match previously described CCR4⁺ Treg cells and CXCR3⁺ Treg cells, respectively (67, 71). We also found three novel subsets at low frequency, within the Th17, Tfh, and Treg cell populations, all sharing coexpression of CXCR3, CCR4, and CCR6 receptors. The large number of subsets that we could find simultaneously by viSNE speaks for the fact that the overlap between baseline populations is biologically relevant.

We initially found a vicinity continuum of CD3⁺CD4⁺ T cell phenotypes, as indicated by the undivided population of cells displayed by viSNE. However, by looking at each Th, Treg, and Tr1 cell population individually, it became apparent that the CD3⁺CD4⁺ T cell viSNE map is still organized into multiple distinct subsets. Along with other single-cell analyses of T cells in the blood, the described similarity based on shared marker expression agrees with the idea that Th and Treg cell populations may have flexible terminal differentiation (12). In this study, we did not evaluate the differentiation of these T cells, but, by phenotype alone, we show the capability of monitoring the type and frequency of cells sharing multiple Th, Treg, and Tr1 cell defining characteristics. By unsupervised clustering analysis, we confirm the heterogeneity within each baseline population and the overlap among some subsets. We found that Th1 and Th2 cells have no similarities, Tfh cells have similarity with Th1 and Th17 cells, and Th17 cells are similar to Th1, Th2, and Tfh cells. Interestingly, Treg cells were predominantly similar to Th2 cells. This Treg-Th2 cell similarity is in line with the large majority of CCR4⁺CXCR3⁻ Treg cell phenotypes previously identified by CyTOF (35). Likewise, the high degree of similarity between Th17 cells, as well as Tfh cells, and other Th populations, but not with Treg cells, is in line with previous descriptions of their heterogeneity (59, 69, 70). By simultaneously examining all Th, Treg, and Tr1 cell populations using unsupervised analyses, we reaffirmed the overlap and similarity between subsets observed by viSNE.

An important part of this investigation is the ability to provide a tool that measures any or all CD3⁺CD4⁺ T cell subsets and allows addition of markers on demand. Indeed, using our clustering approach and measurement by similarity, the interrelationships among Th1, Th2, Th17, Tfh, and Treg cell populations remain fairly the same between the first and second set of donor samples. The same was also true after the addition of LAG3 and CD49b to our unsupervised clustering analysis, showing no overall change in similarity among baseline populations. The Tr1 cell subset identified using viSNE initially shows no overlap with other Th and Treg cell subsets, confirming that coexpression of LAG3 and CD49b distinguishes Tr1 cells within CD3⁺CD4⁺ T cells (26). By unsupervised analysis, we confirmed that Tr1 cells had no similarity to Treg cells, further supporting the differences previously reported between these two populations (29). However, Tr1 cells did share some similarity with Th1 cells, which suggest that Tr1 cells may share some properties with Th1 cells. Indeed, Tr1 cells produce IFN-γ (6) and have killing capacity against myeloid cells (27). Furthermore, the similarity between Tfh and Tr1 phenotypes that we observed suggests they may share a similar role in regulating T cell–dependent B cell activation, which has only been described for mouse Tr1 cells (72).

In conclusion, this comprehensive snapshot analysis of the CD3⁺CD4⁺ T cell compartment provides an expanded image of Th, Treg, and Tr1 cell diversity, heterogeneity, and interrelationships. We identified novel heterogeneity within each baseline population, as well as phenotypic overlap between some subsets. In addition, this expanded single-cell analysis by CyTOF allows a better definition of CD3⁺CD4⁺ Treg cell populations and how they compare with Th cells. This unique opportunity to measure the similarity between each population using unsupervised clustering unveiled specific interrelationships, indicative of a comprehensive organization of CD3⁺CD4⁺ Th cells and Treg cells. The interrelationships may indicate the existence of additional novel subsets with overlapping characteristics or, rather, suggest the presence of T cell subsets in transition between one another, which may represent their differentiation states. An alternative interpretation of our data is the plasticity within the CD3⁺CD4⁺ T cell compartment in healthy conditions, which might represent a readiness within memory CD3⁺CD4⁺ T cells to quickly switch from one functional “program” to another in response to external triggers, as recently suggested (12, 28). We are currently developing a specific CyTOF panel to further characterize the identified subsets by determining their functional properties upon activation.

Overall, these data in healthy human blood provide new insights into the CD3⁺CD4⁺ T cell network and will constitute a reference for future studies aimed at evaluating Th, Treg, and Tr1 cell subset variations in disease conditions and in response to therapies.

We thank E. Hsieh for excellent technical support with the set-up of the mass cytometric assay, B. Carter for quality instrument operation, Anita Khant for assistance with panel design, and S. Wager for statistical advice. We also thank Prof. Jeff Bluestone for helpful scientific discussion and critical review of the manuscript.

The authors have no financial conflicts of interest.