In the application of CCA to transcriptome data, the same genes must be used in both transcriptome dataset matrices. Genes are treated as the “sites” of measurements, and the expression of each gene is considered to occur in each sample, providing gene expression as variable. In other words, gene expression is defined as the amounts of transcripts that occur at each site in the genome (i.e., gene) (rows), producing the gene expression for each sample as variable (column). Thus, explanatory variable represents the inclination to differentiation/activation process at each site in the genome, which corresponds to environmental gradients in the analysis of ecological data by ter Braak (34). The CCA will identify the part of main data that can be interpreted by explanatory variables (constrained space) using linear regression, and subsequently, the regressed data will be analyzed by CA, producing new CCA sample and gene spaces, in which the distance represents similarities between cells (and between genes). Finally, the correlation coefficients between explanatory variables and the new gene space will be determined for the visualization of the immunological processes (Figure 1A).

Mathematically, RNA-seq data X ∈ R^p × m is the measurement of m cell samples for the expression of p genes. The j-th column x_j = (x_1j, x_2j, …, x_pj)^T is the expression data of the j-th sample with the expression of p genes, where T indicates transposed vector. Meanwhile, the explanatory data Z will be obtained to provide k explanatory variables, i.e., Z ∈ R^p × k. Z is scaled and standardized (i.e., mean = 0 and variance = 1). For optimizing the measurement of similarities by distance, X is standardized in the χ² metrics, using the sums of transcripts in each sample (i.e., column sums, c), the sums of transcripts in each gene (i.e., row sums, r), and the grand total of transcripts in X (n), the abundance matrix P is defined: P = 1/n X′ − r c^T, and the standardized matrix S=Dr−1/2 P Dc−1/2, where D_r and D_c are the diagonal matrices of r and c [i.e., Dr = diag(r), Dc = diag(c)]. Next, S is projected onto Z by linear regression while weighting by the sums of transcripts in each gene, using the projection matrix Q=Dr1/2 Z (ZT Dr Z)−1 ZT Dr1/2, and thus, the constrained space S* = Q S (35). Next, CCA performs singular value decomposition (SVD) of S* = U D_α V^T, where U^T U = V^T V = I, and D_α is the diagonal matrix of singular values in descending order (α₁ ≥ α₂ ≥ …). Note that eigenvalue λ_j = α_j² (j = 1, …, J), and J is the minimum value of p-1, m-1, and k. Thus, CCA analyses the constrained space and provides new axes where the dispersion of samples and that of genes are maximized. Sample (cell) scores are defined as Dr−1/2 V Dα or Dr1/2 V (when explanatory variable is only one). Because U becomes a linear combination of explanatory variables (36), weighted average scores (WA scores) G_wa (37) are used as gene score and are obtained by projecting the abundance matrix P onto sample scores while weighting by row sums (i.e., transcript abundance): Gwa=Dr−1 P V Dα−1 (or Dr−1 P V, when explanatory variable is only one).

In order to visualize the relationship between explanatory variables and T cells in the CCA solution, the explanatory variables Z will be linearly regressed to each axis of U. Suppose Z′ = [d₁, …, d_k] is weighted explanatory variables by the transcript abundance (i.e., D_r), and the n-th column of Z′ is a vector with the length p, d_n = (d_1n, d_2n, …, d_pn)^T, the biplot values will be obtained by the correlation coefficient ρ = cov(d_n, u_j)/σ(d_n)σ(u_j)(n = 1, …, k; j = b), where b = min(p-1, m-1, k). In the 1D CCA solution, the single biplot value can either be +1 or −1, indicating the direction (increasing/decreasing) of correlated genes in the explanatory variable against that in the main dataset. In order to use the 1D solution as a scoring system, the CCA score (i.e., Axis 1 score) is multiplied by the single biplot value, which indicates positive or negative correlation to Axis 1, ensuring that cells and genes with high scores have high-positive correlations to the explanatory variable.

The map approach enables the comparison of two or more explanatory variables, while the regression process in CCA allows the analysis across two different experiments (34). Explanatory variables (biplot values) are shown by arrows on the CCA map. CCA provides a map that shows the correlations between samples of interest, explanatory variables, and genes. Highly correlated components are closely positioned on the map. Since CCA is based on linear operations only, the interpretation of CCA solution is relatively straightforward using CCA biplot, which allows to directly compare the cell and gene spaces (34).

Importantly, the CCA in this article is different from a more commonly used multivariate method, canonical correlation analysis, which aims to identify shared correlation structures by maximizing the correlations between the two datasets using cross-covariance matrices, and thereby map samples from two datasets in the same sample space (38).

