MicroRNA dynamics at the onset of primordial germ and somatic cell sex differentiation during mouse embryonic gonad development.

In mammals, commitment and specification of germ cell lines involves complex programs that include sex differentiation, control of proliferation, and meiotic initiation. Regulation of these processes is genetically controlled by fine-tuned mechanisms of gene regulation in which microRNAs (miRNAs) are involved. We have characterized, by small-RNA-seq and bioinformatics analyses, the miRNA expression patterns of male and female mouse primordial germ cells (PGCs) and gonadal somatic cells at embryonic stages E11.5, E12.5, and E13.5. Differential expression analyses revealed differences in the regulation of key miRNA clusters such as miR-199-214, miR-182-183-96, and miR-34c-5p, whose targets have defined roles during gonadal sexual determination in both germ and somatic cells. Extensive analyses of miRNA sequences revealed an increase in noncanonical isoforms on PGCs at E12.5 and dramatic changes of 3' isomiR expression and 3' nontemplate nucleotide additions in female PGCs at E13.5. Additionally, RT-qPCR analyses of genes encoding proteins involved in miRNA biogenesis and 3' nucleotide addition uncovered sexually and developmentally specific expression, characterized by the decay of Drosha, Dgcr8, and Xpo5 expression along gonadal development. These results demonstrate that miRNAs, their isomiRs, and miRNA machinery are differentially regulated and participate actively in gonadal sexual differentiation in both PGCs and gonadal somatic cells.

INTRODUCTION

Primordial germ cells (PGCs) are the precursors of embryonic germ cells and therefore the precursors of gametes (Ohinata et al. 2009). Their specification is mediated by inductive signals such as BMP4, BLIMP1, PRDM14, and TCFAP2C (Lawson et al. 1999; Schemmer et al. 2013; Günesdogan and Surani 2016), and they are detectable in the mouse embryo as a reduced number of cells within the epiblast at embryonic day 6.5 (E6.5). PGCs migrate to the developing hindgut endoderm at E7.75, into the mesentery at E9.5, and colonize the genital ridges at E10.5 (Saitou and Yamaji 2012). From E11.5, during gonad colonization, PGCs start a sex-dependent differentiation. At E13.5, female PGCs enter meiosis (Chuma and Nakatsuji 2001; Nakatsuji and Chuma 2001), while in males PGCs enter mitotic arrest to become prospermatogonia until puberty (Western et al. 2008; Griswold 2016).

MiRNAs are small noncoding RNAs (sncRNAs) ∼22 nucleotides (nt) in length that act as post-transcriptional regulators of gene expression through binding to mRNAs by sequence complementarity (Jonas and Izaurralde 2015). This interaction is determined by a region (called the “seed region”) of 6–8 nt in the 5′ region of the miRNA sequence (usually from the 2nd to the 8–10th nucleotide) (Bartel 2009). Primary miRNA transcripts are processed to pre-miRNAs in the nucleus by a dsRNA nuclear type III endoribonuclease (DROSHA) and DGCR8 and transported to the cytoplasm by XPO5. Then, DICER cleaves the pre-miRNA followed by duplex separation that generates functionally mature miRNA molecules. It is known that DICER and/or DROSHA process miRNA precursors imprecisely, generating miRNA variants with several plus/minus nucleotides at the 5′ and/or 3′ end called isomiRs (Morin et al. 2008; Neilsen et al. 2012). IsomiRs with different 5′ ends (i.e., different seed sequence) have different target repertories than their corresponding canonical miRNAs (Cloonan et al. 2011). This increases the possibility of interaction with different potential targets from the same miRNA precursor and thus extends their functionality. In addition, miRNAs can also be post-transcriptionally adenylated or uridylated, which could alter miRNA targeting properties and their stability (Burroughs et al. 2010; Marzi et al. 2016).

MiRNAs have been described as regulators in most developmental processes including embryogenesis and pluripotency (Pauli et al. 2011; Wright and Ciosk 2013). For instance, miRNA clusters such as miR-290 and miR-17-92 are known to be essential for germ cell development (Hayashi et al. 2008; Medeiros et al. 2011). However, little is known about the role of miRNAs and their isomiRs in mouse gonadal sex determination (E11.5–E13.5) in both PGCs and supporting somatic cells. Some previous studies did not differentiate between PGCs and gonadal somatic cells (Rakoczy et al. 2013; Bhin et al. 2015) or between males and females at E11.5 (Hayashi et al. 2008), and neither characterized the isomiR population and the regulation of genes involved in miRNA biogenesis. Consequently, it is crucial to elucidate the potential participation of specific miRNAs and their isomiRs in both PGCs and gonadal somatic cells during this key developmental window. To achieve this, we isolated PGCs and somatic cells from male and female embryos at E11.5, E12.5, and E13.5 to perform NGS of the sncRNA population. Using molecular and bioinformatics approaches, we have identified and characterized specific sexual and developmental expression patterns of miRNAs/isomiRs and genes involved in miRNA biogenesis. Differential expression analyses identified several miRNAs with targets that have critical roles in gonadal sexual fate and development. Analyses of isomiR sequences and 3′ nontemplate nucleotide additions (3′ NTA) revealed dramatic differences in E13.5 female PGCs, which could be potentially associated with their meiotic entry. Finally, the analyses performed by RT-qPCR of miRNA biogenesis machinery and 3′ terminal uridylyl transferases (Tut4 and Tut7) showed a decrease in the expression of key modulators of pri-miRNA processing and pre-miRNA transport and up-regulation of Tut4 during PGC development.

RESULTS

MiRNAs from PGCs vs. somatic gonadal cells, sex, and development show differential expression

Using our bioinformatic pipeline (Supplemental Fig. S1), we identified between 916 and 721 different miRNAs, which corresponded to a total of between 17,386 and 4,530 miRNA sequences, considering all diverse isomiRs, in the different samples analyzed (Table 1). Previous studies on complete gonads, but using older versions of miRBase, were able to detect only 331 different miRNAs (Rakoczy et al. 2013).

TABLE 1.

Summary of small RNA-seq

Despite the attributed critical role of miRNAs in developing PGCs between E11.5 and E13.5 (Hayashi et al. 2008), significantly higher populations of miRNAs were detected in somatic cells in both sexes at the different stages of development when compared to PGCs (Fig. 1A,B). Interestingly, in both cell types, PGCs and somatic cells, the highest percentage of reads associated to miRNAs was detected in E11.5 female gonads (Table 1). Another surprising finding was the significant increase of abnormally short miRNA reads (16 nt in length) in E12.5 male and female PGCs (Fig. 1D). Interestingly, these samples also showed the lowest percentage of reads associated to miRNAs and detected miRNA sequences (Table 1). The potential roles of these specific variations are yet unknown.

FIGURE 1.

