Identification and Profiling of microRNAs in Goat Endometrium during Embryo Implantation

Background MicroRNAs (miRNAs) are short, highly conserved small noncoding RNAs that had fundamental roles in post-transcriptional gene expression, and they are crucial for proper control of biological processes and known to participate in embryo implantation. However, miRNA expression profiles in the pre-receptive and receptive phases of the goat endometrium during embryo implantation are unknown. Results A total of 1,069 and 847 miRNAs were expressed in receptive (R) and pre-receptive (P) goat endometrium, and 632 miRNAs were co-expressed in both phases. We identified 545 (50.98%) known miRNAs in the R library and 522 (61.63%) in the P library. There were 110 up-expressed miRNAs and 33 down-expressed miRNAs in receptive endometrium compared with the pre-receptive endometrium meeting the criteria of P-values< 0.05. Moreover, GO and KEGG analysis of the target genes of the differentially expressed miRNAs revealed some candidate miRNAs, genes and pathways that may involve in the formation of the receptive endometrium. Based on stem-loop RT-qPCR, 15 miRNAs were detected and the results suggested that the majority of the miRNA expression data measured by Solexa deep sequencing could represent actual miRNA expression levels. Conclusions Our data revealed the first miRNA profile related to the biology of the goat receptive endometrium during embryo implantation, and the results suggested that a subset of miRNAs might play important roles in the formation of endometrial receptivity. Thus, elucidating the physiological roles of endometrial miRNAs will help us better understand the genetic control of embryo implantation in goats.


Introduction
Embryo implantation is a crucial step in mammalian reproduction [1]. Endometrium sensitivity to embryo implantation is classified as a pre-receptive phase and a receptive phase in mammals [2,3]. The development of a receptive endometrium is a complex process and appropriate endometrial preparation with stromal proliferation and epithelial differentiation is stimulated prior to embryo-uterine interactions [4], which is a spatial and temporal phenomenon known as the 'window of implantation' [5,6]. During this restricted period, a folded endometrial epithelial bilayer develops and the endometrium acquires the adhesive properties that allow embryo adhesion and its subsequent invasion [7].
Some factors act directly on the endometrium, and regulated endometrial functions promote the interactions between endometrium and embryo [8]. Accumulated evidence suggests that hormones [9], homeobox proteins [10], morphogens [11] and other factors participate in the preparation of endometrial receptivity [12]. Moreover, molecular studies have extensively investigated the expression and regulation of different factors connected with endometrial receptivity, with the list of possible genes involved in the establishment of a receptive endometrium increasing exponentially in rats, pigs, and humans [13,14]. Gaining a solid understanding of the molecular mechanisms underlying endometrium remodeling is important for investigation into reproduction of female goats, which is the most economically important component in commercial goat breeding [15,16].
Endometrial receptivity is a dynamic process that was believed to be under post-transcriptional regulation because several genes expressed in the endometrium were identified as being epigenetically regulated [17]. MicroRNAs (miRNAs), with an average length between 18 nt and 26 nt [18], are a new class of small non-coding RNAs that function as negative regulators by repressing gene expression, inhibiting translation or regulating mRNA degradation [19]. miRNA recognition is based on either incomplete or complete complementarity to target sequences within the 3' untranslated region (3'UTR) at the post-transcriptional level [7,20]. Between 1 and 5% of genes are predicted to encode miRNA, and these miRNA may regulate the expression of as many as 30% of mRNA [21,22]. Mature miRNA are synthesized through a multi-step process that concludes with the cleavage of the pre-miRNA stem-loop by the RNase III enzyme Dicer1. A single miRNA can target hundreds of mRNAs [23]. Thus, miRNA are potential regulators of endometrial receptivity at the post-transcriptional level via targeting of mRNAs for degradation or translational repression [19,20]. Differentially expressed miRNAs have been identified in the mouse, rat and human uterus [14,24], as well as the porcine placenta [25,26].
Nevertheless, little information was available on miRNA expressed in the goat endometrium, and differences between the pre-receptive and receptive phases during embryo implantation were unknown. Focusing on these points, we aimed to gain more understanding of the complex regulation of endometrial receptivity in healthy multiparous dairy goats by analyzing miRNA expression in the endometrium at gestational day 5 (pre-receptive endometrium phase) and gestational day 15 (receptive endometrium phase) using Illumina Solexa technology. And then we investigated the known miRNAs, potential novel miRNAs, differentially expressed miRNAs and their target genes. Given the biological functions of miRNAs at the posttranscriptional level and genes usually interact with each other to play their roles, GO enrichment and KEGG pathway analysis were also made in this study. The results would provide further understanding of the miRNAs role in the regulation mechanism of the endometrial receptivity.

