Temporal Control of the TGF-β Signaling Network by Mouse ESC MicroRNA Targets of Different Affinities.

Although microRNAs (miRNAs) function in the control of embryonic stem cell (ESC) pluripotency, a systems-level understanding is still being developed. Through the analysis of progressive Argonaute (Ago)-miRNA depletion and rescue, including stable Ago knockout mouse ESCs, we uncover transforming growth factor beta (TGF-β) pathway activation as a direct and early response to ESC miRNA reduction. Mechanistically, we link the derepression of weaker miRNA targets, including TGF-β receptor 1 (Tgfbr1), to the sensitive TGF-β pathway activation. In contrast, stronger miRNA targets impart a more robust repression, which dampens concurrent transcriptional activation. We verify such dampened induction for TGF-β antagonist Lefty. We find that TGF-β pathway activation contributes to the G1 cell-cycle accumulation of miRNA-deficient ESCs. We propose that miRNA target affinity is a determinant of the temporal response to miRNA changes, which enables the coordination of gene network responses.

INTRODUCTION

A mechanistic understanding of how embryonic stem cells (ESCs) maintain and exit the pluripotent state provides insight into development and directed ESC differentiation for cell-based therapies. In transitions from the pluripotent state, RNA post-transcriptional regulation is important in controlling the rapid changes in the ESC transcriptome (Tiscornia and Izpisúa Belmonte, 2010). However, the systems view of how post-transcriptional and transcriptional regulatory networks contribute to changes in pluripotency is still incomplete.

Argonaute (Ago)-microRNA (miRNA) complexes function in a key layer of post-transcriptional gene regulation by reducing target RNA levels or by suppressing translation (Bartel, 2018; Jonas and Izaurralde, 2015). The miRNA seed region (nucleotides 2–8) determines targets for repression through direct base pairing with complementary sequences in target RNAs (Bartel, 2009). Several factors determine the extent of Ago-miRNA repression, including the cellular concentration of the miRNA seed family and the miRNA-target site binding affinity, which is related to the extent of target complementarity (shorter seed matches represent lower-affinity miRNA target sites than longer seed matches). How Ago-miRNA complexes contribute to dynamic gene network responses temporally has not been explored at the systems level. A genome-wide characterization of transcriptome responses to changes in Ago-miRNA levels may reveal additional principles underlying target organization and how they function together to control gene expression patterns.

An evolutionarily conserved miRNA seed family, referred to here as the mouse miR-294 seed family, is highly expressed in ESCs (Houbaviy et al., 2003) and acts as a master regulatory factor shaping the ESC transcriptome (Suzuki et al., 2017). The miR-294 seed family regulates several ESC processes, including cell-cycle progression and response to cellular signals that determine lineage specification (Greve et al., 2013). More specifically, this seed family represses regulators of the G1- > S checkpoint, such as p21 (Cdkn1A/Cip1/Waf1) (Subramanyam et al., 2011; Zheng et al., 2011), and antagonizes the transforming growth factor beta (TGF-β) pathway (Choi et al., 2007; Rosa et al., 2009). In ESCs, TGF-β signaling functions in both self-renewal and differentiation. TGF-β receptor signaling leads to activated Smad transcription factors that execute a transcriptional response (Massagué, 2012). Upon strong activation of the TGF-β pathway, ESCs transition toward mesendoderm differentiation. Interestingly, Ago-miRNA regulation of TGF-β pathway genes in embryos is conserved, but targets are species specific, regulating both activating receptors and pathway antagonist Lefty proteins. The misregulation of these TGF-β pathway miRNA targets alters lineage specification during differentiation, both in vitro and in vivo, supporting important roles in development (Choi et al., 2007; Rosa et al., 2009). How multiple Ago-miRNA regulatory interactions within the TGF-β signaling pathway can function together to control cell signaling activity is still unclear.

Here, we use a series of Ago genetic mutants, a time course of Ago2 loss, and Ago2 rescue at different levels to analyze the early and direct consequences of Ago-miRNA reduction in mouse ESCs. Through the characterization of the transcriptome, histone modification, and RNA polymerase II (Pol II) changes, we find that the TGF-β signaling network responds sensitively to reduction in ESC Ago-miRNA levels. We also identify miRNA target affinity as an important factor in determining the temporal dynamics in response to changes in Ago-miRNA levels. Overall, our study offers insight into the coordinated transcriptional and post-transcriptional regulatory network dynamics in response to changes in miRNA levels.

RESULTS

Reduced miRNAs and miRNA-Mediated Repression with Ago Protein Reduction

Based on the characterization of a mouse ESC mutant series consisting of AB2.2 (wild type, designated AgoWT here), B9 (Ago1−/−; Ago3−/−; Ago4−/−, designated Ago2only here), and E7 (Ago1−/−; Ago2−/−; Ago3−/−; Ago4−/−; hsAgo2; CreERT2; BsdS, designated Ago2low here) cells (Figure 1A) and a failure to isolate stable Ago1–Ago4 knockouts (AgoKO), Ago-miRNAs were previously proposed to suppress apoptosis in mouse ESCs (Su et al., 2009). Because this AgoKO viability phenotype was distinct from other miRNA biogenesis mutants (Leung et al., 2011; Wang et al., 2007), all of which lack miRNA-mediated repression and are viable, we set to further characterize the Ago mutant genetic series.

Figure 1.

Ago Proteins Are Required for miRNA Stability and miRNA-Mediated Repression

(A) Overview of Ago-deficient mouse ESCs described by Su et al. (2009) and AgoKO cells derived in this study.

(B) Western blot of Ago proteins in wild type (AgoWT), Ago-deficient cell lines (Ago2only and Ago2low), and AgoKO ESC clones (AgoKO_C1 and AgoKO_C2). hsAgo2 is a human-specific Ago2 antibody to detect the transgene.

(C) Northern blot for three miRNAs in AgoWT and mutant mouse ESCs.

(D) Quantification of log2 fold changes of mature miRNAs (biological replicates = 2) by small RNA sequencing (p value indicated for comparison with no fold change, one-sample t test; see STAR Methods).

(E) Population miRNA fold repression from the single-cell, two-color fluorescent reporter assay of a perfectly complementary synthetic miR-20a target in indicated cell lines. Error bars indicate SD of biological replicates, n = 3.

(F) Cumulative fraction plots of log2 fold changes of miR-294 family targets of different types (8-mer targets, red; 7-mer targets, orange; 6-mer targets, green) and non-targets (black) in Ago2only/AgoWT (left panel), Ago2low/AgoWT (middle panel), and AgoKO/AgoWT (right panel). The number of genes in each set is indicated in parentheses in the left panel (p values indicate comparison with non-targets, Wilcoxon rank-sum test). Biological replicates = 2.

See also Figure S1.

Because miRNA-deficient cells are more susceptible to cellular stress (Zheng et al., 2011), we hypothesized that prolonged 4-hydroxytamoxifen (4OH-T) treatment stressed cells after loss of the human Ago2 (hsAgo2) transgene in Ago2low cells, resulting in cell death. We found that shortened 4OH-T treatment yielded high-frequency (21 of 76, or 27.6%) AgoKO clones (data not shown). Of the 21 AgoKO clones isolated, 2 were chosen at random for detailed characterization (designated AgoKO_C1 and AgoKO_C2). We confirmed hsAgo2/BsdS transgene cassette loss in AgoKO clones by testing blasticidin sensitivity (data not shown), PCR screening (Figure S1A), and western blot analysis (Figure 1B). We performed growth curves for Ago2low and AgoKO cells and found that 4OH-T inhibited growth (Figure S1B) and probably interfered with previous attempts at isolation of AgoKO cells. These AgoKO cells provided a unique additional cell line to characterize Ago-miRNA regulation.

We next determined the relationship between Ago protein and cellular miRNA levels. By northern blotting for miR-295, miR-20a, and miR-16, we found AgoKO cells lack mature miRNAs (Figure 1C). An intermediate miRNA reduction was found in Ago2low cells. Specifically, a cross-reactive mouse and hsAgo2 antibody, previously determined to recognize each at similar efficiencies in western blots (Zamudio et al., 2014), indicated the Ago2low hsAgo2 transgene is expressed at lower levels than endogenous mouse Ago2 in AgoWT cells (Figure 1B, Ago2 antibody). The reduced Ago2 levels in Ago2low cells occurred with an average 3.8-fold reduction of mature miRNAs compared with AgoWT cells, as determined by northern blotting (Figure 1C). We confirmed that miRNA levels changed globally in Ago2low using small RNA sequencing (Figure 1D, p = 7.1e–13; Figure S1C), which indicated a global average log2 decrease by 2.7-fold in Ago2low compared with AgoWT cells. From these two measures of miRNA levels, we conclude Ago2low cells express between 15% and 26% of wild-type Ago-miRNA complexes. Conversely, no significant miRNA expression difference was detected between Ago2only and AgoWT cells (Figures 1C and 1D, p = 0.15; Figure S1C), suggesting either that Ago2 predominantly stabilizes miRNAs in the presence of the other three Ago homologs or that Ago2 compensates for loss of Ago1, Ago3, and Ago4.

We next used a single-cell fluorescence-based miRNA reporter plasmid that simultaneously expressed mCherry transcripts with perfectly complementary miR-20 binding sites and eYFP transcripts for expression normalization to quantify Ago-miRNA repression at reduced Ago-miRNA levels (Figures S1D and S1E; Mukherji et al., 2011). With this assay, we confirmed that Ago2 alone was sufficient to functionally maintain miRNA-target repression in Ago2only cells (Figure 1E). The repression in Ago2low cells was ~36% of that in AgoWT cells, consistent with a reduced number of Ago2-miRNA complexes in this cell line. We observed that AgoKO cells lacked the ability to repress miRNA targets similar to Dicer1 knockout mouse ESCs; however, in contrast to Dicer1 knockouts (Figure S1F), repression in AgoKO was not rescued by a miRNA mimic (Figure S1G). This indicated a full loss of miRNA repression in AgoKO cells and that Ago proteins are required for repression.

