Research ArticleMALARIA

Single-cell RNA-seq and computational analysis using temporal mixture modeling resolves TH1/TFH fate bifurcation in malaria

See allHide authors and affiliations

Science Immunology  03 Mar 2017:
Vol. 2, Issue 9, eaal2192
DOI: 10.1126/sciimmunol.aal2192

Fork in the road for immune cells

Immune cell differentiation along T helper pathways can profoundly influence the nature of the immune response—from promoting allergy to enhancing inflammation. Lönnberg et al. use single-cell transcriptomics and computational modeling to delineate the molecular cues that guide this decision. During blood-stage Plasmodium infection in mice, the authors track TH1/TFH bifurcation at both the population and single-clone levels and identify genes associated with each path. They demonstrate roles of particular cell types in shaping this decision. This approach provides a broad framework for modeling immune cell differentiation in vivo.


Differentiation of naïve CD4+ T cells into functionally distinct T helper (TH) subsets is crucial for the orchestration of immune responses. Because of extensive heterogeneity and multiple overlapping transcriptional programs in differentiating T cell populations, this process has remained a challenge for systematic dissection in vivo. By using single-cell transcriptomics and computational analysis with a temporal mixtures of Gaussian processes model, termed GPfates, we reconstructed the developmental trajectories of TH1 and TFH (T follicular helper) cells during blood-stage Plasmodium infection in mice. By tracking clonality using endogenous T cell receptor sequences, we first demonstrated that TH1/TFH bifurcation had occurred at both population and single-clone levels. Next, we identified genes whose expression was associated with TH1 or TFH fates and demonstrated a T cell–intrinsic role for Galectin-1 in supporting TH1 differentiation. We also revealed the close molecular relationship between TH1 and interleukin-10–producing Tr1 cells in this infection. TH1 and TFH fates emerged from a highly proliferative precursor that up-regulated aerobic glycolysis and accelerated cell cycling as cytokine expression began. Dynamic gene expression of chemokine receptors around bifurcation predicted roles for cell-cell interaction in driving TH1/TFH fates. In particular, we found that precursor TH cells were coached toward a TH1 but not a TFH fate by inflammatory monocytes. Thus, by integrating genomic and computational approaches, our study has provided two unique resources: a database,, which facilitates discovery of novel factors controlling TH1/TFH fate commitment, and, more generally, GPfates, a modeling framework for characterizing cell differentiation toward multiple fates.


CD4+ T cells are key instructors of the immune system. They can display extensive phenotypic and functional diversity by differentiating into a range of T helper (TH) subsets, including TH1, TH2, TH17, TFH (T follicular helper), TH22, Treg (T regulatory), and TH9 cells, that are distinguished mainly by cytokine and transcription factor expression profiles. Because TH cells can control infections and drive immune-mediated diseases, there remains tremendous interest in the molecular mechanisms that mediate their in vivo differentiation.

Malaria, caused by the protozoan parasite Plasmodium, afflicted 212 million humans in 2015 (1). Both TH1 responses (2) and TFH-dependent antibody responses (3) can independently protect against malaria and are elicited simultaneously in malaria-infected individuals (4), as well as in mice challenged with rodent-infective strains, such as Plasmodium chabaudi chabaudi AS (PcAS) (5). However, the molecular relationships between TH1 and TFH cells remain unclear during Plasmodium infection and, more generally, during any immune challenge. A recent study has demonstrated that the unique T cell receptor (TCR) of a naïve CD4+ T cell imparted a strong preference for either a TH1 or a TFH fate (6). Nevertheless, for many clones, both fates could still emerge, implying that other mechanisms, such as internal stochasticity and cell-extrinsic factors, also govern fate choices in vivo. Transcription factors including T-bet, Gata3, RORγT, and Bcl6 have been reported to drive and stabilize TH fates, leading to their characterization as “lineage-defining” molecules. This has tended to present TH differentiation as a choice between mutually exclusive linear pathways. However, transient coexpression of these transcription factors (for example, of Bcl6/T-bet and Foxp3/RORγT) suggests that overlapping intermediate TH states also exist in vivo. Moreover, substantial heterogeneity occurs in the kinetics of CD4+ T cell responses, resulting in a complex mixture of intermediate states during differentiation, which is not easily resolved via assessment of a small number of molecules.

Conventional dendritic cells (cDCs) are the dominant initial source of antigenic signaling to naïve CD4+ T cells in secondary lymphoid tissues, for example, in the spleens of Plasmodium-infected mice (7). In other models, it was shown that cDCs made long-lasting stable contacts with naïve CD4+ T cells to initiate priming (8). Once activated, CD4+ T cells regained motility, permitting further cellular interactions. Consistent with this observation, activated CD4+ T cells required further antigenic stimulation to optimize clonal expansion and TH differentiation (9); cDCs were considered the most likely candidates (8, 10), with other cell types remaining less explored. Studies of mice with altered monocytic responses suggested roles for these cells in CD4+ T cell priming, specifically in tissues with few cDCs (11). Other reports used cDC deficiency to illustrate that monocytes could activate naive CD4+ T cells (12). However, few in vivo studies have explored roles for monocytes in TH differentiation, where cDC responses remain intact.