Explanatory variables for CCA were prepared as follows. Differentially expressed genes were selected by a moderated t-test result using the Bioconductor package, limma. The top-ranked differentially expressed genes (according to their p-values) were used for making the explanatory variables. The T cell activation explanatory variable (Tact) was defined by the difference in gene expression between anti-CD3/CD28-stimulated (17 h) CD4⁺ T cells and untreated naïve CD4⁺ T cells from GSE42276 (39). Precisely, genes were selected by FDR < 0.01 and log2 fold change (> 1 or < −1) in the comparison of the gene expression profile of the activated and resting T cells. For the 1D CCA of T cell populations, the expression data of GSE15907 (40) were regressed onto gene values in Tact representing the change in gene expression following T cell activation, and CA was performed for the regressed data and the correlation analysis was done between the new axis and the explanatory variable. For the 2D CCA of T cell populations the expression data of GSE15907 (40) were regressed onto Tact, in combination with Foxp3 and Runx1 explanatory variables representing the effects of Foxp3 and Runx1 expression on CD4⁺ T cells [GSE6939; (41)]. Foxp3 explanatory variable is the log2 fold change of Foxp3-transduced naïve CD4⁺ T cells and empty vector-transduced CD4⁺ T cells. Runx1 explanatory variable is the log2 fold change of Runx1-transduced naïve CD4⁺ T cells and empty vector-transduced CD4⁺ T cells. Subsequently, CA was applied to the regressed data and the correlation analysis was performed between the new axes and the explanatory variables.

RNA-seq expression data of GSE72056 were obtained from single-cell suspension of tumor cells with unknown activation and differentiation statuses (42). Genes with low variances and low maximal values were excluded. In order to identify CD4⁺ T cells, single-cell data were filtered by the expression of CD4 and CD3E to obtain only the CD4⁺ CD3E⁺ single cells, and also by k-means clustering of PCA gene plot to exclude outlier cells (30) for subsequent analysis.

Single-cell combinatorial CCA is a composite approach to understand non-annotated single-cell data by identifying distinct populations of cells and the differentiation processes that are correlated with these populations. Since a T cell population is usually identified by a single lineage specification factor, in the application of SC4A, such a factor will represent the cell population and their differentiation process (Figure 1B). The advantage of this approach is that SC4A uses the gene expression data of a part of the dataset analyzed, and thus the regression analysis of CCA becomes more efficient because of the absence of between-experimental variation, which is usually significant in cross-dataset analysis (34).

Single-cell combinatorial CCA is performed by the following three steps: (1) identification of putative cell populations and candidate genes for explanatory variables by conventional CCA, (2) combinatorial CCA to identify the top-ranked genes to be used as explanatory variables, and (3) the final CCA solution using the selected genes as explanatory variables.

Importantly, SC4A uses the same single cells as samples (rows), rather than genes, and analyzes the expression of genes as variables (columns), in order to cross-analyze main data and the explanatory variables (i.e., selected genes). In other words, gene expression in each single cell is cross-analyzed between two sets of genes in SC4A, while the expression of each gene is cross-analyzed between two sets of cell samples in conventional CCA.

Single-cell combinatorial CCA aims to identify a set of genes that make the dispersion of cell populations maximum in the CCA solution. To achieve this, all the combinations of genes will be used as explanatory variables and tested for discriminating each two populations using CCA. During each combinatorial cycle, two genes are chosen from the total selected genes for all defined single-cell populations in the main dataset and tested for their correlations to one defined cell population vs all other T cells. In the analysis of tumor-infiltrating lymphocytes, the following two cell populations were analyzed by the combinatorial CCA: (1) activated T cells vs resting T cells; (2) FOXP3⁺ cells vs FOXP3⁻ cells; (3) BCL6⁺ cells [as T follicular helper (Tfh)-like T cells] vs BCL6⁻ cells. The most correlated gene to each population (activated T cells, resting T cells, FOXP3⁺ cells, or BCL6⁺ cells) was identified, and these four genes were used as explanatory variables in the final output of SC4A.

All microarray datasets were downloaded from GEO site, and normalized, where appropriate using the Bioconductor package Affy. Data were arranged into an expression matrix where each row corresponds with gene expression for each gene and each column corresponds with cell phonotype (sample). Data were log2-transformed and values above log2(10) were used for analysis. Differentially expressed genes (DEG) the TCR KO dataset and the aTreg dataset were identified by moderated t-statistics. DEG for activated CD44^hi and resting CD44^lo Treg were combined. The CRAN package vegan was used for the computation of CCA. Gene scores used the wa scores of the CCA output by vegan. The Bioconductor package limma was used to perform a moderated t-test. RNA-seq data were preprocessed, normalized, and log-transformed using standard techniques (34).

