An integrative computational systems biology approach identifies differentially regulated dynamic transcriptome signatures which drive the initiation of human T helper cell differentiation

T cells are key regulators of the adaptive immune system and have a central role in defense against pathogens and cancer as well as protection from autoimmune diseases. CD4⁺ T lymphocytes can differentiate to functionally distinct effector subtypes, including T helper 1 (Th1), T helper 2 (Th2) and more recently described T helper 17 (Th17) cells [1]. Th1 cells secrete effector cytokine IFN-γ and regulate cell-mediated immunity and play a role in the pathogenesis of autoimmune diseases, such as multiple sclerosis. Th2 cells in turn produce IL-4, IL-5, and IL-13 cytokines, and mediate immunity against extracellular pathogens and allergic reactions. Th17 cells, characterized by the production of a proinflammatory cytokine IL-17, regulate inflammatory responses on the mucosal surfaces. For the overall health in humans and animals, the proper balance between different effector T cell types and T regulatory cells is crucial [2, 3]. Aberrant activation of Th1 and Th17, or Th2 cells can trigger inflammatory autoimmune diseases as well as asthma and allergy. Previous studies utilizing genome-wide expression data and computational modeling have aimed at revealing the master regulators and regulatory networks within the differentiating Th1 and Th2 cells [4–9]. However, studies in human have been less extensive than in mouse due to the difficulty in collecting sufficient amount of samples to comprehensively profile T cell differentiation over time. In addition, lack of appropriate computational methods suitable for analyzing large-scale experimental data from multiple lineages over several time points spanning the lineage commitment process has limited the progress on revealing dynamics and molecular mechanisms underlying multiple lineage commitment.

A number of different time-series analysis approaches have been proposed to solve large-scale lineage commitment analysis problems. The general purpose F-test [10] can be used to test the difference between time-series data sets, but it does not extend to simultaneous comparison of multiple lineages and fails to take into account the correlation between the measurements at different time points. More recent approaches to analyze time-series data, including regression, differential expression, discriminant and clustering methods, are reviewed by Coffey and Hinde [11]. Methods for differential expression analysis include e.g. spline-based methods, generalized F-tests and hierarchical error and empirical Bayes models. Spline-based EDGE method by Storey et al. [12] is relevant for our problem because it provides comparisons for multiple conditions (lineage commitment profiles). Although EDGE computes a p-value for differential expression, it does not quantify the differential expression for all lineage comparisons, such as reciprocal genes (i.e., all lineages behave differently). ANOVA-based TANOVA method is based on the approach where different ANOVA structures are defined and the optimal one is found by evaluating the effects and significancies of the factors [13]. Recently, Stegle et al.[14] proposed an approach based on Gaussian processes (GP) to determine the time interval when a gene is differentially expressed. The methodology of Stegle et al. (2010) was limited to analyzing only two conditions. Moreover, it is often observed at transcriptional level that immediately after a treatment, such as activation of T cells by engagement of T cell receptor and CD28, genes are highly dynamic for some time but activity of gene expression decreases at later time points [15, 16]. Thus, an ideal computational method − that does not exist at the moment − should take into account the temporal correlation, handle a non-uniform measurement grid, cope with non-stationary processes, and be able to do a well-defined analysis of multiple conditions.

Here we developed a computational methodology, LIGAP (Li neage commitment using Ga ussian p rocesses) which analyzes experimental data from any number of lineage commitment time-course profiles and analyzed genome-wide gene expression profiles of human umbilical cord blood T helper cells (Thp) activated through their CD3 and CD28 receptors and cultured in absence (Th0) or presence of cytokines promoting Th1 or Th2 differentiation. The results give insight into differences of the three lineages in the expression landscape and provide marker genes for lineage commitment identification. Key lineage specific, that is, differentially regulated, genes discovered computationally were validated either experimentally at protein level or based on the published literature. Using a module-based analysis, we identified known and putative regulatory control mechanisms by overlaying highly coherent lineage profile clusters with genome-wide transcription factor (TF) binding predictions and pathway information. Consistent with the previously published results on IL-4/STAT6-mediated control of a large fraction of genes in Th2 program [17], our analysis revealed a comparable up-regulated and down-regulated modules, which are suggested to be controlled by STAT6 and other TFs. Interestingly, we also found that the genes which behave differently between all the lineages studied exhibit a consistent characteristic pattern, i.e., they are up-regulated in Th1 polarizing cells, down-regulated in Th2 polarizing cells, and in activated cells (Th0) the expression levels are between Th1 and Th2 cells. In addition, our analysis revealed a large set of novel genes, which are specific for different T cell subsets in human. All the gene expression data and differentially regulated genes as well as software implementing our computational analysis are made publicly available.