Here, we used single-cell RNA sequencing (scRNA-seq) to study Plasmodium-specific TCR transgenic CD4+ T (PbTII) cells during blood-stage PcAS infection in mice. We then used a computational modeling strategy to reconstruct the molecular trajectories of TH1 and TFH cells. Last, we investigated cell-cell interactions based on dynamic expression of chemokines and their receptors, and examined roles for inflammatory monocytes in supporting activated CD4+ T cells toward a TH fate.


scRNA-seq resolves TH1 and TFH cell fates during Plasmodium infection in mice

We used scRNA-seq to elucidate the development and heterogeneity of TH1 and TFH cells during PcAS infection (Fig. 1A and fig. S1). We transferred naïve, proliferative dye-labeled PbTII cells into congenic wild-type mice and recovered them at days 2, 3, 4, and 7 after infection by cell-sorting those expressing the early activation marker CD69 or displaying dilution of the proliferative dye (fig. S2). Flow cytometric measurements of the canonical TH1 markers T-bet (coded by Tbx21) and interferon-γ (IFNγ) and TFH markers CXCR5 and Bcl6 indicated that these subsets emerged in parallel by day 7 after infection (Fig. 1, B to D, and fig. S3) (13, 14). Notably, markers of TH2, TH17, or Treg subsets were not up-regulated by PbTII cells (fig. S4).

Fig. 1 Single-cell mRNA sequencing of PbTII cells.

(A) PbTII cells were transferred from a single donor to multiple recipients. The numbers denote single cells from which mRNA-sequencing (mRNA-seq) data were successfully recorded. Numbers in parentheses refer to the replicate experiment presented in fig. S12. (B and C) Representative fluorescence-activated cell sorting (FACS) plots showing bifurcation of splenic TH1 (T-bet+IFNγ+) and TFH (Bcl6+CXCR5+) PbTII CD4+ T cells at day 7 post-infection (p.i.) with PcAS. (D) Flow cytometry data indicate concurrent differentiation of TH1 (IFNγ+) and TFH (CXCR5+) PbTII CD4+ T cells within the spleen of PcAS-infected mice (n = 4). Index expression is the product of mean fluorescence intensity and proportion IFNγ+ or CXCR5+. Data are representative of two independent experiments. (E) PCA of single PbTII cells at day 7 after infection with PcAS. The arrows represent Pearson correlation with PC1 and PC2. Cell size refers to the number of detected genes. “TH1 signature” and “TFH signature” refer to cumulative expression of genes associated with TH1 or TFH phenotypes [total transcripts per million (TPM) of all genes in the set] (15). (F) Expression levels of the leading 50 genes with the largest PC2 loadings at day 7 (D). Genes were annotated as either TH1- or TFH-associated on the basis of public data sets (15, 37, 44, 47). *Cdk2ap2 appears twice because two alternative genomic annotations exist.

Initially, we used principal components analysis (PCA) to explore the overall transcriptomic landscape of the PbTII cells (fig. S5A). The top principal components were strongly associated with the number of detected unique transcripts [reflective of mRNA content and proliferative status (fig. S5B)] and differentiation (figs. S5C to S7 and table S1). As expected, the variability related to previously established TH1 and TFH gene expression signatures became more prominent with time, separating two subpopulations at day 7 (Fig. 1, E and F) (15). Together, these results suggested a progressive commitment to TH1 and TFH fates, and indicated that single-cell transcriptomes could be used for estimating both proliferative states and degrees of differentiation of individual cells.

Delineation of TH1 and TFH trajectories using a Mixture of Gaussian Processes model

The results from the PCA suggested that variation in PbTII transcriptomes could be used to reconstruct the transcriptional programs that are underlying the TH1 and TFH differentiation. To more explicitly model the temporal dynamics of this differentiation process, we developed GPfates, a temporal mixture model that builds on the Gaussian Process Latent Variable Model (16) and Overlapping Mixtures of Gaussian Processes (OMGP) (17). Briefly, this approach is based on first reconstructing the differentiation trajectory from the observed data (“pseudotime,” Fig. 2, A and B), thereby establishing an order for the cells. Although our model uses the sample time as prior information, the inferred orderings did not strictly adhere to the experimental time points (fig. S8). For example, cells from day 4 after infection were mixed with some of the cells from days 3 and 7 at either end of the day 4 pseudotime distribution. This result is consistent with the idea that bulk assessments of cells at specific time points fail to account for the heterogeneity and differential kinetics of responses made by single cells. To assess the robustness of the established ordering, we repeated this analysis without supplying the experimental sampling times to the model, finding overall consistent results (Comp. Supp. Fig. 8).

Fig. 2 GPfates modeling of bifurcation processes using scRNA-seq data.

(A) Overview of the analysis workflow that underlies GPfates, consisting of dimensionality reduction of high-dimensional single-cell transcriptomes (left), inference of a pseudotemporal ordering of the cells (middle), and the reconstruction of trajectories using temporal mixture modeling (right). These individual steps build on models derived using the Gaussian process framework. Once fitted, GPfates enables for different downstream analyses, including cell orderings, bifurcation time point estimates, and inference of the genes that drive bifurcation events. (B) Illustration of intermediate results obtained from GPfates. Left: A low-dimensional representation, as well as a pseudotemporal ordering of the cells, is inferred using a nonlinear dimensionality reduction (Gaussian Process Latent Variable Model). Temporal trajectories and bifurcations are then reconstructed using a temporal mixture model (Overlapping Mixture of Gaussian Processes), with data-trend assignments per cell. B-GPLVM, Bayesian Gaussian Process Latent Variable Model; 3D, three-dimensional. (C) Low-dimensional representation (2D) of the complete data sets (408 single-cell transcriptomes). The blue line depicts the inferred progression of pseudotime. Text labels illustrate features typical of cells in the corresponding pseudotime region. (D) Inference of two simultaneous trends based on the pseudotime using the temporal mixture model.