Heatmaps were generated using the CRAN package gplots. Venn diagram was generated using the R code, overLapper.R, which was downloaded from the Girke lab at Institute for Integrative Genome Biology (http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/overLapper.R). Gene lists were compared for enriched pathways in the REACTOME pathway database using the Bioconductor packages ReactomePA and clusterProfiler. Violin plots show kernel density plots (outside) and the median and interquartile range (inside) of the original gene expression data and were generated by the Bioconductor package ggplot2. The lineage curve was constructed by clustering SC4A/CCA sample scores using an expectation–maximization algorithm (43), and the nodes of these clusters were identified by constructing a minimum spanning tree using the Bioconductor package Slingshot (44).

Firstly, we investigated how T cell activation-related genes are differentially regulated in rTreg and other CD4⁺ T cell populations including Tmem and Teff. To address this multidimensional problem, we applied CCA to the microarray dataset of various CD4⁺ T cells using the explanatory variable for the T cell activation process, which was obtained from the microarray dataset that analyzed resting and activated conventional T cells (“T cell subset data” and “T cell activation data” in Table 1). Thus, we aimed to visualize the cross-level relationships between genes, the T cell populations, and the T cell activation process (Figure 1A). Using the single explanatory variable, the T cell activation process, the solution of CCA is 1D, and the cell sample scores of CCA (represented by Axis 1) provide “T cell activation score” (see Materials and Methods), indicating the level of activation in each cell population relative to the prototype signature of T cell activation, as defined by the explanatory variable Tact. All the naïve T cell populations had low Axis 1 values [i.e., Foxp3⁻ T naïve cells (Tnaive); Tnaive, and non-draining lymph node T cells from BDC TCR transgenic (Tg) mice, which develop type I diabetes]. In contrast, Foxp3⁺ Treg, Tmem, and tissue-infiltrating Teff in the pancreas from BDC Tg mice (i.e., with inflammation in the islets) had high scores (Figure 2). These CCA results indicate that Treg are as “activated” as Tmem and tissue-infiltrating activated Teff at the transcriptomic level.

Next, we addressed whether the highly “activated” status of Treg is dependent on Foxp3. Since Foxp3 suppresses Runx1-mediated transcriptional activities (41), we investigated the same T cell population dataset using the following three explanatory variables: T cell activation (Tact), retroviral Foxp3 transduction (Foxp3), and Runx1 transduction (Runx1) (see Materials and Methods). The CCA solution was 3D, while the first two axes explained the majority of variance (98.8%, Figure 3A). As expected, Tmem, tissue-infiltrating Teff, and Treg had high negative values and showed high correlations to T cell activation (Tact) in Axis 1, whereas only Treg had high correlations with the Foxp3 variable in Axis 2, while Tmem and Teff were correlated with the Runx1 variable in Axis 2 (Figure 3A). By analyzing the gene space of the CCA solution, genes in the lower left quadrant (i.e., negative in both Axes 1 and 2) were enriched with the genes that are involved in T cell activation, effector functions, and Tfh cells, including Cxcr5, Pdcd1 (PD-1) Il21, Ifng, Tbx21 (T-bet), and Mki67 (Ki-67) (Figure 3B). On the other hand, genes in the upper left quadrant (i.e., negative in Axis 1 and positive in Axis 2) were enriched with Treg-associated genes including Ctla4, Il2ra (CD25), Itgae (CD103), Tnfrsf9 (4-1BB), and Tnfrsf4 (OX40) (Figure 3B). These results indicate that a set of activation genes are operating in all the three non-naïve T cell populations (i.e., Treg, Teff, and Tmem), while some of them are more specific to Treg.

Next, we investigated the modules of genes that are differentially regulated between Treg and Tmem, in order to understand the multidimensional identity of Treg and Tmem transcriptomes (i.e., how these populations can be defined in comparison to all relevant populations). Specifically, we asked if Axis 2 captured the differential transcriptional regulations between Tmem and Treg. Importantly, Axis 2 represents Foxp3-driven and Runx1-driven transcriptional effects, which are correlated with Treg and Tmem/Teff, respectively (Figure 4A). This suggests that Axis 2 provides a “scoring system” for regulatory vs effector functions. Thus, the genes in Axis 1-low (precisely, genes above the 25th percentile for positive correlations with Tact) were identified as Tact genes. These genes were subsequently classified into Axis 2-positive (i.e., positive correlations with Foxp3 and Treg) (designated as “Tact-Foxp3 genes”; top left quadrant of CCA gene space in Figure 3B) and Axis 2-negative genes (i.e., positive correlations with Runx1 and Tmem/Teff) (designated as “Tact-Runx1 genes”; bottom left quadrant of CCA gene space in Figure 3B) (Figure 4A). Tact-Runx1 genes contain genes linked to T cell activation (e.g., Mki67), effector functions (e.g., Tbx21), and Tfh differentiation (e.g., Bcl6, Pdcd1), while Tact-Foxp3 genes contain “Treg markers” such as Il2ra (CD25) and Tnfrsf18 (GITR) (Figure 3B).