We also found reduced miRNA repression in our quantification of endogenous miRNA targets by mRNA sequencing (mRNA-seq). Our target classification was based on Ago2 individual-nucleotide resolution UV crosslinking and sequencing (iCLIP-seq) in this genetic background (Bosson et al., 2014; Figure S1H; see STAR Methods). As shown for miR-294 seed family targets, Ago2low and AgoKO cells displayed significant target derepression compared with AgoWT cells (Figure 1F). Specifically, the log2 median fold changes in comparison to AgoWT for higher-affinity 8-mer, medium-affinity 7-mer, and lower-affinity 6-mer miR-294 target genes were 0.33, 0.20, and 0.12, respectively, in Ago2low cells and 0.52, 0.26, and 0.18 in AgoKO cells. Considering the average median fold changes of these targets in Ago2low compared with full derepression in AgoKO cells, Ago2low cells maintained ~31% of AgoWT repression. miR-294 targets of all target site affinities measured did not change significantly in Ago2only cells compared with AgoWT, consistent with measures of Ago-miRNA levels. Overall, these data support that decreasing levels of Ago-miRNA complexes result in progressive loss of target repression.

The TGF-β/Smad Target Network Is Activated with Partial Ago-miRNA Reduction

The characterization of Ago2low cells indicated an intermediate state controlled by partial Ago-miRNA function. We reasoned that genomic comparisons of Ago2only, Ago2low, and AgoKO cells may reveal an initial gene network signature that accompanies Ago-miRNA reduction. Generally, the overall number of significant (adjusted p < 0.05) differentially expressed genes increased with decreased Ago-miRNA complexes: Ago2only/AgoWT (21 up, 38 down), Ago2low/AgoWT (156 up, 186 down), to AgoKO/AgoWT (550 up, 361 down) (Figure S1I). Because of the sustained miRNA-mediated repression, Ago2only cells had minimal transcriptome changes compared with AgoWT cells, and pathway analysis did not indicate any specific network response (data not shown). We therefore focused on the Ago2low/Ago2only comparison as an indicator of the response to a partial Ago-miRNA decrease and AgoKO/Ago2low for the subsequent effect of complete Ago-miRNA loss.

Overall, we observed a larger range of mRNA fold changes in the AgoKO/Ago2low comparison than in the Ago2low/Ago2only comparison (SD of log2 fold changes, 0.70 and 0.53, respectively; Figure 2A). Interestingly, we found well-characterized TGF-β pathway targets displayed large changes in the Ago2low/Ago2only comparison. Two TGF-β pathway targets, goosecoid (Gsc) (Blum et al., 1992) and Lefty2, were among the highest upregulated genes in the Ago2low/Ago2only comparison (Figure 2A). Although Gsc, a known direct TGF-β target activated in mesendoderm differentiation, is upregulated in Ago2low cells, neither Ago2low nor AgoKO cells indicated differentiation with overall expression patterns similar to the naive pluripotent state (Chen et al., 2018; Figure S1J). This observation of maintained pluripotency is consistent with previous studies indicating miRNAs are required to silence the pluripotency program (Wang et al., 2007).

Figure 2.

TGF-β/Smad3 Network Is Upregulated with Partially Reduced Ago-miRNA Levels

(A) Density plot of log2 fold changes in Ago2low/Ago2only (x axis) and AgoKO/Ago2low (y axis) shown as gray levels. Contour lines marking the 50 and 5 gene density levels of 8-mer, 7-mer, and 6-mer targets of active miRNA families are shown in red, orange, and green, respectively. Circles of highlighted genes are colored according to their longest miRNA seed match or in black for non-targets. Biological replicates = 2.

(B) Table of significantly changing target sets of active miRNA families (p < 0.001, Wilcoxon rank-sum test; see STAR Methods).

(C) Heatmaps of average log2 fold changes (left) and associated p values (right) of indicated Hallmark and Smad3 target gene sets that changed significantly (p < 0.001, Wilcoxon rank-sum test, outlined in orange; see STAR Methods) in a comparison indicated on the x axis. The total number of significantly up- and downregulated gene sets in a comparison is indicated above the heatmap.

(D) Proportions of genes in a gene set that are classified as partial loss up (red), partial loss down (blue), complete loss up (light red), and complete loss down (light blue) according to their relative log2 fold change at Ago2low/Ago2only (for partial loss) or AgoKO/Ago2low (for complete loss) compared with their maximum fold change (see STAR Methods).

To further investigate this transcriptome response, we analyzed the behaviors of diverse gene sets that included 13 Ago2 iCLIP-defined miRNA-regulated networks (see STAR Methods), 50 Hallmark gene sets, which are based on multiple experimental data to represent important molecular signatures (Liberzon et al., 2015), and as direct TGF-β pathway targets, two gene sets defined by proximal Smad3 binding either before or after activin treatment in mouse ESCs (Mullen et al., 2011; see STAR Methods). In the Ago2low/Ago2only comparison, we found all 13 miRNA-regulated networks were significantly (p < 0.001) upregulated (Figure 2B). In the AgoKO/Ago2low comparison, only 2 of these 13 miRNA-regulated networks changed significantly. These 2 target sets were the highly expressed ESC miR-294 and miR-292 seed families, suggesting that in Ago2low cells, these seed families maintained some repression on their target pools. Only 8 miRNA target gene sets were classified as derepressed significantly in the AgoKO/Ago2only comparison. This may be a consequence of a concurrent transcriptional response in AgoKO, which overlaps the miRNA derepression signature (described later). We found 11 of 52 Hallmark and TGF-β/Smad3 target gene sets responded significantly at a stringent cutoff (p < 0.001) in Ago2low/Ago2only and AgoKO/Ago2low (Figure 2C). Strikingly, the two TGF-β/Smad3 target sets represented the most significantly differentially expressed sets in the Ago2low/Ago2only comparison (p < 1e–45 for both TGF-β/Smad3 target sets). TGF-β/Smad3 set upregulation was not significant in the AgoKO/Ago2low comparison (p > 0.02) (Figure 2C). This gene set behavior indicated initial TGF-β pathway activation following partial reduction of Ago-miRNA complexes that was maintained in the activated state after complete Ago protein loss. In contrast, the Myc target network, represented by two Hallmark gene sets, Myc targets v1 and v2 (defined by different experiments), was significantly decreased in both Ago2low/Ago2only and AgoKO/Ago2low comparisons (p = 2.6e–28 and 5.9e–9, respectively, for Myc target v1), indicating progressive downregulation with Ago-miRNA loss. This is in agreement with Myc being among the most strongly downregulated genes in the Ago2low/Ago2only and AgoKO/Ago2low comparisons (Figure 2A) and being part of a previously reported transcriptional response to loss of miR-294 seed family regulation (Melton et al., 2010). We identified other previously reported pathway responses to loss of ESC miRNAs; for example, a switch from glycolytic metabolism to oxidative phosphorylation that is linked to ESC differentiation (Cao et al., 2015) occurred only in the AgoKO/Ago2low response, possibly controlled less sensitively or indirectly by the miR-294 seed family.

To further characterize the gene network responses, genes within differentially expressed gene sets were classified based on how they responded to Ago-miRNA loss: up/down with partial Ago-miRNA loss (in Ago2low/Ago2only) or up/down with complete loss (in AgoKO/Ago2low) (Figure 2D; see STAR Methods). After ranking these pathways based on the fraction of genes that were upregulated with partial Ago-miRNA loss, we found TGF-β/Smad3 targets were top ranked, with almost 40% upregulated genes. This skew toward TGF-β/Smad3 target upregulation was not detected with complete loss. In contrast, Myc target sets demonstrated sustained downregulation upon partial and complete Ago-miRNA loss. We concluded that widespread changes to Ago-miRNA loss include transcriptional responses with an initial activation of the TGF-β/Smad3 target network and an initial and progressive downregulation of the Myc network.

Transcriptome, Pol II, and Histone Modification Profiling Support Transcriptional Responses of TGF-β/Smad3 and Myc Targets to Ago-miRNA Loss

To confirm transcriptional changes with Ago-miRNA loss, we calculated gene fold changes separately from exon and intron read coverage, which are established quantitative estimators for the overall (log2 exon fold change), transcriptional (log2 intron fold change), and post-transcriptional (log2 exon fold change – log2 intron fold change) changes (Gaidatzis et al., 2015). We first examined these changes in an assembled panel of well-characterized ESC-miRNA responsive genes, including direct targets and secondarily responsive genes (Figure 3A). The genes selected are involved in cellular processes that promote ESC differentiation by regulating the cell cycle (p21, Lats2, and Rbl2) (Pauklin and Vallier, 2013), the TGF-β pathway (Lefty1, Lefty2, Tgfbr1, and Acvr2b), the epithelial-to-mesenchymal transition (Cdh1 and Fn1), transcription regulation (Myc, Mdb2, Dnmt3a, and Dnmt3b), and cellular metabolism (Ldha and Pkm) (Cao et al., 2015). We observed that transcriptional (intronic) and post-transcriptional (exon-intron) changes contribute to the overall exon changes, suggesting a coordinated response at genes, including well-characterized miR-294 targets p21 and Lefty mRNAs. We next analyzed transcriptional changes network-wide and observed upregulated intron levels for TGF-β/Smad3 targets (p = 4.6e–10) and downregulated levels for Myc targets (p = 6.6e–3) in Ago2low/Ago2only, supporting a transcriptional response (Figure 3B). Overall, we found 5 of 11 pathways had a significant underlying transcriptional response in Ago2low/Ago2only, indicating a substantial impact on transcription (Figure S2A).

Figure 3.

TGF-β/Smad3 Network Activation with Ago-miRNA Loss Occurs with Transcriptional, Pol II Occupancy, and Histone Modification Changes

(A) Heatmap of log2 fold changes in exon levels (1st panel), intron levels (2nd panel), and exon-intron levels (3rd panel) for selected genes. miRNA target-type information and Smad3 binding are indicated by color code in the last panel. Numbers next to gene names indicate the number of 8-mer:7-mer:6-mer seed matches in the 3′ UTR (see STAR Methods).