Characterization of miRNA expression in male and female PGCs and gonadal somatic cells. (A,B) Read length distribution of the aligned reads to the pre-miRNA database. Reads within 16 and 30 nt are represented. (A) Somatic cells; (B) primordial germ cells. (C,D) Percentage of trimmed reads aligned to the pre-miRNA database. (C) Somatic cells; (D) primordial germ cells. (E) Heatmap representing log2 normalized reads of miRNAs corresponding to sequences classified as “No Change” with unsupervised hierarchical clustering. The range of colors goes from blue (minimum expression) to red (maximum expression). These reads correspond to canonical miRNA sequences and sequences with the same seed (3′ isomiRs) without allowing mismatches.

An unsupervised hierarchical clustering (Fig. 1E) was performed to analyze the effects of sex, developmental stage, and cell type on global miRNA expression. This revealed that cell type was the main factor that determined miRNA expression, as the dendogram grouped the samples into two clearly differentiated groups: one integrating all PGCs, except females at E11.5, and another grouping all the somatic cells. Interestingly, developmental stage was more determinant of global miRNA expression than sex. Previous studies also showed that cell type and developmental stage were the main factors to determine the transcriptome expression of coding genes in these same samples (Jameson et al. 2012). Another interesting fact was that, in contrast to E12.5 and E13.5, at E11.5 males and females were not in the same clade, showing a completely different miRNA expression pattern. E11.5 is key in the genetic regulation of sexual differentiation of embryonic germ cells (Feng et al. 2014). Our expression analyses strongly suggested that differentiation at the level of gene expression to enter either meiotic process in females or arrest of division in males could be mediated by differentially expressed miRNAs.

Differential production of isomiRs during gonadal development

IsomiRs, produced by alternative cleavage within the pre-miRNA, enrich miRNA regulatory networks (Cloonan et al. 2011). We have classified the different miRNA variants in PGCs and gonadal somatic cells to characterize the potential role of isomiRs in the developing gonad of both sexes between E11.5 and E13.5 (Fig. 2).

FIGURE 2.

Characterization of miRNA isoforms during gonad development. “Noncanonical processing” refers to 5′ isomiRs (trimming and additions), which have a different seed sequence than the canonical miRNA; “No Change” refers to sequences with the same seed sequence than the canonical miRNA (3′ isomiRs and canonical miRNAs) without sequence variations. Canonical miRNAs and 3′ isomiRs with sequence variations are classified as “first nucleotide,” “outseed,” or “inseed” depending on where those variations occur. In cases where different variations occurred, the priority was first nucleotide > inseed > outseed. (A,B) Classification of miRNA isoforms as percentage of total reads. Somatic cells and PGCs, respectively. (C,D) Classification of miRNA isoforms as percentage of the total sequences. Somatic cells and PGCs, respectively. (E) Heatmap representing for each miRNA the percentage of reads corresponding to isomiRs (5′ and 3′) with respect to total reads of that miRNA. MiRNAs with low expression (canonical sequence or sum of isomiRs sequences being less than 50 counts) are not represented. Green (0%–49% of reads corresponding to isomiRs), black (50% of reads corresponding to isomiRs), red (51%–100% of reads corresponding to isomiRs). (F) Read coverage over mmu-miR-20a-5p and mmu-miR-106a-5p in E13.5 PGCs.

First, we classified miRNA sequences based on their seed region, since it mainly determines their targeting capabilities (Lewis et al. 2005; Agarwal et al. 2015), and represented them in relation to the total number of reads (Fig. 2A,B) and total number of different sequences (Fig. 2C,D). In all samples, sequences with the same seed region as canonical miRNAs and without mismatches (classified as “No Change”) represented a small fraction of the total sequences (Fig. 2C,D) but accumulated most of the total reads (Fig. 2A,B). That is, sequences with the same expected targets as canonical sequences seemed to be positively selected over sequences with different seed regions (Fig. 2C,D). Additionally, variations outside the seed region (“outseed”) with respect to variations inside (“Inseed”) were also positively selected (Fig. 2A–D). These results suggested the existence of a putative selection mechanism biasing sequences with the same seed region as the canonical miRNA and consequently the same expected targets. Another interesting finding was that at E12.5 there was an increased accumulation of 5′ isomiR variants (classified as “noncanonical processing”), which can target different mRNAs. These changes in the regulation of miRNA processing in PGCs at E12.5 were associated with a lower number of detected miRNAs and total miRNA reads, as shown in Table 1, just after the initiation of morphological sex determination in germ cells (Bowles et al. 2006).

Then, we evaluated the proportion of isomiRs corresponding to each miRNA (Fig. 2E). Interestingly, each miRNA showed differential isomiR accumulation. The number of miRNAs with most of their reads corresponding to canonical sequences (Fig. 2E, red dendrogram) was lower than those miRNAs with isomiRs detected at low or high abundance (Fig. 2E, green and blue dendrogram, respectively). These patterns were very similar among all samples except for E13.5 female PGCs. This suggested that the generation of isomiRs is miRNA-dependent rather than cell type-dependent except in critical differentiating circumstances such as the cytologically detectable entry into meiosis in female PGCs at E13.5. The differences in the accumulation of isomiRs for some miRNAs in E13.5 female PGCs compared to the other samples were due to a drastic increase or decrease in the number of 3′ isomiRs with respect to the canonical miRNA (data not shown), as illustrated in Figure 2F by clear examples of two relevant miRNAs. Variations in the 3′ end of miRNAs have been related with differences in miRNA stability and transport (Hwang et al. 2007; Bail et al. 2010).

Sexual dimorphism in E13.5 PGCs of nontemplate additions

Recent studies have shown that nontemplate nucleotide additions (NTA) to the 3′ end of miRNAs, which occur on a genome-wide scale (Burroughs et al. 2010), are demonstrated as relevant modifications that contribute to the complexity of miRNAs (Wyman et al. 2011). Most of them are uridylations and adenylations (Landgraf et al. 2007). On miRNAs, uridine additions are mediated by 3′ terminal uridylyl transferases (TUT4 and TUT7, also known as ZCCHC11 and ZCCHC6) and adenine additions by PAPD4 (also known as GLD-2) (Katoh et al. 2009; Thornton et al. 2014). We have looked for differences in 3′ NTA during gonad development in both PGC and somatic cells by classifying all NTA by the number of poly-additions (not allowing combinations of different nucleotides) and monoadditions. Since the number of polynucleotide additions was very low (data not shown), just mononucleotide additions were represented in Figure 3. When all miRNAs were analyzed together, monouridylation was the most abundant 3′ NTA followed by cytosine addition (Fig. 3A). However, E13.5 female PGCs showed an impressive amount of adenine additions and a dramatic decrease in cytosines. These data were then confirmed when NTAs at individual miRNA levels (Fig. 3B,C; Supplemental Table S1) were analyzed. Interestingly, among the top miRNAs with adenine additions, we found many members from miRNA clusters 290–295 and 17–92, already known for their importance in PGC development and spermatogenesis (Hayashi et al. 2008; Medeiros et al. 2011). It is remarkable that this dramatic increase in 3′ adenylation of miRNAs in E13.5 female PGCs overlaps with variations in 3′ isomiR production in these cells (Fig. 2E). In fact, some miRNAs that presented a change in the predominant miRNA sequence, for example mmu-miR-20a-5p and mmu-miR-106a-5p (Fig. 2F), also suffered an increase in monoadenylation (Fig. 3C).