Intriguingly, heatmap analysis showed that both Treg and Tmem expressed Tact-Foxp3 genes at high levels, compared to naïve and effector T cells (Figure 4B). On the other hand, Tact-Runx1 genes were selectively downregulated in Treg, while their expressions were sustained in Tmem (Figure 4C). In other words, the repression of Tact-Runx1 genes was the major feature of Treg in comparison to Tmem, and Tact-Foxp3 genes are the activation genes, the expression of which is induced by T cell activation in both Treg and Tmem, and is sustained or enhanced even in the presence of Foxp3. Interestingly, comparable selective downregulation of Tact-Runx1 genes was observed in Teff as well (Figure 4C). This suggests that the set of activation genes operating in Teff is different from the ones in Tmem, and that Tmem and Treg share more activation genes than Treg–Teff and Tmem–Teff (Figures 4B,C). These results collectively indicate that the Treg-ness is composed of the induction of the Treg–Tmem shared activation genes (i.e., Tact-Foxp3 genes) and the Foxp3-mediated repression of Tmem-specific genes (i.e., Tact-Runx1 genes), defining the multidimensional identity of Treg.

While the overall activation levels of Treg and Tmem are similar to the ones of the tissue-infiltrating Teff at transcriptional level (Figure 2), when explained by the prototype signature of activation in CD4⁺ T cells (i.e., the explanatory variable Tact), the compositions of the activation genes are different between Treg, Tmem, and Teff (Figures 4B,C). Importantly, many of these activation genes are shared between Treg and Tmem, but not with Teff. The closer similarity between rTreg and Tmem, compared to Teff, is not surprising, considering that both rTreg and Tmem are at the resting state, while Teff are more recently activated and executing effector functions. In addition, the distinct features of Teff may also include their capacity of tissue infiltration and the effects of the microenvironment. These features were not captured by standard t-test analysis (Figure S1 in Supplementary Material).

Tact-Foxp3 genes included the transcription factors Nfat5, Runx2, and Ahr, which were expressed by most of Tmem cells as well (Figure 4D). The Treg-associated markers, Il2ra (CD25), Itgae (CD103), and Tnfrsf18 (GITR) were expressed not only by Treg but also by Tmem at moderate to high levels. Notably, the expression of Ctla4, Ccr4, and Lag3 was high in Treg and Tmem cells, but it was repressed in Teff (Figure 4D). This suggests that Treg and Tmem are in later stages of T cell activation, when the expression of CTLA-4 is induced as a negative feedback mechanism (46), while it is not induced in tissue-infiltrating Teff, presumably because they are more recently activated and actively proliferating.

Tact-Runx1 genes included many cell cycle-related genes (e.g., Ccna2, Cdca2, and Chek2), suggesting that these cells are in cell cycle and proliferating (Figure 4E). The higher expression of Mki67 and Fos in Tmem suggests that these cells had been activated by TCR signals in vivo before the analysis. Tact-Runx1 genes also included the transcription factors Tbx21, Maf, Hif1a, and Bcl6, which have roles in Th1, Th2, Th17, and Tfh differentiation, respectively (47–49). In accordance with this, the Tfh markers Cxcr5 and Pdcd1 were specifically expressed by Tmem, suggesting that Tmem are heterogeneous populations and composed of Th and Tfh cells. These results are compatible with the model that Treg and Tmem constitute the self-reactive T cell population that have a constitutive activation status (7), and that the major function of Foxp3 is to modify the constitutive activation processes by repressing a part of the activation gene modules (i.e., Tact-Runx1 genes) (Figure 5).

We next asked whether the constitutively “activated” status of Treg is dependent on TCR signals. We applied CCA to the microarray data of CD44^hiCD62L^lo aTreg (CD44^hi aTreg) and CD44^loCD62L^hi naïve-like Treg (CD44^lo naïve Treg) from inducible Tcra KO or WT mice (TCR KO data, Table 1; Figure 6A) using the T cell activation variable as explanatory variable. The CCA result showed that CD44^hi aTreg from WT mice only showed high activation scores, compared with all the other groups. Interestingly, Tcra KO CD44^lo naïve Treg showed the lowest scores, and these scores were lower than WT CD44^lo naïve Treg (Figure 6B). These results indicate that TCR signaling is required for the constitutive activation status of Treg, especially CD44^hi aTreg, and suggest that these aTreg are more enriched with the cells that received TCR signals recently, compared to CD44^lo naïve Treg.