Overview of sequencing data
To systematically identify small RNAs and changes in the expression level of miRNAs between the pre-receptive and receptive phases of the endometrium in Xinong Saanen dairy goats, we purified and sequenced small RNAs from the goat endometrium. A total of 7,717,603 and 6,307,668 raw reads were obtained from the R (receptive) and P (pre-receptive) endometrium, respectively. To assess the efficiency of Solexa sequencing and the quality of the sequences, all reads were annotated and classified by aligning against the miRBase 20.0 database (ftp:// mirbase.org/pub/mirbase/CURRENT/), the Rfam database (http://rfam.janelia.org), the Repbase database (http://www.girinst.org/repbase), the Genome database (http://goat.kiz.ac. cn/GGD/download.htm) and the mRNA database (http://goat.kiz.ac.cn/GGD/download.htm). Small RNAs were classified into different categories according to their annotations. Next, 3' ADT and length filter, junk reads, Rfam, mRNA, repeats, rRNA, tRNA, snRNA and snoRNA sequences, as well as other Rfam RNA sequences, were separated out and discarded. As a result, 4,233,874 clean reads representing 284,622 unique sequences from the R library were obtained, and 5,386,427 clean reads representing 139,040 unique sequences from the P library (Table 1). The lengths of the majority of the clean reads ranged between 18 and 26 nt. The most abundant class size in the small RNA sequence distribution was 22 nt, which accounted for 38.54% and 51.43% of reads in the R and P libraries, respectively, followed by 21 nt (17.03% and 18.27%), 23 nt (13.65% and 11.33%) and 20 nt (8.81% and 11.68%), suggesting that they were typical small RNA (Fig 1, Table 2).
The clean reads were subjected to advanced bioinformatic analyses and divided into four groups. (1) Group 1b, containing 14 miRNAs corresponding to 8 pre-miRNAs and miRNAs in miRBase 20.0; the pre-miRNAs further map to the genome and EST. (2) Group 2a, containing 215 miRNAs corresponding to 174 pre-miRNAs and miRNAs in miRBase 20.0; the mapped pre-miRNAs do not map to the genome, but the reads (including the miRNAs of the pre-miR-NAs) map to the genome. The extended genome sequences from the genome loci may form  (Table 3). From this analysis, 1,284 miRNAs were identified, and 632 of these were co-expressed in both libraries, 1,069 miRNAs were detected in the R library, 847 were detected in the P library (Fig 2).

Identification of Known miRNAs
To identify known miRNAs in goat endometrium, the dataset was compared with the known mammalian miRNAs (miRNA precursors and mature miRNAs) in miRBase 20.0 (ftp://mirbase. org/pub/mirbase/CURRENT/). Sequences with a perfect match or one mismatch were retained in the alignment. We identified 545 known miRNAs in the R library and 522 in the P library (S1 Table). Considering the conservation of mature miRNAs among various species, the sequences of the existing miRNAs in goat were aligned and conservatively analyzed to investigate their evolutionary relationships (S2 Table). The expression levels of oar-miR-10a_R+1_1ss12TA (a variant of oar-miR-10a as described previously in the same way [27]), oar-miR-10b_L+1R-1, oar-miR-26a, oar-miR-7a and aja-miR-   143_1ss22GT were predominately higher, with more than 100,000 reads in both the R and P libraries (Table 4). These miRNAs constituted the most abundant read in the receptive endometrium, suggesting they potentially play an important role in the formation of receptive endometrium in goats. In addition, five known miRNAs (miR-449a, miR-182, miR-187-3p_R+1, miR-183-_L-1, miR-200a-5p) were detected by stem-loop qRT-PCR because they highly and differently expressed between P and R libraries in this study.

Identification of potential novel miRNAs
Sequencing reads that did not match any of the known miRNA were further analyzed to discover potential novel miRNAs. To determine whether these un-annotated small RNA reads were genuine miRNA, their hairpin structures, dicer cleavage sites and minimal free energies were explored using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi/). A total of 643 miRNA candidates were identified as having the typical stem-loop secondary structure, with 488 candidates identified in the R library and 298 in the P library (S3 Table). There were 33 mature miRNAs with NE numbers (Normalized expression) greater than 10 in the libraries ( Table 5). The length of these miRNA sequences ranged from 18 nt to 24 nt. Moreover, 11 of 33 potential novel miRNA candidates were new 3p-derived sequences, and 22 were new 5p-derived sequences.

Differential expression of miRNAs in receptive and pre-receptive endometrium
Pair-wise comparisons revealed important differentially expressed miRNAs between each developmental phase. Although some of their expression quantities were equivalent, there were some miRNAs expressed differently between the pre-receptive and receptive phases. We focused our attention on miRNAs meeting our designated criteria of P-values< 0.05 and |log2 (fold change)|> 1, and a total of 143 differentially expressed miRNAs were selected from the P and R libraries ( Table 6). There were 33 differentially expressed miRNAs that were down-regulated in the receptive endometrium compared with pre-receptive endometrium in goats (fold change 0.5), 231 miRNAs with fold changes were greater than 0.5 and less than 2 (0.5 <fold change 2), and 110 miRNAs with fold changes greater than 2 (fold change> 2) (Fig 4). It deserved to note that there were some miRNAs expressed specifically, such as bta-miR-431_R-3 was detected only in the receptive endometrium and bta-miR-216a_R-2 in the pre-receptive endometrium ( Table 6).
The most differentially expressed miRNA was hsa-miR-449a, which had a 113.2-fold (NE> 1,000) increase in the receptive endometrium compared with the pre-receptive endometrium. Other miRNA that were differentially expressed between our two libraries included bta-miR- Note: rep_miR ID was the representative known miRNA ID. R-NE represented the normalized expression level of miRNAs in the small RNA library generated from receptive endometrium of Xinong Saanen dairy goats. P-NE represented the normalized expression level of miRNAs in the small RNA library generated from pre-receptive endometrium of Xinong Saanen dairy goats.

GO enrichment and KEGG pathway analysis
We used a well-established miRNA-target dataset to investigate the possible function of these miRNAs with the aid of gene function annotation methods, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) progress. GO enrichment analysis using cellular components showed that 490 genes (13.64%) mapped to GO terms for component topology in the database. The analysis of biological processes showed that 2,087 genes (58.10%) were involved in cellular or metabolic processes, and the analysis of molecular function showed that 1,015 genes (28.26%) were assigned different functions based on gene background (Fig 5 and S5 Table). GO biological process analysis based    on the predicted targets showed that the differentially expressed miRNAs were involved in remodeling of the receptive endometrium, including induction of apoptosis by extracellular signals (GO: 0008624) (Fig 6). Other miRNA-gene networks of interest included cell cycle checkpoints (GO: 0000075), vesicle-mediated transport (GO: 0016192), intracellular protein transport (GO: 0006886) and maternal processes involved in female pregnancy (GO: 0060135). KEGG pathway annotation showed 1,886 target genes that were annotated for 220 biological processes (S6 Table). Moreover, there were 25 KEGG pathways with P values less than 0.05 (Fig 7). KEGG pathway analysis based on predicted targets revealed that differentially expressed miRNAs were involved in several pathways affecting the development of endometrium, including ubiquitin mediated proteolysis (84 annotated genes), T cell receptor signaling pathway (53 annotated genes), Fc gamma R-mediated phagocytosis (52 annotated genes) and apoptosis (52 annotated genes) (S6 Table).