(B) Comparison of log2 intron fold changes in Ago2low/Ago2only for Smad3 targets (purple), Myc targets (light blue), and all genes (black) (p values indicated from Wilcoxon rank-sum test; see STAR Methods). Biological replicates = 2.

(C) Cumulative fractions of log2 fold changes in Pol II binding (left) and histone modifications (H3K27me3, middle panel; H3K27ac, right panel) at Smad3-bound (purple) and Myc-bound (light blue) regions for Ago2low/Ago2only (top panels) and AgoKO/Ago2low (bottom panels) comparisons (see STAR Methods). The number of included regions is indicated in parentheses. The p values are indicated from comparison with all ChIP-seq peaks (black; Wilcoxon rank-sum test).

(D) Genome browser view of normalized (see STAR Methods) ChIP-seq read distributions for H3K27ac and H3K27me3 modifications and Pol II occupancy in different Ago mutants at the Lefty2/Lefty1 genomic locus. The positions of statistically significant Smad3 ChIP-seq peaks are indicated in purple at the top.

See also Figure S2.

To further characterize these transcriptional responses, we generated genome-wide profiles of Pol II occupancy and histone H3 lysine 27 acetylation (H3K27ac), histone H3 lysine 27 trimethylation (H3K27me3), and histone H3 lysine 9 trimethylation (H3K9me3) modifications in Ago2only, Ago2low, and AgoKO cells and confirmed the expected correlations with RNA expression levels (Figures S2B and S2C). In the Ago2low/Ago2only comparison, we identified TGF-β/Smad3 targets as the most significantly enriched sets among genes associated with changes in activating H3K27ac (up, p < 1e–23; down, p < 1e–25), as well as increased Pol II binding (up, p < 1e–9; Figure S2D). In addition, targets with decreased repressive H3K27me3 were enriched for proximal Smad3 binding (p < 6e–4), indicating a subset of TGF-β/Smad3 targets are co-regulated by the Polycomb complex (Figure S2D). Overall, this unbiased analysis supported an underlying transcriptional response for gene sets identified based on intron-level changes.

We next examined Pol II, H3K27ac, and H3K27me3 coverage directly at genomic binding sites of Smad3 and Myc (Figure 3C; see STAR Methods). Pol II occupancy was increased at Smad3-bound regions in the Ago2low/Ago2only and AgoKO/Ago2low comparisons (average log2 fold changes, 0.13 and 0.23, respectively; p value compared with all Pol II-bound regions, 7e–27 and 1e–13, respectively), while at Myc-bound regions, a significant decrease in Pol II occupancy was detected in the Ago2low/Ago2only comparison. In support of a subset of Smad3-bound regions losing the repressive H3K27me3 modification, we found a decreased H3K27me3 signal in the Ago2low/Ago2only comparison (p = 2e–3) (Figure 3C, top middle). For H3K27ac changes, Smad3-bound regions displayed a wider range of changes (SD 0.57) than Myc regions (SD 0.32) in the Ago2low/Ago2only comparison, which may reflect Smad3-mediated transcriptional activation and repression. The H3K27ac levels at Smad3 regions were more skewed toward upregulation in the AgoKO/Ago2low comparison (Figure 3C, right bottom).

A genome browser view of the Lefty gene locus provides an example for increased H3K27ac, increased Pol II, and decreased H3K27me3 at TGF-β/Smad3 targets with progressive Ago-miRNA loss (Figure 3D). At each likely intragenic enhancer indicated by the Smad3 and H3K27ac chromatin immunoprecipitation sequencing (ChIP-seq) signal upstream of Lefty genes, there was an accumulation of H3K27ac and Pol II with decreased cellular Ago-miRNA levels. The repressive H3K27me3 mark was found within the Lefty genes and was decreased with reduced cellular Ago-miRNA levels. The profiled factors were not altered at the Pycr2 gene located between Lefty2 and Lefty1 genes, supporting specificity of the transcriptional response (Figure 3D).

TGF-β/Smad3 Targets Are Enriched for miR-294 Seed Family-Directed Post-transcriptional Repression

From our assembled panel of Ago-miRNA responsive genes (Figure 3A), we noted two additional features. First, Smad3 binding was frequently detected proximal to well-characterized miR-294 responsive genes, including p21 and Myc (Figure 3A; Figure S2E). Second, post-transcriptional derepression was partitioned differently between partial and complete Ago-miRNA loss. Specifically, most strong 8-mer miR-294 targets showed a larger subsequent derepression with complete Ago-miRNA loss (Rbl2, p21, Lefty1, and Lats2), while the miR-294 7-mer target TGF-β receptor 1 (Tgfbr1), an upstream regulator of the TGF-β signaling pathway, showed comparable post-transcriptional derepression to partial and complete Ago-miRNA loss. These observations suggested a coordinated regulation by TGF-β/Smad3 and miR-294 and distinctive miRNA target responses, possibly based on target affinity.

To examine the direct Ago-miRNA post-transcriptional control on TGF-β network activation, we analyzed Ago2 iCLIP targets with miR-294 seed matches, including Lefty1, Lefty2, and Tgfbr1 (Figure 4A). Interestingly, we found that TGF-β/Smad3 targets were significantly enriched in miR-294 targets, and ~15% had an Ago2 iCLIP signal at their respective miR-294 seed matches (p = 1.3e–6) (Figure 4B, top). Only three additional pathways were enriched in direct miR-294 seed family targeting (UV response, p = 3.3e–3; G2M checkpoint, p = 1.4e–5; E2F targets, p = 4.2e–3), likely related to the miR-294 seed family regulation of the ESC cell cycle (Figure 4B, bottom). To verify functional regulation, we identified pathways that exhibited significant post-transcriptional (log2 exon fold change – log2 intron fold change) fold changes. TGF-β/Smad3 target sets were significantly post-transcriptionally derepressed (p < 1e–17) in Ago2low/Ago2only (Figure 4C), but not in AgoKO/Ago2low. Thus, Ago-miRNAs post-transcriptionally repress the TGF-β/Smad3 target network, and derepression upon partial miR-294 seed family loss correlates with the TGF-β/Smad3 network transcriptional activation. The transcriptional and post-transcriptional responses of TGF-β/Smad3 target genes were also detected upon knockout of miRNA-biogenesis factors Dicer1 and Dgcr8 (Leung et al., 2011; Wang et al., 2007), supporting the identification of a miRNA-mediated control of the TGF-β network response (Figure S3).

Figure 4.

TGF-β/Smad3 Targets Are Enriched for miR-294 Targeting and Are Post-transcriptionally Derepressed with Partial Ago Loss

(A) Genome browser views of Ago2 iCLIP-seq (TTFHAgo2-iCLIP, first profile in each panel) and background (TTAgo2-iCLIP, second profile in each panel) read distributions across Lefty1 (top), Lefty2 (middle), and Tgfbr1 (bottom) mRNAs. The positions of miRNA target sites are indicated and colored according to the target site type (8-mers, red; 7-mers, orange; 6-mers, green).

(B) Top: overlap between direct TGF-β/Smad3 target genes identified in activin-stimulated cells and miR-294 family targets (p value indicated from hypergeometric test). Bottom: fraction of genes in a gene set that are iCLIP-supported miR-294 family targets separated by target type. Gene sets significantly overlapping (hypergeometric test) with miR-294 family targets are marked with asterisks.

(C) Heatmaps of log2 exon-intron fold changes (left) and associated p values (right) for indicated Hallmark and Smad3 target gene sets that were significantly changed (p < 0.001, Wilcoxon rank-sum test, outlined in orange; see STAR Methods) in a comparison indicated on the x axis. Biological replicates = 2. See also Figure S3.

TGF-β Network Response Occurs Early during Progressive Loss of Ago-miRNA Complexes, Concurrent with a Greater Derepression of Weaker miRNA Targets

We next profiled gene expression changes during progressive Ago-miRNA reduction using doxycycline (dox)-inducible FLAG-hemagglutinin (HA)-tagged Ago2 (TTFHAgo2) cells (Zamudio et al., 2014) to establish a temporal relationship between network responses. TTFHAgo2 cells allow controlled Ago2 protein expression ranging from wild type to levels nearly undetectable by western blot (Figure S4A). For the Ago2 removal time course, TTFHAgo2 cells were induced for wild-type Ago2 expression (2.5 μg/mL of dox for 48 h) and subsequently dox starved, resulting in progressively decreasing Ago2 levels. RNA was collected 0, 24, 48, and 96 h after dox removal. As quantified by northern blot, the average reduction in miRNA levels was ~5-fold and ~20-fold at 24 and 48 h, respectively (Figure 5A). miRNA reduction was confirmed globally by small RNA sequencing (Figure S4B). Through mRNA-seq, we found that overall, the time series agreed well with the Ago mutant genetic series, with expression fold changes after 96 h closest to changes in AgoKO cells (Figure S4C).

Figure 5.

TGF-β Signaling-Related Gene Sets Are Enriched in Early-Responding and Lower-Affinity Ago-miRNA Targets

(A) Northern blot for miRNAs in the Ago-loss time course (++, 2.5 μg/mL of Dox for wild-type induction; −, no Dox).

(B) Boxplots of relative (left panels) and actual (right panels) log2 fold changes in 24/0 h (left in each panel), 48/24 h (middle in each panel), and 96/48 h (right in each panel). Relative fold changes were calculated relative to the maximum log2 fold change of a gene (see STAR Methods) for miR-294 family targets of different affinities. The median relative fold change (%) is indicated on top of each box. Significant differences (Wilcoxon rank-sum test) in relative fold changes between miRNA target affinity groups are indicated by asterisks.

(C) Left: classification of active miRNA family targets into late-, intermediate-, and early-responding genes according to their relative post-transcriptional derepression at 24 h (see STAR Methods). The number of genes in each class is indicated in parentheses. Right: fraction of 8-mer, 7-mer, and 6-mer targets among early-, intermediate-, and late-responding genes. The p value (Fisher’s exact test) indicates a significantly larger proportion of 8-mer targets among late-responding than among early-responding miRNA target genes.