In a second step, GPfates uses the inferred temporal orders as input for a nonparametric time series mixture model [OMGP (17)]. This approach revealed two simultaneous trends emerging during pseudotime (Fig. 2, C and D), which separated from each other, indicating that a developmental bifurcation occurred.

In a third step, GPfates uses a change point model (section 4.2 of Supplementary Computational Methods), thereby facilitating to annotate pseudotime after bifurcation. The cell fate split appeared to initiate among early day 4 post-infection cells (in pseudotime; Fig. 2, C and D), an inference that was robust when using bootstrapped subsets of cells (section 6.2 of Supplementary Computational Methods).

We found that differentially expressed genes between the identified trajectories agreed with known TH1/TFH signature genes (Fig. 3, A and B, and fig. S9) (15), strongly suggesting that the fitted mixture components corresponded to cells with TH1 and TFH phenotypes. Notably, these bifurcation trends could not be identified by other published methods for reconstructing bifurcating single-cell trajectories (Comp. Supp. Fig. 14) (1822). We also successfully applied GPfates to resolve bifurcation events in other published data sets (Comp. Supp. Figs. 11 and 12) (23, 24), suggesting that our approach is more generally applicable for studying cellular differentiation using scRNA-seq data.

Fig. 3 The relationship of known TH1 and TFH transcriptomic signatures and the GPfates trajectories.

(A) TH1 and TFH assignment probabilities of individual cells. For differential expression analysis (B), TH1 and TFH were defined as cells with assignment probability of ≥0.8 for the respective trend. (B) Differential expression patterns between cells assigned to TH1 and TFH states. Fold differences (x axis) and the corresponding adjusted P value (y axis) of differential expression for expressed genes (in at least 20% of cells) are shown. Statistical significance was determined using Wilcoxon rank sum test, with Benjamini and Hochberg correction for multiple testing. The horizontal and vertical dashed lines denote adjusted P value of 0.05 and twofold change, respectively. (C) Parallel TH1 and TFH differentiation within cells of a single CD4+ T cell clone. Colors correspond to individual clones determined by sequence analysis of endogenous TCR genes (tables S2 and S3). (D) Identification of genes associated with TH1 and TFH trajectories. For each gene, the expression correlation with pseudotime (x axis) versus the correlation with the TFH trend assignment (y axis) is shown. Gene relevance was determined using the bifurcation statistic (Materials and Methods, fig. S9C). The top 248 bifurcating genes, with bifurcation statistic >49, are represented in colors according to the functional classification of the genes (Supplementary Materials and Methods and table S4). (E) Genes with the strongest association with TH1 or TFH differentiation, filtered using the bifurcation score as in (D). The genes are ranked in descending order of association with the respective trend. Cdk2ap2 appears twice because of alternative genomic annotations. (F) Web address for GPfates database, where the expression kinetics of genes of interest can be visualized. Examples illustrate the top-ranking bifurcating genes from (E).

Lineage barcoding using endogenous TCR sequences reveals TH1/TFH bifurcation from single CD4+ T cells

Although the TCR transgenic approach used in this study minimized the influence of TCR sequence variability on cell fate determination (6), the strain was Rag-sufficient, thus retaining potential for expression of diverse endogenous TCR chains, in addition to the transgenic TCR. Sequence analysis of TCR transcripts in single PbTII cells confirmed universal expression of the PbTII Vα2 and Vβ12 chains (tables S2 and S3), as well as highly diverse, though lower, levels of expression of endogenous TCRα chains in many cells (fig. S10). Rag-sufficient PbTII cells differentiated as effectively as Rag1−/− PbTII cells into both TH1 and TFH cells (fig. S11), indicating that endogenous TCR sequences had not influenced TH fate bifurcation.

Given the vast combinatorial diversity of the endogenous TCR sequences, we used these as unique molecular barcodes to identify ancestrally related PbTII clones. We identified six clones comprising multiple sibling cells. Of these, two consisted of sibling cells that mapped close to the bifurcation point. For the remaining four clones, siblings exhibited highly diverging patterns of differentiation, with three sibling pairs falling at the extremities of the TH1-TFH phenotype spectrum (Fig. 3C). These results demonstrated that TH1/TFH bifurcation had occurred at both population and single-clone levels in our system, with the progeny of a single cell populating both TH1 and TFH compartments.

Transcriptional signatures associated with bifurcation of TH1 and TFH fates

Next, we sought to identify genes whose expression differed between the TH1 and TFH branches. We derived a bifurcation statistic to estimate the concordance with bifurcation for individual genes (see section 4.2 of Supplementary Computational Methods for details and Fig. 3D). Among the highest-ranking genes, the most common pattern was up-regulation along the TH1 branch (Fig. 3D). This suggested that TFH cells were developmentally closer to the shared progenitor state than TH1 cells, because the TH1 fate involved up-regulation of numerous genes not expressed in either the progenitor or TFH states.

