Morpholinos for splice modificatio

Morpholinos for splice modification

Advertisement

miR-335 promotes mesendodermal lineage segregation and shapes a transcription factor gradient in the endoderm
Dapeng Yang, Dominik Lutter, Ingo Burtscher, Lena Uetzmann, Fabian J. Theis, Heiko Lickert

Abstract

Transcription factors (TFs) pattern developing tissues and determine cell fates; however, how spatio-temporal TF gradients are generated is ill defined. Here we show that miR-335 fine-tunes TF gradients in the endoderm and promotes mesendodermal lineage segregation. Initially, we identified miR-335 as a regulated intronic miRNA in differentiating embryonic stem cells (ESCs). miR-335 is encoded in the mesoderm-specific transcript (Mest) and targets the 3′-UTRs of the endoderm-determining TFs Foxa2 and Sox17. Mest and miR-335 are co-expressed and highly accumulate in the mesoderm, but are transiently expressed in endoderm progenitors. Overexpression of miR-335 does not affect initial mesendoderm induction, but blocks Foxa2- and Sox17-mediated endoderm differentiation in ESCs and ESC-derived embryos. Conversely, inhibition of miR-335 activity leads to increased Foxa2 and Sox17 protein accumulation and endoderm formation. Mathematical modeling predicts that transient miR-335 expression in endoderm progenitors shapes a TF gradient in the endoderm, which we confirm by functional studies in vivo. Taken together, our results suggest that miR-335 targets endoderm TFs for spatio-temporal gradient formation in the endoderm and to stabilize lineage decisions during mesendoderm formation.

INTRODUCTION

The first lineage decision during mouse development occurs when the morula develops to the blastocyst stage at embryonic day (E) 2.5-3.5 (Rossant and Tam, 2009). During this time, the inner cell mass (ICM) of the blastocyst segregates from the trophectoderm (TE) that will form the placenta. The ICM further develops to the epiblast epithelium that will give rise to all differentiated cell types in the mammalian body. At E6.5, gastrulation starts and pluripotent Oct4+ epiblast cells undergo epithelial-mesenchymal transition (EMT) to ingress into the posterior primitive streak (PS) region to form mesoderm and definitive endoderm (DE), whereas the remaining epiblast cells will form the ectoderm (Beddington and Robertson, 1999; Tam and Loebel, 2007; Zorn and Wells, 2009). Both EMT and mesendoderm differentiation are induced by Wnt/β-catenin and Nodal/TGFβ signaling in the mouse embryo (Tam et al., 2006; Arnold and Robertson, 2009). Wnt/β-catenin signaling leads to the activation of target genes in the epiblast, such as the Forkhead transcription factor a2 (Foxa2) (Sasaki and Hogan, 1993; Monaghan et al., 1993; Sawada et al., 2005) and the T-box transcription factor Brachyury (T; Herrmann, 1991; Yamaguchi et al., 1999; Arnold et al., 2000), which mark distinct mesendodermal progenitor cell populations in the posterior epiblast (Burtscher and Lickert, 2009). Genetic lineage tracing experiments revealed that Foxa2+ and T+ mesendoderm progenitors give rise to anterior and posterior mesendoderm populations, respectively (Uetzmann et al., 2008; Horn et al., 2012; Kumar et al., 2007; Verheyden et al., 2005); however, how lineage-inappropriate TF expression is prevented after lineage segregation into mesoderm and endoderm occurs is currently not known.

While delaminating from the epiblast epithelium and migrating into the PS, the mesendoderm progenitors upregulate cell type-specific molecular programs. Foxa2 and the SRY-related HMG box transcription factor 17 (Sox17) become strongly upregulated in the nascent DE during intercalation into the outside visceral endoderm (VE) layer (Burtscher et al., 2012). Knockout studies have revealed that both TFs are crucial cell-fate determinants and master regulators of DE differentiation (Ang and Rossantt, 1994; Kanai-Azuma et al., 2002; Burtscher and Lickert, 2009). Similarly, T highly accumulates in the nascent mesoderm population and is crucially important for mesoderm formation (Wilkinson et al., 1990; Burtscher and Lickert, 2009). At the end of gastrulation (E7.5), all DE cells are recruited and form an epithelial sheet of ∼500 cells on the outside of the mouse embryo (Wells and Melton, 1999). Already at this early stage of development, the endoderm is patterned along the anterior-posterior (A-P) axis with high levels of Foxa2 and Sox17 protein accumulating in the anterior endoderm (Burtscher et al., 2012). The first endoderm recruited from the epiblast progenitors migrates in anterior direction to overlie the forming headfold and is by definition older than posterior endoderm (Lawson and Pedersen, 1987; Thomas et al., 1998). Whether this allows these cells to accumulate more TF protein over time and how spatio-temporal TF gradients are established are far from being understood. Nowadays it is thought that morphogen gradients differentially regulate gene regulatory networks along the major body axes to set up TF gradients and specify distinct cell fates (Tabata and Takei, 2004; Kutejova et al., 2009; Rogers and Schier, 2011; Buechling and Boutros, 2011). However, emerging evidence indicates that micro-RNAs (miRs) are prime candidates to fine-tune signaling pathways and TF gradients in a dose-sensitive manner for pattern formation (reviewed by Inui et al., 2010). In this study, we specifically addressed the function of miRs in Foxa2- and Sox17-mediated A-P endoderm patterning.

Initially, miRs were discovered in regulating the process of developmental timing in Caenorhabditis elegans (Lee et al., 1993; Wightman et al., 1993; Reinhart et al., 2000). The heterochronic miRs lin-4/miR-125 and the let-7 family directly regulate lineage- and larval stage-specific genes, and therefore ensure correct developmental timing (Slack and Ruvkun, 1997; Ambros, 2000; Rougvie, 2001). These miRs target specifically mRNAs that code for regulatory proteins, to switch between symmetric versus asymmetric cell divisions or to induce terminal differentiation, thus ensuring developmental timing in the worm (Ambros, 2011). In vertebrates, miRs and processing enzymes control the developmental progression during oocyte maturation (Murchison et al., 2007) and ESC differentiation (Wang et al., 2007; Xu et al., 2009; Gangaraju and Lin, 2009). This suggests that vertebrate miRs function analogously to regulate spatio-temporal transitions in gene expression programs during cell-fate acquisition. Indeed, this is directly illustrated by miR-15 and miR-16, which target the Nodal receptor Activin receptor type II and are negatively regulated by Wnt/β-catenin signaling in Xenopus laevis to establish a dorso-ventral Nodal signaling gradient (Martello et al., 2007). Moreover, miR-430 targets both the Nodal agonist squint and the Nodal/TGFβ antagonist Lefty to control the availability of Nodal ligands in the extracellular space and establish morphogen gradients in zebrafish embryos (Choi et al., 2007). Thus, several different miRs have co-evolved with a handful of signaling cascades and TFs to ensure developmental timing and establish pattern formation in the early embryo.