(D) Significantly (p < 0.05, Fisher’s exact test) enriched Hallmark, KEGG, and Reactome pathways and Smad3 target genes among early-responding compared with late-responding genes (top) and among late-responding compared with early-responding genes (bottom).

(E) Cumulative fraction plots of log2 fold changes of TGF-β/Smad3 (purple) and Myc (light blue) direct target genes and all genes (black) in the Ago2 loss time course. The number of genes in each set is indicated in parentheses (p values indicated from Wilcoxon rank-sum test; see STAR Methods). Biological replicates = 2.

See also Figure S4.

To track the progressive responses of miRNA targets during the time series, we calculated relative mRNA fold changes for each gene by dividing the log2 fold change at each time point by the maximal derepression in the series, which was measured between wild-type Ago2 expression levels and any subsequent time point (see STAR Methods). We distinguished between weak and strong miRNA targets based on their longest seed match (see STAR Methods). Although it is known that many factors contribute to miRNA targeting, such as secondary structure at the miRNA target site, concentrations of miRNA and target mRNA, or multiplicity of target sites in a 3′ UTR (Agarwal et al., 2015; Denzler et al., 2016), we reasoned that the behavior of 6-mer, 7-mer, or 8-mer miRNA targets generally reflected those of weaker to stronger miRNA targets.

As illustrated for the miR-294 seed family targets (Figure 5B), derepression occurred progressively in the time series, with more repression of 8-mer targets remaining 48 h after Ago2 removal compared with 7-mer and 6-mer targets (median relative derepression between 48 and 96 h: 20%, 12%, and 9% for 8-mer, 7-mer, and 6-mer targets, respectively; p values: 2.5e–3 and 3.6e–4 for miR-294 8-mer versus 7-mer targets and 8-mer versus 6-mer targets, respectively). The miR-294 family targets of different strength also demonstrated differential responses in the time series. Specifically, weaker targets were derepressed by a larger proportion of the maximum fold change at 24 h Ago2 removal (median relative fold change after 24 h: 35%, 45%, and 51% for 8-mer, 7-mer, and 6-mer targets, respectively; p values: 8.3e–3 and 1.1e–5 for 8-mer versus 7-mer targets and 8-mer versus 6-mer targets, respectively) (Figure 5B). This result suggested that in response to decreased Ago-miRNA levels, miR-294 family targets under weaker repression are more responsive and are responsive earlier than stronger targets. This difference was not observed among targets of less abundant miRNA seed families, such as the miR-92 seed family (Figure S4D), possibly because we missed the intermediate target response with our time point selection.

To unbiasedly assess differences between earlier- and later-responding miRNA targets, we classified targets based on their percentage of maximum post-transcriptional derepression at 24 h after dox removal (see STAR Methods), which was large (≥90%), moderate (≥30% and <90%), or small (<30%) corresponding to predominantly early (within 24 h after dox removal), intermediate (within 24 h after dox removal and later), or late (mostly later than 24 h after dox removal) changing genes, respectively (Figure 5C, left). We found that early-responding genes had a significantly smaller proportion of high-affinity 8-mer targets than late-responding genes (Figure 5C, right, p = 6.3e–3). Furthermore, early-responding genes were enriched for TGF-β signaling gene sets, indicating that weaker, more responsive miRNA targets may be responsible for activating the TGF-β pathway response (Figure 5D). TGF-β/Smad3 direct targets were enriched among later-responding genes, suggesting sustained repression during transcriptional activation of some stronger miRNA targets (described later). Overall, the different responses of weak (earlier) and strong (later) miRNA targets to Ago-miRNA loss suggest a function in temporally controlling transcriptional responses.

Importantly, our identification of an initial TGF-β pathway response in the Ago mutant genetic series was confirmed in the Ago-loss time series. TGF-β/Smad3 network transcriptional activation was initiated significantly at 24 h Ago-miRNA depletion (Figures 5E and S4E) with concurrent post-transcriptional fold changes at 24 h Ago-miRNA depletion, indicative of TGF-β/Smad3 target derepression (p < 5e–3) (Figure S4F). This 24 h activation of the TGF-β/Smad3 network was more prominent than the Myc network downregulation, which was strongest in the 48/24 h comparison (Figure 5E), potentially placing TGF-β/Smad3 activation upstream of Myc target downregulation.

We concluded that the Ago-loss time series validated our interpretation of the genetic series response and implicated miRNA target strength as a factor in determining the temporal response to Ago-miRNA loss. Specifically, it revealed that weaker miR-294 targets are proportionally more derepressed earlier during reduction of Ago-miRNAs. These weaker Ago-miRNA targets include those controlling the TGF-β pathway transcriptional activation, such as the Tgfbr1 target (Figures S4G and S4H). Stronger Ago-miRNA targets, such as Lefty proteins, are derepressed more gradually, suggesting a function for miRNAs in dampening concurrent transcriptional induction.

TGF-β Pathway Activation Is Reversible and Requires Wild-Type Levels of Ago-miRNA Complexes for Full Suppression

To test whether the TGF-β pathway response is reversible, we performed Ago2 rescue experiments. TTFHAgo2 cells were dox starved for 96 h and rescued with either low (0.1 μg/mL) or high (2.5 μg/mL, wild-type-level induction) dox levels for 48 h. Ago2 rescue was confirmed by western blot (Figure S5A). We estimated by miRNA northern blotting that miRNA levels in the low-dox condition were 60% of those in wild-type-level induction (Figure 6A), and high-dox treatment yielded an additional nearly uniform increase in miRNA levels (Figures 6A and S5B). Overall, we found that changes of most gene sets significantly responsive to Ago-miRNA reduction were reversible, supporting a direct or secondary control by Ago-miRNA complexes (Figures S5C and S5D).

Figure 6.

Ago2 Rescue Shows Reversibility of TGF-β/Smad3 Response and Reveals a Disproportional Transcriptional Deactivation at Low Ago Induction Levels

(A) Northern blot for three highly expressed microRNAs in the Ago rescue series (++, 2.5 μg/mL [high] of Dox for wild-type induction; +, 0.1 μg/mL [low] of Dox for an intermediate induction; −, no Dox).

(B) Scatterplots of log2 fold changes in exon levels (x axis) and intron levels (y axis) of significantly (adjusted p < 0.01, calculated using DESeq; Anders and Huber, 2010) changing genes in 0.1 mg/mL of Dox versus no Dox (left) and 2.5 mg/mL of Dox versus no Dox (right). Genes circles are colored according to their response: significant exon changes, blue; significant intron change, red; both significant, green. Selected genes are indicated. Biological replicates = 2.

(C) Top: cumulative fraction of log2 fold changes in intron levels of Smad3 (purple) and Myc (light blue) target genes and all genes (black) for 0.1 μg/mL of Dox versus no Dox (left) and 2.5 μg/mL of Dox versus no Dox (right). The number of genes in each set is indicated in parentheses (p values indicated from Wilcoxon rank-sum test; see STAR Methods). Bottom: heatmaps of log2 fold changes (left) and associated p values (right) for selected gene sets in the Ago2 loss time course (24/0, 48/24, and 48/96 h) and Ago2 rescue series (low Dox/no Dox and high Dox/no Dox; p values from Wilcoxon rank-sum test; see STAR Methods).

(D) Top: boxplots of relative post-transcriptional (exon-intron) fold changes at low Ago2 induction for 6-mer, 7-mer, and 8-mer target genes of active miRNA families (see STAR Methods). The p value (Wilcoxon rank-sum test) indicates a significantly larger relative post-transcriptional fold change for 8-mer target genes compared with 6-mer targets at low Ago2 induction. Bottom: boxplots of relative transcriptional (intron) fold changes of selected TGF-β signaling-related gene sets after low Ago2 induction (see STAR Methods).

(E) Model for the responses of stronger and weaker miRNA targets to increasing Ago-miRNA (Ago loss, red arrow) or decreasing Ago-miRNA levels (Ago rescue, blue arrow).

See also Figure S5.

The low-dox rescue condition did not proportionately suppress the TGF-β transcriptional response (Figure 6B). Specifically, after low-dox induction, the medians of log2 intron fold changes of TGF-β/Smad3 direct target sets were ~41% of the changes observed after Ago2 wild-type induction (Figure 6C, left), and larger transcriptional changes occurred between low and full rescue (Figure 6C, right). Myc network suppression was also reversible, displaying a disproportionate response similar to that of TGF-β/Smad3, with low dox activating 45%–49% of the log2 intron fold changes observed after Ago2 wild-type induction. To refine these disproportional Ago-miRNA and transcriptional gene network changes, we specifically examined the relative post-transcriptional changes of miRNA targets and relative transcriptional changes of TGF-β-related gene sets (Figure 6D). For miRNA targets of different affinities, we observed a reverse response from the time series with regard to sensitivity. In this case, stronger targets were repressed to a larger proportion of the maximum repression in the low-Ago2 rescue condition than weaker targets (median relative fold change of 8-mer targets of 75% versus 56% for 6-mer targets, p = 1.1e–5; Figure 6D, top). By comparing the proportions of post-transcriptional repression and transcriptional changes at low Ago2 induction, we found that the relative changes of TGF-β signaling-related networks, which were reverted by 35%–42% in the low-dox condition (Figure 6D, bottom), were closer to those of weaker miRNA targets. Altogether, we conclude that miRNA target site strength contributes to observed gene network responses to Ago-miRNA levels (Figure 6E). During decreases in Ago-miRNA levels, weaker Ago-miRNA interactions contribute to an earlier target derepression, resulting in sensitive responses, while stronger interactions control more gradual derepression, persisting later in depletion. During increases in Ago-miRNA levels, weak Ago-miRNA targets show later responses, while stronger targets are proportionally more suppressed earlier.

A Sensitive miR-294 7-mer Target Site in the Tgfbr1 3′ UTR Controls TGF-β Network Activation, and Lefty Proteins Show Dampened Induction