FIGURE 3.

Analyses of 3′ nontemplate additions during gonad development. (A) Percentage of reads, calculated with all miRNAs with mono 3′ NTAs (A in green, U in red, C in blue, G in yellow). (B) Abundance of 3′ addition of monoadenine and monocytosine in E13.5 PGCs for individual miRNAs. (C) Abundance of mono 3′ NTA with respect to the total number of reads in E13.5 PGCs for individual miRNAs from miR-290–295 and 17–92 clusters. (D) Relative mRNA expression of Tut4 and Tut7 in somatic cells and PGCs during gonad development. Data were obtained by RT-qPCR using Ppia and U6 as housekeeping genes to normalize the data, which is represented as the mean of three technical replicates ±SD.

To find a possible explanation for this interesting phenomenon, we assessed the expression of genes involved in miRNA 3′ NTA, Tut4, Tut7, and Papd4 (Fig. 3D). Surprisingly, we were not able to detect Papd4 expression in any sample (data not shown). This indicated that the increase in 3′ miRNA monoadenylation observed in E13.5 female PGCs is independent of Papd4 expression. Low levels of Tut7 were observed in all samples, except by E11.5 female somatic cells. Tut4 was the most expressed gene across all samples and presented a similar expression across the somatic cell population. On the other hand, in E11.5 PGCs, its expression was barely detected and increased steadily toward E13.5. At that point (E13.5), female PGCs presented a higher expression of Tut4 with respect to males (1.4 log2 fold change, P-value <0.05). Those differences in Tut4 expression were not reflected in 3′ uridylation of miRNAs on a global scale (Fig. 3A) or miRNA level (Supplemental Fig. S2). Tut4 and Tut7 also participate in the uridylation of mRNAs (Lim et al. 2014), which could explain why a correlation between their expression and miRNA uridylation is not observed in our data.

Regulation of miRNA biogenesis during mouse gonadal development

Due to the conspicuous differences in miRNA populations among the cellular samples (Fig. 1B), we analyzed by RT-qPCR the regulation of expression of genes encoding key elements of the miRNA biogenesis during gonadal sex determination.

The microprocessor complex Drosha/Dgcr8 together with Ago2 were the most expressed genes across all samples (Fig. 4). On the other hand, Ago3 and Ago4 showed lower accumulation than the other Ago transcripts during all developmental stages, as has also been detected in germ and somatic testicular cells (González-González et al. 2008) and early embryos (García-López and del Mazo 2012). With PGC development, the expression of Drosha, Dgcr8, Xpo5, and Ago2 decayed (levels at E11.5 were higher than at E13.5, except for Ago2 in E13.5 male PGCs) (Fig. 4B), indicating that PGCs differentiation was concomitant with a decline in the expression of miRNA biogenesis machinery. On top of that, we found significant differences in the expression of Ago2 and Dicer1 transcripts between male and female E13.5 PGCs, which was much higher in males (Fig. 4C). We also found seven miRNAs up-regulated in E13.5 female PGCs with respect to males with Ago2 and/or Dicer1 as validated targets (Fig. 4D), which may help to explain the differences in the expression of these genes.

FIGURE 4.

Relative expression of genes encoding proteins involved in miRNA biogenesis. (A) Expression of genes encoding proteins involved in miRNA biogenesis in somatic cells and (B) in PGCs. Data were obtained by RT-qPCR using Ppia and U6 as housekeeping genes to normalize the data, which is represented as the mean of three technical replicates ±SD. (C) Expression of Dicer1 and Ago2 in E13.5 PGCs. (*) P < 0.05, one-way ANOVA followed by Tukey's HSD post hoc tests. (D) Network of miRNAs up-regulated in E13.5 female PGC with respect to male that target Ago2 and/or Dicer1.

Differential miRNA expression in gonads considering: sex, PGCs vs. somatic cells and development

After data normalization, differential expression analyses (DE) were performed using the R/Bioconductor package DESeq (Anders and Huber 2010). To avoid false positives, we only considered as units for the analyses those sequences without mismatches and the same seed region as the canonical miRNA (sequences classified as “no change”). Also, we validated NGS data by performing RT-qPCR (Pearson correlation; R = 0.74; P < 0.001; Supplemental Fig. S4) of some representative miRNAs (see Materials and Methods). The DE analyses were designed as follows: Data from female cells were compared with male cells of the same age (i.e., PGC11F versus PGC11M); PGCs were compared with somatic cells of the same age (i.e., PGC11F versus SC11F), and finally each cell type was analyzed following their temporal development from E11.5 to E13.5. The thresholds to consider a miRNA differentially expressed were minimum fold change of two and 100 counts. Venn diagrams were created using the R/Bioconductor package VennDiagram (Chen and Boutros 2011).

From the large number of validated mRNA targets of the miRNAs that we found to be differentially expressed, we selected some miRNAs that could have particular relevance in the process of differentiation of gonad cells in this period, which are compiled in Table 2 and mentioned in this section. Our DE analyses demonstrated that the dynamics of miRNAs are tightly regulated during gonadal development. We detected remarkable differences in relation to sex, PGCs versus somatic cells, as well as between the different developmental stages (Fig. 5A–C).

TABLE 2.

Summary table of the most representative DE miRNAs and their targets

FIGURE 5.

Differential miRNA expression during gonad development. (A) Results of females versus males differential expression analyses. In the first column are represented the analyses of PGCs, and in the second, of somatic cells. The first row of pie charts corresponds to females versus males at E11.5, the second at E12.5, and the third at E13.5. (B) Results of PGCs versus somatic cells differential expression analyses. The first column corresponds to female PGCs versus female somatic cells and the second to male analyses. The rows correspond to the same developmental stages as in A. (C) Results of differential expression analyses based on gonad development. The first row corresponds to E11.5 versus E12.5 analyses and the second to E12.5 versus E13.5. The first two columns correspond to PGCs (females and males, respectively) and the last two to somatic cells. Normalized expression of miRNA reads corresponding to sequences classified as “no change.” (Yellow) Female PGCs; (orange) male PGCs; (blue) female somatic cells; (aquamarine) male somatic cells. (D) Up-regulated miRNAs in E13.5 male somatic cells with respect to females. (E–I) MiRNA families that are differentially regulated across the different samples analyzed.