miRs are small non-coding RNAs that either cause the degradation of the targeted mRNAs or block mRNA translation (Bartel, 2004; Guo et al., 2010; Mourelatos et al., 2002; Pratt and MacRae, 2009). Many miR genes are transcribed by RNA polymerase II or are encoded within protein coding host genes as intragenic miRs. Intragenic miRs located in introns are termed intronic miRs and are often co-expressed with their host genes. The primary transcript or pri-miR is processed by Drosha to form pre-miR hairpin loops. The pre-miRs are transported from the nucleus into the cytoplasm, where Dicer forms a mature miR complex that is then incorporated into the RNA-induced silencing complex (RISC) (Hammond et al., 2000; Kim and Kim, 2007). Importantly, both knockouts of Dicer and Dgcr8 prevent the differentiation of ESCs, because they fail to efficiently downregulate the pluripotency network and to upregulate differentiation factors (Kanellopoulou et al., 2005; Murchison et al., 2005; Wang et al., 2007). Global inhibition of the miR pathway can serve as a starting point to investigate miR function in a biological process (Konopka et al., 2010), but these experiments will not reveal the specific function of certain miRs in ESC and embryonic differentiation.

Here, we have specifically addressed the issue of whether miRs could potentially regulate mesendoderm lineage segregation and endoderm pattern formation. Therefore, we concentrated on miRs that are differentially expressed during mesendoderm and endoderm differentiation, and specifically target the 3′-UTRs of the endoderm TFs Foxa2 and Sox17. This identified the intronic miR-335 encoded in the mesoderm-specific transcript Mest that specifically targets Foxa2 and Sox17. Strikingly, we found one of the first examples of a miR that is important to shape a TF gradient in the endoderm germ layer. Additionally, the high expression of miR-335 (Mir335 - Mouse Genome Informatics) in mesoderm progenitors and lineage strongly suggests that miR-335 promotes lineage segregation by suppressing lineage inappropriate expression of endoderm TFs.

RESULTS

Identification of miR-335 as a potential regulator of mesendoderm development

We have recently noticed that the endoderm-specific TFs Foxa2 and Sox17 accumulate in an A-P gradient in the endoderm germ layer (Burtscher et al., 2012). How this TF gradient in the endoderm is established is not clear. To test whether miRs are involved in setting up this gradient, we focused on differentially expressed miRs during mesendoderm formation that have predicted target recognition sites in the 3′-UTRs of Foxa2 and Sox17. For this purpose, we performed a four-step filtering approach to identify miRs that potentially target one or both of the selected TFs (Fig. 1A). Initially, 37 miRs potentially targeting Foxa2 and Sox17 were selected by using miRecords (Xiao et al., 2009), an integrated resource for miRNA target prediction that uses several miRNA target prediction tools (including PITA, miRanda, RNAhybrid) and several hand-curated miRNA target interactions. To further select miRs that were potentially expressed during mesendoderm differentiation, we filtered for all intragenic miRs by assuming correlated expression of host gene and miR (Baskerville and Bartel, 2005), which resulted in 22 host gene-miR pairs. Intragenic miRs were identified using mirBase (Kozomara and Griffiths-Jones, 2011). Differential expression of host genes during embryonic body-induced ESC differentiation was analyzed using previously published mRNA expression profiles (Sene et al., 2007) (see Materials and Methods). This reduced our list to six host gene-miR pairs (supplementary material Table S1). We have previously shown that functional relationships between intronic miRs and their targets are often associated with correlated expression of the host and target gene (Lutter et al., 2010). Hence, in a last step, we decided to choose the miR that most putatively targets the two TFs by calculating correlation between host gene and Foxa2 and Sox17 mRNA, respectively (Fig. 1B,C). We found significant negative correlation between the host gene Mest and its intronic miR-335 to all targeted TF mRNAs (P<10-7). The miR-335 is located in the second intron of Mest and well conserved within mammals (Fig. 1D).

Fig. 1.

Identification of miR-335 as a possible regulator for endoderm-specific TFs. (A) Workflow for prediction of miRs that possibly target Foxa2 and Sox17. (B) To select host genes showing reciprocal expression patterns to Foxa2 and Sox17, we used a correlation test for each host-TF pair. The boxplot shows the distribution of all host gene-TF pairs. (C) Fold change in Foxa2 and Sox17 mRNA during mesendodermal development stages. Whereas the endoderm-specific TFs show a slight upregulation during the mesendodermal stage, Mest is downregulated when cells pass the mesendodermal stage. (D) Mest locus showing the genomic sequence encoding miR-335, which is located in the second intron of Mest. For miR-335, the sequence, the stem-loop and conservation score for each base is shown.

Spatio-temporal expression of miR-335 during mesendoderm development

To obtain further evidence for a potential miR-335 function during mesendoderm lineage formation, we first investigated the spatio-temporal expression pattern of miR-335 and Mest in differentiating ESCs and mouse embryos. miR-335 is an intronic miR embedded in intron 2 of the Mest mRNA transcript (Fig. 1D) and according to previous findings should be co-expressed with its host gene (Ronchetti et al., 2008). Using quantitative PCR (qPCR), we confirmed that the relative miR expression level correlates with the host gene expression during mouse development at embryonic day (E) 6.5-8.5 (Fig. 2A). Furthermore, we used whole-mount in situ hybridization with antisense probes against Mest and miR-335 to analyze the tissue-specific expression of host gene Mest and miR-335. We used the heart-specific miR-1 as a positive control (Kloosterman et al., 2006) and a scrambled miR-335 as negative control (Fig. 2B). The results revealed that miR-335 and Mest mRNA are co-expressed in a tissue-specific manner in the brain, branchial arches, heart, limb buds, and somites at E9.5.

Fig. 2.