Based on its behavior in the Ago genetic mutant series (Figure 3A) and the Ago2 loss time series (Figures S4G and S4H), Tgfbr1 appeared to be a strong candidate for an upstream regulator of the TGF-β network, because it was bound by Ago2 at a 7-mer miR-294 seed family target site (Figure 4A) and responded early and predominantly post-transcriptionally (Figures 3A, S4G, and S4H).

To test a role for TGF-β receptor-mediated activation in pathway induction, we simultaneously performed Ago2-miRNA depletion in TTFHAgo2 cells, together with TGF-β pathway inhibition using SB-431542, a small molecule that is a selective receptor inhibitor of TGF-β signaling (Inman and Hill, 2002). We monitored TGF-β-activated Lefty mRNA levels in which early induction was predominately transcriptional (Figure S4G). We confirmed Lefty mRNA expression was induced early and increased progressively with loss of Ago-miRNAs (Figure 7A, fold change at 24 h: 1.7, p = 0.02) and expression was lost with SB-431542 treatment (all time points, p < 0.001). This result indicated that Ago-miRNA complexes control Lefty transcriptional activation in a manner dependent on TGF-β receptor signaling. The receptor-mediated induction was additionally supported by increased Smad2 phosphorylation during the Ago-loss time series (Figure S6A).

Figure 7.

miR-294 Site in TGF-β Receptor 1 Suppresses Pathway Activation, and Lefty Proteins Are Suppressed Early with Ago Loss

(A) Fold changes quantified by qPCR in the Ago2 loss time series in the presence of TGF-β pathway inhibitor SB-431542 (SB) or carrier-only control (DMSO). Fold changes are in comparison to the wild-type-level Dox treatment condition (2.5 μg/mL of Dox) and normalized to Gapdh fold changes. Asterisks indicate statistical significance of the response to TGF-β pathway inhibition (**p < 0.01, ***p < 0.001; two-tailed t test). Biological replicates = 3.

(B) Locations of the sgRNA sites to delete the miRNA target site within the Tgfbr1 3′ UTR. The 7-mer target site for miR-294 is indicated by the yellow bar.

(C) qPCR quantification of expression fold changes for mature Tgfbr1 mRNA, mature Lefty1 mRNA (exon-exon), and Lefty1 pre-mRNA (exon-intron) in Tgfbr1 miRNA-target site knockout ESCs (mutant) compared with wild-type ESCs. Bars and error bars represent mean and SD of three biological replicates from different mutant clones. The p values (two-tailed t test) are indicated by asterisks (**p < 0.01, *p < 0.05).

(D) FACS quantification of Lefty and Ago2 protein levels in single cells using intracellular staining after Ago2-miRNA loss. The percentage of cells in each quadrant is indicated.

(E) Model scheme for the coordinated control of the TGF-β pathway response to Ago-miRNA loss by miR-294 targets of different affinities (weaker target Tgfbr1, orange; strong targets, red). The dashed line indicates weaker repression and solid thick lines indicate stronger repression by Ago-miRNAs.

(F) Boxplots of percentages of cells in the G1 phase, quantified using DNA staining. The percentage of G1 cells increased significantly with Ago2 loss, among cells treated with DMSO (p < 0.01 at 24, 48, and 96 h compared with wild type), and SB (p < 0.05 at 24, 48, and 96 h compared with wild type). At 96 h, the percentage of G1 cells among SB-treated cells is significantly smaller (p < 0.05) than among DMSO-treated cells. The p values were calculated using two-tailed t test. Biological replicates = 3.

See also Figure S6.

To determine whether direct miR-294 targeting of Tgfbr1 mRNA suppresses transcriptional activation of the TGF-β pathway, we used Cas9 to delete the miR-294 family binding site in the Tgfbr1 3′ UTR (Figure 7B). Homozygous deletion mutants were validated by PCR screening (Figure S6B) and Sanger DNA sequencing (data not shown). Using qPCR to profile gene expression, we found that Tgfbr1 miR-294 family binding site mutants had 1.8-fold increased Tgfbr1 levels compared with non-edited controls, supporting functional repression by this binding site (Figure 7C, p < 0.01). This increase in Tgfbr1 was sufficient to increase Lefty1 mRNA as measured by qPCR probe sets for mature mRNA (exon-exon, 2.4-fold, p < 0.05) and nascent RNA (exon-intron, 3.0-fold, p < 0.05). These results supported the role of the miR-294 binding site within the Tgfbr1 mRNA in mediating the sensitive response of the TGF-β pathway to Ago-miRNA depletion.

To verify Lefty upregulation with Ago-miRNA loss at the protein level and in single cells, we performed Lefty and Ago2 intracellular staining and fluorescence-activated cell sorting (FACS) quantification with Ago2-miRNA depletion. Protein changes confirmed the population response, in which a decrease in Ago2 protein increased intracellular Lefty protein levels overall (Figure 7D). However, at 24 h after Ago2 loss, only a small subset of single cells (<3%) displayed increased Lefty protein levels. These unchanged Lefty protein levels in most cells at 24 h were in contrast to the Lefty mRNA increase observed at this time, possibly indicating some contribution of translational repression. Subsequent reductions in intracellular Ago2 protein levels occurred with increased Lefty protein levels at 48 and 96 h after Ago2 loss.

Thus, differential responses to Ago-miRNA loss based on miRNA target strengths allow both a sensitive activation of the TGF-β pathway and an extended repression of certain transcriptionally induced, strong miRNA targets to dampen their output (Figure 7E).

TGF-β Pathway Transcriptional Response Contributes to the Cell-Cycle Defect of miRNA-Deficient ESCs

Because the mouse miR-294 family functions in the control of the ESC cell cycle, we next explored the consequences of TGF-β pathway activation on cell-cycle progression. We focused on the G1-S checkpoint cell-cycle regulator p21, which has proximal Smad3 binding sites (Figure S2F) and responded to Ago-miRNA loss with concurrent transcriptional activation and post-transcriptional derepression in our transcriptome profiling (Figures 3A, S4G, and S4H). Similar to Lefty genes, p21 mRNA has a well-characterized 8-mer miR-294 family target site. We first confirmed that in response to Ago-miRNA depletion, p21 is transcriptionally induced (Figure S6C) and activation occurred independent of signatures of p53 activation (Figure S6D). We found that p21 induction was suppressed by SB-431542 treatment to 63% of the 96 h levels (Figure 7A, p < 0.01), indicating a role for TGF-β-signaling in p21 activation. To determine whether TGF-β signaling contributed to the cell-cycle defect with Ago2-miRNA loss, we profiled cell-cycle stages under Ago-miRNA depletion and simultaneous TGF-β pathway inhibition using DNA staining (Figure S6E). We observed that accumulation of cells in G1 phase continued progressively during the depletion of Ago-miRNA complexes and that suppression of TGF-β signaling decreased the G1 accumulation significantly by ~30% of that observed under DMSO treatment between 0 and 96 h (Figures 7F and S6E, p < 0.05). These results indicated that TGF-β signaling contributes to regulation of the ESC cell cycle, possibly through direct p21 regulation.

DISCUSSION

Using a genetic series of Ago mutants, a time course of Ago-miRNA loss, and an Ago-miRNA rescue series, we provide insight into gene network responses to Ago-miRNA-level changes. Our characterization of these responses indicates that miRNA target affinity is a factor in determining temporal pathway activity. We found that although weak (lower affinity) miRNA targets may confer less overall repression, they displayed relatively more derepression earlier during the Ago-miRNA reduction time course. Such behavior was apparent for upstream TGF-β pathway regulators, including the receptor protein Tgfbr1. The loss of the miR-294 target site in the Tgfbr1 3′ UTR was sufficient to cause upregulation of Lefty mRNA, suggesting an important role of this miRNA family in the temporal regulation of TGF-β pathway activation. In contrast to weaker and more sensitive target sites, the robust repression of strong targets provided repression later during the Ago-miRNA depletion time course. We characterized such repression on the TGF-β pathway downstream target Lefty. Because Lefty proteins are antagonists of TGF-β receptor signaling, a longer window of miRNA-mediated repression may act in dampening their transcriptional activation to possibly delay the negative feedback loop.

Although miRNA regulation of agonist receptor complexes and antagonist TGF-β pathway proteins was previously characterized in development and stem cell differentiation (Choi et al., 2007; Rosa et al., 2009), our characterization indicates that these regulatory interactions can have variable impact at various ESC-miRNA concentrations, and thus different times, during ESC-miRNA decrease. In zebrafish, the TGF-β Nodal agonist Ndr1/squint was shown to have a weaker miR-430 target site than the Lefty antagonist (Choi et al., 2007), and translational repression of Lefty by miR-430 was reported to extend the temporal window for TGF-β signaling to impart spatial patterning information (van Boxtel et al., 2015). This suggests a conserved arrangement of a sensitive TGF-β pathway activation and a delayed induction of Lefty proteins. In addition, TGF-β signaling may alter differentiation through a prolonged ESC G1 cell-cycle stage at which they are more responsive to differentiation cues (Pauklin and Vallier, 2013). Upon increasing Ago-miRNA concentration, weak miRNA targets were proportionally repressed less at lower miRNA levels. This behavior would translate into later responses of weak miRNA targets during ESC-miRNA increases. Such temporal control may function in the final stages of ESC reprogramming.

The role of binding affinity has been proposed in temporal responses in other processes. For example, the iron stress response in a cyanobacteria induces a small RNA proposed to coordinate the temporal response of target genes independent of target site affinity (Legewie et al., 2008). Similar observations have been made for transcription factor induction in other bacteria in response to changes in phosphate levels (Gao and Stock, 2015). There, strong transcription factor targets are activated earlier and impart a specific initial response. In mammalian development, a combination of weak enhancers controls the temporal activation of the transcriptional factor Pax6, supporting a developmental connection between regulatory factor, target affinity, and temporal control of activation (Rowan et al., 2010). In respect to mammalian RNA regulatory factors, the relationship between RNA binding factor concentration and target site affinity has been explored in in vitro binding assays for alternative splicing regulators, where weaker binding sites become bound upon higher regulatory factor concentration, suggesting differential patterns of regulation dependent on concentration of regulatory proteins (Lambert et al., 2014). Here, we propose that the same basic biochemical principle functions with Ago-miRNA levels to determine the temporal responses of target pathways during cell-state transitions. In this context, highly expressed miRNA seed families, such as miRNAs associated with super-enhancers (Suzuki et al., 2017), may mediate more dynamic responses because they confer a larger repression even on weaker targets, as opposed to lower-expressed miRNA families. How combinations of miRNA targets sites within the same transcript may impart specialized responses would be important to characterize in future experiments.