In order to further address whether the TCR signal-dependent activation signature of Treg is constitutively maintained or specifically induced by in vivo activation events [presumably as tonic TCR signals (7)], we analyzed the RNA-seq dataset of in vivo aTreg (25) (Table 1). The dataset was generated by depleting a part of Treg by diphtheria toxin (DT) using bone marrow chimera of Foxp3^GFPCreERT2:Rosa26YFP and Foxp3^{GFP DTR} (25). The DT treatment depletes DT receptor (DTR)-expressing Treg from Foxp3^{GFP DTR}, and thus induces a transient inflammation through the reduction of Treg. Foxp3^GFPCreERT2 allows to label Foxp3-expressing cells by YFP at the moment of tamoxifen administration. van der Veeken et al. thus analyzed rTreg from untreated mice (rTreg), aTreg from mice with recent depletion (11 days before the analysis) in an inflammatory condition (aTreg), and ‘memory’ Treg from mice with a distant depletion (60 days before the analysis) (Figure 6C). As expected, the CCA analysis using the T cell activation variable showed that aTreg had higher activation scores than both rTreg and mTreg (Figure 6D). This indicates that the activation mechanisms are more actively operating in aTreg in an inflammatory environment.

In order to further dissect the activation signature of Treg, we obtained the lists of differentially expressed genes (DEG) between WT Treg vs Tcra KO Treg (designated as TCR-dependent genes), and between aTreg and rTreg (designated as aTreg-specific genes, see Materials and Methods). Interestingly, 94/286 genes of Tact-Runx1 genes (Tmem-specific activation genes, repressed in rTreg) are also used during the activation of Treg (Figure 7A), while only 8/119 of Tact-Foxp3 genes (used by Tmem and rTreg) are induced during the activation of Treg (Figure 7B). This indicates that the activation of Treg does not enhance the genes that are used in rTreg, but induces the expression of the Tmem-specific genes that are suppressed in rTreg. On the other hand, 51/286 of Tact-Runx1 and 19/119 of Tact-Foxp3 genes are regulated by TCR signaling (Figures 7A,B), suggesting that the activation status of rTreg and Tmem may be sustained by TCR signals. Pathway analysis showed that Tact-Runx1 and aTreg-specific genes were enriched for cell cycle-related pathways. In contrast, Tact-Foxp3 genes were enriched for pathways related to signal transduction only (Figure 7C). Collectively, the results above suggest that rTreg are maintained by TCR and cytokine signaling, and that the activation of Treg induces the transcriptional activities of a part of Tact-Runx1 genes, which promote proliferation and cell division.

The analyses above showed that Treg are on average more activated than naïve T cells and that the activation status of Treg can be variable. However, it is still unclear whether individual Treg are more activated than any naïve T cells at the single-cell level. The alternative hypothesis is that Treg are enriched with the T cells that have recognized their cognate antigens and have been activated. In order to determine this and thereby understand the dynamics of T cell regulation in vivo, we investigated the single-cell RNA-seq data of tumor-infiltrating T cells from human patients (42) (Table 1), and further enquired how the activation mechanisms are operating in Treg at the single-cell level.

First, we in silico-sorted FOXP3⁺ and FOXP3⁻ CD4⁺ CD3⁺ T cells from unannotated single-cell data from tumors, the tissues of which were dispersed and CD45⁺ cells were sorted by flow cytometry without the use of any other lymphocyte markers (GSE72056, Table 1). Thus, the identities of individual single cells were needed to be identified in a data-oriented manner, and Treg and non-Treg cells in these tumor tissues had unknown activation and differentiation statuses. Thus, we applied CCA to the in silico-sorted single-cell T cell data using the explanatory variables of human activated conventional CD4⁺ T cells (Tact) and resting T cells (Trest; GSE15390, Table 1), aiming to define individual single cells according to their level of activation by their correlations to these two variables (Figure 8A). Here, we used these two variables, Tact and Trest, in order to generate a 2D CCA solution, instead of a single explanatory variable that represents T cell activation by the log2 fold change in gene expression between activated and resting CD4⁺ T cells (c.f. Figures 2 and 6), which produces a 1D CCA solution visualized as a single axis, because we aimed to identify any additional major differentiation process(es) in the Axis 2. The explanatory variables Tact and Trest are both captured by Axis 1 because they represent two poles of one continuum—the spectrum of activation—ranging from “resting” to “activated” cell state. Thus, the CCA aimed to sort the single cells according to their individual levels of activation along the spectrum of activation, capturing the heterogeneity in activation levels in single cells. Remarkably, in the cell sample space of the CCA solution, the majority of FOXP3⁺ T cells were positively correlated with the T cell activation variable Tact (Figure 8B), and thus had negative scores in the Axis 1 (Figure 8B). Here, CCA Axis 1 × (−1) score is designated as the T cell activation score. Thus, using the activation score and FOXP3 expression, the following four subpopulations were defined: “activated FOXP3⁺,” “resting FOXP3⁺,” “activated FOXP3⁻,” and “resting FOXP3⁻” (Figure 8B).