Co-expression of Mest/miR-335 in mouse ESCs and embryos. (A) The temporal expression levels of miR-335 correlate with those of its host coding-gene Mest in mouse embryos at E6.5-8.5. (B) Whole-mount in situ hybridization shows co-expression of Mest and miR-335 in the brain (b), heart (h), limb bud (lb) and branchial arches (ba) at E9.5. The heart- and somite-specific expression of miR-1 served as a positive control, whereas a scrambled miR-335 served as a negative control. s, somite. (C) The level of miR-335 in the Foxa2+ mesendoderm and endoderm lineage is transiently upregulated during endoderm differentiation (left graph), whereas miR-335 accumulates at high levels in the T+ mesendoderm and mesoderm lineage during mesoderm differentiation (right graph). (D) Differentiation paradigms for mesoderm and endoderm differentiation. ESCs cultured in Wnt3a+/ActA++ (high ActA) media differentiate into endoderm that expresses Foxa2 and Sox17, whereas cultures in Wnt3a+/ActA+ (low ActA) differentiate into mesoderm that expresses T. We monitored the following linage markers: ESC, Oct4+; mesendoderm, T+ and Fox2a+; endoderm, Foxa2+ and Sox17+; mesoderm, T+.

As we have previously shown that miRs regulate pathways in a tissue-specific manner (Kowarsch et al., 2011), we further analyzed cell type-specific Mest and miR-335 expression during ESC differentiation (Fig. 2C). We monitored the progression from pluripotent ESC (Oct4high) to mesendoderm (Foxa2+, T+) and further to mesoderm (T+) and endoderm (Foxa2+, Sox17+) using an established ESC differentiation system and fluorescent protein (FP) reporter knock-in ESC lines. For this purpose, we used a transcriptional Brachyury (T)-GFP (Fehling et al., 2003) and two directly in-frame TF-FP fusion reporter ESC lines: Foxa2-Venus fusion (FVF) and Sox17-mCherry fusion (SCF) (Burtscher et al., 2012; Burtscher et al., 2013). Triggering Wnt/β-catenin and Nodal/ActivinA (ActA) signaling, which induce gastrulation in the mouse embryo, leads to differentiation of adherent ESC colonies into the mesendoderm lineage under serum-free conditions (Yasunaga et al., 2005). Further differentiation into the endoderm lineage requires high levels of ActA, whereas differentiation into the mesoderm lineage requires low levels of ActA (Kubo et al., 2004; Tada et al., 2005; Fig. 2D). Using these two different ESC differentiation paradigms, we isolated the Foxa2+ mesendoderm and endoderm lineage using the FVF ESC line under endoderm conditions, and isolated the T+ mesendoderm and mesoderm populations using the T-GFP ESC line under mesoderm conditions. This was carried out with a time-course experiment that used fluorescent-activated cell sorting (FACS). miR-335 expression levels were analyzed by qPCR (Fig. 2C; supplementary material Fig. S1). This revealed that miR-335 is transiently upregulated ∼40-fold in the Foxa2+ mesendoderm progenitors after 24 hours of ESC differentiation. By contrast, miR-335 levels steadily increase over 4 days by 1000-fold or more in the T+ mesendoderm and mesoderm lineage, consistent with the strong mesoderm-specific expression in heart and somites, but not in the gut endoderm in vivo (Fig. 2B). Together, the spatio-temporal expression of miR-335 suggests that it functions transiently in the Foxa2+ endoderm progenitors and later during mesoderm formation.

Functional analysis of miR-335 during mesendoderm development

As miR-335 has predicted target recognition sites in the 3′-UTRs of the endoderm-specific TF mRNAs Foxa2 and Sox17 (Fig. 1A), we first tested whether it directly targets these TFs. Therefore, we generated a position weight matrix for the miR-335-5p binding motif (Fig. 3A; supplementary material Fig. S2 and Table S2) using the previously described miR-335 target mRNA-binding sites (Tavazoie et al., 2008; Miranda et al., 2006). Scanning the 3′-UTRs identified one and two miR-335-5p-binding motifs for the Sox17 and Foxa2 mRNA, respectively, which we mutated at the highly conserved positions (Fig. 3B). Using a dual Renilla-firefly luciferase system, we analyzed the effect of transient miR-335 precursor transfection in reporter gene assays on the wild-type or mutant 3′-UTRs (Fig. 3C). This revealed that both Foxa2 and Sox17 3′-UTR are significantly downregulated (P<0.01) by miR-335, which can be rescued to different degrees following mutation of the miR-335-5p-binding site (Fig. 3D). In summary, full rescue of the Sox17 3′UTR reporter activity by mutation of the miR-335-5p-binding motif suggests direct regulation, whereas additional miR-335-3p or non-consensus binding motives in the Foxa2 3′UTR might be targeted by miR-335-3p (supplementary material Fig. S2).

Fig. 3.

Foxa2 and Sox17 are targets of miR-335. (A) miR-335-5p target site consensus sequence generated using six known target sites from mouse and human. (B) Summary of miR-335-5p target sites in the 3′-UTR of Foxa2 and Sox17. Nucleotides shown in red indicate changes in the mutant 3′-UTR. (C) Scheme of the reporter constructs used for target validation. huR-Luc, humanized Renilla luciferase; huF-Luc, humanized firefly lucferase; pA, polyadenylation signal; WT, wild type. (D) miR-335 directly represses its targets in a luciferase assay in HEK293T cells. Renilla luciferase activity was assayed 40 hours after transfection and the values were normalized to the activity of firefly luciferase encoded in the same vector. One-way ANOVA, n=3.