Finally, we refined the phenotype of AgoKO cells and determined Ago null cells are viable and susceptible to cellular stress. The observation that Ago2 alone was sufficient to maintain target repression suggested that either Ago homologs are redundant or Ago2 is primarily responsible for target repression. To address this, additional Ago mutants are required. Our characterization of stable AgoKO cells provides a resource for future studies, permitting direct comparisons between different miRNA pathway mutants and offering a more complete view of miRNA function in mouse ESCs.

STAR★METHODS

LEAD CONTACT AND MATERIALS AVAILABILITY

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Jesse Zamudio (jesse.zamudio@ucla.edu). Mouse lines generated in this study are available on request without restrictions.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

The AB2.2, B9 (Ago1−/−; Ago3−/−; Ago4−/−, designated Ago2only here), E7 (Ago1−/−; Ago2−/−; Ago3−/−; Ago4−/−; hsAgo2; BsdS, designated Ago2low here) cell lines are described in Su et al. (2009) and were provided by the Wang lab. The V6.5 and DGCR8 (DGCR8−/−) knockout cells are described in Wang et al. (2007) and were provided by the Blelloch lab. The Dicer (Dicer−/−) knockout cells are described in Leung et al. (2011) and were provided by the Leung lab. TTFHAgo2 (Ago1−/−; Ago2−/−; Ago3−/−; Ago4−/−; CreERT2; pSLIK_TTFHAgo2) cells were characterized in Zamudio et al. (2014) and are available upon request from the Zamudio lab. Mouse ESCs were grown on gelatinized tissue culture plates in Dulbecco’s Modified Essential Media supplemented with HEPES (Thermo Fisher), 15% fetal bovine serum (Hyclone/Thermo Fisher), 1000U/mL leukemia inhibitory factor (Millipore), 0.1mM non-essential amino acids, 0.1mM L-glutamine, 0.1mM Pen/Strep and 0.11mM β-mercaptoethanol. All cells were cultured in 5% CO2 and ambient oxygen. For growth measurements, AgoKO clones and Ago2low cells were trypsinized at indicated times and counted in Trypan Blue stain. TTFHAgo2 cells were maintained at 0.1 μg/ml of dox.

METHOD DETAILS

Western Blotting

Mouse ESCs were lysed with RIPA buffer (1% NP-40, 0.1% SDS, 0.5% Deoxycholate in PBS pH 7.0) containing 1x complete protease inhibitors (Roche). Protein was quantified using the Pierce BCA protein assay (Thermo Scientific), resolved on Novex gradient denaturing PAGE gels (Life Technologies) and transferred in 20% methanol to PVDF membrane (Millipore). Membranes were blocked with 5% skim milk in TBS-T and incubated with primary antibodies overnight. Blots were washed in TBS-T. Secondary HRP-conjugated antibodies were purchased from GE Healthcare.

RNA Isolation and qRT-PCR

RNA was isolated using Trizol (Life Technologies) and treated with Turbo DNase (Thermo Fisher). cDNA was synthesized using random hexamers in the QuantiTect Reverse Transcription Kit (QIAGEN), SuperScript III First-Strand Synthesis System (Thermo Fisher) or M-MLV Reverse Transcriptase (Thermo Fisher). cDNA was measured using SybrGreen (Applied Biosystems) on the Applied Biosystems 7500 Real Time PCR System or LightCycler 480 (Roche).

miRNA Reporter Assays

For fluorescence reporter assays, cells were plated at 2e5 cells/well in a 12-well plate 24h prior to transfection. Cells were transfected using Lipofectamine 2000 (Life Technologies) with 500ng combined of the two-color fluorescent reporter and rtTA plasmid (2:1 mass ratio). The day following transfection cells were treated with 1 μg/mL doxycycline overnight and collected the next day for FACS quantification on a LSRII (BD Biosciences). Cells were binned for YFP expression and mCherry mean and standard deviation calculated for each bin. For luciferase reporters, cells were plated in 24-well plates 24h prior to transfection at 2e5 cells/well. Cells were transfected using Lipofectamine 2000 (Life Technologies) with pRL-CMV containing either wild-type or seed mutant 3′UTR derived from the Slc31a1 gene (Leung et al., 2011), along with PGL3 untargeted control and pWhiteScript. The miR-295 small RNA mimic was ordered from Dharmacon as an siRNA duplex containing 2-3′ nt overhangs (Table S1). At 24h after transfection, cells were lysed in 1x passive lysis buffer (Promega). Renilla and Firefly luciferase were quantified using Promega Dual Luciferase Reporter Assay System. Values represent targeted Renilla luciferase relative to untargeted Firefly luciferase as a transfection control.

Northern Blotting

20 μg of total RNA isolated with Trizol Reagent (Life Technologies) was separated on a 12% polyacrylamide/Urea/TBE gel and transferred to HyBond N+ membrane (GE Healthcare) using the TransBlot SD. Semi-dry Transfer System (Bio-Rad). DNA oligo probes were γ-32P end labeled using T4 polynucleotide kinase (New England Biolabs) and purified using Illustra G-25 MicroSpin columns (GE Healthcare). Blots were blocked in UltraHyb Oligo (Life Technologies) at 42°C for 30 min and hybridized with 5′ labeled DNA probe overnight at 42°C. Blots were washed twice with pre-warmed 2xSSC and 0.5% SDS at 42°C for 30 minutes. RNA was visualized with phosphor screens from GE-HealthCare.

Cell Cycle Profiling and Intracellular Staining

For cell cycle analysis, TTFHAgo2 cells were induced to wild-type Ago2 levels (2.5 μg/ml) for 48h and plated into media without dox in the presence of 10nM SB-431542 (Abcam) or DMSO (Stem Cell) carrier only. Cells were collected at indicated time points and fixed in 2% paraformaldehyde for 20 min. The fixed cells were washed twice in 1x HBS and stained with 5 μM DRAQ5 (Life Technologies). For intracellular staining, fixed cells were permeabilized in BD PermWash solutions (BD Biosciences) and blocked in 5% BSA/1xPBS-0.1% Triton X-100. Cells were incubated with Lefty or HA epitope antibody followed by fluorescent secondary antibodies from Life Technologies. All samples were analyzed on the LSRFortessa (BD Biosciences).

Nuclear Run On Assay

Mouse ESCs were washed 1x with cold PBS, trypsinized and collected, centrifuged, washed once with cold 1x PBS and counted. Equal numbers of cells were centrifuged and then lysed in lysis buffer (10mM Tris-Cl pH 7.5, 2mM MgCl2, 3mM CaCl2, 10% Glycerol, 0.5% NP-40) containing 1x complete protease inhibitors (Roche). The lysate was centrifuged and the nuclear pellet was retained. Nuclei were washed once with 1x Lysis Buffer and then resuspended in 100 μL freezing buffer (50mM Tris-Cl pH 8.0, 40% Glycerol, 5mM MgCl2 0.1mM EDTA) containing 1x complete protease inhibitors (Roche) per 20 million nuclei. Nuclei was mixed with equal volume reaction buffer (10mM Tris pH 8.0, 5mM MgCl2, 300mM KCl, 1mM DTT, 1.8mM ATP, 0.5mM CTP, 0.5mM GTP, 0.375mM UTP, 0.125mM Br-UTP) containing 1x Complete protease inhibitors (Roche) and 0.5U/μL SUPERase-In (Life Technologies). The transcription reaction was incubated 30 min at room temperature and halted by the addition of an equal volume acid phenol chloroform (Life Technologies). BrdU antibody (Sigma) was conjugated to protein G coated Dynabeads (Life Technologies) for 1h in 0.1M Na-Phosphate pH 8.0 buffer at 4°C with rotation. RNA was immunoprecipitated in 250 μL total volume for 1h at 4°C with rotation in RSB-100 buffer (10mM Tris-Cl pH 7.4, 100mM NaCl, 2.5mM MgCl2, 0.4% Triton X-100) with 40U SUPERase-In (Thermo Fisher) and 2.5 μg yeast tRNA. Beads were washed three times with RSB-100 buffer, RNA was extracted with acid phenol chloroform (Life Technologies) and precipitated. After resuspension, RNA was treated with Turbo DNase (Thermo Fisher) and cDNA was synthesized using QuantiTect Reverse Transcription Kit (QIAGEN). cDNA was measured using PowerSybr (Applied Biosystems) on the Applied Biosystems 7500 Real Time PCR System.

RNA-seq and small-RNA cloning

mRNA-seq libraries were prepared from a UTP-based protocol (Parkhomchuk et al., 2009) for barcoded stranded illumina library preparations. Paired-end sequencing reads were obtained on a HiSeq 2000. For small-RNA cloning, small-RNA fractions were size selected on a 12% polyacrylamide/Urea/TBE gel for RNAs from 18-75nt. After gel extraction, the isolated small-RNA was ligated with a 3′ adaptor and cloned in a 5′ phosphate independent manner using circularization of cDNA. Illumina sequencing primer binding sites and sample barcodes were incorporated during PCR amplification of circularized cDNA template. All sample libraries were prepared with two biological replicates.

ChIP protocol and library generation

