In silico discovery of small molecules for efficient stem cell differentiation into definitive endoderm

Abstract

Improving methods for human embryonic stem cell differentiation represents a challenge in modern regenerative medicine research. Using drug repurposing approaches, we discover small molecules that regulate the formation of definitive endoderm. Among them are inhibitors of known processes involved in endoderm differentiation (mTOR, PI3K, and JNK pathways) and a new compound, with an unknown mechanism of action, capable of inducing endoderm formation in the absence of growth factors in the media. Optimization of the classical protocol by inclusion of this compound achieves the same differentiation efficiency with a 90% cost reduction. The presented in silico procedure for candidate molecule selection has broad potential for improving stem cell differentiation protocols.

Keywords: bioinformatics; definitive endoderm; differentiation; drug repurposing; growth factor; pathway analysis; stem cell; transcription.

Graphical abstract

Figure 1

Analysis of public DE differentiation datasets (A) The expression trend of key genes involved in pluripotency and mesoderm/endoderm development; shows consistency among public datasets. (B) Heatmaps with activation values of pathways and transcription factors along differentiation timeline respectively; clustering is based on GEO: GSE75748 values.

Figure 2

Single-cell RNA-seq analysis (A) t-Stochastic neighbor embedding (tSNE) plot of the analyzed dataset; major cell clusters are labeled with the corresponding cell type names. (B) Distribution of the expression values of the marker genes, which were used to annotate the clusters. (C) Pseudotime ordering of the cells shows the expected bifurcation point with DE and ME branches. (D) Expression trends of the key mesoderm markers along the differentiation axis. The absence of the “spike” that is present in bulk data (A) can be explained by low depth and resolution of the single-cell data. (E) The regulon specificity score heatmap of the most informative TFs for each of the cell types. (F) The distribution of GSEA scores per single cell for the key developmental pathways. The asterisks (^∗∗∗∗) highlight the significance of the change with respect to hESCs (Wilcoxon test).

Figure 3

Drug repurposing pipeline (A) Different sources of input for drug repurposing; all are based on the RNA-seq results of the public datasets. (B) ssCMAP result for endoderm transcriptomic profile (intersection of GEO: GSE109658 and GSE75748). (C) Distribution of top 40,000 compound signatures ranked by GSEA enrichment score. The threshold of 0.35 is highlighted with a dashed line. (D and E) Spearman (D) and Pearson (E) correlation of LINCS signatures with transcriptomic profiles from the public datasets (GEO: GSE75748 - x axis; GSE109658 - y-axis); correlation threshold of 0.15 for both datasets is highlighted with transparent rectangles. (F) Pipeline for the identification of potential TGF-β/SMAD inducers. MCF-7 drug profiles are filtered based on dose and duration of treatment. The filtering is followed by GSEA with Hallmark and KEGG pathways and TF targets from RegNetwork. (G) Overlap between three different approaches.

Figure 4

SOX17-NG system for molecule screening (A) Schematic representation of mNeonGreen cassette used in this study. (B) Schematic representation of hESC differentiation toward DE with AA protocol. (C) mNeonGreen-2A-SOX17 knockin H1 (WA01) cells (homozygous mutant) differentiated by planar culture for 3 days. (D) Suspension culture with standard protocol. Days 2 and 3. (E) DE-induction efficiency oriented by AA. AC-A-A, activin A and CHIR99021 on day 1 and AA on day 2 and 3; AC-N-N, AA and CHIR on day 1 and no morphogens on day 2 and 3, etc. (n = 5 independent biological replicates). Data are presented as the mean ± SEM. (F) DE-induction efficiency by small molecules. Every treatment includes CHIR on day 1, and AA is replaced by each molecule described on y axis. LoAA, treatment with CHIR on day 1 and AA by 1:10 concentration for 3 days; LoAA w/BRD-K42644990, treatment with CHIR on day 1 and BRD and AA simultaneously by 1:10 concentration for 3 days (n = 3 independent biological replicates). Data are presented as the mean ± SEM. (G) DE-induction efficiency by small molecules, MCF7 LINCS profile of those not enriched for TGF-β or SMAD targets. Every treatment includes CHIR on day 1, and AA is replaced by each molecule described on y axis (n = 3 independent biological replicates). Data are presented as the mean ± SEM. (H) DE-induction efficiency by small molecules that are listed by all three drug repurposing approaches (n = 3 independent biological replicates). ^∗p < 0.05 (compared with C-N-N in F and compared with C-N-N in E). Data are presented as the mean ± SEM.

Figure 5

Characterization of DE cells induced by small molecules (A) FACS data of DE cells differentiated from ESCs for 72 h by each treatment. DE markers (CXCR4 and mNeonGreen reporter for SOX17) are analyzed. (B) FACS data of induced DE cells for surface markers CD177/CD275 indicating DE heterogeneity. (C) Quantification analysis of induced DE cells for mNeonGreen, CXCR4, CD177, and CD275. ESC, non-treated ESCs as negative control; AA, activin A and CHIR99021 on day 1 and AA on days 2 and 3 as standard treatment; JNJ, JNJ and CHIR on day 1 and JNJ-only treatment on days 2 and 3; BRD, BRD and CHIR on day 1 and BRD-only treatment on days 2 and 3; BRD+1:10AA, treatment with CHIR on day 1 and BRD and AA simultaneously by 1:10 concentration for 3 days. ^∗p < 0.05; n = 4 independent biological replicates. Data are presented as the mean ± SEM.

Figure 6

Characterization of pancreatic progenitor cells induced by small molecules (A) mRNA expression levels of a stem cell marker (Nanog) and pancreatic progenitor markers (PDX1, NKX6-1, and NEUROD1) for pancreatic progenitor cells differentiated by small molecules for 15 days are quantified by nCounter assay (Nanostring) (top panels). Each mRNA expression level is adjusted by the ratio of mNeonGreen cells in all cells on differentiation day 3 (bottom panels). ^∗p < 0.05 compared with the other three groups; n = 3 independent biological replicates. Data are presented as the mean ± SEM. (B) PDX1 protein+ cells are quantified for day 15 pancreatic progenitor (PP) cells by FACS. (C) Immunostaining for PDX1 and NKX6-1 proteins on AA- and BRD-induced PP cells. Scale bar, 50 μm. (D) BRD− or BRD+ AA-induced DE can support differentiation of insulin-expressing β-like cells during a 22 day differentiation protocol as assessed by flow cytometry using an INS-2A-GFP knockin add-on hESC line (n = 6 independent biological replicates for AA, n = 3 independent biological replicates for the other samples). Data are presented as the mean ± SEM.

Figure 7

Effect of different compounds on hESC differentiation Heatmap of marker gene expression and pathway activities (left). Scaled expression of the marker genes for different differentiation protocols. (Right) Single-sample gene set enrichment scores of selected Hallmark pathways.