Matched miRNA and mRNA signatures from an hESC-based in vitro model of pancreatic differentiation reveal novel regulatory interactions

Type I diabetes, also known as juvenile diabetes, is a chronic condition characterized by autoimmune destruction of the insulin-producing beta islet cells in the pancreas. The death of these cells leads to insulin deficiency, hyperglycemia and dependence on insulin replacement therapy. Whole-pancreas transplantation or infusion of islet cells into the portal vein are effective treatments, but limited by the scarce supply of materials and the necessity for life-long immunosuppression (Guo and Hebrok, 2009). The generation of pancreatic beta-like cells from human pluripotent stem cells (hPSCs), which include embryonic stem cells (hESCs) and induced PSCs (hiPSCs), is a potential alternative source of insulin-producing cells for cell therapy (D'Amour et al., 2006). Other options may involve direct reprogramming of other somatic cell types to beta islet-like cells by forced expression of beta cell-specific transcription factors, which has been achieved in an in vivo mouse model (Zhou et al., 2008), but has not been successfully performed in vitro.

Using defined growth factors and small molecules, beta islet-like cells can be obtained from hESCs using a five-stage protocol: (1) definitive endoderm (DE); (2) primitive gut tube; (3) posterior foregut; (4) pancreatic progenitor; and (5) pancreatic endoderm (D'Amour et al., 2006; Kroon et al., 2008; Van Hoof et al., 2009). This method mimics the in vivo process of pancreatic development, and has been utilized in experiments designed to identify cell-surface markers for the isolation of pancreatic progenitor cells (Kelly et al., 2011). It is important to note that the hormone-producing cells generated in this manner express only low levels of insulin, glucagon, or both (D'Amour et al., 2006), which suggests that these cells are either immature or have followed an aberrant differentiation trajectory. However, transplantation of hESC-derived pancreatic endoderm cells into a mouse model of type I diabetes results in generation of beta-islet like cells that secreted high levels of insulin in a glucose-responsive manner (Kroon et al., 2008; Xie et al., 2013). These results suggest that the hESC-derived pancreatic endoderm population contains cells that can mature into fully functional beta cells, supporting the notion that in vitro differentiation of hESCs can be used as a model system to study early pancreatic differentiation. Furthermore, we postulate that comparisons made between hESC-derived pancreatic endoderm cells and bona fide beta islet cells from the human pancreas can be used to instruct the development of methods for in vitro maturation of hESC-derived pancreatic endoderm cells into functional beta cells.

MicroRNA (miRNA) repression of mRNAs is an important mechanism for regulation of expression during cell fate specification, apoptosis, and metabolism (Bartel, 2004). These small (∼22 nt) non-coding RNAs bind to partially complementary sequences on target mRNAs (most commonly in their 3′UTRs), and lead to post-transcriptional gene repression. miRNAs are classified in two ways: (1) by families, which share similar target recognition motifs, or seed sequences, located at positions 2–8 from the 5′-end of mature miRNAs (Bartel, 2009; Doench and Sharp, 2004; Lewis et al., 2003); and (2) by clusters, which are miRNAs encoded by sequences found in close proximity in the genome. High complementarity between miRNAs and cognate sequences in the target mRNAs, commonly seen in plants, generally leads to direct mRNA cleavage. In contrast, partial complementarity, the most common case in metazoans, results in alteration of mRNA stability through de-adenylation, translational repression, or other mechanisms (Bartel, 2004). Work from our lab and others has shown that hPSCs have a unique miRNA signature, which plays an important role in regulation of the pluripotent state, and that cell fate can be changed by altering the miRNA content of cells (Ivey et al., 2008; Lin et al., 2008). In fact, it has recently been shown that overexpression of a few miRNA(s) can directly reprogram fibroblasts into iPSCs (Liao et al., 2011), or help other transcription factors to directly reprogram fibroblasts into functional neurons (Ambasudhan et al., 2011) or cardiomyocytes (Srivastava and Ieda, 2012), indicating that miRNAs are powerful regulators of cell lineage specification.

As is the case in other lineages, miRNAs play an important role in pancreatic organogenesis and functional maturation. In a mouse study, conditional deletion of Dicer1, a critical enzyme for miRNA and siRNA biogenesis, from the pancreas led to gross defects in all pancreatic lineages, especially beta islet cells, with the endocrine defects associated with persistent overexpression of progenitor genes such as Neurog3 (Lynn et al., 2007). In other studies, it has been shown that miR-375, one of the most abundant miRNAs in the endocrine pancreas, regulates insulin secretion by regulating myotrophin expression (Poy et al., 2004), while miR-30d regulates the response of beta cells to glucose and inhibits the epithelial–mesenchymal transition (EMT) in human beta islet cells during prolonged in vitro culture (Joglekar et al., 2009; Ozcan, 2009). The miR-200 family, which inhibits EMT by directly targeting ZEB1 and ZEB2, two transcription factors that repress E-cadherin (CDH1) expression during EMT (Gregory et al., 2008), is abundantly expressed in both the mouse and human pancreas and has been implicated in pancreatic carcinogenesis (Brabletz et al., 2011). While these studies have revealed some of the roles that miRNAs play in pancreatic development and differentiation, the discrepancy between the catastrophic effect of pancreatic lineage-specific loss of Dicer1 and the more subtle effects of perturbations in miR-375, miR-200 and miR-30d expression suggests that many critical miRNA/target mRNA interactions regulating this process remain to be uncovered.

In this study, WA09 hESCs were differentiated into pancreatic endoderm cells following a published five-stage differentiation protocol (D'Amour et al., 2006). Differentiated cells at each stage were collected and compared to undifferentiated hESCs, fetal pancreatic tissue, and purified adult islets, which consist of ∼70% beta cells, ∼20% alpha cells, ∼10% delta cells, and small percentages of pancreatic polypeptide producing (PP) and epsilon cells (Brissova et al., 2005; Cabrera et al., 2006). These samples were subjected to genome-wide microarray-based miRNA and mRNA expression profiling, and the resulting data were bioinformatically analyzed. We focused on identifying transcripts and miRNAs with two expression patterns: those with differential expression among the five stages of hESC differentiation, and therefore are more likely to have roles in early pancreatic differentiation; and those that were differentially expressed between the hESC-derived pancreatic endoderm and adult islet samples, which are candidates for beta cell maturation factors. We found that cells in the early and late stages of hESC differentiation possessed distinct miRNA signatures. In addition, we used integrated miRNA and mRNA analyses to identify critical miRNA–mRNA interactions occurring during the differentiation process. qRT-PCR and dual luciferase experiments were used to validate functional roles for miR-200a in regulating both the EMT process and definitive endoderm specification during early stages of differentiation. In addition, we verified that miR-30d and let-7e, two miRNAs that are induced during differentiation, regulate the expression of RFX6, a pancreatic progenitor (Smith et al., 2010). Our findings suggest precise control of target mRNA expression by miRNA contributes to proper lineage specification during the process of pancreatic differentiation.

The WA09 cell line was cultured without feeder cells in Stempro (Invitrogen) medium on Geltrex (Invitrogen) and was differentiated into pancreatic beta islet-like cells using a modification of a published five-stage differentiation protocol (D'Amour et al., 2006) (Fig. 1A). Biological replicates were collected at each stage, with efficient differentiation confirmed by immunocytochemical staining and quantitative real-time PCR (qRT-PCR) with stage-specific differentiation markers (Fig. 1B–D).

Twenty-three samples were included in the miRNA expression analysis, including two to five replicates of: undifferentiated hESCs; hESCs differentiated toward beta islet-like cells (stages 1–5); fetal pancreatic tissue; and adult islets (supplementary material Table S1). To characterize changes in miRNA expression during differentiation, we first divided the 16 undifferentiated and differentiated hESC samples into six groups according to their differentiation status (undifferentiated and Stages 1–5) for ANOVA, followed by hierarchical clustering. The initial ANOVA resulted in a total of 220 differentially expressed miRNA probes (Pearson correlation coefficient as the distance measure, FDR<0.05; Fig. 2A,B). Bi-directional hierarchical clustering indicated that Stage 1 and 2 samples clustered together with, and shared some expression patterns with undifferentiated hESCs (see Fig. 2B, lower portion of miRNA Group 1 and also Group 4), while samples from Stages 3, 4 and 5 grouped together in a separate cluster (Fig. 2B). We thus considered Stages 1 and 2 to represent early stage differentiation, and Stages 3, 4 and 5 to represent late stage differentiation. The 220 probes were separated into four major groups by the hierarchical clustering analysis (Fig. 2B, labeled 1–4 on left side of figure): hESC-specific miRNAs that decreased in expression with differentiation (Group 1); miRNAs that were induced in the early stages of differentiation (Group 2); miRNAs enriched in both undifferentiated cells and late stage differentiated cells, but downregulated at early stages of differentiation (Group 3); and miRNAs that were induced only at late stages of differentiation (Group 4).

The 108 hESC-specific miRNA probes in Group 1 were further divided into two subgroups, which were rapidly (Fig. 2B, top) or gradually (Fig. 2B, bottom) downregulated upon differentiation. Inspection of the genomic positions of these miRNAs revealed that more than 40 miRNAs in the rapid downregulation group were located on chromosome (chr) 19q13.42 (Fig. 2C). This particular miRNA cluster is known as the chromosome 19 miRNA cluster (C19MC), which contains 54 miRNAs and spans 96 kb. miRNAs from this cluster have been demonstrated to be highly expressed in both hESCs and placenta (Laurent et al., 2008; Miura et al., 2010; Ren et al., 2009; Zhu et al., 2009). In addition, consistent with a previous report (Hinton et al., 2010), expression of the hESC-specific miR-371/372/373 family, which is located 100 kb distal to the C19MC on chr19q13, gradually decreased during differentiation (Fig. 2B,C). Another hESC-specific miRNA family, the miR-302 family located on chr4q25, showed an expression pattern similar to the miR-371/372/373 cluster, raising the possibility that these two miRNA families may be regulated by a common mechanism that differs from the regulation of C19MC (Fig. 2C). Interestingly, all three of these miRNA families are highly expressed in hESCs and share similar seed sequences, and hence are likely to regulate a common set of downstream mRNAs (Laurent et al., 2008). In addition, members of the miR-302 and miR-371/372/373 families (or their homologues in mouse) have been shown to play an important role in regulation of the ESC cell cycle (Wang et al., 2008a) and promote reprogramming of somatic cells to pluripotency (Judson et al., 2009). Taken in the context of these prior findings, our results suggest that miRNA regulation of differentiation may be more complex than previously suspected, and that the multiple pluripotency-associated miRNAs containing similar, but not identical, seed sequences may have different roles, which are influenced by the precise timing of downregulation during differentiation and subtle differences in target specificity.

Twenty miRNA probes (Group 2) were significantly induced at early stages of differentiation, Stage 1 in particular (Fig. 2B). Among these miRNAs, induction of miR-489 appeared to be endoderm specific, while miR-375 was immediately elevated in Stage 1 definitive endoderm cells and remained at a high level of expression throughout the entire differentiation process (Fig. 2B).

Thirty miRNA probes (Group 3), including miR-200a/429, miR-30a and miR-29a/b, had a more complex pattern of expression, and were highly expressed in undifferentiated cells, downregulated at early stages of differentiation, and then upregulated at later stages of differentiation (Fig. 2B).

Sixty-two miRNA probes (Group 4) showed increased expression specifically during late stages of differentiation, and contained several miRNA clusters than have been previously recognized for their roles in differentiation of other lineages (Fig. 2B–D). Among these were the let-7e/miR-99b/125a-5p cluster on chr19 (2MB upstream of the C19MC), which is intergenic and has been shown to be involved in the regulation of hematopoietic differentiation (Cabrera et al., 2006; Gerrits et al., 2012; Guo et al., 2010), the miR-23b/24/27b cluster on chr9, which is intronic and involved in hepatocellular differentiation (Rogler et al., 2009), and the miR-532/501/362/500/660 cluster on chrX, which is intronic and as yet has no known functions. We observed discordant expression patterns among members of several miRNA families, including the miR-30 family and let-7 family, with only miR-30c/d and let-7e/f specifically induced at late stages of differentiation (Fig. 2D). In addition, unlike miR-200a, which showed a transient dip in expression in the early stages of differentiation, miR-200c was upregulated at the beginning of differentiation and maintained its increased expression throughout the entire differentiation process, suggesting (as in the case with the C19MC, miR-302, and miR-371/372/373 miRNAs) that miRNAs with similar seed sequences (in this case miR-200c and miR-200a/429) are regulated by different mechanisms during pancreatic differentiation.

To confirm the dynamic changes of miRNA expression in differentiated samples, we performed quantitative RT-PCR of selected miRNAs in each group. Our data confirmed the expression patterns of miR-372 (Group 1), miR-375 (Group 2), miR-200a (Group 3), miR-30d and let-7e (Group 4) that correlates well with the microarray expression pattern (supplementary material Fig. S1).

To identify miRNAs induced during the process of differentiation that are likely to be important for pancreatic development, we studied the 92 miRNA probes that were highly expressed in the late stages of hESC differentiation (Groups 3 and 4) in the context of bona fide pancreas tissue, including fetal pancreas and adult beta islet samples. Hierarchical clustering showed that, among the 92 probes, 53 miRNA probes, including the miR-200a/c/429 (miR-200 family), miR-30b/c/d (the miR-30 family), let-7e/99b/125a and miR-23b/24/27b, were enriched in both the late stage differentiated derivatives and adult islet cells (Fig. 2D; Table 1), indicating that these miRNAs are likely to be bona fide beta islet cell-associated miRNAs that are functionally involved in the differentiation or maturation of human islet cells. The other 39 miRNA probes, including the miR-532/501/362/500/660 cluster on chrX, were not enriched in the human adult islet samples, but a large subset was highly expressed in the fetal pancreatic tissue samples, indicating that they may be transiently expressed during differentiation.

To identify gene expression changes involved in differentiation of human beta islet cells, genome-wide mRNA expression profiling was performed using Illumina HT-12v4 expression BeadChips on the same samples as the miRNA analysis. Using ANOVA, we identified 10,063 probes for genes differentially expressed among the differentiated hESC samples (FDR<0.05, supplementary material Fig. S2A). As seen in the miRNA analysis, bi-directional hierarchical clustering separated the hESC samples into two clusters, with the Stage 1 and 2 samples representing early stages of differentiation that clustered together with undifferentiated hESCs, while the Stage 3, 4, and 5 samples were clustered together and represented late stage differentiation. Clustering of the probes resulted in four distinct categories (supplementary material Fig. S2B), with 3601 probes upregulated during late stage differentiation. We then studied these 3601 probes in the context of the same fetal pancreas and adult beta islet samples as were used in the miRNA analysis, and identified 1524 mRNA probes that were highly expressed both in the late stage hESC differentiation and human pancreatic tissue samples (supplementary material Fig. S2A,C), suggesting that the corresponding transcripts may have functional roles in differentiation to pancreatic beta islet cells. Gene ontology (GO) analysis revealed that this set of 1524 genes is enriched for genes in several relevant biological processes, including the regulation of insulin secretion (supplementary material Fig. S2D).

Post-hoc ANOVA pair-wise comparisons were performed using Tukey's test to identify miRNAs and mRNAs that were uniquely upregulated and downregulated at each stage of differentiation of hESCs into beta islet-like cells (Tukey P-value <0.05 and expression fold-change >1.5; Fig. 3A). As expected, compared to the undifferentiated hESCs, the derivatives from Stage 5 had the highest number of significantly upregulated and downregulated miRNA and mRNA probes (Fig. 3B). To discover factors that regulate the passage from one stage to the next during islet cell differentiation, we also compared each differentiation stage to the following stage. As shown in Fig. 3B, there were 98 miRNA probes and 982 mRNA that were differentially expressed between undifferentiated hESCs and Stage 1 cells. The top 50 miRNAs and mRNAs are shown in Fig. 3C, and include mir-200a/429 and the known definitive endoderm-associated genes SOX17 and FOXA2. Comparison of Stage 2 to Stage 3 revealed upregulation of critical pancreatic progenitor genes, such as PDX1, ONECUT2 and RFX6 in the Stage 3 samples, consistent with a transition between early stage and late stage differentiation. In this way, we were able to identify miRNAs and mRNAs associated with each transition during pancreatic differentiation. Detailed information on all significantly differentially expressed miRNA and mRNA probes is provided in supplementary material Table S2 and Table S3.

To identify potential regulatory effects of miRNAs on gene expression at each stage of differentiation, we performed integrative analysis of the miRNA and mRNA expression data in GALAXY, using bioinformatically predicted miRNA–mRNA interactions from TargetScan 6.1 to link the two datasets. For each predicted miRNA–mRNA interaction, the miRNA and mRNA expression levels measured by the microarray analyses were used to calculate a correlation coefficient (supplementary material Table S4). We were most interested in two transitions: between the pluripotent and definitive endoderm stages; and between early and late differentiation. Given the interesting expression patterns of the miR-200 and miR-30 miRNA families as noted above, we focused on identifying the potential mRNA targets of miR-200a/429 in the undifferentiated to Stage 1 transition, and mRNA targets of miR-30d in the Stage 2 to Stage 3 transition. Depending on the overall intracellular milieu, which may contain other regulatory factors, and the specific mechanism of miRNA repression at work for a particular miRNA–mRNA pair, it is possible for miRNAs and their direct mRNA targets to show anti-correlated, correlated, or uncorrelated expression. We reasoned that anti-correlation might indicate a predominant effect of the miRNA on regulation of the cognate mRNA at a specific stage of differentiation, while correlation might indicate either an indirect effect or ‘priming’ of the levels of miRNA to enable rapid repression of the cognate mRNA at a later stage of differentiation. Therefore we chose some anti-correlated and some correlated miRNA–mRNA pairs for further study, and used luciferase reporter assays to distinguish between direct and indirect interactions.

It is well known that the EMT process is required for germ layer formation during early embryonic development (Arnold et al., 2008). To investigate the potential role of the miR-200 family, which has previously been described to regulate EMT in other systems (Gregory et al., 2008), in endoderm specification via EMT, we explored the expression of E-cadherin (CDH1), a cell adhesion protein that is expressed on epithelial cells, during differentiation from the pluripotent state to definitive endoderm, using expression of the endoderm-specific transcription factor SOX17 to verify proper differentiation to definitive endoderm. As shown in Fig. 4A, the expression of CDH1 was significantly downregulated in the definitive endoderm cells. Moreover, the expression of SOX17 and CDH1 in the cells was anti-correlated, suggesting that EMT may be required for SOX17 expression and endoderm-lineage specification (Fig. 4B). Interestingly, the expression of CDH1 rebounded in the Stage 2 samples, with mRNA levels (Fig. 4C) increasing more dramatically than protein levels (Fig. 4B). We then examined the expression of ZEB1 and ZEB2, two CDH1 repressors that were previously shown to be negatively regulated by the miR-200 family, in the early stage differentiation samples. qRT-PCR analysis revealed that the expression of ZEB1 increased at Stage 1 but then decreased at Stage 2, while ZEB2 was increased at Stage 2 (Fig. 4C), suggesting that the previously described regulatory system consisting of miR-200a/429, ZEB1/ZEB2 and CDH1 might be involved in the regulation of EMT during endoderm specification.

Since the 3′UTRs of the endoderm-associated transcription factors SOX17 and FOXA2 contained the miR-200a binding sequence and had expression levels that were anti-correlated with the miR-200a levels, we performed luciferase-based reporter assays (Fig. 5A) to functionally test the regulatory effect of miR-200a on SOX17 and FOXA2, and used a ZEB2 3′UTR construct as a positive control. Synthetic miR-200a mimics or inhibitors were transiently co-transfected into human neonatal fibroblasts or undifferentiated WA09 cells with reporter plasmids containing the Renilla luciferase cDNA linked to the 3′UTRs of ZEB2, FOXA2, SOX17 or a mutated SOX17 3′UTR (Fig. 5A). Constitutive Firefly luciferase expression, encoded on the same plasmid, was used to normalize transfection efficiency among samples. As shown in Fig. 5B, the normalized luciferase activities in the cells transfected with the SOX17 and ZEB2 3′UTR reporter plasmids were significantly decreased by treatment with miR-200a mimics and increased by treatment with miR-200a inhibitors, indicating direct repression of these transcripts by miR-200a. In contrast, the luciferase activities in cells transfected with the FOXA2 or the mutated SOX17 3′UTR reporter plasmid were not affected by miR-200a mimics or inhibitors (Fig. 5B). We further tested how the endogenous miRNA content of cells undergoing differentiation changed the normalized luciferase activities in cells transfected with the 3′UTRs of SOX17, FOXA2 and other pancreatic differentiation-related genes (supplementary material Fig. S3A). As shown in Fig. 5C, the Stage 1 cells containing the SOX17 3′UTR reporter plasmid showed an increase in normalized luciferase activity compared to the undifferentiated cells, with specificity demonstrated by significant blockage of this effect by co-transfection of a synthetic miR-200a mimic (Fig. 5D). In contrast, the mutated SOX17 3′UTR reporter plasmid showed a significantly lower increase in normalized luciferase activity compared to wild-type SOX17 3′UTR (Fig. 5C). The fact that the mutated SOX17 3′UTR reporter plasmid showed some increase in luciferase activity with differentiation suggests that there are other sequences (outside the miR-200a binding site) regulating SOX17 expression. This finding is consistent with the observation in Fig. 5D, showing that the addition of miR-200a mimic only partially blocks the induction of the wild-type SOX17 reporter with differentiation. Cells receiving the FOXA2 3′UTR reporter plasmid also showed increased normalized luciferase activity on days 2 and 3 of differentiation, but this increase could not be blocked by miR-200a treatment, further suggesting that miR-200a directly interacts with SOX17, but not FOXA2, and that there is likely a different miRNA or other repressive factor that interacts with the FOXA2 3′UTR (supplementary material Fig. S3B,C). Moreover, overexpression of miR-200a during early stage differentiation resulted in a population of CDH1+/SOX17− cells, which is not normally seen during the differentiation process (Fig. 5E). Real-time PCR and western blotting further confirmed the opposite regulation of SOX17 and CDH1 after miR-200a treatment (Fig. 5F,G). These results suggest that miR-200a is critically involved in the normal regulation of EMT and endoderm lineage commitment through direct repression of ZEB2 and SOX17.

To understand miRNA–mRNA interactions during the transition between early and late stage differentiation, we focused on miR-30d because it is one of the miRNAs highly expressed in both late stage differentiation and human islet cells (Fig. 2D; Table 1). Its rapid induction during the transition from Stage 2 to Stage 3 was paradoxically positively correlated with RFX6 expression, a pancreatic beta cell progenitor gene that has two predicted binding sites for miR-30d in its 3′UTR (Fig. 6A). RFX6 is also a predicted target of the let-7 family, and we demonstrated in our microarray profiling that let-7e and let-7f were significantly induced during late stage differentiation (Fig. 2B). To validate the predicted miR-30 and let-7 binding sites in the RFX6 3′UTR, we performed luciferase assays using the RFX6 3′UTR-containing reporter plasmid in undifferentiated WA09 cells, MIN6 cells and fibroblasts. Our results demonstrated robust and dose-dependent repression of the normalized luciferase activity in all cells transfected with the RFX6 3′UTR-containing reporter by the miR-30d and let-7e mimics, and significant induction of luciferase activity for the same reporters by the miR-30d inhibitor (data in MIN6 cells are shown in Fig. 6B,C). RT-PCR confirmed the inhibitory effects of miR-30d on endogenous RFX6 mRNA in both MIN6 and human islet cells (data in MIN6 cells are shown in Fig. 6D). We further demonstrated that treatment with miR-30d mimics at Stage 3 of differentiation muted the stage-specific induction of the RFX6 transcript normally seen during in vitro differentiation of hESCs to beta islet-like cells by 25% (Fig. 6E).

In this study, we differentiated WA09 hESCs into pancreatic endoderm cells in a staged manner as a model system to identify the miRNA and mRNA expression signatures specific to each stage of differentiation, and to identify the functional elements of these signatures by comparing them to data from fetal and adult pancreas samples.

Although previous studies have reported on the miRNA (Poy et al., 2004; Rosero et al., 2010) or mRNA (Webb et al., 2000) profiles of beta islet cells and fetal pancreas, or examined the effects of global miRNA deficiency (Lynn et al., 2007; Melkman-Zehavi et al., 2011) or specific miRNAs (Baroukh et al., 2007; El Ouaamari et al., 2008; Melkman-Zehavi et al., 2011) on pancreatic development and expression of beta-cell specific proteins, this report is the first to use matched miRNA and mRNA data from an in vitro model of beta cell differentiation in the context of bona fide fetal pancreas and adult islet samples. The analysis of these data enabled us to identify several patterns of miRNA expression in the differentiated hESCs, including loss of expression of pluripotency-associated miRNAs and gain of expression of some, but not all, of the pancreas-specific miRNAs. In terms of the pluripotency-associated miRNAs, we noted that they were downregulated in response to differentiation cues at variable rates. In contrast to the rapid downregulation of certain miRNAs (e.g. C19MC miRNA cluster), the expression of the miR-371/372/373 and miR-302 families showed a more gradual trajectory. The sustained expression of these pluripotency-associated families during early stage differentiation (definitive endoderm in particular) is consistent with findings from a previous study (Hinton et al., 2010), and suggests that endoderm formation may actually require persistent expression of certain pluripotency-associated miRNAs.

We observed that hESCs that had undergone differentiation to the later stages of the protocol (Stages 3, 4, and 5) acquired a portion of the fetal pancreas-specific signatures, as well as a portion of the signatures shared by the fetal pancreas and adult islet samples, suggesting that these differentiated hESCs were following a trajectory toward the beta islet cell lineage, but were at an immature stage. The identification of miRNAs that were highly expressed in late stage differentiated samples and fetal pancreas but absent in adult beta islet cells suggests that the differentiated hESC may be at an earlier stage than the cells present in the adult islet samples. These miRNAs may need to be further suppressed in driving islet cell precursors to maturity.

We reasoned that the set of miRNAs that were highly expressed in the late stage differentiated hESCs as well as the fetal pancreas and adult islet samples were likely enriched for miRNAs that are functionally important in beta islet cell differentiation. Supporting this hypothesis, several of these miRNAs have been implicated in biological processes of potential relevance to the development of endoderm cells, including miR-7, the miR-200 family, the miR-30 family, miR-29a/b, miR-199a/b, the miR-23b/24/27b cluster and the let-7e/99b/125a cluster. For example, miR-7 affects the differentiation of endocrine cells in the pancreas through inhibition of PAX6 (Kredo-Russo et al., 2012) and miR-200 has been implicated as a factor that inhibits pancreatic carcinogenesis by inhibiting EMT (Brabletz et al., 2011; Burk et al., 2008; Li et al., 2009; Soubani et al., 2012). In addition, the expression of miR-30d has been shown to help maintain the morphology of beta islet cells and their epithelial characteristics (Joglekar et al., 2009), although the mRNA targets through which it acts were not identified in this previous study. miR-29a/b regulates DNA methylation through repression of the de novo DNA methyltransferases, which have been shown to affect beta islet cell-specific gene silencing (Pullen et al., 2011) and to repress OCT4/POU5F1 and NANOG (Li et al., 2007), while miR-24/27b are involved in the differentiation of liver stem cells and lipid metabolism (Rogler et al., 2009; Vickers et al., 2013).

Among these miRNAs, the miR-200 family is of particular interest to us due to its complex pattern of expression, namely, rapid downregulation of miR-200a during early endoderm differentiation followed by increasing expression through the later stages of differentiation, and a gradual increase in expression of miR-200c throughout the process of differentiation. It is worth noting that although they originated from a common ancestral miRNA sequence, miR-200a and miR-200c have different seed sequences. The miR-200s, particularly miR-200a, regulate the EMT process through repression of ZEB1/ZEB2, which in turn repress CDH1 (Burk et al., 2008), and have been implicated in the development of pancreatic cancer (Yu et al., 2010). Our results have demonstrated for the first time that this regulatory system is intact and functional during hESC differentiation to endoderm, with downregulation of miR-200a in the first stage of differentiation to definitive endoderm being required both for proper downregulation of CDH1, and upregulation of miR-200a in the later stages being associated with upregulation of CDH1. A potential role for miR-200c as a repressor of SOX2 has been proposed (Wellner et al., 2009), and is supported by the presence of a predicted miR-200c binding site in the 3′UTR of SOX2 (TargetScan 6.1), but has not been validated. Integrative mRNA and miRNA expression analysis identified SOX17 and FOXA2 as two additional potential targets of miR-200a, which had both anti-correlated expression with miR-200a and contained miR-200a binding sites in their 3′UTRs. We experimentally demonstrated for the first time that SOX17, but not FOXA2, is directly repressed by miR-200a. Taken together, these results suggest that complex patterns of expression and targeting of the different miR-200 family members appears to promote endoderm lineage differentiation through coordinated repression of a pluripotency factor (SOX2), promotion of EMT, and de-repression of SOX17.

RFX6 contains predicted binding sites for miR-30d and let-7, but unexpectedly had an expression pattern that was positively correlated with those of these miRNAs. Since we were able to experimentally validate direct repression of RFX6 by miR-30d and let-7, we hypothesize that the inhibitory effects of miR-30d and let-7e were counteracted during hESC differentiation from Stage 2 to Stage 3 by other regulatory mechanisms. When cells were further differentiated into Stage 4 pancreatic precursors, the expression of RFX6 and the miRNAs were anti-correlated, mostly due to decreased RFX6 expression. This pattern of RFX6 expression was also observed in a mouse study, in which RFX6 was initially expressed in all pancreatic and gut endoderm cells and then became progressively restricted to the endocrine lineage by an unknown regulatory mechanism (Smith et al., 2010; Soyer et al., 2010). Our results, including the novel findings that RFX6 is a direct target of miR-30d and let-7, and that expression levels of the RFX6 transcript and these miRNAs are correlated in some contexts and anti-correlated in others, suggest that there are multiple competing regulatory factors controlling the expression of RFX6 expression during the course of pancreatic development, and that the interplay among these factors may influence the order and timing of the events in this process.

In conclusion, we performed comprehensive miRNA and mRNA expression profiling to identify the dynamic expressions patterns that are present during in vitro differentiation of hESCs to beta islet-like cells, to assess the relationships between differentiated hESCs and bona fide pancreatic tissue samples (fetal pancreas and adult islets), and to identify novel miRNA–mRNA interactions that regulate the differentiation process. Experimental confirmation of miRNA–mRNA interactions established a critical role of miR-200a in regulating both EMT and definitive endoderm formation, through direct repression of ZEB2 and SOX17, during early stage differentiation. In addition, we found that miR-30d and let-7e regulate the pancreatic progenitor gene RFX6 during late stage differentiation. Our data highlight the importance of precise control of target mRNA expression by miRNAs to ensure proper lineage specification during the differentiation process. The miRNA signatures identified in this study may also be useful for future investigations aimed at characterizing specific miRNAs involved in pancreatic beta islet cell regeneration and pathogenesis.

Undifferentiated WA09 cells were maintained in the absence of feeder cells on Geltrex (1∶200, Life Technologies) at 37°C, 5% CO₂ in DMEM/F12 medium supplemented with 2% Stempro supplement (Life Technologies, Carlsbad, CA), 2% BSA, 4 ng/ml recombinant human FGF2 (R&D Systems), 2-mercaptoethanol (0.55 mM, Life Technologies) and 1% Glutamax. Cells were plated at a density of 60,000–100,000 cells per well in 24-well Geltrex-coated tissue culture plates and fed daily with fresh medium. Differentiation was carried out as previously described (D'Amour et al., 2006) with minor modifications. Briefly, 1% BSA and 1% NEAA were added in the first 3 days of treatment to increase cell viability. The ABI-13D cell line and hiPSCs were cultured and differentiated in the same conditions as WA09. MIN6 was cultured as previously described (Tang et al., 2009). Fetal pancreas samples were collected according to an approved IRB protocol at UCSD (IRB #081510ZX) in collaboration with Planned Parenthood of the Pacific Southwest. Purified human beta islet cells were obtained through subscription to the Integrated Islet Distribution Program (IIDP, http://iidp.coh.org/).

Cultured cells were fixed for 15 minutes at room temperature in 4% paraformaldehyde dissolved in 1× phosphate-buffered saline (PBS). After fixation, cells were washed three times in PBS, and blocked for 1 hour in blocking solution that contains 0.1% Triton X-100, 2% BSA and 2% low fat milk in PBS. Primary antibodies diluted in the same blocking solution were incubated for 1–2 hours at room temperature or overnight at 4°C. Secondary antibodies diluted in blocking solution were incubated for 0.5–1 hour at room temperature. The following antibodies and dilutions were used: primary mouse monoclonal anti-TRA-181 IgG 1∶100 (Santa Cruz Biotechnology, sc-21706); mouse anti-POU5F11∶100 (Santa Cruz Biotechnology, SC-9081); goat anti-SOX17, 1∶100 (R&D system, AF1924); goat anti-HNF1β, 1∶100 (Santa Cruz Biotechnology, SC-7411); goat anti-PDX1, 1∶100 (R&D system, AB3243); rabbit anti-HNF6, 1∶100 (Santa Cruz Biotechnology, SC-13050); mouse anti-NKX2-2, 1∶10 (Developmental Studies Hybridoma Bank, 74.5A5); mouse anti-NKX6-1 (F55A12); guinea pig anti-insulin, 1∶200 (Dako, A0564), mouse anti-E-Cadherin (R&D, MAB1838). Cy3- and Cy5-conjugated donkey antibodies against mouse, rabbit and guinea pig (Jackson ImmunoResearch Laboratories), or Alexa-Fluor-488- and Alexa-Fluor-555-conjugated donkey antibodies against mouse, rat, rabbit, guinea pig or goat (Molecular Probes) were used at 1∶400 dilution.

Total RNA, including miRNA, was purified from all cell types using the mirVANA miRNA Isolation Kit (Ambion, Austin, TX). RNA quantitation was performed using Qubit® RNA Assay Kit (Life Technologies). RNA quality was assessed using the Agilent 2100 Bioanalyzer. For cells cultured in dishes, the RIN (RNA integrity number) was always more than 9.0. For tissue samples, RIN >8.0 was accepted. For real-time quantitative PCR of mRNA, 500 ng RNA was first reverse transcribed into cDNA using QuantiTect reverse transcription kit (Qiagen). PCR reactions were carried out in duplicate with QuantiTect SYBR Green master mix (Qiagen) according to the manufacturer's instructions. Quantified values for each gene of interest were normalized against the input determined by two housekeeping genes (GAPDH and β-actin). For real-time quantitative PCR of miRNA, 10 ng RNA was reverse transcribed using RNA-specific stem-looped RT primers and amplifications were carried out by Taqman probes (Applied Biosystems Inc.). RNU44 were used as endogenous control for real-time quantification of miRNA according to the manufacturer's instructions. After normalization, the samples were plotted relative to control samples in the dataset and the standard deviation of at least three gene expression measurements was reported. Information about primer sequences is provided in supplementary material Table S5.

The miRNA expression data were acquired on the Illumina Human v2 miRNA platform (with 1146 miRNA probes) and the mRNA expression data was acquired on the Illumina Human HT-12 v4 gene expression microarray (>40,000 mRNA transcripts), using 200 ng total RNA per sample each. Gene expression values were first filtered by detection P-value in Genome studio with a cutoff value of <0.01 in at least one sample. The miRNA expression data was then quantile normalized using CLC bio (Muehltal, Germany), while the mRNA data was normalized using robust spline normalization (RSN) with lumi in R 2.13.1 (scripts available upon request). ANOVA and post-hoc ANOVA Tukey's tests were performed using SAS 9.3 software (scripts are available upon request). Hierarchical clustering was performed by Cluster 3.0 with uncentered correlation and complete linkage (Eisen et al., 1998), or by Expander 6.0 CLICK clustering module with homogeneity value set up as 1.0. The heatmaps were produced by Java Treeview (after cluster 3.0) or by Expander 6.0 after CLICK clustering. The integrated miRNA and mRNA data analysis was performed in GALAXY. The miRNA target prediction algorithms were set as TargetScan 6.1. Gene ontology analysis was performed using DAVID Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/).

psiCHECK2 vector was purchased from Promega. To generate reporter plasmids of multiple genes simultaneously, the Multisite Gateway system (Life Technologies) was used following the manufacturer's instructions. Briefly, 10 µg of psiCHECK2 DNA was digested using EcoRI and blunt-ended with Klenow I DNA polymerase. The R1-ccdB-chloramphenicol-R2 (R1R2) cassette was purified from pLenti 6.2V5 (Life Technologies) after being digested with EcoRV. The R1R2 cassette was then ligated into psiCHECK2 to create a dest vector for future use. Genomic DNA from WA09 hESCs was isolated using the Qiagen DNAeasy Blood and Tissue kit, to serve as template for the following PCR reactions. The full length 3′UTR sequences of human genes RFX6, NKX2_2, ONECUT, FOXA2, SOX17, PAX6, NEUROD, ISL1 and ZEB2 were obtained from GenBank and amplified respectively using standard procedures with the primers listed in supplementary material Table S6. BP reactions were set up using the pDONR221 vector (Life Technologies) and the above PCR products containing B1B2 adapters. The products were L1L2 entry vectors for each gene. To generate the mutated SOX17 3′UTR reporter plasmid, we eliminate the binding site of SOX17 3′UTR (cagtgtt) with primer pairs designed to amplify the entry vector p221 L1-Sox17 3′UTR-L2 using Platinum® PCR SuperMix High Fidelity (cat. no. 12532-016 Invitrogen). Finally, these entry vectors were cloned into R1R2 dest vector respectively to assemble the 3′UTR sequences downstream of the Renilla luciferase expression cassette, which allows assessment of the effect of the 3′UTR on transcript stability by Renilla luciferase activity. Firefly luciferase expression is driven off the HSV-TK promoter on the same plasmid to control for well-to-well variation in transfection efficiency. All of the reporter constructs were verified by sequencing.

WA09 cells of passages 40–60 maintained feeder-free on Geltrex-coated plates were fed daily with Stempro medium (Life Technologies) and split ∼1∶6 when confluent using Accutase (Life Technologies). About 60,000–100,000 cells were plated in each well of a 24-well tissue culture plate coated with Geltrex. On day 1, cells were transfected with 4 µg reporter plasmid with or without 50–100 nM miRNA mimics or hairpin inhibitors (Dharmacon). The transfection reagent used for hESCs, including WA09, ABI-13D and iPSCs was FuGene HD (Roche) at 2–4 µl/well in 24-well plates. Dy547-labeled miRNA inhibitor (Dharmacon) was used as a transfection control. Synthetic miRNA mimics or hairpin inhibitors were all purchased from Dharmacon.

Adherent cells transfected with miRNA reporter plasmids with or without miRNA mimics or inhibitors were processed using the Promega Dual Luciferase Reporter Assay Kit according to the manufacturer's instructions. The assay was performed on a Synergy HT Multi-Mode Microplate Reader (Biotek) using injectors preprogrammed according to the manufacturer's recommendations (Gen5). Data was first normalized to the constitutively expressed firefly luciferase gene and then compared to scramble control (for mimics) or control inhibitors (for hairpin inhibitors). Data was reported as a percentage with control being 100.

Differentiating WA09 cells were transfected with miR-200a mimics (200 nM) on day 0 and day 2 of differentiation and harvested on day 3 in RIPA buffer containing 1% NP40, 0.5% sodium deoxycholate, 0.1% SDS, 0.5 mM phenylmethylsulfonyl fluoride, and 1% Halt™ protease inhibitor cocktail (Thermo Scientific). Equivalent amounts of proteins (∼10 ng) were resolved by SDS-PAGE and then transferred to nitrocellulose membranes for immunoblotting as described previously (Wang et al., 2008b). The primary antibodies (SOX17 and CDH1) used here were the same as used for immunocytochemistry. ACTIN antibody was purchased from MP Biomedicals. The intensity of western blot bands was quantified with ImageJ software.

For experiments with three independent replicates, a paired two-tailed t-test was used to test for significance between the controls and experimental samples. In all figures, the mean ± standard deviation was plotted. A single asterisk indicates P≤0.05 and a double asterisk indicates P≤0.01 in the t-test.

The authors thank Candace Lynch, Ha Tran, Trevor R. Leonardo, Sergio Mora Castilla and Francesca S. Boscolo for technical assistance, and Robert Morey, Rathi Thiagarajan and Ivka Afrikanova for critical reading of this manuscript. The fetal tissues were obtained under UCSD IRB approval through a collaboration with Planned Parenthood of the Pacific Southwest, with the assistance of Mana Parast and Sandra Leon-Garcia. The adult islets were obtained from the NIH-supported Islet Cell Resource.