To analyze the function of miR-335 in more detail, we generated several stable ESC lines that constitutively express miR-335. For this purpose, we used the RNA polymerase II-driven human ubiquitin C promoter to co-express miR-335 from a modified intronic miR-155 precursor followed by an exon coding for histone 2B-cyan fluorescent reporter protein (H2B-CFP) (Chung et al., 2006; supplementary material Fig. S3). Using this bicistronic vector expression system, we could directly correlate miR-335 expression with H2B-CFP fluorescent reporter activity (Chung et al., 2006). Two independent ESC lines with medium (miR-335 #1) and high (miR-335 #2) H2B-CFP reporter activity were used to analyze a dose-dependent effect of miR-335 on ESC differentiation (see Materials and Methods). We compared the miR-335 overexpression (miR-335 #2) with endogenous miR-335 levels measured in embryos to assure that the overexpression is in a physiological range (supplementary material Fig. S1C). Qualitative (Fig. 4A) and quantitative (Fig. 4B; supplementary material Fig. S4 and Table S3) analyses revealed that Foxa2- and Sox17-mediated endoderm differentiation was blocked by miR-335 gain of function in a concentration-dependent manner. Upon miR-335 gain of function, Oct4+ pluripotent ESC colonies remained round and only a few flattened Foxa2 and Sox17 double-positive DE cells appeared at the edge of the colonies (Fig. 4A). Western blot analysis confirmed these results and revealed that Foxa2-mediated mesendoderm induction at 48 hours occurred normally, whereas further differentiation into the Foxa2+ and Sox17+ DE lineage was strongly reduced (Fig. 4C). Even after 96 hours, Oct4 protein levels remained high, suggesting that differentiation was blocked at the Foxa2+ mesendoderm progenitor cell stage (Fig. 4C). Finally, we analyzed the effect of miR-335 gain of function in completely ESC-derived mouse embryos in vivo (Nagy et al., 1993; Tam and Rossant, 2003). Tetraploid (4n) embryo⟷miR-335#2 ESC aggregation chimera implanted normally and proceeded through development until early gastrulation stage, confirming that miR-335 gain of function has no generic and deleterious effect on cell biology per se (Fig. 4D). Strikingly, whereas control chimera normally generated Foxa2+ mesendoderm and Foxa2+/Sox17+ double-positive DE, DE induction was almost completely blocked in miR-335-overexpressing chimeras at the Foxa2+ mesendoderm stage (Fig. 4D). These in vivo results directly reflect our ESC differentiation results (Fig. 4A-C) and suggest that miR-335 gain of function blocks DE differentiation at the Foxa2+ mesendoderm progenitor stage in vitro and in vivo.

Fig. 4.

miR-335 overexpression represses endoderm differentiation. (A-C) During endoderm differentiation, both Foxa2 and Sox17 were significantly downregulated by overexpression of miR-335 at medium and high levels, when compared with a control clone, as shown by immunostaining (A), quantitative analysis of immunofluorescence (B) and western blot analysis of Foxa2, Sox17 and Oct4 in differentiating mESCs (C). (D) Immunostaining of completely ES cell-derived mouse embryos at gastrulation stage at E7.5. Foxa2 and Sox17 expression was found only in newly formed endoderm cells in the anterior primitive streak region, whereas Foxa2 and Sox17 were strongly suppressed in older endoderm cells in the anterior and lateral regions of miR-335-overexpressing embryos when compared with the wild-type controls. Anterior is towards the left, distal is towards the bottom.

To gain further insight into the physiological function of miR-335 during mesendoderm formation, we performed loss-of-function experiments in ESC differentiation cultures using competitive inhibition by overexpression of a sponge construct. For this purpose, we used the CMV-enhancer β-actin promoter (CAG) to express either H2B-BFP (H2B-blue fluorescent protein) followed by a polyadenylation sequence (pA) alone or H2B-BFP followed by a 3′-UTR with 18 miR-335-5p-binding motifs (sponge-3P). For live-cell analysis at the single-cell level, we used a Foxa2Venus/+; Sox17Cherry/+ ESC line (Burtscher et al., 2013) to generated three different ESC subclones to constitutively express the control or sponge-3P construct (Fig. 5A; supplementary material Fig. S5). Both FVF and SCF knock-in reporter ESC lines use the endogenous 3′-UTR of the Foxa2 and Sox17 mRNAs, and act as miR sensors for direct analysis of TF levels. The analyses of FVF and SCF reporter activity by quantitative single cell FACS analysis using three independent ESC clones expressing either sponge-3P construct or corresponding control in a time-course experiment suggested that miR-335 loss of function led to an increase in FVF and SCF reporter-positive cells after 96 hours of Wnt3a/ActA-mediated DE induction (Fig. 5B; supplementary material Fig. S5). Together, these results suggest that endogenous miR-335-5p activity blocks Foxa2 and Sox17 translation, which can be specifically released by competitive inhibition for enhanced endoderm formation.

Fig. 5.

Knockdown of miR-335 leads to the increase in endoderm formation. (A) The control and sponge-3P expression vectors. (B) FACS analysis of a time-course experiment over 96 hours of endoderm differentiation shows an increase in the Foxa2 and Sox17 double-positive population in sponge-3P-expressing clones. Pooled FACS data of three control and three sponge clones were analyzed.

Mathematical modeling predicts miR-335 function in TF gradient formation

To analyze and describe the effect of miR-335 on Foxa2- and Sox17-mediated mesendoderm and DE differentiation in a continuously quantitatively resolved fashion, we generated an in silico molecular mathematical model based on ordinary differential equations (ODEs; Fig. 6A). The model is based on the generic models previously used to analyze miR-mediated effects on protein expression (Levine et al., 2007; Mukherji et al., 2011). In contrast to these models, where miR concentrations were assumed to be constant, we introduced dynamic miR turnover rates, as deduced from qPCR data. For all RNA molecules, distinct transcription rates, as well as the free-mRNA and protein degradation rates, were unknown and therefore estimated from the data (supplementary material Appendix S1). To test for parameter identifiability and to estimate confidence intervals, we exploited the profile likelihood estimation (PLE) (Raue et al., 2009). We assume that degradation of the complex leads to a partly degradation of the miR; thus, a fraction of the active miR is recycled. Furthermore, we assume that miR half-lives are much longer compared with mRNA half-lives (Krol et al., 2010); thus, the model does not include a separate degradation rate for the free miR. Referring to the Foxa2 expression data, we assume a constant transcription kf rate for the target mRNA and, as miR expression decreases after 24 hours, we model the miR transcription using a bell-shaped function with an estimated maximum at time=0. The transcription rate km thus refers to the maximal miR transcription rate.

Fig. 6.

Modeling miRNA-mediated protein expression dynamics. (A) Model describing the dynamic behavior of the free mRNA, the free miRNA, the complex (miR bound to mRNA) and the translated protein. (B) Model parameter estimation for Foxa2-positive cells. Experimental data (black dots) were measured from FACS-sorted Foxa2+ cells. Foxa2 mRNA and miR-335 were measured using qPCR; Foxa2 protein expression was measured using western blot. Pulsed miR-335 expression was modeled using a Gaussian expression function. Gray areas denote the confidence intervals for the parameter estimation. (C) Model prediction for different miR expression rates (km). Protein expression was predicted using different miR-335 expression rates. The green solid line displays protein dynamics for the estimated miR-335 expression rate (km=1.8) for pulsed expression, as estimated from experimental data. Protein dynamics for a simulated complete miR-335 knockdown (km=0) is shown by the broken orange line. Linear miR-335 (km=10) expression as estimated for the mesodermal T+ cells is shown by the broken blue line. The vertical broken red line indicates t=48 hours. (D) Predicted protein gradient for different miR-335 transcription rates (km). The simulated protein gradient at time=48 hours is shown as a function of miR-335 transcription rates (km log10). The gradient increases for increasing transcription rates until it drops when a miR-mediated knockdown is reached. The red line indicates the estimated transcription rate for Foxa2-positive cells. The curve to the left of the red line shows simulated loss-of-function behavior; the curve to the right of the red line shows gain-of-function behavior.

To calibrate the model parameter, we used qPCR-measured mRNA and miR data from FACS-isolated Foxa2 lineage-positive cells (Fig. 2C; supplementary material Fig. S1B, Appendix S1), and FVF and SCF reporter activities measured in cell culture at 0, 48 and 96 hours after Wnt3a/ActA-mediated differentiation (Fig. 4B, Fig. 6B; see Materials and Methods). To study the model dynamics, we start simulations with very low RNA and protein concentrations, as predicted from the experimental data. After calibrating the model parameters, we studied model dynamics for the protein readout depending on different miR transcription rates. To measure miR-dependent effects on the protein dynamics, we calculated the time-dependent protein gradient at time=48 hours, which corresponds to half of simulated differentiation time (Fig. 6C; see Materials and Methods). For the estimated physiological miR transcription rate (km) of 9.8 [miR/time a.u.] we observed a sigmoid-like expression curve. We then analyzed miR loss of function and gain of function by simulating the model with km between 0 [miR/time a.u.] and 25 [miR/time a.u.] (Fig. 6C). Overexpression of miR-335 completely blocks protein expression, whereas loss of miR-335 function allows for a faster accumulation of the target protein, thus decreasing the gradient at time=48 hours. We proved the confidence of the model by comparing gain- and loss-of-function trajectories for estimated parameter below the 95% PLE threshold (supplementary material Appendix S1). The continuous miR-335-dependent protein (Foxa2) gradient for time=48 hours is displayed in Fig. 6D.

miR-335 shapes a TF gradient in the endoderm

In differentiating ESCs, miR-335 is transiently expressed in the Foxa2-mesendoderm progenitors, but quickly downregulated in the DE (Fig. 2C). Moreover, miR-335 gain of function blocks DE formation after mesendoderm induction (Fig. 4), whereas miR-335 loss of function increases DE formation and Foxa2/Sox17 protein accumulation (Fig. 5). Therefore, we assume that miR-335 acts at the level of the Foxa2+ mesendoderm progenitor to shape a gradient of Foxa2 and Sox17 in the anterior-posterior patterned DE, as predicted by the mathematical model. To validate this model experimentally, we analyzed TF gradient formation in the gastrula-stage mouse embryo, where different morphogen activities translate into an A-P gradient of Foxa2 and Sox17 (Fig. 7A; Burtscher et al., 2012). Foxa2+ mesendoderm epiblast progenitor cells are recruited at the posterior side of the embryo, where miR-335 levels are still high, and differentiating Foxa2+/Sox17+ DE cells intercalate and migrate to the anterior side of the embryo (Burtscher et al., 2012), where according to our ESC data the miR-335 levels decreases (Fig. 2C). This is reflected in the accumulation of Foxa2 and Sox17 protein in an A-P gradient, as revealed by IHC and LSM analysis of wild-type embryos at E7.5 (Fig. 7A). We developed an automated image quantification method to determine the protein amounts in the DE along the A-P axis from fluorescent images (see Materials and Methods; supplementary material Figs S6-S8). Quantification of the Foxa2 and Sox17 protein levels revealed a spatial and temporal protein gradient along the A-P axis in single embryos (Fig. 7C) and in pooled embryo groups at gastrulation (Fig. 7E), which we quantified explicitly (Fig. 7C; supplementary material Fig. S8). As DE cells migrate over time along the A-P axis, we quantify a spatial gradient that reflects the time-dependent protein accumulation. From these data, we predicted that the miR-335 loss of function should lead to an increase of Foxa2 and Sox17 protein levels in developmentally younger cells and therefore should accumulate faster at the posterior side of the embryo. To confirm this, we generated completely ESC-derived sponge-3P expressing embryos and analyzed the TF gradient (Fig. 7B; supplementary material Fig. S8). Analysis of Foxa2 and Sox17 protein accumulation (Fig. 7D-E; supplementary material Fig. S9) and gradient formation (Fig. 7F) confirmed our model predictions and revealed that miR-335 functions to shape a TF gradient in the endoderm in vivo. Finally, we tested whether the change in the TF gradient at gastrulation stage leads to patterning defects in the gut tube at E8.5. Whole-mount in situ hybridization revealed that the foregut marker Pyy was normally expressed (Fig. 7G,H), whereas the midgut marker Nepn was strongly reduced in miR-335 sponge-expressing embryos when compared with controls (Fig. 7I,J). Taken together, these results suggest that miR-335 functions to shape a TF gradient in the endoderm that translates into gut tube patterning.

Fig. 7.

miR-335 shapes transcription factor gradients in the endoderm. (A,B) Upper row: whole-mount immunofluorescence staining for FVF and SCF fusion proteins in control and sponge-3P ESC-derived embryos showing levels of proteins at E7.5. Lower row: normalized intensities of TF gradient. (C,D) Plots display the normalized mean TF levels along the A-P axis for FVF and SCF in control and sponge-3P embryos. Loss of miR-335 function leads to increased FVF and SCF protein at the posterior, developmentally younger side. (E) Summary plot of all TF gradients in multiple embryos (colored lines). The mean (thick blue line) and standard deviation (blue area) is shown for FVF of control embryos (upper) and miR-335 loss-of-function (Sponge-3P) embryos. (F) Estimated protein gradient for control (n=7) and sponge-3P (n=4) embryos. The Foxa2 and Sox17 protein gradient was estimated in the medial range (dashed black lines in C-E) along the A-P axis. For both proteins, a less distinctive gradient was observed for loss of miR-335 function, consistent with the modeling results in Fig. 6. (G-J) Whole-mount in situ hybridization (WISH) with indicated probes at E8.5. CD1 control (G,I) and sponge-3P embryos (H,J) showed similar foregut (fg) expression of Pyy, whereas the midgut (mg) expression of Nepn is strongly reduced in sponge-3P embryos. Scale bars: 500 μm in G-J.

DISCUSSION

After the identification of hundreds of miRs, the challenge is now to understand their specific biological function and to define their molecular targets. Here, we have identified miR-335 as an intragenic miR embedded in the second intron of the mesoderm-specific transcript Mest. miR-335 is transiently expressed in Foxa2+ endoderm progenitors and highly accumulating in the T+ mesoderm lineage. It specifically targets the 3′-UTR of Foxa2 and Sox17, and gain of function blocks DE differentiation, whereas loss of function enhances DE formation. Quantitative mathematical modeling predicted miR-335 to have a function in endoderm TF gradient formation, which we confirmed experimentally in developing embryos. Taken together, our results suggest two physiological functions of miR-335: first, low miR-335 expression levels dampen Foxa2 and Sox17 protein levels in nascent Foxa2+ DE progenitors to establish a TF gradient along the A-P axis in the endoderm germ layer; second, high miR-335 expression levels promote mesendoderm lineage segregation and prevent lineage inappropriate expression of endoderm TFs in the T+ mesoderm lineage by default repression.

How morphogen gradients are established and how they pattern developing tissues in a dose-dependent manner is a long-standing question in developmental biology (Rogers and Schier, 2011). Both duration and concentration of morphogens, as well as the state of target cells determines cell fate in a spatio-temporal manner. It was long thought that the expression of miRs dampen rather than silence the expression of mRNA targets, which makes them prime candidates to fine-tune morphogen and TF gradients. This idea is based on the fact that genome-wide computational and transcriptome analyses showed that miR-mRNA target pairs correlate more positively than negatively in their tissue expression (Martinez et al., 2008). Prime examples are the regulation of Nodal signaling by Wnt/β-catenin-inhibited miR-15 and miR-16 during dorso-ventral patterning in Xenopus laevis (Martello et al., 2007) and the control of extracellular Nodal morphogen availability by miR-430 (Choi et al., 2007). Here, we provide another example of a miR that establishes a gradient not by targeting morphogens, but rather by directly targeting the mRNAs of the endoderm-determining TFs Foxa2 and Sox17. As TF gradients were often viewed as the integration and net result of several synergistic- and antagonistic-acting morphogen gradients, this example suggests that miRs dose-dependently regulate TF accumulation at the post-transcriptional level by dampening rather than silencing the target mRNAs, which adds an additional layer of regulation to fine-tune morphogen gradients.

But how can a miR dampen rather than silence its target mRNAs to generate a spatio-temporal gradient? Through quantitative mathematical modeling based on ordinary differential equations, we predicted that the Mest-miR-335 expression level directly correlates with the dose-dependent degree of TF accumulation in the DE. These predictions were confirmed by gain- and loss-of-function experiments in ESC differentiation and developing mouse embryos. The exact concentration and duration of miR-335 levels along the spatio-temporal axis in vivo is currently difficult to determine. However, our gain- and loss-of-function results demonstrate how different thresholds of miR expression either completely block or increase Foxa2 and Sox17 protein accumulation, thus explaining how low levels of miR-335 in endoderm progenitors lead to a delayed accumulation of TFs along the spatio-temporal axis for gradient formation.

Another intriguing question is how gene expression profiles are stabilized and lineage-inappropriate expression is suppressed after binary lineage decisions have occurred. In contrast to the co-expression of low levels of miR-335 with TF mRNA in the endoderm, our data are also an obvious example of mutually exclusive tissue-specific expression of miR-mRNA target pairs. Pioneering work in the fruit fly suggested that miRs function to repress their mRNA targets in tissues where they should not be expressed to confer robustness of gene expression (Stark et al., 2005). By embedding miR-335 in the second intron of the mesoderm-specific transcript Mest and under the control of mesoderm-specific enhancers, miR-335 is highly expressed in the heart, somites, limb bud and branchial arch mesenchyme. Interestingly, the same tissues are Foxa2+ lineage positive (Uetzmann et al., 2008; Horn et al., 2012), but express neither Foxa2 nor Sox17 mRNA at E9.5 (Kanai-Azuma et al., 2002; Sasaki and Hogan, 1993). This indicates that high levels of miR-335 suppress an endoderm TF program by default repression in mesoderm-derived tissues.

Translation of embryonic principles to ESC culture might eventually lead to the generation of functional cell types for tissue replacement and regenerative medicine approaches. Therefore, developmental design principals have to be translated to the culture dish for the generation of functional cell types, such as liver hepatocytes or pancreatic β cells. Triggering Wnt/β-catenin and Nodal/TGFβ signaling that induce gastrulation in the mouse embryo leads to differentiation of pluripotent ESCs into the DE lineage. Modulating miR-335 levels after DE induction may help to guide ESC differentiation to generate anterior versus posterior endoderm populations that can give rise to lung, liver, pancreas and gastro-intestinal tract along the A-P axis. Specified cell types, such as insulin-producing β-cells or GLP-1-producing L-cells of the gut are of great therapeutic interest to treat metabolic disease. Thus, our SCF and FVF miR sensor ESC lines might be valuable tools for understand the relationship between endoderm TF expression and cell-fate specification. Future analysis might allow us to model the effect of different levels of miR-335 on endoderm cell-fate specification in ESC differentiation culture.

Furthermore, several cancers, such as breast, brain, lung, pancreas and gastric cancers, are associated with differential miR-335 expression (Xu et al., 2012; Polytarchou et al., 2012; Yan et al., 2012; Dohi et al., 2013; Lynch et al., 2012; Vickers et al., 2012). The function of miR-335 in tumor initiation and progression is controversial. For example, high levels of hsa-miR-335 correlated with a high frequency of recurrence and poor survival in individuals with gastric cancer (Yan et al., 2012) and epigenetic silencing of miR-335 and its host gene MEST is associated with hepatocellular carcinoma (Dohi et al., 2013). Our investigation of the miR-335 function during normal development indicates that miR-335 targets key endoderm TFs downstream of Nodal/TGFβ and Wnt/β-catenin signaling. Interestingly, both signaling pathways are implicated in tumor formation and progression and miR-335 is differentially regulated during both processes (Polytarchou et al., 2012; Lynch et al., 2012). Our results suggest that overexpression of miR-335 in tumors might silence the epithelialization factor Foxa2 (Burtscher and Lickert, 2009) and contribute to EMT and metastasis formation. Furthermore, the loss of endodermal organ-specific expression of Foxa2 might lead to de-differentiation of mature cell types and acquisition of a more naïve proliferative cellular status. By contrast, silencing miR-335 and its host gene MEST by methylation might cause inappropriate upregulation of its target genes, such as Foxa2, Sox17 and Sox4, and promote tumor growth in several cancer forms (Dohi et al., 2013; Huang et al., 2012). Overall, our findings allow for new perspectives in several disease related aspects, both for cell replacement and cancer therapy.

MATERIALS AND METHODS

Mouse ESC culture and differentiation

Mouse ESCs were maintained and passaged on gelatin-coated plates with MEF feeders in ESC medium (Burtscher et al., 2012). Prior to differentiation, Wnt3a-expressing feeders (Kispert et al., 1998) were seeded at a density of ∼3×104 cells/cm2. ESCs were passaged onto gelatin-coated plates for 30 minutes twice to remove feeders. After washing with PBS, ESCs were transferred to the Wnt3a-expressing feeder plate with a seeding density of 7×104 cells/cm2. SFO-3 medium supplemented with low (3 ng/ml) and high (12 ng/ml) ActA (R&D Systems) to induce the endoderm differentiation and medium was changed daily.

Generation of expression vectors, fluorescent reporter ESC lines and chimeras

The generation of miRNA expression vector is described in supplementary material Fig. S2. For the sponge vector, the oligos contained each three miR-335-3p binding motives were ligated into the control pCAG-H2B-BFP-2A-IRES-Puro polyA vector. Subcloning of this oligos into the NsiI site was repeated six times, resulting in a total of 18× mir335-5p-binding sites in the 3′-UTR of the sponge-3P vector. These vectors were stably introduced into the Foxa2Venus/+; Sox17Cherry/+ double heterozygous ESC line (Burtscher et al., 2013). Electroporation was performed as previously described (Burtscher and Lickert, 2009). Cells were selected with 1 μg/ml puromycin, and resistant clones were screened for uniform and ubiquitous reporter expression by quantitative live fluorescent imaging. Tetraploid chimeras were generated according to standard protocols (Nagy et al., 1993).

Immunofluorescence imaging

Cells were fixed with 4% paraformaldehyde and then permeabilized in cold methanol. Cells were blocked with PBS containing 3% donkey serum, 10% fetal bovine serum (FBS), 0.1% BSA, 0.2% Tween-20 and incubated with primary antibodies overnight at 4°C. Then secondary antibodies were incubated for 1 hour at room temperature after washing with PBST. Images were acquired with a Zeiss Axiovision inverted microscope with a 20× objective. Immunofluorescence whole-mount staining was performed as previously described (Burtscher and Lickert, 2009). The following antibodies and dilutions were used: mouse anti-Oct3/4 (1:1000, Santa Cruz), rabbit anti-RFP (1:1000, Biotrend) and goat anti-Foxa2 (1:1000, Santa Cruz).

FACS and flow cytometry analysis

For FACS sorting, cells were dissociated using 0.05% trypsin-EDTA, suspended in PBS and filtered through nylon. Cells were analyzed and sorted using a Becton Dickinson FACS Aria IIIU equipped with a 488 nm laser.

For FACS analysis, cells were stained with primary antibodies (list above) diluted in blocking medium, and detected using fluorescent-conjugated secondary antibodies raised in donkeys. The following lasers were used to analyze the frequency of the cells: 405 nm, 488 nm, 561 nm and 633 nm.

RNA isolation and real-time PCR analysis

Total RNA was extracted from cells with the miRNeasy mini kit (Qiagen). cDNA synthesis was performed with a high-capacity RNA-to-cDNA kit (Applied Biosystems). The mouse miR-335, Foxa2 and Brachyury transcripts were quantified by real-time RT-PCR using the corresponding TaqMan Gene Expression Assays (Applied Biosystems) in an ABI PRISM 7900HT sequence detection System (Applied Biosystems) according to the manufacturer’s instructions. U6 and β-actin were used as endogenous normalization controls for miRNA and protein-coding genes, respectively. Endogenous mRNA levels of Mest and Gapdh were measured using SYBR Green PCR Master Mix (Applied Biosystems). The primer list is provided in supplementary material Tables S4, S5.

Protein extraction and western blot analysis

Cells were lysed in 1×RIPA lysis buffer in the presence of protease inhibitor. Protein samples were resolved by 10% SDS-PAGE and gels transferred to nitrocellulose membranes (Bio-Rad). After blocking with 5% skim milk in TBS/0.2% Tween-20 (TBST), the membrane was incubated overnight at 4°C with specific primary antibodies at the following dilutions: Oct3/4 (1:1000, Santa Cruz); anti-RFP (1:1000, Biotrend); Foxa2 (1:1000, Santa Cruz); Brachyury (1:200, Santa Cruz); Sox17 (1:5000, Acris/Novus) and Actin Ab-5 (1:5000, Bio-Rad). The membrane was washed several times in TBST, incubated with secondary antibodies for 1 hour at room temperature, and developed by chemiluminescence HRP Substrate (Millipore).

Luciferase reporter transfection and dual luciferase assay

The full-length 3′-UTRs and mutated 3′-UTRs of miR-335-5p target genes Foxa2 and Sox17 were amplified and individually cloned into the psiCHECK-3 (Promega) dual luciferase reporter vector (for primer sequence see supplementary material Fig. S3). HEK293T cells were co-transfected with each reporter construct and synthetic pre-miR-335 (Ambion) at the following concentration: 100 ng of the UTR reporter, 30 pmol miRNA precursor molecules and 25 μl of PEI (polyethylenimine, Polysciences) per 24-well plate. Cells were lysed 40 hours after transfection and the ratio of Renilla to firefly luciferase was measured with the dual luciferase assay (Promega). Pre-miR negative control#2 (Ambion) and pre-miR-132 (Ambion) were used as control. Significance was estimated by performing a multiple comparison test using one-way ANOVA.

Whole-mount in situ hybridization

Whole-mount in situ hybridization for Mest and miR-335 were performed as previously described (Lickert et al., 2002). All embryos were stained using BM Purple (Roche) according to the manufacturer’s instruction. The miR-1 probe was used as a positive control and scrambled miR probe was used as a negative control. Embryos were photographed using a Zeiss Stereo Lumar V12 microscope.

Array data preprocessing and analysis

To study large-scale effects of intronic miRs in ESC development, a previously published (Hailesellasse Sene et al., 2007) study of three cell lines differentiated into embryoid bodies (EB) was used. The datasets were downloaded from the GEO (Barrett et al., 2011) database (GSE2972, GSE3749, GSE3231). For expression profiling, total RNA was measured at eleven time points and three replicates at t=0 hours, 6 hours, 12 hours, 18 hours, 24 hours, 36 hours, 48 hours, 4 days, 7 days, 9 days and 14 days. Affymetrix raw data was preprocessed and expression values were calculated using MATLAB’s RMA (Irizarry et al., 2003) implementation. Probe sets with a total expression below 100 in all data points and probe sets with a variance less than the 10th percentile were removed from the data. We calculated Pearson correlation score based on single probe sets for all three ESC lines. The fold change of Mest, Foxa2 and Sox17 mRNA was calculated for the R1 ESC line (GSE2972) only.

IF image analysis

IF images were taken from in vitro differentiated ESCs under endoderm conditions after 0, 48 and 96 hours. Images were preprocessed to improve cell identification. First, raw images were filtered using a two-dimensional digital FIR filter averaging in a 10×10 window. Image contrast was increased by adjusting pixel values such that 1% of the data is saturated at low and high intensities (supplementary material Fig. S4). Cells and boundaries were identified using image segmentation based on the maximally stable extremal regions (MSER) algorithm (Matas et al., 2002). For each identified cell, we calculated mean intensity by averaging over the raw image IF intensities, as adjusted intensities were not comparable between distinct images. Significance was estimated using a two-sample Kolmogorov-Smirnov test.

Modeling the gradient shaping effect of miR-335

In order to study the miR-mediated effect on the TF gradient during embryogenesis, we developed a dynamic mathematical model based on ordinary differential equations (ODE). The core model is derived from the miR-mRNA binding models proposed by Levine (Levine et al., 2007) and Mukherji (Mukherji et al., 2011). In contrast to these models, we here allow for miR turnover to be of similar scale as gene expression. The model describes the dynamics of the free target mRNA (foxm) the formation and dynamics of the complex (Cmlx) formed by the target mRNA and the miR and the free miR itself. To differentiate between mesodermal and endodermal differentiation, we generated two slightly different models. Both models were identical except for miR transcription: for the mesodermal model, RNA molecules are transcribed from DNA with constant rates and functions assuming linear expression for Foxa2 and miR-335, as deduced from the data. For endodermal differentiation, we also assume linear mRNA expression but a temporal expression for the miR. Thus, miR expression is modeled using a bell-shaped (Gaussian) expression function with an estimated maximum at time=0. The model consists of four ordinary differential equations: Embedded Image Embedded Image Embedded Image Embedded Image Embedded Image

Furthermore, the model assumes that only free mRNA can be translated into protein (foxP). Furthermore, we assume independent degradation rates for the free mRNA (gf) and for the protein (gC), and that miR-mediated degradation of the complex leads to a partial degradation of the miR. The system of ODEs was solved numerically by applying the CVODES solver (Hindmarsh et al., 2005) compiled as C-executable for MATLAB together with the ODE equations. As biochemical rates are not available, we performed parameter estimation by maximum likelihood estimation for three observables (mRNA and miR concentrations measured by qPCR and protein concentrations measured by fluorescence imaging) and a total of 193 data points (supplementary material Appendix S1). The profile likelihood was calculated applying the MATLAB implementation of the trust-region method (lsqnonlin). We used a log-normal error model to calculate the likelihood, including an individual noise parameter for each experimental measurement (FACS/PCR, Image Analysis). For further model details, estimated parameter, profile likelihood estimation and confidence intervals see supplementary material Tables MM1 and MM3, and Figs MM1 and MM2 in Appendix S1.

Foxa2 and Sox17 gradient estimation in tetraploid embryos

Transcription factor gradients were estimated form whole-mount immunofluorescence stainings of gastrulating embryos at E7.5. For each image, we manually selected a polygon mask covering the entire embryo (Fig. 7A, middle row). Images were normalized by subtracting the overall mean and dividing each pixel by the median intensity. To estimate the gradient for each embryo, we calculated the mean IF intensity along the A-P axis (right to left image axis) (supplementary material Figs S6, S7). Gradients were normalized to the interval [0, 1] and binned into five bins along to the A-P axis. For each bin, mean intensity was calculated. We then calculated the gradient for the middle bin referring to the gradient at the middle of the A-P axis (supplementary material Fig. S8). Significant differences between the gradients were estimated using ANOVA on binned intensities (supplementary material Fig. S9).

Acknowledgments

We thank Wenke Barkey, Anne Theis and Heide Oller for excellent technical support, and Andreas Raue as well as Sabrina Hock for fruitful discussion and helpful suggestions on the modeling part. We are also thankful to David L. Turner and Gordon Keller for generously providing us with plasmid constructs.

Footnotes

  • Competing interests

    The authors declare no competing financial interests.

  • Author contributions

    H.L., I.B., D.Y. and L.U. conceived and designed the experiments. D.Y., I.B. and LU performed the experiments. D.L. analyzed the data and performed bioinformatic analysis. D.L. and F.J.T. designed and tested the mathematical model. H.L. and D.L. wrote the paper. F.J.T. and H.L. co-directed the study.

  • Funding

    This work was supported by the Initiative and Networking Fund of the Helmholtz Association within the Helmholtz Alliance on Systems Biology (project ‘CoReNe’), the DFG priority program SPP1356, ‘Pluripotency and Cellular Reprogramming’ and the European Union with the ERC starting grants LatentCauses and CiliaryDisease dedicated to F.J.T. and H.L., respectively. H.L. thanks the Helmholtz Society, German Research Foundation and German Center for Diabetes Research (DZD e.V.) for financial support. This work was funded (in part) by the Helmholtz Alliance ICEMED - Imaging and Curing Environmental Metabolic Diseases, through the Initiative and Networking Fund of the Helmholtz Association.

  • Supplementary material

    Supplementary material available online at http://dev.biologists.org/lookup/suppl/doi:10.1242/dev.104232/-/DC1

  • Received September 26, 2013.
  • Accepted November 11, 2013.

References

View Abstract