Identification of the key T helper cell regulators provides possible targets for modulation of immune response. To reveal T cell subset specific genes and their often subtle differences in expression, we developed a novel computational method, LIGAP. Traditional ways of identifying differentially expressed genes, such as the t-test, are problematic in studying time-series data since there is a need to carry out hypothesis tests on individual time points. On the other hand, commonly used statistical tests for whole time-course, including e.g. F-test, do not account for the inherent correlation between measurement time points. LIGAP overcomes many problems that have previously prevented quantitative comparisons of multiple differentiation profiles, with or without replicates. Among several beneficial features, LIGAP models correlation between time points and can cope with non-stationarities and non-uniform measurement grid. Other methods, such as EDGE, uses splines to estimate smooth time-course profiles but does not quantify the differential expression for all lineage comparisons. TANOVA uses standard regression framework and lacks explicit correlation structure between time points. Our study highlights the validity of the method by identifying known and novel differentially regulated genes and their kinetic differences during T helper cell differentiation. In addition, the non-parametric computational analysis automatically provides informative illustrations of time-course profiles together with associated uncertainty.

LIGAP calculated Th0 specific gene set contains only 18 genes and Th1 specific 49 genes compared to 466 genes that are specific to Th2 conditions. Activation of Thp cells through TCR and CD28 results in induction of IFNγ, which in turn leads to activation of Th1 signature genes. Addition of IL-12, however, results in enhanced induction of these genes and Th1 programming. Consistent with our previous results genes differentially regulated in response to Th1 programming are much more limited than those detected in response to initiation of Th2 response [16, 94].

Most of the Th1 specific genes encode well-known Th1 signature molecules. However, also genes new in this context were discovered. Interestingly, we identified RORC as one of the Th1 specific genes. Up-regulation of RORC in Th1 cells and existence of Th17/Th1 cells, however, remain conflicting as the master regulator of Th1 differentiation, T-bet, is known to inhibit transcription of RORC through RUNX1 [95], and expression of IL12Rβ2 is down-regulated by IL-17 [96]. It has been suggested that the high concentration of TGFβ required for in vitro Th17 polarization would inhibit IFNγ production [97], hence, it remains an open question whether some conditions would drive the differentiation of IL-17 and IFNγ producing cells from same naïve precursor T cell. Notably, ex vivo Th17 cells could be induced to develop further into Th17/Th1 cells by the combined actions of IFNγ and IL-12, and such conditions resulted in permissive chromatin remodeling at the IL12RB2 locus and loss of repressive histone modification at the TBX21 locus [29, 98].

As an example of previously uncharacterized differentially regulated genes, we validated the expression of Th2-associated phosphatases DUSP6 and PPP1R14A on protein level. PPP1R14A was shown in human pancreatic and melanoma tumor cell lines to positively regulate Ras/MAPK signaling [99], which are also involved in IL-4 induced signaling cascades. In T cells, the ERKs are activated though TCR stimulation and a TCR-mediated activation of Ras/MAPK signaling is required in differentiating murine Th2 but not in Th1 cells [100]. Furthermore, the Ras/MAPK cascade was shown to enhance the stability of GATA3 protein [101] as well as STAT6 independent CD3 and CD28 induced initial IL4 production [102]. DUSP6 on other hand is known to negatively regulate members of the mitogen-activated protein (MAP) kinase superfamily associated with cellular proliferation and differentiation [103]. More specifically, DUSP6 expression was shown to be induced by ERK1/2 signaling in differentiating mouse embryonic cell line and in human retinal pigment epithelial cells [104, 105] and it was hypothesized that DUSP6 is an essential part of a negative feedback loop of ERK1/2 signaling [106]. However, the T cell associated functions of both PPP1R14A and DUSP6 are completely unknown. Therefore, their significance in the signaling cascades of differentiating Th2 cells remains a highly interesting area of future research.