Next, we aimed to determine whether individual activated FOXP3⁺ Treg are more activated than activated FOXP3⁻ non-Treg at the single-cell level. According to the T cell activation score established by the CCA solution in Figure 8B, FOXP3⁺ Treg had significantly higher T cell activation scores than FOXP3⁻ non-Treg on average, as indicated by the higher median in the violin plots and greater density of samples with higher T cell activation scores (Figure 8C), confirming the results by bulk cell analysis (Figure 2). Using the CCA definition of activated and rTreg and non-Treg established in Figure 8B, the T cell activation score neatly captured the activated status of single cells, allocating high positive and negative scores to activated and resting cells, respectively (Figure 8D). Importantly, there was no significant difference between activated FOXP3⁺ and activated FOXP3⁻ cells and between resting FOXP3⁺ and resting FOXP3⁻ cells (Figure 8D), indicating that in the tumor microenvironment, Treg cells are as activated as non-Treg CD4⁺ T cells, which may be enriched with Teff. Strikingly, 32.5% of activated T cells expressed FOXP3, while only 8.2% of resting T cells expressed FOXP3 in Figure 8B. In other words, FOXP3 expression occurred more frequently in activated T cells. Given that the activation signature of Treg is dependent on TCR signals (Figure 6), these results suggest that FOXP3 expression occurs predominantly in the activated T cells that have recognized the tumor antigens and received TCR signals, as a negative feedback mechanism to suppress the effector response against the tumor antigens (7). Alternatively, but not exclusively, FOXP3⁺ T cells may have high-affinity TCRs to self-MHC and/or tumor antigens and be more prone to activation (10).

In the gene space of the CCA solution, genes with strong correlations to activated FOXP3⁺ T cells included FOXP3 itself and common Treg markers such as CTLA4 and IL2RA (CD25), which were found in the upper left quadrant (Axis 1-negative Axis 2-positive). Interestingly, the lower left quadrant (Axis 1-negative Axis 2-negative) contained more Tfh-like or effector-like molecules PDCD1 (PD-1), BCL6, IL21, and IFNG. The chemokine receptor CCR2 had negative scores in Axis 1 (i.e., correlated with Tact), while CCR7 had a high positive score in Axis 1 (i.e., correlated with Trest) (Figure 8E).

Next, we aimed to identify major differentiation and activation processes in the single-cell transcriptomes above. To this end, we have developed a new CCA approach for single-cell analysis (SC4A), which aims to visualize major differentiation/activation processes and the underlying gene regulations (Figure 9A, see Materials and Methods; Figure 1B). First, we classified single cells into the four populations (activated and resting cells, and FOXP3⁺ Treg and FOXP3⁻ non-Treg; Figure 8B), and thereby identified the following four processes as putative differentiation and activation processes in the dataset: T cell activation (activated cells), and naïve-ness (resting cells), FOXP3-driven process (activated FOXP3⁺), and Tfh-like process (activated FOXP3⁻) (Figure 8). Second, based on their high scores in the CCA solution (i.e., either high positive or high negative scores in either Axis 1 or 2 in Figure 8E) and abundant expressions in FOXP3⁺ and FOXP3⁻ cells (data not shown), we selected 12 genes (CCR7, CCR5, CCR4, IL2RA, IL2RB, CTLA4, ICOS, TNFRSF4, TNFRSF9, FOXP3, BCL6, and PDCD1) as the candidate genes for the four processes. From these genes, we identified the most positively correlated gene to each of the four processes using the combinatorial CCA, which tests all the combinations of variables by CCA and obtains the most correlated gene for each population (see Materials and Methods). Thus, PDCD1, FOXP3, CTLA4, and CCR7 were identified as the most correlated genes for activated FOXP3⁻, activated FOXP3⁺, activated T cells, and resting T cells, respectively (Figure S2 in Supplementary Material), which represent the four immunological processes (see above). Finally, using these four genes as explanatory variables, we applied CCA to the single-cell transcriptomes, obtaining the solution of the SC4A approach.