Note that 5′ isomiRs (“noncanonical processing”) were considered separately. DE analyses of 5′ isomiRs did not reveal relevant differences between the different samples, even in E12.5 PGCs, where we found an increase of 5′ isomiRs. Since these sequences correspond to <20% of the total reads and most of them have 100 or less reads, they have little relevance compared to 3′ isomiR and mature miRNA sequences.

Males vs. females

A total of 198 different miRNAs were up-regulated in female PGCs compared to males, while only 47 miRNAs were down-regulated (Supplemental Table S2A). Interestingly, somatic cells showed a similar proportion of differentially expressed miRNAs, 257 up-regulated in female somatic cells, and 65 down-regulated (Supplemental Table S2A). Particularly, at E11.5 we found that female gonads showed an impressive number of up-regulated miRNAs compared to males (192 in PGCs, 248 in somatic cells, Fig. 5A), just when sexual determination occurs (Bowles et al. 2006). We also found a switch in the expression of some particular miRNAs between female and male PGCs at E12.5 (five miRNAs up-regulated in females, 38 down-regulated) and E13.5 (37 miRNAs up-regulated in females, two down-regulated) (Fig. 5A). This suggests that the meiotic process, initiated in females E13.5, is potentially regulated by an increased presence of miRNAs, which would eliminate mRNAs interfering with the onset of meiosis.

Some of the DE miRNAs have validated targets with key roles in gonad sexual determination. miR-103-3p, let-7g-5p, miR-107-3p, and miR-26a-5p, whose targets Cyp26b1 and Fgf9 are key in male germ cell differentiation (Bowles et al. 2006, 2010), were down-regulated in E11.5 male somatic cells with respect to females. We also found 18 miRNAs (mmu-miR-19a-3, mmu-miR-22-3p, mmu-miR-30a-5p, mmu-miR-30d-5p, mmu-miR-30e-5p, mmu-miR-125a-3p, mmu-miR-139-5p, mmu-miR-140-5p, mmu-miR-149-5p, mmu-miR-185-5p, mmu-miR-202-3p, mmu-miR-204-5p, mmu-miR-214-3p, mmu-miR-500-3p, mmu-miR-532-3p, mmu-miR-532-5p, mmu-miR-665-3p, mmu-miR-667-3p) that were up-regulated in E13.5 male somatic cells with respect to females, which target different elements of the activin pathway, such as Smad3/4 and Acvr1b (Pauklin and Vallier 2015). Since Cyp26b1 expression is inhibited by activin, a member of the Tgf-β superfamily (Kipp et al. 2011), inhibition of this pathway would allow its expression and thus contribute to the degradation of RA in male germ cells. In addition, at E13.5 we found down-regulated miRNAs in female somatic cells with respect to males (Fig. 5D,E), which have been validated to target key regulators of female and male gonad differentiation such as Ctnbb1 (Liu et al. 2009) (miR-34c-5p and miR-455-3p), Six1 and Six4 (Fujimoto et al. 2013). Also, miR-30 family and miR-140-3p, up-regulated in E13.5 male somatic cells, have been reported to be key in the adhesion between Sertoli cells and spermatids (miR-30) (Nicholls et al. 2011), as well as the modulation of Leydig cell numbers in testis development (miR-140-3p) (Rakoczy et al. 2013).

With respect to PGCs, we also found two miRNA clusters differentially expressed between males and females. The cluster miR-199-214, up-regulated in female PGCs, has been described to be up-regulated during in vitro differentiation induced by RA (Le et al. 2009) to participate in RA-induced differentiation (Juan et al. 2009; Laursen et al. 2013) and in the regulation of cell proliferation (Wang et al. 2012). On the other hand, the clusters miR-182-183-96 and miR-17-92 were down-regulated in female PGCs at E11.5 and E12.5. Interestingly, one 5′ isomiR of miR-183-5p (first 5′ nucleotide trimmed, miR-183-5p_T1) was also down-regulated in females but at E11.5 its mature form or 3′ isomiRs was not. miR-183 targets the 5mC-specific dioxygenase Tet1. This gene plays an important role in mouse oocyte meiosis activation (Yamaguchi et al. 2012) and would explain the relative decrease of this cluster in female PGCs. Then, miR-183-3p_T1 was predicted to target Rxra, a receptor of the meiotic inductor RA. miR-17-92 is known to be the key regulator in spermatogenesis (Tong et al. 2012; Xie et al. 2016), being down-regulated in the presence of RA (Beveridge et al. 2009). In addition, it participates in the regulation of pluripotency (for review, see Wang et al. 2013). These data indicate that these miRNAs are key for the development and sexual differentiation of PGCs.

PGCs vs. gonadal somatic cells

As expected by the differences in miRNA populations between PGCs and somatic cells (Fig. 1A), the number of miRNAs down-regulated in PGCs with respect to somatic cells was very high (332 different miRNAs across all comparisons of PGC versus somatic cells) (Supplemental Table S2B).

We identified two DE clusters of miRNAs between PGC and somatic cells. One up-regulated in somatic cells against PGCs, let-7 family, and one down-regulated, cluster miR-290-295. These clusters participate in pluripotency (miR-290-295) and differentiation (miR-let-7) (Wright and Ciosk 2013). Interestingly, the expression of miR-290-295 increased from E11.5 to E13.5 in PGCs and let-7 family in somatic cells (Fig. 5F,G). Surprisingly, the levels of the cluster miR-290-295 were very low in E11.5 female PGCs with respect to males, which could be interpreted as a consequence of differentiation in females and the loss of its relative pluripotency (Jouneau et al. 2012).

miR-182-183-96 cluster was also down-regulated in somatic cells. miR-182-5p and miR-183-5p_T1 target Foxo3 and Foxo1, respectively, which are regulators of the cell cycle (Schmidt et al. 2002). This cluster has already been reported to be highly expressed in E13.5 male PGCs (García-López et al. 2015) and may be necessary for the proliferation of PGCs, which takes place once they colonize the genital ridges until their cell cycle is arrested.

Comparing developmental stages

Differential expression analyses by developmental stage showed a negative regulation of miRNAs from E11.5 to E12.5 in all samples, except for male somatic cells that showed the opposite pattern (Fig. 5C), followed by an up-regulation and partial recovery (Supplemental Fig. S3) of miRNAs from E12.5 to E13.5. This may suggest that a strong regulation of miRNAs between E11.5 and E13.5 is critical for the development and sexual determination of PGCs.