ChIP was performed as described previously (Marson et al., 2008). Approximately 30 million cells were crosslinked for 15 min at room temperature at 1% final formaldehyde (Sigma) solution by adding directly to the growth media. Excess formaldehyde was quenched for 5 min with 125 mM glycine. Cells were scraped off the plates, washed twice with PBS, and then liquid removed and the cell pellet flash frozen with liquid nitrogen and stored at −80°C. For immunoprecipitations, 100 μL Protein G Dynabeads (Thermo Fisher) were blocked with 0.5% BSA (w/v) in PBS and conjugated with 40 μL of anti-H3K27me3 (Abcam, ab6002), anti-H3K9me3 (Abcam, ab8898), anti-H3K27ac (Abcam, ab4729) or anti-Pol II (Santa Cruz, sc-899) serum. To prepare extracts from cell pellets, nuclei were isolated and sonicated in lysis buffer (20mM Tris-HCl, pH 8.0, 150mM NaCl, 2mM EDTA, 0.1% SDS, 1% Triton X-100) in a Bioruptor (Diagenode) sonication system with 60 cycles at 30 s each on ice with 30 s on ice between cycles. Sonicated lysates were cleared once by centrifugation and incubated overnight at 4°C with magnetic beads bound with antibody to enrich for DNA fragments bound by the indicated factor. Beads were washed twice in Sonication buffer (20mM Tris-HCl, pH 8.0, 150mM NaCl, 2mM EDTA, 0.1% SDS, 1% Triton X-100), twice in Sonication buffer + 500mM NaCl (20mM Tris-HCl, pH 8.0, 500mM NaCl, 2mM EDTA, 0.1% SDS, 1% Triton X-100), LiCl wash (10mM Tris-HCl, pH 8.0, 250mM LiCl, 1mM EDTA, 1% NP-40), and TE + 50mM NaCl (10mM Tris-HCl, pH 8.0, 1mM EDTA, 1% NP-40). DNA was eluted, cross-links were reversed and DNA was purified with phenol chloroform extraction and ethanol precipitation. Libraries for Illumina sequencing were prepared using a protocol based on the Illumina TruSeq DNA Sample preparation kit.

Isolation of Ago Knockout Cells

Ago2low cells were grown on gelatin-coated plates at ~50-100 cells per 10cm dish for isolation of colonies derived from single cells. These cells were plated directly in media containing 1 μM 4OH-T (Sigma) for 24h or 48h, after which the media was replaced with ESC media without 4OH-T. Colonies from single cells were trypsinized in sterile cloning cylinders (Bel-Art Products), and transferred to a single well of a 12-well plate. These cells were then assayed for Ago expression by western blotting and transgene deletion by PCR from gDNA isolated with the GenElute Mammalian Genomic DNA Miniprep kit (Sigma). Cells were also scored for resistance to Blasticidin S (Invitrogen) treatment. Those clones that were annotated as true AgoKO clones were susceptible to Blasticidin treatment whereas Ago2low cells were resistant to treatment.

Tgfbr1 miRNA-binding site mutants

Tgfbr1 miRNA-binding site mutants were created in TTFHAgo2 cells using two separate CRISPR-Cas9 cassette (px330) plasmids with sgRNA sequences that targeted flanking regions surrounding the Tgfbr1 3′-UTR miR-294 site. These two plasmids and a plasmid conferring puromycin resistance were mixed in a total of 2.5 μg of DNA for transfection into cells using Lipofectamine 3000 reagents (Life Technologies). The day after transfections cells were selected with 1 μg/ml puromycin (Life Technologies). After 48h selection, puromycin-resistant cells were diluted into 15cm2 plates for growth of colonies from single cells. The colonies from single cells were manually picked and expanded in isolation. Clonal lines were lysed and gRNA isolated with QuantaBio extraction reagents for PCR screening using primers that flanked the Tgfbr1 miRNA deletion site. Deletion of homozygote mutants was validated using DNA sequencing. The oligonucleotides used for sgRNA and genotyping are described in Table S1.

Identification of active miRNA families

Active miRNA families were identified by a significantly larger fraction (p value ≤ 1e-10, Fisher’s exact test) of Ago2-iCLIP-supported 7-mer seed matches compared to other 7 nucleotide-long sites at random positions in the same 3′UTRs. Ago2-iCLIP-supported 7mers were defined by requiring at least 2 TTFHAgo2-iCLIP reads within ± 10 nt around the 7-mer site. Additionally, active miRNA families were required to be highly expressed (>= 500 small RNA-seq reads in wild-type mESCs). Ago2-iCLIP data is available at Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/; accession code GEO: GSE61347).

Identification of active miRNA targets and non-target genes

Target genes of active miRNA families were categorized as 6-mer, 7-mer or 8-mer targets if the longest Ago2-iCLIP-supported miRNA target site in the 3′UTR of a gene was a 6-mer, 7-mer or 7-mer seed match followed by an A nucleotide, respectively, according to target site definitions used by TargetScan (Lewis et al., 2005). miRNA target sites were defined as Ago2-iCLIP-supported if at least 2 Ago2-iCLIP reads were found within ± 10 nucleotides around the miRNA target site.

Non-targets were defined separately for each active miRNA family and were genes with no miRNA seed match in their 3′UTR, but with similar expression levels, 3′UTR lengths and GC contents as active miRNA family targets (in particular, within the second and third quartiles of those values for active miRNA family targets).

Quantification of miRNA levels from small RNA-seq