The single-cell space of the SC4A solution showed that activated and resting T cells had negative and positive scores, respectively (Figure 9B). This indicates that Axis 1 represents T cell activation vs naïve-ness. Single cells were successfully clustered into activated FOXP3⁺ Treg, activated FOXP3⁻ non-Treg, and resting T cells. Resting FOXP3⁺ Treg and resting FOXP3⁻ T cells were mostly overlapped (Figure 9B), indicating that the major features in the dataset dominated the difference between these two resting T cell groups. Importantly, the explanatory variable CTLA4, which represents the T cell activation process, was highly correlated with both activated FOXP3⁺ Treg and activated FOXP3⁻ non-Treg at the middle, indicating its neutral position in terms of Tfh and Treg activation processes. As expected, the variable CCR7, which represents naïve-ness, was correlated with both resting FOXP3⁺ Treg and resting FOXP3⁻ T cells. The explanatory variable PDCD1, which represents the Tfh-like process, was highly correlated with activated FOXP3⁻ non-Treg cells, while the variable FOXP3 was correlated with activated FOXP3⁺ Treg. Thus, the single-cell transcriptomes were modeled by the correlations between gene expression, single cells, and the expression of the four key genes, which represent the four immunological processes (Figures 9B,C). PCA and t-distributed stochastic neighbor embedding (t-SNE) did not provide insights into such cross-level relationships or clear separations of the populations (Figure S3 in Supplementary Material).

Next, in order to understand the relationship between the T cell activation signature and FOXP3-driven and Tfh-like processes (Figures 9B), we aimed to identify and characterize genes with high correlations to these processes, which were represented by CTLA4, FOXP3, and PDCD1 explanatory variables, by analyzing the gene space of the final output of SC4A (Figure 9C; see Materials and Methods). As expected, the Tfh genes, IL21 and BCL6 (50), were highly correlated with PDCD1 explanatory variable. IL2RA (CD25) is a Treg marker (51) and was highly correlated with FOXP3 explanatory variable. IL7R and BACH2 are known to be associated with naïve T cells (52, 53) and were positively correlated with CCR7 explanatory variable, which represents the naïve-ness (Figure 9C). Thus, we defined FOXP3-driven Treg genes (magenta circles) and Tfh-like genes (blue circles) according to their high correlation to the FOXP3 and the PDCD1 explanatory variables, respectively, while we designated as activation genes (red circles) the genes that have high correlations with the CTLA4 variable, including LAG3 and CCR5, which were positioned around 0 in Axis 2 (Figure 9C).

The analyses above strongly suggested that there are two major differentiation pathways for those tumor-infiltrating T cells, which are regulated by FOXP3-driven and Tfh-like processes. In order to identify these lineages, we applied an unsupervised clustering algorithm to the sample space of the SC4A/CCA result (Figure 9B), and identified six clusters, to which a pseudotime method (54) was applied, constructing “lineage curves” (Figure 9D; see Materials and Methods). Importantly, the lineage curves had a bifurcation point at Cluster II, which leads to the two distinct differentiation pathways, Tfh-like and FOXP3-driven differentiation. Since cells may change and mature their phenotypes in different dynamics between these two lineages, we designated Tfh-like-associated and FOXP3-associated pseudotime as Tfh-pseudotime and FOXP3-pseudotime (Figure 9D).

In fact, the expression of activation genes was progressively increased in the shared clusters (i.e., Clusters I and II) for the two pseudotimes, and throughout the rest of the FOXP3-pseudotime and the early phase of Tfh-like differentiation (i.e., Cluster III) in Tfh-pseudotime, while it was suppressed toward the end of Tfh-like differentiation (Cluster IV; Figure 9E) in Tfh-pseudotime. Given that Tfh-pseudotime is correlated with PDCD1 expression (Figure 9C), this suggests that PDCD1 expression and the Tfh-effector process are induced during the earlier phases of effector T cell activity, and that the activation processes in PDCD1^high T cells are suppressed, presumably through PD1–PDL1 interactions in the tumor environment (55). Interestingly, FOXP3-driven genes had similar dynamics to activation genes in both FOXP3-pseudotime and Tfh-pseudotime (Figure 9F). By contrast, Tfh-like genes were mostly suppressed throughout FOXP3-pseudotime, while they were progressively induced throughout Tfh-pseudotime (Figure 9G). These differential regulations of two gene modules resonate with those of Tact-Foxp3 genes (which are expressed by both Treg and Tmem) and Tact-Runx1 genes (which are expressed specifically in Tmem, and repressed in Treg) (Figure 4). In fact, FOXP3 expression is weakly induced in some cells in the bifurcating Cluster II and the early phase of Tfh-like differentiation (Cluster III) in Tfh-pseudotime and is progressively increased at and beyond Cluster V in FOXP3-pseudotime (Figure 9H).

RUNX1 is highly expressed in the common Clusters I and II and is downregulated in the transition from Clusters II to III in Tfh-pseudotime, and from Clusters II to V in FOXP3-pseudotime (Figure 9I), which is compatible with the known dynamics of RUNX1 expression: Runx1 is downregulated when naïve CD4⁺ T cells differentiate into activated/effector cells following TCR signaling (56). By analyzing other key genes used as CCA explanatory variables, CTLA4 was induced at the bifurcating point, Cluster II, and onward in both of the lineages at equivalent expression levels (Figure 9J), reflecting the activated status of both effector Tfh-type cells and Treg. Importantly, CTLA4 is a marker of Treg as well as activated effector T cells, where it acts as a negative regulator of T cell proliferation (57).