SPINT2 was recently identified as a direct STAT6 target in differentiating human Th2 cells [17] and in this study we are the first to show that SPINT2 is upregulated in Th2 cells at protein level as compared to other Th cell subsets. We found SPINT2 to be specifically expressed on Th2 cell surface as well as secreted into the culture medium, suggesting presence of a multiple transcripts of which some may lack the anchoring transmembrane domain [47]. Human SPINT2 (HAI-2) is a physiological inhibitor of matrix cleaving proteases and decreased expression of SPINT2 has been linked to progression of several cancers [107–109]. Up-regulated expression of extracellular proteases is crucial for pro-cancerous pathways as this enables efficient remodeling of the extracellular matrix as well as cleavage and activation of growth factors and their receptors. Interestingly, a truncated and secreted SPINT2 may act as an inhibitor for the activator of hepatocyte growth factor (HGF) and HGF is prominently expressed in lung tissue and is linked to reduced expression of Th2 cytokines and TGFβ, reduction of allergic airway inflammation, airway hyperresponsiveness and remodeling as well as reduced recruitment of eosinophils to the site of allergic inflammation in vivo[110, 111]. This suggests that SPINT2 might enhance Th2 response in allergic airway inflammation by inhibiting HGF signaling.

The LIGAP method elegantly identified the reciprocally regulated genes within the Th0, Th1 and Th2 conditions. Essentially, the list included genes encoding the hallmark Th1 specific transcription factor T-bet and cytokine IFNγ as well as the transmembrane receptor for IL-12. This list also included few cytoskeleton associated proteins, such as dystrophin (DMD), and palladin (PALLD), of which there is no current knowledge for their function in differentiating T helper cells. The observation suggests differences in cellular structures or putatively in the interaction of APC with the Th cell subsets as rearrangement of the cytoskeleton in T cells plays an important role in the organization of the immunological synapse (IS) and Th1 and Th2 cells are known to form morphologically distinct ISs [112, 113]. In addition to MAP3K8, molecules that participate in phosphorylation signaling cascades e.g. P2RY14, LPAR3, PPP1R14A, and PTPRO suggest their potential role for initiation or regulation of differentiation cascades. Importantly, the results presented here enable opportunities for further data mining and follow-up studies addressing the functions and importance of the novel Th subset specific genes.

The identification of STAT6 as the most significant TF regulating Th2 specific enhancement of transcription by the TF binding analysis is well in line with our previous STAT6 ChIP results [17]. Furthermore, the analysis between the predicted STAT6 target gene promoters and experimentally observed promoter associated binding sites showed statistically significant correlation. Interestingly, the overlapping STAT6 targets included INO80, which has been identifies as a part of a chromatin remodeling complex [114] and may hence, be involved in Th2 specific epigenetical regulation of Th cell differentiation. STAT6 specific regulation of Mannosyl (alpha-1,3-)-glycoprotein beta-1,2-N-acetylglucosaminyltransferase (MGAT1), a N-glycan-processing enzyme [115], may on one hand be involved in modifying the Th2 cell specific surface glycoprotein structures [116]. The overlapping target sites included also the promoter for SPINT2. The number of predicted STAT6 binding sites, however, was much larger than the experimentally observed binding sites, which may reflect the typically observed high false positive rate of computational binding predictions and the cell type specific state of chromatin as well as other competing factors affecting binding in vitro. The data created here also further suggests novel control mechanisms involving GATA3 regulated NKX3A as well as chromatin modification associated CDP. Only less than 10% of the Th2 down-regulated genes were reported to be direct targets of STAT6 by Elo et al., ([2010]) suggesting other major regulatory mechanisms play role among the IL-4 induced down-regulated genes. We found enrichment of IRF family and ISGF3 binding motifs in promoter regions of genes that are repressed in Th2 polarizing conditions, indicating that these TFs may play a significant role in the suppressing undesired gene expression in differentiating Th2 cells. Indeed, several IRF family members have been identified as differentially expressed during Th cell differentiation and necessary for both Th1 and Th2 polarization. As the IRF family proteins, excluding IRF1, share the same binding specificity model in TRANSFAC, the individual regulatory role for these factors is, however, difficult to postulate based on in silico TF binding site analysis.