To validate the robustness of these gene signatures and the timing of the bifurcation, we repeated the infection and, at days 0, 4, and 7, sequenced additional PbTII cells using the Smart-seq2 protocol (Fig. 1A and fig. S12A). A nonlinear dimensionality reduction indicated that the single cells from both experiments populated similar transcriptional landscapes (fig. S12B) and that the subset characteristic coexpression patterns of the bifurcating genes identified by GPfates emerged by day 7 (fig. S12C). Notably, the day 7 cells from each mouse could be separated into distinct TH1 and TFH subpopulations using the top bifurcating genes (fig. S12D). These results indicated that the bifurcation-associated gene expression patterns were reproducible across experiments and sequencing platforms.

The highest-ranking transcription factors for the bifurcation included Tcf7 for the TFH fate and Id2 for the TH1 fate (Fig. 3, D and E). Tcf7 is required for T cell development and has been recently shown to be instrumental for TFH differentiation (25, 26). It also represented one of the rare genes defined by a decrease in expression when moving toward the TH1 fate. Id2 is an antagonist of Tcf7 and was recently identified as a key driver of TH1 responses (27). As expected, the hallmark TFH transcription factor Bcl6 was also strongly associated with the TFH fate. In TH1 cells, many bifurcating genes encoded immune-related receptors (Fig. 3, D and E), such as Cxcr6 (fig. S13, A and B), Ifngr1, and S1pr1, which mediate egress from secondary lymphoid organs. This was consistent with the notion that TH1 cells can migrate to peripheral tissues and remain receptive to external signals. In contrast, the only bifurcating chemokine receptor associated with a TFH fate was Cxcr5, which is important for trafficking into B cell follicles (28).