In these differential expression analyses, we identified a very interesting cluster of miRNAs: miR-880-881-741-470-465. These miRNAs were almost absent at E11.5 in male and in female PGCs while at E13.5 they had 1500–3000 reads (Fig. 5H). In somatic cells its levels also increased, but the expression was lower than in PGCs. Some of the validated targets of this cluster are Oct4, Nanog, and Lin28b, which participate in the maintenance of PGC pluripotency in early stages (Irie et al. 2014) and begin to be down-regulated after E13.5 (Yoshimizu et al. 1999; Yamaguchi et al. 2005). In the same context, miR-200a-3p was up-regulated in E13.5 PGCs with respect to E11.5 and E12.5. Other members of this miRNA family also suffered an increase in expression, but less than twofold. Interestingly, while the levels of miR-200a-3p were increased in PGCs, levels in somatic cells were down-regulated (Fig. 5I). miR-200a-3p miRNA is up-regulated in differentiating neural stem cells, where an inverse correlation between miR-200a-3p levels and Oct4 and Nanog (among others) was observed (Pandey et al. 2015).

Gene Ontology (GO) enrichment analyses

Analysis of GO revealed interestingly enriched terms from targets of miRNAs up-regulated in female gonads (PGCs and somatic cells) at E13.5 with respect to males that were related to “male sex differentiation” (Supplemental Table S3). This indicated that the expression of male differentiating factors could be being inhibited in female gonads of E13.5 by relative overexpression of the miRNAs targeting “male expressed” mRNAs. In contrast, GO enrichments of targets from DE miRNAs comparing male and female PGCs at E11.5 and E12.5 revealed that in males there were enriched terms related to “chromatin modification” and “reproductive system development” (Supplemental Tables S4, S5), while in females, the most representative terms were related to “ubiquitin-dependent protein catabolic process,” and “steroid metabolic process” (Supplemental Tables S6, S7).

Finally, GO enrichment analyses of targets of miRNAs down-regulated from E11.5 to E13.5 in PGCs revealed a selective regulation of processes involved in morphogenesis and differentiation, such as “regulation of cell morphogenesis involved in differentiation” and “positive regulation of neuron differentiation” (Supplemental Table S8). Interestingly, the most enriched processes regulated by miRNAs that were up-regulated at E13.5 with respect to E11.5 were “chromatin modification,” “response to growth factor,” and “cell-type specific apoptotic process” (Supplemental Table S9).

Altogether, these data are indicative of the key and solid participation of specific miRNAs in both PGCs and in accompanying somatic cells in the crucial period between E11.5 to E13.5 of male and female mammalian gonad differentiation.

DISCUSSION

In mice, the developmental window corresponding to E11.5–E13.5 is crucial in the early stages of differentiation and fate of the germline. Until now, most of the studies carried out during this developmental window have not distinguished between the somatic cell population and PGCs (Rakoczy et al. 2013; Bhin et al. 2015). Most likely, the main reasons for such a relative lack of studies in this field may be the difficulties in obtaining high enough amounts of purified samples to perform molecular and gene expression analysis by NGS.

The miRNA accumulation in types and number of molecules was significantly higher in somatic cells than in PGCs (Fig. 1A). In addition to other roles, miRNAs in developing PGCs are known for their participation in the maintenance of naïve pluripotency and germ cell lineage (Hayashi et al. 2008; Medeiros et al. 2011). However, their roles in developing gonadal somatic cells at this time (E11.5, E12.5, and E13.5) have not been assessed. MiRNAs are known to be necessary for the development of Sertoli and granulosa cells (Nagaraja et al. 2008; Kim et al. 2010), in which male and female somatic stromal cells of the gonads are, respectively, differentiated during development.

The results of the hierarchical clustering showed that E11.5 male and female miRNAome (in PGCs and somatic cells) were dramatically different compared to E12.5 and E13.5 (Fig. 1E). At this stage (E11.5), female PGCs initiate the meiotic process to oocytes, which is regulated by their surrounding somatic cells (Bowles and Koopman 2007). This fate decision implies a regulation of a variety of signaling pathways, where miRNAs could play important fine-tuned regulatory roles, explaining the differences detected in this study between male and female miRNAome.

Another interesting fact was that in PGCs at E12.5, the relative proportion of miRNA levels dropped in both males and females. This fast-flux of miRNAs in E12.5 PGCs could possibly be explained by rapid processes of miRNA degradation or turnover (Guo et al. 2015; Sanei and Chen 2015). However, some members of the cluster miR-290-295, known for their role in PGC development (Medeiros et al. 2011), maintain their expression from E11.5 to E12.5 (Fig. 5F), suggesting that changes in miRNA dynamics seem to be selective for a group of miRNAs rather than random. This down-regulation of miRNA expression at E12.5 in PGCs coincides with a relative increase in shorter sequences (Fig. 1D) and “noncanonical” miRNA isoforms (Fig. 2C,D). Thus, the patterns detected by comparatively analyzing the expression in PGCs at E11.5 versus E12.5 strongly suggest a reprogramming of PGCs concerning the sex fate determination and the entry of PGCs into meiosis in females. It is interesting to note that the levels of mRNAs encoding key elements of biogenesis and function of miRNAs were also decreased in E12.5 with respect to E11.5.

Differentially expressed miRNAs with validated targets involved in gonad cell differentiation

MiRNAs participating in meiotic entry regulation

In comparing female somatic cells with males at E11.5, we found that three miRNAs (miR-103-3p, let-7g-5p, and miR-107-3p) targeting Cyp26b1 were up-regulated in females. Blocking the expression of Cyp26b1 (which degrades RA, the inductor of meiotic entry) in pre-granulosa cells (female somatic cells of the gonad) allows female PGCs to enter into meiosis (Bowles et al. 2006). Thus, these miRNAs could be a cornerstone in the down-regulation of Cyp26b1 in gonadal somatic cells to allow the entry into meiosis of female PGCs. On the other hand, 18 miRNAs targeting the activin signaling pathway (such as Smad3/4 and Acvr1b) were up-regulated in E13.5 male somatic cells in contrast with females. Activin has been reported to be a strong inhibitor of Cyp26b1 (Kipp et al. 2011). Interestingly, fetal gonads at E12.5 start to express activin and activin receptors in a sex-dependent manner (Feijen et al. 1994) and fetal ovaries at the same developmental stage strongly down-regulate Cyp26b1 expression (Bowles et al. 2006). Moreover, the addition of activin to cultured E12.5 ovaries enhanced PGC progression throughout meiotic prophase I (Liang et al. 2015). These data suggest that inhibition of some members of the activin pathway by miRNAs in male somatic cells may potentiate the expression of Cyp26b1 and promote male fate.

Another miRNA, miR-26a-5p, which is down-regulated in E11.5 male somatic cells in contrast with females, targets Fgf9. This gene acts as a meiosis inhibitor factor in male gonads by preventing the expression of Stra8, which is male specific at E11.5 in gonadal somatic cells, and is slowly down-regulated toward E13.5 (Bowles et al. 2010). Interestingly, miR-26a-5p is differentially expressed at low levels in E11.5, but then increases dramatically in male somatic cells (threefold from E11.5 to E13.5). These data indicate that miR-26a-5p may participate in the regulation of Fgf9 expression in male somatic cells of the gonad during sexual determination.