Small RNA-seq reads were trimmed for adaptor sequences and mapped to the mouse mm10 genome assembly with Bowtie 1 (Langmead et al., 2009) allowing up to two mismatches. The mapped reads were considered as being derived from a miRNA if the read start matched exactly the annotated start position of the mature miRNA, and ended within 2 nucleotides of the annotated mature miRNA end. Annotations of mature miRNAs were downloaded from miRbase v21 (http://www.mirbase.org) [CITATION?]. To compare miRNA levels in different Ago conditions, normalization factors were determined for each sample using a DESeq2-like normalization approach (Love et al., 2014) considering eight miRNAs that clustered by correlation with read counts of miR-2137, which appeared to be an Ago-independent miRNA. Only miRNAs with at least 10 reads in all conditions were included in the comparisons.

Exon-intron split analysis (EISA)

EISA was performed as described in Gaidatzis et al. (2015). mRNA-seq reads were mapped to the mm10 genome assembly using Tophat 2.0.14 (Trapnell et al., 2012). RefSeq exon annotations were downloaded from UCSC genome browser (https://genome.ucsc.edu). Intronic regions of genes were defined as those parts of annotated introns that did not overlap any exon annotations. Additionally, 10 nucleotides were trimmed at intron starts and ends to eliminate inaccuracies in splice site usage or annotation. mRNA-seq reads were counted that overlapped exons or intronic parts using custom scripts and bedtools multicov using the -split option. Read pairs were considered as two separate reads.

Differential gene expression analysis

Differential expression analysis was performed using DESeq (Anders and Huber, 2010). Reads mapping to exons or introns were analyzed separately.

Identification of significantly changing gene sets

Gene sets from the molecular signature database (mSigDB) were downloaded from the GSEA (Subramanian et al., 2005) webpage (http://software.broadinstitute.org/gsea). Hallmark gene sets (Liberzon et al., 2015) were supplemented with target gene sets for Smad3 and Myc transcription factors and with 8-mer, 7-mer and 6-mer target gene sets of active miRNA families.

Smad3 ChIP-seq peaks, obtained before and after stimulation with activin (Mullen et al., 2011), were downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/; accession codes GSM539541 and GSM539542). c-Myc ChIP-seq peaks (Chen et al., 2008) were downloaded from Gene Transcription Regulation Database (http://gtrd.biouml.org; experiment EXP032706). Smad3 and Myc ChIP-seq target genes were identified as RefSeq genes closest to a ChIP-seq peak, within a maximum distance of 100,000 nucleotides.

Significantly changing gene sets were identified using Wilcoxon rank sum test (p value < 0.001) to compare fold changes of genes contained in a gene set with fold changes of all other genes not contained in this set. A minimum normalized read count of 10 reads in each condition being compared was required for a gene to be included.

Identification of genes changing with partial or complete loss of Ago-miRNAs (Figure 2D)

Genes upregulated (downregulated) with partial Ago-loss were defined as those genes that increased (decreased) by more than 50% of their maximum positive (minimum negative) fold change in the Ago2low/Ago2only comparison. Genes upregulated (downregulated) with complete Ago-loss were defined as those genes that increased (decreased) by more than 50% of their maximum positive (minimum negative) fold change in the AgoKO/Ago2only comparison. The maximum positive or minimum negative fold change was determined among the Ago2low/Ago2only and AgoKO/Ago2only comparisons. Only genes with a maximum positive or minimum negative fold change < −0.05 were considered as being upregulated or downregulated, respectively. Each gene was assigned to only one category, and it was evaluated in the following order: partial loss up, partial loss down, complete loss up, complete loss down.

Analysis of histone modification and Pol II occupancy ChIP-seq data

ChIP-seq peak identification

ChIP-seq reads were mapped to the mm10 genome assembly with Bowtie 1 (Langmead et al., 2009) allowing up to two mismatches. Reads with a multiplicity > 1 were randomly assigned to a single genomic region. Peaks of histone modifications or Pol II binding were identified using MACS (Zhang et al., 2008) using ChIP-seq reads of whole cell extract (WCE) as a background sample. Reads in the larger dataset of signal or background were down-sampled to the number of reads in the smaller dataset. Peaks with a p value < 0.01 and a log2-enrichment > = 2 compared to WCE were considered for further analysis.

Normalization of ChIP-seq peak reads for each condition

Positions of peaks that were significant in any of the conditions (Ago2only, Ago2low, AgoKO_C1, AgoKO_C2) were merged. Then reads overlapping these peak positions were counted for each condition, after extending reads by 200 nucleotides.

The read count of each peak in each condition was normalized to the sum of reads overlapping those peaks that were among the top 30% of peaks in every condition ranked by the number of reads overlapping a peak.

Identification of proximal genes

To compare changes at ChIP-seq peaks with intron changes of proximal genes, ChIP-seq peaks were assigned to the closest RefSeq gene. In case several ChIP-seq peaks were assigned to the same gene, the average fold change at all ChIP-seq peaks assigned to that gene was compared with its intron fold change.

Identification of gene sets associated with strongly changing ChIP-seq peaks (Figure S2E)

Strongly changing ChIP-seq peaks between two conditions were identified by a z-score of log2 fold changes > 2 or < −2 for strongly increased or decreased peaks, respectively. Genes were assigned to strongly changing peaks (see above), and a set of reference genes was defined consisting of genes that were assigned to not strongly-changing peaks (z-scores > −2 and < 2) and that were not among the genes assigned to strongly changing peaks. To identify gene sets enriched with genes assigned to strongly changing ChIP-seq peaks, the proportion of genes contained in a gene sets among genes assigned to strongly changing peaks was compared to the proportion of genes contained in that gene sets among genes of the reference set. Fisher’s exact test was used to identify gene sets with a significant (p value < 0.05) difference in the proportion of genes.

Comparison of ChIP-seq reads changes at Smad3 and Myc ChIP-seq peaks (Figure 3D)

To compare histone modification and Pol II occupancy changes at binding sites of Smad3 and Myc transcription factors, Smad3- and Myc-bound regions were defined as a 2kb region around the center of the Smad3- or Myc ChIP-seq peak. Smad3-bound regions that overlapped Myc-bound regions were excluded. As a reference, 2kb regions were defined around the centers of histone modification or Pol II occupancy ChIP-seq peaks. Then histone modification and Pol II occupancy ChIP-seq reads overlapping these regions were counted for each condition. Read counts were normalized using the normalization factors described above. Regions were considered if in each condition of a comparison more than five ChIP-seq reads overlapped that region.

Significance of overlap between two gene sets

The significance of the overlap of genes between two gene sets was calculated using a hypergeometric test, considering all quantifiable genes with more than 10 reads in every condition.

Relative fold change of miRNA targets (Figures 5B and 6D; Figure S4D) and transcriptional networks (Figure 6D)

The relative fold change of miRNA targets in the Ago2 loss time course was calculated as the ratio of the log2 fold change between two conditions (24h/0h, 48h/24h, or 96h/48h) to the maximum positive log2 fold change occurring between 0h and any subsequent time point (24h/0h, 48h/0h, or 96h/0h). In the case of the Ago2 rescue series it was calculated as the ratio of the log2 fold change between two conditions (0.1/0 dox or 2.5/0.1 dox) to the minimum negative log2 fold change occurring between 0 dox and any of the two Ago2 induction levels (0.1/0 dox or 2.5/0 dox).

Only genes with at least 10 reads and with a maximum positive log2 fold change > 0.1 (or with a minimum negative log2 fold change < −0.1) were included in the Ago2 loss time course (or in the Ago2 rescue series).

The relative post-transcriptional decrease with 0.1 dox induction was calculated for target genes of active miRNA families with at least 10 exonic and 10 intronic reads.

In the case of relative intron level changes of genes contained in transcriptional networks, the dominant direction of change was first determined by comparing the maximum positive and absolute minimum negative log2 fold changes and selecting the larger of the two. Then the relative fold change of intron levels was calculated by dividing the log2 fold change at 0.1/0 dox by the dominant fold change, which was either the maximum or minimum of 0.1/0 dox and 2.5/0.1 dox log2 fold changes of introns.

Pathways enriched among early or late changing genes (Figure 5D)

miRNA target genes that responded early or late to Ago-miRNA loss were defined as explained in the main text based on post-transcriptional (exon-intron) fold changes, and only genes with a positive post-transcriptional fold change in the 24h/0h comparison were included. Fisher’s exact test was used to identify gene sets with a significant (p value < 0.05) difference in the proportion of genes contained in a gene set among early and late responding genes. Hallmark, KEGG and REACTOME gene sets (contained in mSigDB and downloaded from http://software.broadinstitute.org/gsea) and Smad3 target gene sets were included in the analysis.

QUANTIFICATION AND STATISTICAL ANALYSIS

For small RNA seq comparisons, the statistical significance for differences in distributions was evaluated using a one sample t test. For the comparison of cumulative fraction curves, the p values indicated are from a Wilcoxon Rank-Sum test. The significance of gene set overlap was determined by hypergeometric and Fisher’s exact tests. For qPCR data, a two-tailed t test was used and error bars indicate the standard deviation between biological replicates. A Spearman correlation was used to determine the relationship between proximal ChIP signal and RNA expression. All statistical analysis for qPCR data was done with Microsoft Excel. All other statistical analysis was done with custom python and R scripts. Details of the statistical analysis can be found in the figure legends and under Method details.

DATA AND CODE AVAILABILITY

The datasets (RNA-seq, ChIP-seq and small RNA-seq) generated during this study are available at GEO (https://www.ncbi.nlm.nih.gov/geo/) under the accession code: GEO: GSE102175. Original data for northerns, westerns and FACS in the paper is available in Mendeley Data under https://doi.org/10.17632/9r28fw3x3s.1

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
HA	Roche	11095200; RRID:AB_390918
Ago2	Cell Signaling	2897; RRID:AB_2096291
human Ago2	Wako	011-22033; RRID:AB_1106836
Ago1	Millipore	04-083; RRID:AB_11215803
Tubulin	GenScript	A01410-40; RRID:AB_1968943
Lefty	Abcam	ab204283
Smad2	Cell Signaling	5339S; RRID:AB_10626777
pSmad2	Cell Signaling	3108S; RRID:AB_490941
H3K27me3	Abcam	ab6002; RRID:AB_305237
H3K9me3	Abcam	ab8898; RRID:AB_306848
H3K27ac	Abcam	ab4729; RRID:AB_2118291
Pol II	Santa Cruz	sc-899; RRID:AB_632359
ECL Anti-Rabbit IgG, Horseradish peroxidase	GE Healthcare	AB_772206
ECL Anti-Rat IgG, Horseradish peroxidase	GE Healthcare	AB_772207
ECL Anti-Mouse IgG, Horseradish peroxidase	GE Healthcare	AB_772210
AlexaFluor 568 Goat anti-rat	Life Technologies	A11077
AlexaFluor 488 Goat anti-rabbit IgG	Life Technologies	A11008
Chemicals, Peptides, and Recombinant Proteins
Doxycycline	Sigma-Aldrich	D9891
4-Hydroxy-Tamoxifen	Sigma-Aldrich	H7904-5MG
Blasticidin S	Life Technologies	R21001
Puromycin	Life Technologies	A1113802
DRAQ5	Life Technologies	62252
SB-431542	Abcam	ab146590
Lipofectamine 2000	Life Technologies	11668-019
Lipofectamine 3000	Life Technologies	L3000008
Critical Commercial Assays
TRIzol Reagent	Fisher Scientific	15596018
DNase Turbo	Fisher Scientific	AM2238
SuperScript III First-Strand Synthesis System	Fisher Scientific	18080051
QuantiTect Reverse Transcription Kit	QIAGEN	205311
Sybr green mix	Applied Biosystems	4364344
CircLigase II	Epicenter	CL9021K
T4 RNA Ligase 2, truncated K227Q	NEB	M0351L
Deposited Data
RNA-seq, smallRNA-seq, Chip-seq	This paper	GEO: GSE102175
Raw data for northerns, westerns and FACS	This paper	Mendeley https://doi.org/10.17632/9r28fw3x3s.1
Exonic read counts of naive mESCs and primed epiblast cells	(Chen et al., 2018)	GEO: GSE112749
Ago2-iCLIP-Seq	(Bosson et al., 2014)	GEO: GSE61347
c-Myc ChIP-Seq peaks	(Chen et al., 2008)	EXP032706
Smad3 and Smad3_activin ChIP-Seq peaks	(Mullen et al., 2011)	GEO: GSM539541 and GEO: GSM539542
Experimental Models: Cell Lines
Mouse: AgoWT (AB2.2): WT	Su et al., 2009	NA
Mouse: Ago2only (B9): Ago1−/−; Ago3−/−; Ago4−/−	Su et al., 2009	NA
Mouse: Ago2low (B9): Ago1−/−; Ago2−/−; Ago3−/−; Ago4−/−; hsAgo2floxed; CreERT2; BsdS	Su et al., 2009	NA
Mouse: TTFHAgo2: Ago1−/−; Ago2−/−; Ago3−/−; Ago4−/−; CreERT2; pSLIK hsAgo2 Hygro	Zamudio et al., 2014	NA
Mouse: AgoKO (C1,C2): Ago1−/−; Ago2−/−; Ago3−/−; Ago4−/−; CreERT2;	This paper	NA
Mouse: TTFHAgo2_Tgfbr1miRNAKO: Ago1−/−; Ago2−/−; Ago3−/−; Ago4−/−; CreERT2; pSLIK hsAgo2 Hygro; Tgfbr1_miR295site−/−	This paper	NA
Mouse: V6.5 ESC WT	Wang et al., 2007	NA
Mouse: V6.5 Dgcr8−/−	Wang et al., 2007	NA
Mouse: 2fc ESC WT	Leung et al., 2011	NA
Mouse: 2fc1 ESC Dicer1−/−	Leung et al., 2011	NA
Oligonucleotides
miR-295 small RNA mimic	This paper	Table S1
sgRNA primers	This paper	Table S1
Northern and RT PCR primers	This paper	Table S1
Recombinant DNA
pX330	Addgene	#42230
pTRETightBI-RY-1pf	Addgene	#31467
pRL-CMV-Slc31a1	Leung et al., 2011	NA
Software and Algorithms
Bowtie1	Langmead et al., 2009	http://bowtie-bio.sourceforge.net/index.shtml
Tophat 2.0.14	Trapnell et al., 2012	http://ccb.jhu.edu/software/tophat/downloads/
MACS	Zhang et al., 2008	http://liulab.dfci.harvard.edu/MACS/Download.html
DESeq	Anders and Huber, 2010	https://bioconductor.org/packages/release/bioc/html/DESeq.html
Bedtools v2.26.0	Quinlan and Hall, 2010	https://bedtools.readthedocs.io/en/stable/index.html
Other
HALLMARK, REACTOME and KEGG gene sets	mSigDB	http://software.broadinstitute.org/gsea/msigdb/collections.jsp