The proposed LIGAP method can quantify a well-defined probabilistic specificity score for each gene and for each condition promoting a certain lineage commitment. In addition to grouping and ranking genes based on their dynamics, LIGAP summarizes all time-course measurements, together with the associated uncertainty, in an illustrative summary plot for visualization and manual assessment purposes. While here we have demonstrated the utility of LIGAP in analysis of gene expression dynamics, the LIGAP method is widely applicable to many types of datasets including quantitative time-course experiments and generalizes to any number of conditions.

Human CD4+ T cell purification and culturing. The human naïve umbilical cord blood CD4+ T cells were isolated as previously described [17]. Briefly, umbilical cord blood was collected from healthy neonates born in Turku University Hospital, Finland. Mononuclear cells were separated with Ficoll-Paque gradient centrifugation (#GEHE17-1440-3, Amersham Biosciences) and CD4+ T cells were then isolated with magnetic beads (Dynal CD4 Positive Isolation Kit, #113-31D, Invitrogen). After isolation the CD4+ cells were pooled to prepare cell cultures consisting cells from several neonates. The same pooled cells as utilized for Th0 (activated) and Th2 (activated and IL-4 stimulated) culture conditions by Elo et al. ([2010]) were used parallel for Th1 polarizing cultures. For activation, the cells were treated with plate-bound anti-CD3 (500 ng/24-well culture plate well, #IM1304, Immunotech) and soluble anti-CD28 (500 ng/ml, #IM1376, Immunotech) in density of 2-4 × 10⁶ cells/ml of Yssel’s medium (Iscove modified Dulbecco medium, #31980-048, Invitrogen) supplemented with Yssel medium concentrate [117], 1% human AB serum (#C11-011, PAA) and 100 U/ml Penicillin and 100 μg/ml Streptomycin (#P0781, Sigma) at 37°C in 5% CO2. For induction of Th1 cell polarization, IL-12 (2.5 ng/ml, # 219-IL, R&D Systems) was added to the cultures. At 48h after activation, IL-2 was added (17 ng/ml, #202-IL, R&D Systems) to all the cells and the polarizing conditions were maintained throughout the culture. The polarizing Th cells were harvested at time points 0, 12, 24, 48 hours in three replicates and at 72 hours in two replicates.

All the data included in this manuscript has been acquired under the permission from the Ethics Committee of the Hospital District of Southwest Finland approving the anonymous collection of cord blood samples after a parental consent, and the permission being in compliance with the Helsinki Declaration

Microarray studies. The preparation of samples for microarray detections was done as described in [17]. Essentially, total RNA (RNeasy Mini Kit, Qiagen) was extracted from the cultured cells and cRNA hybridized on Affymetrix GeneChip HG-U133 Plus 2.0 arrays (Affymetrix, Santa Clara, USA). All the microarray samples included in this study have been prepared at Finnish DNA Microarray Centre, Turku. The raw microarray data were processed using robust multi-array average normalization and log2-transformed in R (version 2.12.0) using the Bioconductor affy package (version 1.28.0).

Flow cytometry. The Th0, Th1 and Th2 condition cells at 24 hours were stained for SPINT2 expression studies. Purified anti-SPINT2 (8.7 μg/ml, #HPA 011101-100UL, Sigma-Aldrich) was used as primary antibody followed by staining with FITC-conjugated F(ab’)2 anti-rabbit IgG secondary antibody (1:1000 dilution, #11-4839-81, eBioscience). The stainings were analyzed with LSR II flow cytometer (BD Biosciences) and Flowing Software (http://www.flowingsoftware.com).

ELISA. The cell culture supernatants (at 48 hours) from Th0, Th1 and Th2 conditions were assayed for SPINT2/HAI-2 secretion by ELISA (# DY1106, R&D) according to the manufacturer instructions.

LIGAP. We construct our model-based lineage commitment comparison and visualization methodology, called LIGAP, using non-parametric GP regression similar to that in [14], extend the methodology to any number of conditions and propose to use a non-stationary neural network (NN) covariance function k(x^p,x^q) = σ*asin(x^p '*diag(l^-2)*x^q / sqrt[(1+x^p'*diag(l^-2)*x^p)*(1+x^q'*diag(l^-2)*x^q)]). The vectors x^p and x^q are augmented by an extra bias unit value entry and the parameter l defines the length-scale and σ controls the signal variance [118]. A non-stationary covariance function is chosen because often after cell activation or other stimulation the effects on temporal behavior of gene expression are very active and dynamic right after the stimulation but they mellow down over time and, thus, the observed behavior is non-stationary. For each gene at a time, LIGAP makes all comparisons between different cell subsets over the whole time-course data sets. In our application, the multiple hypotheses H_j are defined by the different partitions of the cell lineages. For example, if there are only two different lineages, then there are two different partitions (or hypothesis): H₁ denotes that lineages are similar and H₂ denotes that lineages are different. In our application consisting of three lineages, Th0, Th1 and Th2, we have 5 alternative hypotheses; (i) “Th0, Th1, Th2 time-course profiles are all similar” (hypothesis H₁), (ii) “Th0 and Th1 are similar and Th2 is different” (hypothesis H₂), (iii) “Th0 and Th2 are similar and Th1 is different” (hypothesis H₃), (iv) “Th1 and Th2 are similar and Th0 is different” (hypothesis H₄), and (v) “Th0, Th1, and Th2 are different from each other” (hypothesis H₅). LIGAP comparisons and quantifications are illustrated in Figure 1. In general, the total number of different partitions of N lineages is known in literature as the Bell number B_n (e.g., B₁ = 1, B₂ =2, B₃ = 5, B₄ = 15, etc.) [119].

Bayes factor is commonly used to see the evidence of the two alternative hypotheses; differentially expressed or not within a given time interval. To extend this to multiple lineages, we use the marginal likelihood p(D_i | H_j) to define the posterior probabilities of the different hypotheses H_j. For each of the hypothesis H_j, the data D_i for the i^th gene is split according to the partitioning. For example, for our application containing three lineages, hypothesis H₁ corresponds to grouping data from all lineages, hypothesis H₂ corresponds to splitting the data so that Th0 and Th1 time-course profiles are grouped together and time-course profiles from Th2 forms its own subset of data, hypothesis H₃ corresponds to splitting the data so that Th0 and Th2 time-course profiles are grouped together and Th1 forms its own subset of data, etc.

For each hypothesis, non-parametric regression is carried out separately for each subset of the data. For example, for the hypothesis H₃ we fit a GP to the combination of Th0 and Th2 time-course profiles and another GP to the Th1 time-course profiles. Following the standard GP regression methodology [118], the marginalization is done over the latent regression function and the hyperparameters are estimated using type II maximum likelihood estimation with a conjugate gradient based optimization algorithm initiated with ten randomly chosen hyperparameter values. Under the assumption of Gaussian likelihood and noise, the marginal likelihood can be written out analytically, and thus its value can be easily evaluated [118]. The marginal likelihood of a certain hypothesis (i.e., partitioning) is the product of the marginal likelihood of the separate subsets. The key idea behind the modeling is to find the marginal likelihood of the data under different hypotheses and thus have a probabilistic score to objectively compare different hypotheses.

Using the Bayes’ theorem and assuming unbiased, equal prior probabilities for different hypotheses (i.e., P(H_k) = P(H_l) for all k and l), we can write the posterior probabilities for the i^th gene as P(H_j | D_i) = P(D_i | H_j)P(H_j)/C, where C = Σ_j P(D_i | H_j)P(H_j) is a normalizing constant. Finally, these quantities can be combined to quantify the score of differential regulation for each gene. For example, the probability of the i^th gene being differentially regulated in Th2 lineage can be quantified as P(“Gene is differentially regulated in Th2” | D_i) = P(H₂ | D_i) + P(H₅ | D_i) .

ProbTF . ProbTF [76] method is used to make TF binding predictions on promoters of all RefSeq genes. Sequence specificities of TFs are taken from the TRANSFAC database [77] version 2009.3. All non-redundant PSFMs associated to human were taken, totaling 248 matrices. Promoters are defined as the [−1000,500] bp region around TSS. To assess statistical significance, we construct a TF specific null distribution by randomly sampling 50000 genomic locations of size 1501 nucleotides, against which the p-values of TF binding are computed.

Hierarchical clustering. The hierarchical clustering in Figure 5 was done using complete linkage and Euclidean distance metric.

Data access. The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [120] and are accessible through GEO Series accession number GSE 32959 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32959).