We have also found miRNAs that could potentially regulate the entry into meiosis in female PGCs, such as miR-214-3p and miR-183-5p. miR-214 (from the cluster miR-199-214) is up-regulated during in vitro differentiation of SH-SY5Y cells induced by RA (Le et al. 2009). It is tempting to hypothesize that the elements of this miRNA cluster may be involved in differentiation processes mediated by RA such as sexual differentiation. Also, TWIST1, a transcription factor that regulates miR-199-214 transcription (Gu and Chan 2012) and inhibits SOX9 transactivation (Gu et al. 2012), is a downstream effector of the WNT signaling pathway (Reinhold et al. 2006), which is key in the determination of female fate and development (Vainio et al. 1999). miR-214 has been described as an accelerator of differentiation by targeting the polycomb repressive complex (Prc2) (Juan et al. 2009), which attenuates RA signaling by methylation of its target promoters (Laursen et al. 2013). Since this miRNA is up-regulated in female PGCs, depletion of PRC2 may enhance the response to RA in females, while in males it would be inhibited. The other member of this cluster, miR-199a, has been described as an inhibitor of pluripotency by enhancing the expression of p53 (Wang et al. 2013). Interestingly, Nanog, a pluripotency associated gene, is naturally down-regulated at early stages of development in female PGCs (Yamaguchi et al. 2005), which could explain these differences in the regulation of pluripotency between males and females during this developmental window.

Opposite to the miR-199-214 cluster, the miR-182-183-96 cluster is up-regulated in male PGCs compared to female PGCs (E11.5 and E12.5). These miRNAs are enhancers of cell proliferation by targeting FOXO proteins (Gebremedhn et al. 2016) while the inhibition of the miR-182-183-96 cluster has been linked to a rise in p53 expression (Tang et al. 2013). Another validated target is Tet1, which regulates meiotic gene expression in mouse oocytes (Yamaguchi et al. 2012). In this sense, reducing the expression of these miRNAs allows higher levels of Tet1 in female PGCs and hence the correct entry into meiosis. Furthermore, in recent studies we have reported a reduction in the expression levels of this cluster in the transition from male PGCs at E13.5 to spermatogonia (García-López et al. 2015). It seems that a reduction in the expression levels of the cluster miR-182-183-96 may be related to differentiating processes in germline.

Pluripotency and differentiation of PGCs is regulated by miRNAs

While comparing the miRNAome with PGCs and their respective gonadal somatic cells, we found a group of differentially expressed miRNAs with mRNA targets related to pluripotency, proliferation, or differentiation. Many members of the cluster miR-290-295, which we detected as one of the most expressed in PGCs, have been associated with pluripotency and are indispensable for PGC development (Hayashi et al. 2008; Medeiros et al. 2011). This correlated with a higher expression of this cluster in PGCs with respect to somatic cells. However, an unexpectedly low expression of this miRNA cluster was detected in E11.5 female PGCs compared to males, which was recovered at E12.5 and E13.5 in most of its members, but not all of them, such as miR-293-3p (Fig. 5F). As mentioned, the miR-290-295 inhibits some regulators of differentiation pathways (Gruber et al. 2014). As female PGCs start the meiotic process at E11.5, the reduction in the levels of expression of this cluster could be expected. After this, its levels are gradually recovered to promote the maintenance of the germline.

In contrast, the let-7 family is up-regulated in somatic cells with respect to PGCs. The let-7 family induces a more differentiated (less pluripotent) state by targeting transcripts such as N-myc and C-myc, while cluster miR-290-295 supports a pluripotent state by targeting p27 and Last2 transcripts (Wang et al. 2013).

These results are in agreement with the state of both cell types: PGCs are more pluripotent with respect to somatic cells of the gonads during this developmental window.

Lastly, we pointed out that the cluster miR-880-881-741-470-465, up-regulated in PGCs at E13.5 when compared to E11.5, targets Oct4; Nanog and Lin28b. OCT4 and NANOG proteins are transcription factors of the miR-290-295 cluster (Marson et al. 2008), and Lin28 is an inhibitor of the miR-let-7 family (Heo et al. 2008). miR-200a-3p, also up-regulated in PGCs at E13.5 when compared to E11.5, is up-regulated during neuronal differentiation and its expression is inversely proportional to Oct4 and Nanog expression. In consequence, the miR-880-465 cluster and miR-200a-3p could be participating in the shutdown of PGC pluripotency from E13.5 in two different ways: promoting the expression of the miR-let-7 family by targeting Lin28b and reducing the expression of the miR-290-295 cluster by down-regulating its upstream regulators (Oct4 and Nanog).

Sex differential miRNA expression

Differential expression analyses revealed miRNAs up-regulated in male gonadal somatic cells at E13.5, when compared to females, thus targeting genes involved in sexual determination and germline development. The miR-30 family was up-regulated from E11.5 to E13.5 in male somatic cells, while in females it was down-regulated in the same period (Fig. 5E). This miRNA family has been reported to be highly expressed in testis (Mishima et al. 2008) and is thought to be key in the adhesion between Sertoli cells and spermatids (Nicholls et al. 2011), which are indispensable in maintaining the equilibrium during postnatal spermatogenesis (Griswold 1995).

Another interesting group of miRNAs that showed a dramatic increase in expression from E12.5 to E13.5, specifically in male somatic cells (Fig. 5D), was miR-202-5p, miR-140-3p, miR-34c-5p, and miR-455. One of them, miR-202-5p, is transcriptionally regulated by Sox9/Sf1, which promotes male gonad differentiation (Kanai et al. 2005; Wainwright et al. 2013). In this sense, this miRNA may be considered a potential candidate as a regulator of male gonad differentiation, being already reported to be enriched in E13.5 gonads (Rakoczy et al. 2013). miR-34c-5p and miR-455-3p have been validated to target the mRNA encoding β-catenin (Ctnbb1). β-Catenin repress the expression of Sox9, which promotes Sertoli cell differentiation in developing gonads (Kanai et al. 2005) and via Wnt4 allows the expression of female markers (Boyer et al. 2012). SRY and SOX9 also inhibit Ctnbb1, as studied in Xenopus and humans (Akiyama et al. 2004; Bernard et al. 2008). In this sense, these miRNAs could be participating together with SOX9in the inhibition of the WNT pathway to promote the male fate of the gonad. In addition, Sry expression could be down-regulated by several miRNAs, specifically in E13.5 male somatic cells since Sry expression reaches its maximum at E11.5 and is then silenced (Kashimada and Koopman 2010). Members of the miR-30 family, miR-22-3p, miR-19a-3p, miR-540-3p, and miR-665-3p, have also been validated to target the transcription factor Six4 (miR-22-3p also has Six1 as target), which promotes the expression of Sry (Fujimoto et al. 2013). It is tempting to consider that these miRNAs could potentially participate in the inhibition of SRY in E13.5 male somatic cells after the initiation of male sex determination and preserving the stability of the male germline.