Many of the bifurcating genes had no known role in TH differentiation. For example, lgals1 (encoding Galectin-1), a molecule generally implicated in cDC (29) and Treg function (30), was unexpectedly up-regulated in PbTII cells around bifurcation and maintained along the TH1 but not the TFH trajectory (fig. S14A). This observation was confirmed at the protein level (fig. S14B). Next, comparison of TH1/TFH fates in cotransferred wild-type and lgals1−/− PbTII cells during PcAS infection (fig. S14C) revealed a specific role for Galectin-1 in supporting TH1 but not TFH fate (fig. S14D). Together, these data illustrate the potential for the GPfates model to enable identification of factors controlling TH1 and TFH fates. Further examination of bifurcating genes is facilitated by an online database ( accompanying this paper (Fig. 3F).

Coinciding with TH1/TFH differentiation, we also noted up-regulation of Il10 particularly in the TH1 branch (fig. S15A). Most of the Il10-expressing cells also expressed Ifng at equal or higher levels as those expressing Ifng alone (fig. S15, B and C). These data revealed the development of Tr1 cells, defined as interleukin-10 (IL-10)/IFNγ–coexpressing CD4+ T cells. Given that Il10 expression was associated with the TH1 branch, this suggested that Tr1 cells were developmentally related to TH1 cells. Unexpectedly, we found that aside from Il10, only two genes, Trib2 and BC017643, were differentially expressed between Il10/Ifng-coexpressing Tr1 cells and Ifng-expressing TH1 cells (fig. S15D). Furthermore, a comparison of gene expression frequencies between Tr1 and TH1 cells revealed a substantial degree of similarity across the genome (fig. S15E). Together, these data strongly suggest that Tr1 cells derive directly from TH1 cells during blood-stage Plasmodium infection.

Pseudotemporal relationships between cell cycling, aerobic glycolysis, and cytokine expression

Clonal expansion, increased aerobic glycolysis, and cytokine expression are hallmarks of TH cell development whose temporal relationships with each other remain to be fully resolved in vivo. We noted that PbTII cells became highly proliferative around bifurcation, as shown by the up-regulation of Mki67 (Fig. 4, A and B, and fig. S16A) and other known proliferation marker genes (fig. S16B) (31). This correlated with cell cycle activity, as inferred from the scRNA-seq data using the Cyclone tool, and confirmed by flow cytometric measurements of DNA content and cell size (Fig. 4, C and D, and fig. S16C). On day 4 after infection, the cells also increased expression of genes associated with aerobic glycolysis but not oxidative phosphorylation (Fig. 4F), an indication of increased metabolic requirements being met by glucose metabolism and increased mammalian target of rapamycin complex 1 (mTORC1) activity. Consistent with this was the observed elevated levels of ribosomal protein S6 phosphorylation on day 4 after infection (Fig. 4E).

Fig. 4 The bifurcation of T cell fates is accompanied by changes in transcription, proliferation, and metabolism.

(A) Expression kinetics of Mki67, encoding the proliferative marker Ki-67, as a function of pseudotime. (B) Representative FACS plots showing kinetics of CellTrace Violet (CTV) dilution and Ki-67, IFNγ, or CXCR5 expression, with summary graphs showing % of PbTII cells expressing these (after 106 PbTII cells transferred) in uninfected (day 0) and PcAS-infected mice at indicated days after infection (n = 4 mice per time point, with data from individual mice shown in summary graphs; solid line in summary graphs indicates temporal trends fit using a third-order polynominal regression). Data are representative of two independent experiments. (C) Relative cell cycle speed of PbTII cells, determined by measuring the fraction of cells in S, G2, or M phase. Results when allocating cells to cell cycle phases using flow cytometry (fig. S16C) or computational assignments on the basis of the scRNA-seq data are shown. (D) Cell size estimation using forward scatter (FSC) measurements of PbTII cells. (E) Cellular metabolic activity of PbTII cells in naive mice (n = 3) and at days 4 and 7 after infection (n = 6), as determined by flow cytometric assessment of ribosomal protein S6 phosphorylation (p-S6). Histogram and proportions are representative of two independent experiments. ***P < 0.001, one-way analysis of variance (ANOVA) and Tukey’s multiple comparisons tests. (F) Expression kinetics of genes associated with the cell cycle [251 genes derived from Cyclebase (48)], glycolysis (41 genes, GO:0006096), and oxidative phosphorylation (30 genes, GO:0006119) during PcAS infection. Cumulative expression levels of genes in the respective categories per single cell are shown. Data from all cells and mice (Fig. 1A) were pooled. (G) Differential expression analysis comparing the experiment-corrected expression of genes associated with cell cycle (P < 10−103), glycolysis (P < 10−4), and oxidative phosphorylation (P < 10−5) (F) in Ifng-positive (≥10 TPM) and Ifng-negative cells (<10 TPM) at day 4 after infection with PcAS from both experiments combined.

By day 4 after infection, PbTII cells had gone through several rounds of cell division with differing kinetics and with some cells expressing IFNγ. By comparing Ifng-expressing and nonexpressing cells on day 4 after infection, we noted that early Ifng-expressing cells cycled faster and expressed aerobic glycolysis genes more highly than non–cytokine-expressing counterparts (Fig. 4G). Together, our data suggest that around bifurcation, PbTII cells exhibited a highly proliferative and metabolically active state, with those cells cycling fastest and exhibiting most glycolytic activity being the first to acquire the capacity to secrete IFNγ.

Gene dynamics identifies potential decision-making molecules

To elucidate how PbTII cells transitioned from the proliferative precursor state to TH1 and TFH fates, we sought to resolve the hierarchy of gene expression before and during cell fate bifurcation. In addition to genes directly following the bifurcation trend, we reasoned that expression of genes encoding key decision-making molecules might also be likely to be dynamic and peak before the bifurcation. First, to identify these, we selected those genes displaying interesting nonlinear trends in their expression patterns over pseudotime by Gaussian Process regression. This was achieved via a D statistic (model likelihood ratio), where each gene’s expression pattern over pseudotime was tested for variation unexplained by random noise (32). On the basis of the D statistic (>50.0, Fig. 5C), we identified 2061 dynamic genes (Fig. 5A). Second, we ordered these genes according to their peak expression time to provide a temporal overview (Fig. 5A) and noted that a substantial fraction of them peaked around bifurcation. These included the TH1-driving genes Tbx21, Il2ra, and Il2rb, supporting our initial hypothesis. Moreover, cells around bifurcation also transcribed the highest number of genes compared with those at all other points in pseudotime (Fig. 5B).

Fig. 5 Temporal gene expression dynamics during PcAS infection.

(A) Expression patterns over pseudotime shown for top 2061 dynamic genes (defined by D statistic >50). Genes are ordered per peak expression time. TH1 and TFH probability trajectories from the GPfates OMGP model presented at the bottom to depict bifurcation and provide temporal context between the gene expressions and cellular fates. Various dynamically expressed immune receptors, transcription factors, and secreted molecules are annotated. (B) Relationship of transcriptional activity and divergence of TH1 and TFH fates. The number of detected genes per cell is shown across pseudotime. The color of the data points represents trend assignment probability (Fig. 2). TH1 and TFH trajectories from the GPfates OMGP model presented to depict relation to the bifurcating behavior. (C) Gene expression dynamics assessed using D statistic and optimal squared exponential kernel length-scales. Genes with a D statistic of >50 selected as displaying nonlinear expression patterns over pseudotime. The optimal length-scales of the squared exponential kernels of the Gaussian Processes plotted on x axis, where small values indicate that some rapid changes in expression over pseudotime occur. (D) Model summarizing the expression patterns of key chemokine receptors and the transcription factors Id2 and Tcf7 during TH1-TFH cell fate determination. The size of the cell represents proliferative capacity (Fig. 4, A to F).

This model also infers the length-scale of the dynamic model, namely, the degree of fast-acting behavior over pseudotime (Fig. 5C). Using this additional feature, we noted roughly equivalent dynamics for Tbx21, Il2ra, and Il2rb. Furthermore, we noted similar dynamics, though with slightly later peak times, for the chemokine receptors Cxcr5 and Cxcr3. Closer examination of all chemokine receptor genes also revealed peak expression around bifurcation for Ccr4 but not others (fig. S17). Given that Cxcr5 and Cxcr3 have been associated with TFH and TH1 cells, respectively (10, 33, 34), and because they exhibited similar dynamics, we hypothesized that these were competing receptors that directly influenced TH1/TFH fate (Fig. 5D). Assessment of Cxcr3/Cxcr5 coexpression around bifurcation revealed a substantial portion of cells expressing both receptors (fig. S18). Thus, our examination of gene expression dynamics revealed large numbers of genes being expressed and peaking around bifurcation, including not only those associated with clonal expansion but also numerous sequentially expressed transcription factors and receptors with potential to influence TH fate.

Monocytes support activated PbTII cells toward a TH1 but not a TFH fate

Given similar dynamics and peak expression times for Cxcr3 and Cxcr5, and peak expression around bifurcation for Ccr4 (fig. S17), we reasoned that cell-cell interactions via these receptors controlled TH1/TFH fate. Hence, we considered cell types that could control TH fate, specifically around bifurcation. Because B cells supported a TFH fate (fig. S19), we hypothesized that coordinated action by myeloid cells provided competing signals to support a TH1 fate.

To study this, we examined splenic cDCs and inflammatory monocytes before PbTII bifurcation. We sorted CD8α+ and CD11b+ cDCs and Ly6Chi monocytes from naïve and infected mice (fig. S20) and performed scRNA-seq. PCA of cDCs distinguished the two naïve cell types along PC2 (Fig. 6A and fig. S21) with an efficiency consistent with recent data (35) and further highlighted a number of expected and previously unknown cDC subset-specific genes (fig. S22). We next compared naïve cDCs with those from infection (Fig. 6A and fig. S21) and separated these along PC6 (Fig. 6A). Analysis of differential gene expression in cDCs due to infection identified 30 genes, 29 of which were up-regulated (Fig. 6B and fig. S23), including transcription factors Stat1 and Irf1 and CXCR3-attractant chemokines Cxcl9 and Cxcl10. Notably, gene expression patterns among individual cDCs varied according to the gene. For example, Stat1 and Irf1 were expressed by several naïve cDCs and further up-regulated during infection (Fig. 6C). This was similar for Cxcl9, which was expressed by CD8α+ cDCs in naive mice, whereas Cxcl10 was induced only upon infection (Fig. 6C). These data suggested interactions between cDCs and uncommitted CXCR3+ PbTII cells, consistent with a recent study (10). Next, PCA of Ly6Chi monocytes from naïve and infected mice distinguished them from each other along PC2 (Fig. 6D and fig. S24). Differential gene expression analysis between naïve and infected groups uncovered ~100 genes, both up-regulated and down-regulated during infection (Fig. 6E and fig. S25). A high proportion (~40%) of genes up-regulated in cDCs were also induced in Ly6Chi monocytes, including Stat1, Irf1, and Cxcl10 (Fig. 6, E and F), suggesting possible overlapping functionality. In addition, monocytes expressed other chemokines, including Cxcl2, Ccl2, and Ccl3 (Fig. 6, E and F). Furthermore, specific examination of all immune cellular interaction genes (fig. S26) revealed emerging variable expression of Tnf, Cd40, Pdl1, Ccl4, Ccl5, Cxcl16, Cxcl9, and Cxcl11 in monocytes, thus suggesting complex cell-cell interactions for Ly6Chi monocytes during infection.

Fig. 6 Myeloid cells influence TH bifurcation in uncommitted PbTII cells.

(A to C) Splenic CD8α+ and CD11b+ CD8α cDCs from a naïve mouse, mixed cDCs from an infected mouse, and (D to F) Ly6Chi monocytes from naïve and infected mice were analyzed by scRNA-seq, with mRNA reads filtered by minimum expression of 100 TPM in at least two cells. (A and D) PCA showing clustering of (A) cDCs or (D) monocytes. (B and E) Fold change and confidence for differentially expressed genes (19) between infected and naive (B) cDCs or (E) monocytes; genes were filtered on expression in >10 cells; genes satisfying q < 0.05 are colored per function. (C and F) Differentially expressed genes (q < 0.05) in (C) cDCs and (F) monocytes, between naive and infected mice: Cells and genes are ordered according to PC score and loading, respectively. Common genes between heat maps are annotated in (F). (G) Representative FACS histograms and proportions of splenic CD8α+ cDCs, CD8α cDCs, and Ly6Chi monocytes expressing CXCL9 in naive and infected mice; data show individual mice with line at mean and are representative of two independent experiments (n = 4 mice per time point per experiment). (H) PbTII cells were transferred into LysMCre × iDTR mice 1 day before infection. At 3 days after infection, mice were treated with diphtheria toxin (DT) or saline. Proportions of TH1 (T-bethi IFNγ+) and TFH (CXCR5+) PbTII cells at 7 days after infection; data pooled from three independent experiments (n = 5 to 6 per experiment). ****P < 0.0001, Mann-Whitney U test; NS, not significant. (I) Summary model proposes that chemokine interactions between nonbifurcated PbTII cells and myeloid cells support a TH1 fate, whereas B cells support a TFH fate.

Because Cxcl9, Cxcl10, Cxcl11, Ccl2, Ccl3, and Ccl5 signal through Cxcr3 or Ccr4, which were expressed by activated PbTII cells, we next hypothesized that Ly6Chi monocytes, in addition to cDCs (10), interacted with PbTII cells and influenced TH1/TFH fate. To test this, we first assessed chemokine expression at the protein level by Ly6Chi monocytes (Fig. 6G). Kinetics of CXCL9 production was similar in cDCs and Ly6Chi monocytes. Next, we used LysMCre × iDTR mice, in which Ly6Chi monocytes were depleted after PbTII cell activation but before bifurcation (Fig. 6H and fig. S27). We noted a modest reduction in CD68+ splenic macrophages using this approach (fig. S27B). Using this approach, we found that monocytes/macrophages supported a TH1 but not a TFH fate (Fig. 6H). Together, our data support a model in which activated PbTII cells are supported toward either a TFH fate by B cells (fig. S19) or a TH1 fate by chemokine-expressing myeloid cells, including Ly6Chi inflammatory monocytes.


By capturing single CD4+ T cell transcriptomes during an experimental malaria infection, and computationally reconstructing the course of events, we have resolved the bifurcation of naïve CD4+ T cells into TH1 and TFH cells at a molecular level. GPfates modeling of scRNA-seq data is not limited to immune cells or single bifurcation events. This model can also be combined with existing computational workflows, including alternative methods to estimate pseudotemporal dynamics (see section 6.2 of Supplementary Computational Methods) (19, 36). The GPfates approach permits analysis of cellular differentiation toward two fates (Comp. Supp. Fig. 11) and, in principle, toward multiple fates (Comp. Supp. Fig. 12). However, GPfates exhibits some limitations. Most notably, the ability to identify and pinpoint bifurcation events is linked to changes in the transcriptome that reflect these cellular decisions. In particular, because scRNA-seq profiles are subject to high levels of noise, this means that changes will only be detectable with some lag time (Supplementary Computational Methods). The processed expression data and the GPfates model presented in this study can be accessed at, where users can visualize their genes of interest.

Our data provide the framework for revealing molecular insights into the early stages of TH cell differentiation and describe the sequence of transcriptional events before and after the bifurcation of TH1 and TFH fates. Transcriptomic profiling previously suggested developmental similarities between TFH and TH1 cells (37). However, highly immunogenic viral or bacterial infections induced CD4+ T cells to segregate into Bcl6+ (TFH) or Blimp-1+ (TH1) subpopulations within 2 days (38, 39). In our parasitic model, single CD4+ T cell transcriptomes remained remarkably similar until 4 days of infection. Although it is difficult to directly compare infection models, we speculate that Plasmodium infection in mice may not drive TH bifurcation as early as observed with highly immunogenic viruses or bacteria, particularly given evidence of immunosuppression (40).

IL-10–producing Tr1 cells can suppress immune responses, which could aid the treatment of immune-mediated disorders (41) or be detrimental for chronic infections (42). Despite this, their relationship to TH1 cells is not clear (43). In our model, Tr1 cells emerged from the TH1 trajectory. This observation, coupled with similar transcriptomes for TH1 and Tr1 cells, provides evidence that Tr1 cells are highly related to, and derive directly from, TH1 cells in this model. Thus, our modeling of scRNA-seq data revealed molecular relationships between TH1, Tr1, and TFH cells and showed that a single naïve CD4+ T cell can simultaneously give rise to more than one cell fate during experimental malaria.

Activated CD4+ T cells may experience different microenvironments within secondary lymphoid tissue. The observation that bifurcation toward TH1 and TFH fates was preceded by up-regulation of chemokine receptors prompted us to investigate possible cell-cell interactions with chemokine-expressing myeloid cells. Previous studies have highlighted the potential for cDCs in lymph nodes to produce TH1-associated chemokines (10). Our study, which focused on the spleen, further implicated inflammatory monocytes in TH1 support. However, because our transgenic approach for depleting monocytes also removed a small portion of splenic red pulp macrophages, we cannot fully discount the possibility that they may partly contribute to a TH1 fate. Nevertheless, we propose that splenic monocytes/macrophages influence bifurcation by supporting a TH1 fate during Plasmodium infection. Our studies emphasize that although cDCs are key for initiating CD4+ T cell activation in the spleen, other myeloid cells can also promote a TH1 fate in the presence of cDCs. In contrast, given that CXCR5 was the only chemokine receptor notably associated with bifurcation toward a TFH fate, cellular interaction with B cell follicles may be the primary mechanism for supporting a TFH fate. Our model proposes that activated, uncommitted CD4+ T cells become receptive to competing chemoattractant signals from different zones of the spleen, and suggests intercellular communication as a major driver of bifurcation. However, upstream of these processes, internal stochasticity in uncommitted CD4+ T cells may control the balance of chemokine receptor expression. Future experiments combining our integrated single-cell genomics and computational modeling with in vivo positional and trafficking data may reveal molecular relationships between internal stochasticity, migratory behavior, and TH cell fate.


Study design

The goal of this study was to use scRNA-seq to capture the transcriptomes of individual splenic PbTII cells at various time points during the first week of a blood-stage PcAS infection. Multiple mice were used for most time points to test for possible batch effects, with an independent experimental repeat performed on a different scRNA-seq platform. scRNA-seq data were modeled using Gaussian processes, with statistical testing for significance of both genes and cells associated with the Gaussian processes.

Experimental mice and infections

Wild-type and transgenic inbred mouse strains were housed and used in blood-stage Plasmodium infections, as described in Supplementary Materials and Methods.

Flow cytometry

Splenocytes were isolated and assessed by flow cytometry as described in Supplementary Materials and Methods.

Single-cell mRNA sequencing

Single-cell capture and processing, as well as quality control analysis of scRNA-seq data, were performed as described in Supplementary Materials and Methods.


Statistical analyses were conducted using R, Python, or GraphPad Prism. The types of statistical tests and significance levels are described in respective figure legends.


Materials and Methods

Fig. S1. Enrichment of PbTII cells for adoptive transfer.

Fig. S2. Sorting strategy for PbTII cells.

Fig. S3. Flow cytometric assessment of TH1/TFH responses during PcAS infection.

Fig. S4. Expression of subset-specific marker genes in PbTII cells.

Fig. S5. Heterogeneity of activated PbTII cells and variability associated with cell size and differentiation.

Fig. S6. Heterogeneity of TH1/TFH signature gene expression in activated PbTII cells.

Fig. S7. Heterogeneity of the entire PbTII time series and the contribution of TH1 and TFH signature genes to the overall variability.

Fig. S8. The relationship of pseudotime with time points and with the TH1 assignment probability.

Fig. S9. Correlation of GPfates trends with TH1 and TFH signature genes.

Fig. S10. Expression of transgenic and endogenous TCRs.

Fig. S11. Expression of endogenous TCRs does not influence PbTII cell TH1/TFH differentiation.

Fig. S12. Robustness of top bifurcating genes across experiments.

Fig. S13. Flow cytometric validation of CXCR6 expression in PbTII cells before and after bifurcation.

Fig. S14. T cell–intrinsic Galectin-1 supports TH1 fate commitment.

Fig. S15. IL-10– and IFNγ-coproducing Tr1 cells derive from TH1 cells.

Fig. S16. Proliferative burst of activated PbTII cells.

Fig. S17. Kinetics of chemokine receptor expression during PcAS infection according to the GPfates model.

Fig. S18. Coexpression of chemokine receptors at single-cell level during PcAS infection.

Fig. S19. B cells are essential for TFH responses in PbTII cells during PcAS infection.

Fig. S20. Sorting strategy for myeloid cells.

Fig. S21. PCA of cDCs from naïve and infected mice.

Fig. S22. Differential gene expression between single splenic CD8α+ and CD8α cDCs.

Fig. S23. Differentially expressed genes between single naïve and day 3 PcAS-infected cDCs.

Fig. S24. PCA of Ly6Chi monocytes from naïve and infected mice.

Fig. S25. Differentially expressed genes between single Ly6Chi monocytes from naïve and day 3 PcAS-infected mice.

Fig. S26. Expression of immune signaling genes by cDCs and monocytes.

Fig. S27. Myeloid cell depletion in LysMCre × iDTR mice.

Table S1. The expression data from day 7 after infection with functional annotations for genes (15, 37, 4446) (external file).

Table S2. TraCeR detection statistics for the original data set (external file).

Table S3. TraCeR detection statistics for the replicate data set (external file).

Table S4. Annotation of receptors, cytokines, and transcription factors.

Supplementary Computational Methods—The GPfates model (external file)

References (4964)


Acknowledgments: We thank the Wellcome Trust Sanger Institute Sequencing Facility for performing Illumina sequencing, Wellcome Trust Sanger Institute Single Cell Genomics Core Facility for single-cell sample processing, and the Wellcome Trust Sanger Institute Research Support Facility for care of the mice used in these studies. We thank QIMR Berghofer Flow Cytometry and Animal Facilities for expert advice and care of wild-type and transgenic mice. We acknowledge S. Lorenz, J. Cartwright, and T. Metcalf for expert technical assistance. We thank G. Emerton for constructing the database and the interface for accessing the data. We thank M. Raymond for work in defining cytokines and cell surface receptors. S. Ng is acknowledged for assistance with graphic design. Funding: This work was supported by Wellcome Trust (no. WT098051), European Research Council grant ThSWITCH (no. 260507), Australian National Health and Medical Research Council Project grant (number 1028641), and Career Development Fellowship (no. 1028643), University of Queensland; Australian Infectious Diseases Research Centre grants; and the Lister Institute for Preventive Medicine. K.R.J. was supported by grants from European Molecular Biology Laboratory Australia and OzEMalaR. F.O.B. was supported by the Lundbeck Foundation. M.Z. was supported by the Marie Curie Initial Training Networks grant “Machine Learning for Personalized Medicine” (EU FP7-PEOPLE Project Ref 316861, MLPM2012). Author contributions: T.L. and K.R.J. performed the scRNA-seq experiments. V.S. developed the GPfates model in collaboration with M.Z., N.D.L., O.S., and S.A.T. D.F.-R. and W.R.H. generated the PbTII mouse model. K.R.J., R.M., I.S., M.S.F.S., L.G.F., A.S.N., U.N.L., F.S.-F.-G., P.T.B., and C.R.E. performed the mouse experiments. T.L., V.S., K.R.J., L.-H.L., and F.O.B. analyzed the data and interpreted the results M.J.T.S. performed the TCR clonality analysis. T.L., K.R.J., R.M., O.B., A.H., and S.A.T. designed the experiments. O.S., A.H., and S.A.T. cosupervised the study. T.L., V.S., K.R.J., O.S., A.H., and S.A.T. wrote the manuscript. All authors have read and approved the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The data presented in this paper are publicly available in the ArrayExpress database with accession number E-MTAB-4388.
View Abstract

Navigate This Article