Dissecting the dynamic transcriptional landscape of early T helper cell differentiation into Th1, Th2, and Th1/2 hybrid cells

Abstract

Selective differentiation of CD4+ T helper (Th) cells into specialized subsets such as Th1 and Th2 cells is a key element of the adaptive immune system driving appropriate immune responses. Besides those canonical Th-cell lineages, hybrid phenotypes such as Th1/2 cells arise in vivo, and their generation could be reproduced in vitro. While master-regulator transcription factors like T-bet for Th1 and GATA-3 for Th2 cells drive and maintain differentiation into the canonical lineages, the transcriptional architecture of hybrid phenotypes is less well understood. In particular, it has remained unclear whether a hybrid phenotype implies a mixture of the effects of several canonical lineages for each gene, or rather a bimodal behavior across genes. Th-cell differentiation is a dynamic process in which the regulatory factors are modulated over time, but longitudinal studies of Th-cell differentiation are sparse. Here, we present a dynamic transcriptome analysis following Th-cell differentiation into Th1, Th2, and Th1/2 hybrid cells at 3-h time intervals in the first hours after stimulation. We identified an early bifurcation point in gene expression programs, and we found that only a minority of ~20% of Th cell-specific genes showed mixed effects from both Th1 and Th2 cells on Th1/2 hybrid cells. While most genes followed either Th1- or Th2-cell gene expression, another fraction of ~20% of genes followed a Th1 and Th2 cell-independent transcriptional program associated with the transcription factors STAT1 and STAT4. Overall, our results emphasize the key role of high-resolution longitudinal data for the characterization of cellular phenotypes.

Keywords: T helper cell; cell differentiation; lineage commitment; regression analysis; time-course transcriptomics.

Figure 1

A high-resolution time course of Th-cell differentiation. (A) Experimental setup. Th-cell subsets were induced by polarizing signals in vitro, and gene expression profiles were obtained at 10 time points between 0 and 120 h after activation. (B) Flow-cytometric characterization of Th-cell subsets at day 5 after activation with polarizing conditions as described in (A). Normalized geometric mean indices for T-bet and GATA-3 expression are shown. Geometric mean intensities for IFN-γ- and IL-4-positive cells are indicated in bold. (C) Gene-expression profiles of four groups of genes (top to bottom): Th1-related, Th2-related, Tfh, and Th17-related, and other important Th cell-related genes. (D) Kinetics of master regulator transcription factors and signature cytokines for individual CD4+ T-cell subsets. Shown are normalized expression intensities as fold change relative to the first measured timepoint (0 h). (E) Principal component (PC) analysis of the differentiation time course. Cell subsets are indicated by marker shape. Time of measurement is indicated by color. (F) Numbers and overlap of kinetic genes between cell subsets. (G) Evolution of PC1 over time. Shown is the difference of PC values with respect to the Th0 condition. Genes with high correlation between subsets were removed (bottom) or kept for comparison (top).

Figure 2

Early Th-cell differentiation features three major patterns of kinetic gene expression. (A) Expression heatmap for kinetic genes. (B) Normalized expression kinetics of the three identified kinetic gene expression clusters, shown as averages over all genes and all cell types contained in each cluster. (C) Quantification of the numbers of identified kinetic genes across cell types within each kinetic cluster. (D) Gene classification as non-kinetic or kinetic including cluster association, for the four groups of Th cell-related genes introduced in Figure 1C . (E) Pathway enrichment analysis for genes uniquely assigned to kinetic clusters C1–C3. Pathways were pooled from REACTOME and Msigdb:Hallmark data bases; for a list of all enriched pathways see Table S1 .

Figure 3

A refined selection procedure identifies quantitative and qualitative differences in kinetic gene expression. (A) Workflow illustration. We employed a combination of regression fitting in maSigPro to derive quantitative differentially expressed genes (DEG), followed by a correlation filter to identify qualitative DEG and by an analysis of switching of kinetic clusters between cell types. (B) Correlation volcano plots based on the workflow in (A). Genes are categorized as kinetic (gray), quantitative DEG (black), or qualitative DEG (red). See Methods for details. (C) Numbers of qualitative and quantitative DEG obtained for each comparison of cell types. Brackets indicate the numbers of kinetic cluster switches. (D) Expression heatmap for all qualitative DEG exhibiting kinetic cluster switches in at least one comparison of cell types, as indicated on the left. (E) Venn diagrams of quantitative (left) and qualitative (right) DEG shared between Th1 or Th2 cells and Th1/2 hybrid cells. (F) Pathway enrichment analysis of DEG for all pairwise comparisons between cell types. Pathways were pooled from REACTOME, GO:BP, Msigdb:C2:WikiPathways, and Msigdb:C3:TFT databases. Shown are pathways with significant enrichment in at least two comparisons.

Figure 4

Superposition and independence of genes in Th1/2 hybrid cells. (A) Workflow sketch. Based on the identified qualitative Th1 vs. Th2 DEG (see Figure 3C ), similarity of genes to the expression profile of Th1/2 hybrid cells was assessed by a linear regression model. (B) Time courses of representative genes, and (C) quantification of gene classification into the four different categories. (D) Enrichment analysis using published transcription factor target gene sets (see text). Left: Analysis of upregulated genes in Th1 (Th2) cells, which are taken as Th1 vs. Th2 DEG with expression values higher (lower) in Th1 compared to Th2 cells. Right: Analysis of the Th1/2 hybrid transcriptional profile along the gene categories obtained in (A–C).