Finally, studies in null mice for pre-miR-140 demonstrated that this miRNA was able to modulate Leydig cell numbers in testis development (Rakoczy et al. 2013). Together, these data point to the existence of a coordinated and operative regulatory mechanism that can participate in the sex-dependent organization of embryo gonads by complex networks of interactions mediated by miRNAs.

Role of isomiRs and nontemplate additions in developing gonads

With the rise of NGS, miRNA isoforms started to be reported as isomiRs. These variants can enrich miRNA networks, such as 5′ variations with new seed sequences and consequently different targets (Cloonan et al. 2011). The results presented in this work showed that these miRNA variants were also generated between E11.5 and E13.5 in developing gonads and were differentially regulated across the samples analyzed (Fig. 2). Somatic cells of the gonad presented a more stable profile among the different sexes and developmental stages while PGCs were more variable. This was probably due to the dramatic and fast changes that occur in these cells during the developmental window that we studied. Differences between the percentage of sequences and percentage of reads corresponding to 3′ isomiRs and canonical sequences (50% of total reads and <10% of total sequences) indicated that they are not randomly selected. The increase in “noncanonical processing” isoforms in PGCs at E12.5 also supports the hypothesis of the existence of a putative selective mechanisms in miRNA sequences generated during their biogenesis.

Surprisingly, the signature of the proportion of isomiRs per miRNA detected in E13.5 female PGCs was clearly different from all other samples (Fig. 2E). This behavior in isomiR generation was not biased in the sense of being shorter or longer than in the other samples, but the sequences were different. This unusual fact was also coupled with the fact that adenylation was by far the most frequent modification found in the analysis of nontemplate additions (Fig. 3A). We cannot demonstrate a dependence from PAPD4 from this increase in the 3′ addition of adenines because the expression of its mRNA by RT-qPCR was not detected. Recent studies have shown that the 3′ end of miRNAs was important to determine monoadenylation-mediated stability and sensitivity of specific miRNAs (Katoh et al. 2009; D'Ambrogio et al. 2012). Other studies suggest that variations in the 3′ end of a miRNA can alter its stability due to differences in the interaction with the RISC complex and/or exoribonucleases (Hwang et al. 2007). However, in those experiments, the 3′ end of the miRNA was mutated but not trimmed or elongated as we saw in our results. Also, it has been suggested that 3′ adenylation reduces miRNA uptake by AGO2 and AGO3 proteins (Burroughs et al. 2010). Note that the 3′ ends of miRNAs have also been reported to participate in miRNA–mRNA pairing, acting as an auxiliary base-pairing mostly when seed regions did not match perfectly (Moore et al. 2015). Taken together, these data demonstrate that pre-miRNA cleavage and monoadenylation are differentially regulated in E13.5 female PGCs and could potentially affect miRNA stability and miRNA–mRNA interaction. The function and potential consequences of such extraordinary miRNA modification in such a specific cell type, where the cytological meiosis pattern is initiated, is a question to explore further.

Characterization of miRNA biogenesis in the developing gonad

To understand the dynamics of miRNAs in the developing gonad, we assessed the expression of genes involved in miRNA biogenesis by RT-qPCR. Strikingly, the differences in the expression of these genes among the samples analyzed did not correlate with the differences in the miRNA populations. This indicates that maybe the differences observed in miRNA population among PGCs are due to other biological processes such as turnover or stockpile accumulations as has been already described (Wang et al. 2013); and/or by the presence of other RNAs that may act as miRNA sponges or “competing endogenous RNAs” (ceRNAs) (Kartha and Subramanian 2014). Also, the expression of miRNA biogenesis genes, such as Drosha, Dgcr8, Xpo5, and Ago2 (except for Ago2 in E13.5 male PGCs) decays accordingly with the differentiation state of the cells. We previously described something similar during the transition from zygote to blastocyst (García-López and del Mazo 2012).

With respect to Ago genes, Ago2 was the most expressed Ago gene across all the analyzed samples followed by Ago1. AGO2 is the only Argonaute protein with validated clear cleavage function, since the other family members seem to have lost this property (Rand et al. 2005; Valdmanis et al. 2012). Intriguingly, there was a significant difference in the expression of Ago2 and Dicer1 transcripts between male and female E13.5 PGCs, which was much higher in males. We found seven miRNAs up-regulated in E13.5 female PGCs with respect to males with Ago2 and/or Dicer1 as validated targets. Two of them, miR-27a-3p and miR-19a-3p, had both transcripts as a target. We have not been able to find the biological reason of these differences in the expression of those transcripts. However, it could be related to the differences in miRNA processing that we have found in E13.5 female PGCs with respect to the other samples analyzed in this work.

Our current work presents an in-depth analysis of miRNAs both in germline and somatic cells during the gonadal fate determination window in mice. This work clearly revealed a complex miRNA regulation associated with mRNAs with critical roles in the differentiation of male and female gonads, where meiosis is initiated in mammals. We demonstrated that this regulation operates both in PGCs and stromal somatic cells in the gonads in a coordinated manner. Also, gene encoding proteins involved in miRNA biogenesis and 3′ NTA are differentially regulated in this developmental window, especially during PGC development and between male and female PGCs. Finally, E13.5 was demonstrated to be a critical moment in female PGCs’ miRNAome regulation due to the dramatic differences in the generation of 3′ isomiRs and 3′ monoadenylation. All of these results indicate that miRNAs are critical mediators in the complex regulatory mechanisms of early gonadal sex determination and development in mice.

MATERIALS AND METHODS

Animal procedures

Mus musculus CD1 strain was used as the animal model. All procedures relating to the care and handling of the mice used in the present study were carried out in the CIB-CSIC bioterium under specific pathogen-free (SPF), temperature (22°C ± 1°C), and controlled humidity (50%–55%) conditions. All animals were housed on 12 h light–dark cycles with ad libitum access to food and water.

PGC and somatic cells isolation

PGCs and somatic gonadal cells were isolated from embryos obtained from pregnant females at E11.5, E12.5, and E13.5. Sex identification of gonads at E11.5 and E12.5 was performed by the PCR method. Briefly, we isolated and stored individually each set of gonads and embryos in PBS at 4°C. DNA was obtained from each single embryo after boiling somatic fragment of embryo carcass. Sex identification was performed using primers to the Sry gene and Jarid1d(a single band for XX and a double band for XY) (Supplemental Table S10; Clapcote and Roder 2005) and grouping gonads by sex after identification. Gonads of embryos at E13.5 were sexed based on morphological characteristics.

After sexing the embryos, mesonephros was removed before germ and somatic cell separation. PGCs and somatic cells were sorted using a paramagnetic technology according to Pesce and De Felici (1995). Briefly, groups of about 80 gonads were incubated in 0.25% trypsin–EDTA during 15 min at 37°C, washed and followed by incubation with PGC-specific antibodies (anti-CD15) bound to paramagnetic-microbeads (Miltenyi Biotec). Finally, PGCs were isolated from somatic stromal cells by a magnetic column (miniMAcs, Miltenyi Biote) and both separated fractions were stained with the naphtol AS-MX/ FAST-RED (Sigma-Aldrich) (PGC specific) to determine the level of enrichment of cell populations. In all cases the enrichment of the PGCs was 93%–96%. Similarly, the purity of somatic cells was over 90% in all samples.

RNA isolation and quality control

RNA was isolated using TRIzol Reagent (Invitrogen) according to the manufacturer's instructions. RNA integrity and concentration were measured on a NanoDrop ND-1000 spectrophotometer (NanoDrop) and later in a 2100 Bioanalyzer (Agilent), respectively. All RNA used showed a factor of integrity RIN over 8.

Small noncoding RNA sequencing, quality control, and adapter trimming

Small RNAs were sequenced from 1 μg of purified total RNA from each sample (PGCs and somatic cells from male and females at E11.5, E12.5, and E13.5) by a commercial facility company according to the small RNA-seq Illumina protocols. Briefly, RNAs from each different sample were fractionated by electrophoresis in acrylamide gels and the fractions of about 100 nt were isolated. Adapters were ligated to the RNA molecules and a reverse transcription was performed. MiSeq Sequencing System (Illumina) was used in single-end mode with a read length of 75 bp with an average depth of 10 million reads.

The FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) program was used to assess the quality of the raw and trimmed sequencing libraries. Adapters and low quality bases were trimmed using the wrapper script Trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). The quality Phred threshold considered was 28 and the minimum sequence length was 16 nt. Clean reads were collapsed using the Fastx-toolkit suite (http://hannonlab.cshl.edu/fastx_toolkit/index.html). The bioinformatic miRNA analysis pipeline is described in Supplemental Figure S1.

Identification and classification of miRNAs and their isomiRs

MiRNA identification was performed by mapping collapsed reads to mouse pre-miRNA sequences from the miRBase v. 21 database (Kozomara and Griffiths-Jones 2014) using Bowtie 1.1.2 (Langmead et al. 2009; Ziemann et al. 2016) with the following parameters: three mismatches in v mode (-v 3) and -m 20 --strata --best -y --chunkmbs 256 --norc -S. Mismatches in sequence comparison should be required due to potential miRNA editing and 3′ nucleotide nontemplate additions (Landgraf et al. 2007; García-López et al. 2013).

Identification and classification of isomiRs was performed using a custom Perl script, available, together with other scripts used in this work, at https://github.com/dfernandezperez/miRNA-scripts. Our custom Perl script combines the information of two files: a sam file generated by Bowtie and a str file from miRBase v. 21 containing the positions of the mature miRNAs (3p and 5p) over the pre-miRNA. Comparing the alignment results from the sam file (alignment start/end and CIGAR) (Li et al. 2009) with this miRBase file it is possible to classify all miRNA variants, including nontemplate additions.

Normalization and differential expression analysis and functional annotation of miRNAs

Normalization and differential expression analyses were performed using the Bioconductor package DESeq (Anders and Huber 2010). To perform the normalization step in the present work, for each sample: First, collapsed reads were mapped against mm10 genome with Bowtie and then, the fasta files containing the collapsed aligned reads were transformed into a fasta-tabulated file using the script fasta_ formatter with –t flag from fastx_tools. These new files were joined together in one table that was loaded into DESeq to obtain the coefficients for library normalization using the function estimateSizeFactors. Then, the expression data obtained with our custom Perl script was used to perform the differential expression analyses, which were normalized with the coefficients that we obtained previously. The threshold to consider any miRNA differentially expressed was: twofold change and minimum of 100 counts.

To annotate miRNA targets we used databases with only experimentally validated miRNA–target interactions to avoid false positives. MiRNet database (Fan et al. 2016) and Tarbase7 (Vlachos et al. 2015) were downloaded and joined together. Targets of 5′ isomiRs with different seed sequences than canonical miRNAs were predicted locally with TargetScan7 (Agarwal et al. 2015). After target annotation enrichment of Gene Ontology terms (Blake et al. 2015) applying a hypergeometric test (adjusted P-value <0.05, q-value <0.1) was performed with the R/Biocondcutor package clusterProfiler (Yu et al. 2012).

RT-qPCR of mRNAs

RT-qPCR for genes encoding proteins involved in miRNA biogenesis and postranscriptional miRNA modifications was carried out using 125 ng of total RNA, 2.5 µM of Oligo dT17, 0.5 mM of dNTP mix (0.5 mM each), 5 mM 1× SSIV Buffer (Invitrogen), 0.01 M dithiothreitol (DTT), 2 U of RNase inhibitor (RNAsin Promega), and 200 U of SuperScript IV Reverse Transcriptase (Invitrogen) with RNase-free water up to 20 µL, in each reaction. cDNAs were amplified by quantitative PCR using specific primers (Supplemental Table 10). qPCR reactions contained 5 µL of 2× SYBR Green PCR supermix (Roche), 1 µL of cDNA, 0.0625 µM of each specific primer with DNAse free water up to 10 µL. qPCR reactions were performed in a LightCycler 480 system (Roche). Data obtained from qPCR was normalized using the method ΔΔCq (Livak and Schmittgen 2001) taking as reference two different housekeepings: Ppia and U6. The protocol of the qPCR was: 95°C 5 min and 45 cycles at 95°C 15 sec, and 60°C 1 min.

RT-qPCR of miRNAs

To validate NGS data, expression of some representative miRNAs: miR-20a-5p, miR-let-7a-5p, and miR-293-3p were analyzed using TaqMan probes and stem–loop primers for RT (Thermo Fisher Scientific) (Pearson correlation; R = 0.74; P < 0.001; (Supplemental Fig. S4). RT was carried out following manufacturer's instructions. Program: 16°C 30 min, 42°C 30 min, and 85°C 5 min. qPCR protocol was carried out in a LightCycler 480 system (Roche) following manufacturer's instructions. Program: 95°C 15 sec, 60°C 60 sec, 45 cycles. Data obtained from qPCR were normalized using the method ΔΔCq (Livak and Schmittgen 2001) taking the gene U6 as reference.

DATA DEPOSITION

The data sets supporting the conclusions of this article are available in the Gene Expression Omnibus (GEO) repository, under the accession number GSE98713.The data sets corresponding to GO enrichments of this article are included within the article as additional spreadsheet files. The scripts created to analyze miRNA-seq data are available at https://github.com/dfernandezperez/miRNA-scripts.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.