All of the animal procedures were approved by the Institutional Animal Care and Use Committee (IACUC) at the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, which is where the experiments were conducted. All methods were carried out in accordance with the relevant guidelines and regulations set by the Ministry of Agriculture of the People's Republic of China.

Small Tail Han sheep (Han) and Dorset sheep (Dorset) were bred by the Ao-Te sheep breeding farm in Qingdao (Shandong, China). Blood samples were obtained from the Han ewes to identify the FecB mutation in the BMPR1B gene. Among these ewes, the high-fecundity Han sheep were the adult Han ewes with the genotype BB (HanBB) and the genotype ++ (Han ++), and they were regarded as the two groups. In addition, the adult Dorset ewes were sacrificed to obtain the ovary samples and were considered the low-fecundity group and served as controls. All of the animals were raised under similar conditions and had free access to water and food in natural lighting. Three animals from each group (Dorset, HanBB, and Han++) were used for the experiments described below.

All of the experimental ewes were administered with intravaginal sponges (40 mg, Chronogest, Intervet, Federal District, Mexico) that were impregnated with fluorogesterone acetate. Ten days later the sponges were removed, and then, the estrus synchronization, the ewes received an i.m. injection of pregnant mare serum gonadotropin (Ningbo HormoneCo., Ningbo, China) at a dose of 400 IU (Quintero-Elisea et al., 2011). Estrus was confirmed, and all of the sheep that reached spontaneous estrus after an estrous cycle were sacrificed between 24 and 36 h after the spontaneous estrus was detected. The animals were all between 2.5- and 3-year-old when they were euthanized for ovary harvesting. Whole ovaries were removed, and in order to obtain better ovulation points on the surface of the ovaries, samples were collected. All of the samples (3 animals per group) were immediately snap-frozen in liquid nitrogen and were stored at −80°C until the RNA isolation.

