Construction and validation of a regulatory network for pluripotency and self-renewal of mouse embryonic stem cells

Abstract

A 30-node signed and directed network responsible for self-renewal and pluripotency of mouse embryonic stem cells (mESCs) was extracted from several ChIP-Seq and knockdown followed by expression prior studies. The underlying regulatory logic among network components was then learned using the initial network topology and single cell gene expression measurements from mESCs cultured in serum/LIF or serum-free 2i/LIF conditions. Comparing the learned network regulatory logic derived from cells cultured in serum/LIF vs. 2i/LIF revealed differential roles for Nanog, Oct4/Pou5f1, Sox2, Esrrb and Tcf3. Overall, gene expression in the serum/LIF condition was more variable than in the 2i/LIF but mostly consistent across the two conditions. Expression levels for most genes in single cells were bimodal across the entire population and this motivated a Boolean modeling approach. In silico predictions derived from removal of nodes from the Boolean dynamical model were validated with experimental single and combinatorial RNA interference (RNAi) knockdowns of selected network components. Quantitative post-RNAi expression level measurements of remaining network components showed good agreement with the in silico predictions. Computational removal of nodes from the Boolean network model was also used to predict lineage specification outcomes. In summary, data integration, modeling, and targeted experiments were used to improve our understanding of the regulatory topology that controls mESC fate decisions as well as to develop robust directed lineage specification protocols.

Figure 1. Signed and directed network extracted from ChIP-seq and knockdown or over-expression followed by genome-wide expression.

Edges are established where there is evidence for transcription factor binding to the gene proximal region and also change in expression after knockdown or over-expression. The 15 pluripotency nodes are color coded in light blue, and the lineage markers are color coded for the four major early differentiation lineages. Red diamond-heads denote inhibition and black arrowheads denote activation. Dashed lines connect pluripotency regulators to lineage markers and solid arrows connect pluripotency regulators to other pluripotency regulators.

Figure 2. 30 mRNA gene expression measurements in 96 single mESCs cell using the Fluidgm single cell microfluidic device.

Histograms of gene expression Ct-values from RT-PCR data collected from single mESCs in +serum/LIF (serum/LIF, blue) and -serum/2i/LIF (2i/LIF, red). Ct-values were normalized using the housekeeping gene Gapdh. x-axis represents -ΔCt values and y-axis represents the percentage of total cells. Alpha blending is used to show the results from the other condition for clearer visual comparison.

Figure 3. (A-B) Plots of normalized Ct-values partitioned into two clusters.

(C–D) Hierarchical clustering of gene expression in single mESCs using normalized Ct-values from RT-PCR data measuring levels of 15 pluripotency genes (orange) and 15 lineage-specific genes (purple) in 96 cells cultured in serum/LIF (C) and in 2i/LIF (D). The x-axis represents individual cells and the y-axis represents genes. Yellow represents ‘high’ and blue represents ‘low’ levels of gene expression. (E–F) Association of 30 genes based on gene expression in single mESCs cultured in serum/LIF or 2i/LIF. Numeric values in the color bars represent the distance score calculated as 1 – Pearson-Correlation Coefficient (PCC). Average-linkage, Euclidean-distance-based hierarchical clustering was performed the gene expression data.

Figure 4. Learned Boolean networks.

(A–B) Networks with learned Boolean logic transition functions consisting of 30 genes/proteins. Learning was achieved using the serum/LIF (A) or 2i/LIF (B) single cell data. Light cyan nodes represent genes; gray squares represent learned regulatory logic transition functions. The shadowed inset box exemplifies one such learned transition function upstream of Gata6. The sign represents ‘NOT’, ∧ represents ‘AND’ and ^∨ represents ‘OR’ logic operators. Boolean functions for all genes are available in Supplementary Information. Links from upstream parent nodes appearing in more than 90% of equally well-fitted Boolean functions are colored in green (activation) and red (repression). Other links are not shown. Novel links that resulted from the learning process are highlighted in dark red. (C) Overlap between the networks learned using data from mESCs cultured in serum/LIF or 2i/LIF. Links that are shared in the two conditions are in gray. Non-shared links are solid or dashed based on their source (solid for serum/LIF and dashed for 2i/LIF).

Figure 5. Comparison between computational in silico and experimental knockdowns followed by expression measurements.

(A) Design of in silico and experimental knockdowns. (B) Simulation results of computational knockdowns. Red represents up-regulation, blue represents down-regulation and white represents no change with regards to the unperturbed condition. (C) Results of experimental knockdowns followed by mRNA measurements of all network nodes, where each experiment was repeated twice. Colorbar illustrates log-scaled fold-change. Red represents up-regulation and blue represents down-regulation with regard to the unperturbed condition. (D) Comparison of computational and experimental knockdowns (see Text S1). Colors correspond to discordance score defined in the objective function in the supplementary material. Concordant results are colored in yellow while discordant results are colored in blue with gradients representing the degree of (dis)agreeability.

Figure 6. Computational in silico knockdowns of all possible single and double perturbations linked to predicted lineage differentiation outcomes.

(A, C) Predicted lineage specification based on dynamical model simulations. (B, D) Lineage predictions based on direct effects of knockdowns. Predictions are based on differentially expressed genes from LOF studies and promoter binding of transcription factors to differentially expressed genes based on mESC ChIP-chip/seq studies from the ESCAPE database.