Confirmation and comparison of miRNA data
To date, several miRNAs have been identified in the goat [29,30], but the information of miR-NAs on the development of the receptive endometrium remained quite limited. In the present study, we conducted a detailed overview of experimentally identified miRNAs during the 'window of implantation,' which is the receptive phase when the endometrium was permissive to embryo implantation [3,31]. The analysis resulted in 7,717,603 and 6,307,668 small RNA raw reads obtained from the R (receptive) and P (pre-receptive) endometrium libraries, respectively. After separating out and discarding junk sequences, we obtained 4,233,874 clean reads representing 284,622 unique sequences from the R library and 5,386,427 clean reads representing 139,040 unique sequences from the P library (Table 1). The length distribution showed that more than 90% of the small RNA sequences were primarily distributed in the 19-24 nt range in both libraries (Fig 1), which was consistent with the typical size of mature mammalian miRNA [32]. These results were also consistent with the typical small RNA distribution in mammals, such as cattle [33], sheep [34] and pigs [35]. The size distribution of the small RNAs from receptive and pre-receptive goat endometrium was similar, but there were differences in the length distribution and abundance between the two libraries. Furthermore, the proportion of total rRNA was used for a quality check of the samples. rRNA proportions are usually less than 60% in high quality plant samples [36] and 40% in high quality animal samples [28]. The proportion of total rRNA in our samples was 7.80% and 3.04% in the R and P libraries, respectively, indicating that the collected endometrium samples were of high quality. The above Note: R-NE represented the normalized expression level of miRNAs in the small RNA library generated from receptive endometrium of Xinong Saanen dairy goats. P-NE represented the normalized expression level of miRNAs in the small RNA library generated from pre-receptive endometrium of Xinong Saanen dairy goats. Fold (R /P) and log2 (R/P) indicated the fold change of the miRNAs between samples. P-Value manifested the significance of miRNAs differential expression between two samples. Number of Targets represented putative target genes for a certain miRNA predicted by two software programs (Target Scan v.50 and Mireap v.3.3a).
doi:10.1371/journal.pone.0122202.t007 microRNAs Profiling of Goat Endometrium during Embryo Implantation results suggest that the high-throughput sequencing data were highly enriched for small RNA sequences. Next, we aligned the clean reads to the miRNAs of all known animals in the miRBase 20.0 database (ftp://mirbase.org/pub/mirbase/CURRENT/). A total of 1,284 miRNAs were identified: 632 of these were co-expressed in both libraries, 1,069 miRNAs were detected in the R library and 847 in the P library. Compared with previous research, more conserved miRNAs (1,069 miRNAs in the R library and 847 in the P library) were identified in this study [28,37]. There are two principal explanations for these results. The first explanation is that the detected clean reads were aligned to the latest database of miRBase 20.0 released in June 2013 in this manuscript, and 3,355 new hairpin sequences and 5,393 new mature miRNAs were added to the database based on miRBase 19.0 [28]. The second explanation is that the two libraries were compared with all animals in miRBase 20.0 because few goat miRNA sequences are included, and thus the reference data were not species-specific.

Differential expression of miRNAs in receptive and pre-receptive endometrium of goats
To investigate differences in endometrium activity in the receptive and pre-receptive phases, differentially expressed miRNA were identified in the two constructed libraries. In the present study, 292 miRNAs were up-regulated and 150 miRNAs were down-regulated in the receptive phase compared with the pre-receptive phase (P< 0.05), and there were 110 miRNAs increased White circular nodes represent genes, and red rectangular and pink rounded rectangle nodes represent miRNAs. The size of the nodes represents the power of the interrelation among the nodes, and edges between two nodes represent interactions between genes. The more edges a gene has, the more miRNAs that interact with it, and the more central a role it had within the network. The top five key miRNAs in the network were miR-449a, miR-182, PC-5p-6424_145, PC-5p-2673_284, and miR-138-5p. The top three key mRNAs' serial numbers were JR128745.1 (Caspase 8, apoptosis-related cysteine peptidase (CASP8)), JR133467.1 (Y-linked ubiquitin-specific protease 9 (USP9Y)), and JR132030.1 (CD5 molecule (CD5)) in NCBI (National Center for Biotechnology Information).
doi:10.1371/journal.pone.0122202.g006 microRNAs Profiling of Goat Endometrium during Embryo Implantation at least a 2-fold in expression. Xia reported that 28 miRNAs were up-regulated at least two-fold and 29 miRNAs were down-regulated at least two-fold during the receptive phase compared with the pre-receptive phase in the rat uterus [12]. Chakrabarty reported that 32 miRNAs were significantly up-regulated at least 1.5-fold and five miRNAs were down-regulated at least 1.5-fold on day 4 of gestation (receptive phase) compared with day 1 of gestation (pre-receptive phase) in the mouse uterus [3]. There are three possible explanations for the increase in differentially expressed miRNAs identified in this study. The first explanation is the differences in species, as the experimental animal in this study was the dairy goat, while previous studies focused on mice and humans [12,24], whose endometrium structure was different from goat [38]. The second explanation is that the miRNAs in our two libraries were compared with all animals in the miRBase 20.0, which lacks goat miRNAs, so the reference data were not speciesspecific [28]. Third, in contrast to traditional miRNA identification methods such as direct cloning and computational prediction, Solexa deep sequencing was performed to discover species-specific or poorly expressed miRNAs. This technique has been widely utilized to identify conserved and novel miRNAs in various species [39,40]. Thus, more potential miRNAs were identified with the development of new sequencing technology. The miRNA with the largest differential increase in expression was hsa-miR-449a, which had a 113.2-fold increase in the receptive phase compared with the pre-receptive phase. These results suggest that it might play an important role in regulation of endometrium receptivity in goats. Previously, miR-449a had been identified in various types of cancer tissues where it plays a tumor-suppressive role [41], in part through targeting HDAC1 and activating p27 expression [42]. Lize found miR-449 potently induced apoptosis and up-regulated p53 activity [43]. To date, several targets of miR-449a had been identified, such as CDK6, CDC25a, E2F1, CCND1 and BCL2 [44,45]. Moreover, Paik reported that down-regulation of miR-449a and subsequent up-regulation of CCND1 and BCL2 was a novel mechanism for cell proliferation [46]. miR-449a directly bound to the seed sequence of the LEF-1 (lymphoid enhancer-binding factor-1) 3' UTR, causing effective repression of its expression and ultimately leading to a subsequent reduction in Sox 9 gene expression [47]. Altogether, these results suggested that miR-449a potentially participated in regulating dynamic changes in goat uterine gene expression patterns that occur during the transition from the pre-receptive to the receptive phase. These result need to be further validated under well-controlled conditions in animal models.
In addition, bta-miR-182 aroused our interest for it was the highest expressed miRNA in receptive endometrium (NE = 3518.63) and increased 14.55-fold compare with pre-receptive endometrium. Studies have identified that miR-182 was significantly up-regulated in endometrial Values were "means ± SD" in receptive endometrium (n = 9) that were relative to pre-receptive endometrium (n = 9), whose values all were 1 in this manuscript because of the method of 2 -ΔΔCt , ΔΔ Ct = (CT miRNA − CT 18s ) R − (CT miRNA − CT 18s ) P .
doi:10.1371/journal.pone.0122202.g008 microRNAs Profiling of Goat Endometrium during Embryo Implantation carcinoma tissues (EC) compared with complex atypical hyperplasia, simple hyperplasia and normal endometrial tissues [48][49][50]. Further study suggested that miRNA-182 binds directly to a conserved 8 bp sequence in the 3'-UTR of its target gene transcription elongation factor Alike 7 (TCEAL7), and then promots cell proliferation by targeting the tumor suppressor gene TCEAL7 and modulating the activity of its downstream effectors c-Myc, cyclin D1 and NFκB in EC cell lines compared with normal endometrial epithelial cells [51]. What's more, considering that EC was an estrogen-dependent malignancy and that miRNAs were shown to be regulated by estradiol [52], the association between miR-182 and receptive endometrium needs further investigation.
The specific expressed miRNAs in receptive and pre-receptive goat endometrium Next, we turned our attention to the specific expressed miRNAs. The miRNA with the highest tissue-specific expression was bta-miR-431_R-3 from the R library. After careful analysis, we found that its sequence was consistent with miR-431, which was initially identified as a central nervous system specific miRNA cloned from the brain tissue of mouse embryos [53]. Wu reported that miR-431 expression induced by nerve injury stimulates regenerative axon growth by silencing Kremen1, an antagonist of Wnt/beta-catenin signaling. Both the gain-of-function of miR-431 expression and knockdown of Kremen1 significantly enhanced axon out growth in murine dorsal root ganglion neuronal cultures [54]. This research suggests that miR-431 might participate in the regulation of endometrial receptivity through the nervous system.
The highly tissue-specific miRNA in the P library was bta-miR-216a_R-2, whose sequence was identical to miR-216a with one mismatch. This miRNA regulates Ybx1 post-transcriptionally, which is a key mechanism for enhanced mRNA translation in kidney cells [55]. Another target of miR-216a was shown to be the tumor suppressor in lung cancer-1 gene (TSLC1) mRNA through three target sites in its 3' untranslated region [56]. Moreover, miR-216a plays a relevant role in the pathogenesis of cardiovascular disorders and atherosclerosis [57]. Thus, miR-216a potentially regulates the decidualization of the endometrium in the pre-receptive phase, but this hypothesis requires further study.

Identification of potential novel miRNAs
One of the most important characteristics that distinguishes miRNAs from other endogenous small RNAs is the capacity of the precursor miRNA sequence to fold back into a canonical stem-loop hairpin structure [58]. In this study, a total of 643 miRNA candidates were identified as having the typical miRNA stem-loop secondary structure. Precursor sequences ranged from 50 nt to 123 nt in length and could all be shaped into representative stem-loop structures, with the mature miRNA located at either the 5 0 or the 3 0 end. Additionally, the GC content of pre-miRNA ranged from 23% to 77%, and the free energy (dG) ranged from -80.5 kCal/mol to -17.7 kCal/mol; these factors were considered because unstable pre-miRNA structures were required to produce mature single-strand miRNAs [59]. There were 33 novel miRNA candidates with lengths ranging from 18 nt to 24 nt, with the majority at 22 nt. The numbers of these miRNA were shown to be greater than 10. Other miRNAs were not considered for this study because of their low abundance in the goat endometrium.

GO enrichment and KEGG pathway analyses
Previous studies showed that miRNAs modulated gene expression by inhibiting mRNA translation or regulating mRNA degradation at the posttranscriptional level based on Watson-Crick pairing between the 5' end of the miRNA (2-8 nt, the "seed" region) and the 3' un-translated regions (3' UTR) of target mRNAs [60]. An animal miRNA was predicted to target at least 200 genes on average because of the short length of the seed region [61,62]. However, 7,608 annotated mRNA transcripts were predicted as putative target genes for the 143 differentially expressed miRNAs in the present study, suggesting a miRNA targeted 53 genes on average. The small number of predicted target genes may be due to a lack of relevant data on goat genes compared with other animals (i.e., human and rat).
GO is an international standardized classification system for gene function, supplying a set of controlled vocabulary to comprehensively describe the properties of genes and gene products [28]. GO analysis provides insight into the molecular functions of genes in various biological processes [63]. In the present study, GO enrichment analysis revealed that 58.10% of the genes were involved in biological processes, 13.64% of the genes were annotated to intracellular component ontology, and 28.26% of the genes referred to molecular functions. Moreover, the miRNA-gene network further revealed core genes involved in remodeling of the endometrium, such as induction of apoptosis by extracellular signals, cell cycle checkpoint, vesicle-mediated transport, intracellular protein transport and maternal processes involved in female pregnancy.
In organisms, genes usually interact with each other to play different roles in certain biological functions [28]. KEGG pathway analysis could facilitate the understanding of the biological functions of genes [64,65]. A critical event during embryo implantation is the extensive tissue remodeling at the maternal-fetal interface, which was characterized by cell proliferation, cell migration and cell adhesion [1,66]. In the present study, ubiquitin mediated proteolysis, natural killer cell mediated cytotoxicity, T cell receptor signaling pathway, Fc gamma R-mediated phagocytosis, and apoptosis were involved in the significantly enriched miRNA-associated pathways in receptive goat endometriums. These significantly enriched pathways might imply that the organism was coping with the criteria of the endometrium phase during the 'window of implantation' in goats. In summary, GO and KEGG analyses provide a better understanding of the cellular components, molecular functions and biological processes of target genes [30], and provided a reference for future research

Validation of miRNA expression with Stem-loop RT-qPCR
We performed stem-loop qRT-PCR, which is a reliable method used to detect, measure and validate the expression levels of miRNAs. The expression trends of 8 miRNAs measured by Solexa deep sequencing were consistent with those determined by stem-loop qRT-PCR, but the differences in expression levels between the two methods did not agree well. Moreover, the expression patterns of miR-200a measured by Solexa deep sequencing were opposed to those determined by stem-loop RT-qPCR.
This inconsistency might result from differences in the two methods. The expression levels of miRNAs measured by Solexa deep sequencing required the following seven steps: RNA extraction, small RNA fractionation, 5' adaptor ligation, 3' adaptor ligation, RT-PCR, sequencing, and data analysis; thus, the information from some sequences might be lost or produce lowquality tags that were filtered during data analysis [58]. Moreover, some miRNAs might be hard to sequence because of their physical properties [67]. Stem-loop RT-qPCR was a new method designed to detect and quantify mature miRNAs in a fast, specific, accurate, and reliable manner [58]. It was designed to quantify miRNAs using three steps: RNA extraction, Stem-loop RT (reverse transcription) and RT-qPCR. Moreover, the stem-loop RT primer contained a highly stable stem-loop structure that lengthened the target miRNA, and the RT product was amplified and monitored in real time by use of a specific forward primer and a universal reverse primer [68]. Thus, Solexa deep sequencing was considered inferior to RT-qPCR in terms of miRNA quantification [58]. Nevertheless, the trends of expression of 12 miRNAs detected by the two methods in our study were consistent, suggesting that the majority of the miRNA expression data could represent actual miRNA expression levels.

Conclusion
We obtained high-quality miRNA expression profiles from pre-receptive endometrium and receptive endometrium, and predicted 545 known and 33 novel miRNAs in endometrium of dairy goats for the first time. Target gene predictions for the 143 differentially expressed miRNAs, functional annotations and pathway analyses in GO and KEGG databases could contribute to a better understanding of the miRNA mediated regulation of target genes in the development of the endometrium receptivity.

Ethics statement
All animals in this study were maintained according to the No. 5 proclamation of the Ministry of Agriculture, P. R. China. And animal protocols were approved by the Review Committee for the Use of Animal Subjects of Northwest A&F University.

Study design and tissue collection
A total of 20 healthy 24-month-old multiparous dairy goats (Xinong Saanen) were induced to estrous synchronization for this study. The first day of mating was considered to be Day 0 of pregnancy. Gestational days 5 and 15 were important time points for the embryo implantation in goats [38]. The experimental goats were observed 3 times daily to ascertained estrous sign and mated naturally twice in estrus. And then the goats were euthanized when the goats lost consciousness caused by intravenous injection of barbiturate (30mg/kg) at gestational day 5 (pre-receptive endometrium) and gestational day 15 (receptive endometrium). Endometrium samples from 10 goats at gestational day 5 and 10 goats at gestational day 15 were obtained from the anterior wall of the uterine cavity. All tissue samples were washed briefly with PBS (Phosphate Buffered Saline) and then immediately frozen in liquid nitrogen.

Small RNA library construction and sequencing
Total RNA was extracted from every sample using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure. The total RNA quantity and purity were analysis of Bioanalyzer 2100 (Agilent, CA, USA) and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number> 7.0. Approximately 1 ug of total RNA were used to prepare small RNA library according to protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA). The total RNA with lowest quality was not used for further study from the pre-receptive endometrium samples and receptive endometrium samples, respectively. P library was constructed by mixing the nine samples into three samples with the same concentration, performing the single-end sequencing (36 bp) on an Illumina Hiseq2500 (Illumina, San Diego, USA) thrice at the LC-BIO (Hangzhou, China) following the vendor's recommended protocol, and R library was constructed in the same way.
Data processing followed the procedures described in a previous study [69] by LC Sciences Service. Briefly, the raw reads were subjected to the Illumina pipeline filter (Solexa 0.3), and then the dataset was further processed with an in-house program (ACGT101-miR, LC Sciences, Houston, Texas, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA and snoRNA) and repeats.

Identification of known miRNAs
Unique sequences 18-26 nt in length were mapped to specific species precursors in miRBase 20.0 (ftp://mirbase.org/pub/mirbase/CURRENT/) by BLAST search to identify known miR-NAs and novel 3p' and 5p' derived miRNAs. The remaining sequences were mapped to other selected species precursors (with the exclusion of specific species) in miRBase 20.0 by BLAST search, and the mapped pre-miRNAs were further BLASTed against the specific species genomes to determine their genomic locations. The above two were defined as known miRNAs. On the basis of analysis for detected miRNAs, we further analyzed the conservatism of these miRNAs with selected species on frequency for the number in statistics.

Identification of potential novel miRNAs
Sequencing reads that did not match any known miRNAs were further analyzed to discover novel miRNAs using the characteristic hairpin structure of the miRNA precursor [28,58]. To determine whether these unmapped small RNA reads were genuine miRNA, the unmapped sequences were BLASTed against the specific genomes, and the sequences containing hairpin RNA structures were predicated from the flanking 80 nt sequences using RNA fold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi/). The criteria for secondary structure prediction were described by Mi and Lian [27,70].

Analysis of differential expressed miRNAs
To compare differentially expressed miRNAs in the receptive endometrium (R library) and pre-receptive endometrium (P library) of goats, the expression abundances of miRNAs in the two samples were normalized to obtain the expression of transcripts per 1,000,000. And the formulae was: Normalized expression (NE) = (Actual miRNA count/Total count of clean reads) × 1,000,000.
When the NE of a certain miRNA was 0 in one of the two libraries, we revised the 0 to 0.001 for the comparative, and the miRNA whose NE was less than 1 in both libraries were discarded described in detail previously [28]. The differential expression of miRNA based on normalized counts was analyzed using Fisher's exact test, and the significance threshold was set to be 0.05 for each test. The fold-change and P-value for each miRNA were calculated based on the normalized expression using the formulae as described previously [28] and shown below: Fold change ¼ Log2 ðR À NE=P À NEÞ: P-value formula: Cðy y min jxÞ¼ X y y min y¼0 pðyjxÞ Dðy!y min jxÞ¼ X 1 y!y max pðyjxÞ N1 and X represent the total number of clean reads and normalized expression level of a given miRNA in the small RNA library generated from P library, respectively, and N2 and Y represent the total number of clean reads and normalized expression level of a given miRNA in the small RNA library generated from the R library, respectively.

Target genes prediction of differential expression miRNAs
Two computational target prediction algorithms (TargetScan 5.0 chttp://www.targetscan.org/ cgi-bin/targetscan/data.download.cgi?db=vert.61 and miRanda 3.3a chttp://www.microrna. org/microrna/home. do) were used to predict target genes of differentially expressed miRNAs. Only when the target was identified by both programs, it was considered to be the predicted target for a given miRNA.
GO enrichment and KEGG pathway analysis of target genes GO enrichment analysis was performed on predicted target gene candidates of differentially expressed miRNA. This method mapped all target gene candidates to GO terms in the database (ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz), calculated gene numbers for each term and used hyper-geometric tests to find significantly enriched GO terms in target gene candidates. The results were compared with the reference gene background and to genes corresponding to certain biological functions.
The calculating formula: The N represented the number of GO annotated genes in genome, n represented the number of differentially expressed genes in N. M represented the number of particular GO annotated genes in genome, m represented the number of particular GO annotated genes expressed differentially in M [30]. And then P was less than 0.05 was used as the threshold to judge significant enrichment GO term in this study.
KEGG pathway analysis facilitates the understanding of biological functions of genes [28]. The analysis revealed the main pathways in which the target gene candidates were involved. KEGG pathway analysis (http://www.genome.jp/kegg/) used the same formula as GO analysis to determine target gene candidates. Here, N was the number of all genes with KEGG annotation, n was the number of target gene candidates in N, M was the number of all genes annotated in a certain pathway and m was the number of target gene candidates in M [30].
The miRNA-gene network was constructed using Cytoscape Software [71] (Cytosca-pe_v2.8.1, http://www.cytoscape.org/) to analyze the interactions between miRNA and genes. Firstly, the relationship and the interactions between miRNAs and genes were determined by using their differential expression values in the databases of TargetScan 5.0 and miRanda 3.3a. And then we used the random variance model (RVM) corrected t-test (P< 0.05) to identify the differently expressed genes and miRNAs from the gene expression profiles. The overlap between the genes whose expression was induced and the target genes whose expression was reduced by miRNA were then chosen to construct the network based on the two respective expression values by using the miRNA-mRNA modules via population-based probabilistic learning methods [72,73]. A = [ai, j] was defined as the adjacency matrix of miRNA and genes based on the attribute relationships between target genes and miRNA, represented the relationship weight between gene i and miRNA j [74].

Validation of differential expression miRNAs using Stem-loop RT-qPCR
Stem-loop qPCR was used to validate the known and potential novel miRNAs in present manuscript. Total RNA was extracted using Trizol reagent (TaKaRa, Dalian, China) following the manufacturer's instructions. Total RNA was converted to cDNA with a Stem-loop primer using the Prime Script RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China) according to the manufacturer's instructions. The reverse transcription reaction system were: added 5 μl total RNA (800 ng), 2 μl 5× gDNA Eraser Buffer, 1μl gDNA Eraser and RNase-Free dH 2 O to a final volume of 10 ul, incubated the mixture at 65°C for 3 min, and then added 4μl 5× Prime Script Buffer 2 (for real-time PCR), 1μl PrimeScript RT Enzyme Mix 1, 1μl Stem-loop primer and 4 μl RNase-Free dH 2 O to a final volume of 20 ul, incubated the mixture at 42°C for 15 min, followed by 85°C for 5 sec. The cDNA products were diluted 1:10 (v/v) with sterile water.
RT-qPCR was performed using the Bio-Rad CFX 96 Real Time Detection System and SYBR Green PCR Master Mix (TaKaRa, Dalian, China) in a 20 μl reaction according to the manufacturer's instruction. The PCR mixture included 4 ul cDNA for each miRNA, 2 μl primer mix (10 uM), 10 μl 2× SYBR Green Mix with ROX, and 4 μl RNase-Free dH 2 O. The reaction mixtures were amplified in a 96-well plate at 95°C for 10 min, followed by 40 cycles of 95°C for 10 s, 60°C for 1 min, and all reactions were performed in triplicate. 18S rRNA was used as the reference, and all primers for the Stem-loop RT-qPCR are shown in S7 Table. The relative expression levels of the miRNAs were calculated using the equation N = 2 -ΔΔCt, ΔΔ = (CT miRNA − CT 18s ) R − (CT miRNA − CT 18s ) P . The Ct (cycle threshold) was defined as the cycle number at which the fluorescence intensity passed a predetermined threshold [75].

Statistical analysis
All the data were processed with SPSS 17.0 (SPSS Inc., Chicago, IL, USA). One-way ANOVA was used to compare the differences, and the method of the least significant difference (LSD) was used for further analysis, and the differences were considered significant when P was <0.05 and very significant when P was <0.01.
Supporting Information S1