The total RNA from the ovaries was isolated using the TRIzol reagent (Invitrogen Carlsbad, CA), and the RNA integrity number (RIN) for all of the samples was >8. An additional RNase-free DNase I (Ambion, Inc., Austin, TX, USA) digestion step was conducted to ensure that the samples were not contaminated with genomic DNA. After isolating the total RNA from the samples (3 Dorset, 3 Han BB and 3 Han ++), we constructed whole transcriptome libraries using TruSeq Stranded Total RNA with Ribo-Zero Globin (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. The quality control and quantification of the libraries were performed using the BioAnalyzer 2100 system and qPCR (Kapa Biosystems, Woburn, MA, USA), respectively. Initially, the resulting libraries were sequenced on a HiSeq 2000 instrument, which produced 100 nucleotide paired-end reads. Additionally, small RNAs, between 18 and 30 nt long, were purified from the total RNA from each sheep group. Then, we constructed small RNA libraries that represented each species using the Illumina® TruSeq™ Small RNA Sample Preparation protocol. Briefly, the total RNA was ligated with an RNA 3′ adaptor and 5′ adaptor. Next, the RNA was reverse transcribed with SuperScriptII Reverse Transcriptase (Invitrogen). To obtain the miRNA library, the cDNA was amplified by PCR and was then purified. The Qubit™ dsDNA HS kit and a Qubit® 2.0 Fluorometer was used to assess the quality of the library. The size and purity of the library were assessed using High Sensitivity DNA Chips and the Agilent 2100 system. Clusters were generated using the cBot of Illumina Genome Analyzer IIx according to the manufacturer's instructions. Finally, the RNA libraries were sequenced using the single-read multiplex program in the Shanghai Biotechnology Corporation. The data collection software from Illumina was used to control the sequencing process.

Using the TopHat default settings, the RNA-Seq reads from each sample were aligned to the oar 3.1 sheep reference genome (Trapnell et al., 2009). The transcript annotation of the sheep genome sequence (Release number: Oar3.1) was downloaded (ftp://ftp.ncbi.nlm.nih.gov/genomes/Ovis_aries/).

The gene expression analysis only included the uniquely mapped reads. Based on the previously described rigorous significance test for digital gene expression, the significantly differentially expressed genes (DEGs) were identified using the DESeq package (Anders, 2010). The error rate adjustment in the multiple significance tests was assessed using the false discovery rate (FDR) (Benjamini and Hochberg, 1995). If the fold change was >1.5 or <0.667 and the FDR was <0.05, then the genes and the lncRNAs were considered to be differentially expressed.

Genomic DNA was isolated from each group (3 Han BB, 3 Han ++ and 3 Dorset ovaries) by a phenol-chloroform extraction method. Briefly, the DNA was sheared by sonication on ice to generate random fragments ranging from 100 to 500 bp. The sonicated DNA (5 μl) was used for the immunoprecipitation (IP). Next, the recovered DNA was 5′ and 3′ end blunted, phosphorylated and was repaired by T4 polynucleotide kinase and T4 DNA polymerase (NEB, USA). Then, the DNA was ligated to an Illumina sequencing primer after the addition of an ATP in the 3′ end. Our MeDIP method was the same as the one described by Zhang et al. (2006). After MeDIP, the DNA was PCR amplified with sequencing primers provided by Illumina, and the PCR products were recovered and used for sequencing.

We downloaded the sheep reference sequence and mapping information from the University of California Santa Cruz Genome Bioinformatics Site (http://genome.ucsc.edu). The sequencing reads and the resulting files were aligned to the sheep reference genome. Model-based analysis for ChIP-Seq (MACS) was applied for peak detection and analysis of the IP sequencing data.

The identified methylated DEGs were analyzed with GO (Ashburner et al., 2000), which is an international standardized gene functional classification system. The DEGs were annotated by three structured controlled ontologies, including biological process (BP), molecular function (MF), and cellular component (CC), which were used to maintain a consistent description of the gene products. In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to determine the association of the genes with different pathways (Kanehisa and Goto, 2000). To select the significant GO terms and pathways, a Fisher's exact test followed by the Benjamini-Hochberg (BH) multiple testing correction was performed to calculate the threshold of significance. The GO terms and pathways with a p < 0.05 were considered significant.

The small RNA reads acquired by the Illumina Hiseq 2000 were subjected to several filtering processes, including low quality read filtering and removing the adaptor sequences. Finally, only short trimmed reads sized from 18 to 30 nt were retained. To annotate and classify the small RNAs into different categories and to identify novel miRNAs, the small RNA reads were compared with the known non-coding RNAs deposited in the Rfam database. The rest were subjected to RNAfold, which identifies miRNA candidates according to the hairpin secondary structure (Hofacker et al., 1994). The identified miRNAs were regarded as novel candidate miRNAs if they were absent in the miRBase database (Griffiths-Jones et al., 2006). To study the regulatory relationship between the mRNA and the miRNA, an mRNA-miRNA regulatory network was constructed and was displayed using Cytoscape software (Smoot et al., 2011).

To generate the co-expression network, methylated genes that encoded the DEGs and the differentially expressed lncRNAs were selected. Meanwhile, the differentially expressed lncRNAs with length of 200–800 bp were used to construct a network. The Pearson correlation was calculated for each pair of genes that were analyzed, and the pairs (only lncRNA-mRNA) with significant correlations were selected to construct the network. Only the strongest correlations (0.99 or greater) were included in the network to generate a visual representation. Cytoscape software was applied to draw the co-expression networks (Smoot et al., 2011). Finally, the degree of nodes in each network was calculated.

Meanwhile, a weighted gene co-expression network was constructed to identify the interactions among the methylated genes that encoded the DEGs and the differentially expressed lncRNAs by applying the WGCNA approach (Zhang and Horvath, 2005; Wang et al., 2013). The weighted Pearson correlation matrices, corresponding to the gene or lncRNA expression, were first calculated to generate the co-expression network, which was followed by the standard procedure of the WGCNA (Wang et al., 2009). Briefly, the weighted correlation matrices were transformed into connection strength matrices. Then, the connection strengths were used to calculate the topological overlap, which is a powerful and biologically meaningful measurement of the strength of the relationship of the co-expression of two genes with other genes in the network (Yip and Horvath, 2007). The genes that had highly similar co-expression relationships were grouped using hierarchical clustering based on the topological overlap. The Dynamic Hybrid Tree Cut algorithm was used to generate the hierarchal clustering tree, and the modules were defined as branches from the tree cutting (Langfelder et al., 2008). The minimum spanning tree with a dissimilarity matrix from the WGCNA was used to generate the network for each module. The modules were initially denoted by colors, and the gray color denoted the background genes outside of the modules. The highly correlated genes in each module were summarized by the first principal component, which was referred to as the module eigengene. Modules with a p < 0.1 were merged for further analyses. To identify the genes in the significant modules associated with proliferation, we calculated the intramodular connectivity (k) by summing the connection strengths with other module genes and dividing the number by the maximum intramodular connectivity. The genes with the highest intramodular connectivity were referred to as hub genes.

To validate the identified differentially expressed transcripts, miRNAs and lncRNAs in the three groups, we performed real-time PCR using SYBR Green I Master Mix in a Roche LightCycler 480 II Real Time PCR system. 5S rRNA (for miRNA) and 18S rRNA (for mRNA and lncRNA) were used as internal controls in these experiments. The PCR was performed using the following program: 95°C for 10 min; 40 cycles of 95°C for 15 s with a melting temperature for 30 s; and 72°C for 45 s; and 72°C for 5 min.

The DEG data and the bioinformatics analyses are described above. The data from the qPCR are presented as the mean ± SD. Comparisons between the groups (the Han group vs. the Dorset group) from the qPCR were made using a Student's t-test, and p < 0.05 was considered statistically significant.

The mRNA libraries that were derived from the Han BB, Han ++ and Dorset sheep ovaries were sequenced and analyzed. Functional genetic differences among the three groups of sheep were assessed by identifying the DEGs in the Han ++ vs. Dorset, the Han BB vs. Dorset and the Han BB vs. Han ++ pairwise groups (Tables S1–S3). Compared with the Dorset group, 2,511 genes were differentially expressed in the Han BB group (Table S1), and 1,280 genes were differentially expressed in the Han ++ group (Table S2). Compared with the Han ++ group, 2,241 genes were differentially expressed in the Han BB group (Table S3). To identify candidate lncRNAs that might be related to the proliferation of sheep, we focused on the differentially expressed lncRNAs in the three comparison groups. Compared with the Dorset group, there were 138 lncRNAs that were differentially expressed in both the Han BB and ++ groups. There were 106 lncRNAs that were differentially expressed in the Han BB group compared with the Han ++ group.

Next, we performed methylated DNA immunoprecipitation sequencing (MeDIP-Seq) with an Illumina Hiseq 2000 to reveal the genome-wide DNA methylome of the sheep. The MeDIP-Seq reads were aligned to a reference sequence, and for scanning the methylation peaks, we only used the unique mapped reads. A range of 112,892,114 to 114,963,972 raw reads was generated for the three sheep groups. More than 93% of the reads were mapped, and about 55% of the reads in each group were uniquely mapped to the sheep genome for each sample (Table 1).

In addition, the distribution of the DNA methylation in the regions 2 kb upstream 2 kb downstream of the transcription termination sites, in the CpG island and in the intragenic region were analyzed. In general, higher levels of DNA methylation were revealed in the CpG islands and the gene body regions compared to the upstream and downstream regions. In sheep, the DNA methylation level increased sharply in the CpG islands and the intragenic region (Figure S1). To assess the effect of methylation on transcriptional activity and its relation to sheep reproduction, we analyzed the differentially methylated genes based on the peak coverage within the different elements in the three comparison groups (Figure S2). We found that the number of differentially methylated genes in the three comparison groups was similar.

To identify whether there was a relationship between the methylated genes encoding the DEGs and the differentially expressed lncRNAs, the expression data were compared with the MeDIP data. Compared with the Dorset group, among the DEGs and the differentially expressed lncRNAs, 73 were methylated in the Han BB group, and 129 were methylated in the Han ++ group. When compared with Han ++, there were only 7 methylated genes in the Han BB group.

To highlight important biological processes related to the methylated DEGs among the three groups of sheep, we performed a GO enrichment analysis. The methylated DEGs in the Han BB group compared with Dorset were mainly associated with rhythmic process, lipopolysaccharide metabolic process and metal ion homeostasis (Figure 1A). A comparison of Han ++ group to the Dorset group revealed highly enriched GO terms that were related to development, cell differentiation and signal transduction (Figure 1B). A comparison between the Han BB and Han ++ groups revealed categories associated with biosynthetic process and immune response (Figure 1C). The GO terms for CC and MF were compared among the 3 sheep groups, and the results were similar (Figure 1).

To further explore the pathways enriched by the methylated DEGs, we carried out a KEGG enrichment analysis. A comparison of the genes in the Han BB group with those in the Dorset group, showed pathways that were significantly enriched in secretion, metabolism and protein functions (Table S4). When the Han ++ groups were compared to the Dorset group, the DEGs were mainly related to endocytosis, metabolism and secretion (Table S5). A comparison of the Han BB group with the Han ++ group revealed that the genes were associated with metabolism (Table S6).

To investigate whether miRNAs play a key role in sheep fecundity, the small RNAs were extracted from the three sheep groups and were sequenced using the Illumina platform. A systematic analysis enabled the identification of 132, 98, and 100 differentially expressed and previously unreported miRNAs in the Han BB vs. Dorset, the Han ++ vs. Dorset and the Han BB vs. Han ++ comparison groups, respectively. In addition to identifying novel miRNA, we also identified previously reported miRNAs. A total of 10, 6, and 16 known miRNAs were considered to be differentially expressed in the Han BB vs. Dorset, the Han BB vs. Han ++, and the Han ++ vs. Dorset comparison groups, respectively. Next, the Miranda program was used to predict the putative targets of the differentially expressed miRNAs. To study whether the miRNAs play a regulatory role in sheep fecundity, a miRNA-mRNA regulatory network was constructed (Figure 2).

Methylated coding-non-coding gene co-expression networks for the Han sheep and the control were constructed with all of the genes encoding the methylated DEGs and the differentially expressed lncRNAs identified in this study. The methlyated lncRNA-mRNA pairs with Pearson correlation coefficients that were equal to or >0.99 were selected to construct the co-expression network (Figure 3). The number of nodes was similar in each comparison, but there was different in degree of nodes (Tables S7–S9). In the Han BB vs. Dorset comparison, the lncRNAs tetratricopeptide repeat domain 26 (TTC26) and myosin heavy chain 15 (MYH15) had the largest degree, but the degree of these lncRNAs was different between the Han BB and Dorset groups. Similarly, the lncRNA neuronal tyrosine-phosphorylated phosphoinositide-3-kinase adaptor 1 (NYAP1) had the largest degree in the Han ++ vs. Dorset comparison. However, none of the methylated lncRNAs were co-expressed with either an mRNA or a methylated mRNA in the Han BB vs. Han ++ comparison. These methylated lncRNAs may play a vital regulatory role in sheep breeding.

To further understand the co-expression relationships between the genes that were significantly correlated to sheep breeding at a systems level, we applied the weighted gene co-expression network analysis (WGCNA), which analyzed the RNA-Seq datasets using a gene co-expression network analysis approach that is an unsupervised and unbiased analysis. The WGCNA identified distinct co-expression modules with respect to clusters of the correlated genes or lncRNAs (Figure 4A). Notably, 8 out of the 13 co-expression modules were significant, and these modules consist of genes that correspond to sheep breeding (Figure 4B). These modules were labeled by color. For example, the black and yellow modules correlated with a high significance with the Dorset samples, but the purple module was strongly correlated with Han sheep sample. Additionally, to determine the GO terms that were enriched in the genes of all of the modules, DAVID was applied. Table 2 shows the top 3 GO biological processes of the modules that are associated with sheep breeding. Importantly, the hub genes and the lncRNAs in the eight modules that might be useful for sheep breeding are listed in Table 3.

To validate the RNA-Seq data, 6 methylated genes encoding differentially expressed mRNAs, 2 methylated genes encoding differentially expressed lncRNAs and 4 differentially expressed miRNAs were randomly selected and analyzed via real-time PCR (Figure 5, Table S10). For the methylated genes encoding the differentially expressed mRNAs, the results demonstrated that ankyrin repeat and BTB (POZ) domain containing 1 (ABTB1) were up-regulated and that the glutamate receptor, ionotropic, AMPA 1 (GRIA1), paired box 8 (PAX8), leucine rich colipase-like 1 (LRCOL1), NGFI-A binding protein 2 (NAB2), and LOC101116109 genes were down-regulated in the Han sheep group in comparison to the Dorset group. Additionally, with regard to the lncRNAs, MYH15, and NYAP1 were down-regulated in the Han BB sheep group in comparison to the Dorset group. For the miRNAs, the expression of oar-miR-494-3p, oar-miR-496-5p, oar-miR-376c-3p, and oar-miR-3959-5p showed statistically significant differences between the Han sheep and the Dorset group (p < 0.05). These results validated the RNA-Seq data, providing reassurance that the above findings can be trusted.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.