PDCD1 expression was also induced at the bifurcating point, and throughout Tfh-pseudotime, but specifically suppressed in the early phase of FOXP3-pseudotime (Figure 9K), which is compatible with the known dynamics that PDCD1 is transiently upregulated in activated CD4⁺ T cells as a negative regulatory mechanism to restrain proinflammatory immune responses and maintain peripheral tolerance (58). Further supporting this dynamic perspective, IL2 expression occurs mainly in Cluster II, indicating that these cells are enriched with the T cells that recently recognized antigens (59) (Figure 9L). Consistently, the expression of the naïve T cell marker CCR7 was the highest in the cells with a relatively naïve phenotype in shared Cluster I and was moderately downregulated in the early and late phase of Tfh-pseudotime, and suppressed in most Treg in FOXP3-pseudotime (Figure 9M).

These results collectively support the model that constant activation processes in the tumor microenvironment promote terminal differentiation of the Treg- and Tfh-like lineages in both previously committed and non-committed lineages of T cells. Interestingly, Cluster II is the bifurcation point, in which T cells show moderate activation and together with simultaneous expression of FOXP3 and Tfh-like genes, as well as RUNX1 and PDCD1 expression. These cells are most probably engaged in decision-making about their cell fate and the cell type-specific usage of these genes—whether their transcriptional mechanisms would be used to generate a proinflamamtory or regulatory response. This understanding was possible because SC4A effectively annotated genes and cells, and thereby allowed to identify new cell populations.

Lastly, we aimed to demonstrate the utility of the current approach by discovering exemplary marker genes that distinguish cells in FOXP3⁻ and Tfh-pseudotime (i.e., the FOXP3-driven pathway I–II–V–VI, and the Tfh-like pathway I–II–III–IV) (Figure 10A), and to identify the T cell subpopulations by a flow cytometric visualization of single-cell data. Since activation genes (Figure 9C) are shared by early phases of Tfh-like and FOXP3-driven differentiation (Figure 9E), we took the intersect of these genes and the Tact-Foxp3 genes, which were expressed by both resting Tmem and rTreg in mice (Figure 4). DUSP4 and NFAT5 were such genes and were in fact induced in cells at the activated bifurcating Cluster II and onward in both lineages (Figure 10B). Similarly, in order to identify a marker to distinguish Treg- and Tfh-like cells, firstly, we identified CCR8 and IL2RA in the intersect of FOXP3-driven genes (Figure 9C) and the Tact-Foxp3 genes, which were induced highly and progressively in Treg-lineage cells throughout FOXP3-pseudotime, while mostly suppressed across Tfh-pseudotime (Figure 10C). By contrast, BCL6 and KCNK5 [found in the intersect of Tfh-like genes (Figure 9C) and the Tact-Runx1 genes, which are expressed in resting Tmem but suppressed in rTreg (Figure 4)] were progressively induced across Tfh-pseudotime, while suppressed in FOXP3-pseudotime (Figure 10D).

Lastly, in order to make the newly obtained knowledge easily accessible to experimental immunologists, we showed the expression of NFAT5, IL2RA, CCR8, BCL6, and KCNK5 in the tumor-infiltrating T cells in a flow cytometric format (Figure 10E). The common activation gene NFAT5 in fact captured the majority of Treg-lineage cells (i.e., cells in the Clusters V and VI) and Tfh-like-lineage cells (i.e., cells in the Clusters III and IV). The expression of the Treg-specific genes IL2RA and CCR8 occurred in the majority of FOXP3⁺ Treg-lineage cells, whether NFAT5-positive or negative, but not in most of Tfh-like-lineage cells. By contrast, the Tfh-like-specific genes BCL6 and KCNK5 were expressed by the majority of Tfh-like-lineage cells but were not expressed in Treg-lineage cells (Figure 10E).

Collectively, these results indicate that the SC4A analysis successfully decomposed the gene regulations for T cell activation and Treg and effector T cell differentiation, identifying new cell populations, which include activated cells at the bifurcation point, early and late phases of Treg and Tfh-like differentiation, and their feature genes. In addition, although there must be considerable differences between resting T cells in the secondary lymphoid organs and between humans and mice, our study successfully identified the shared activation processes and the conserved genes that are differentially used between the Treg- and the Teff-lineage cells. We thus identified a shared systems-level mechanism for the differential regulation of activation and differentiation processes in CD4⁺ T cell populations.