Lambs used in experiments originated from the natural mating of mature Merino ewes kept in the Sheep Research Centre, Urmuqi, Xinjiang Uygur Autonomous Region. Ovarian superstimulation was induced three times (two individuals each time) in a total of 6 donor Merino lambs aged 7 weeks following a previously reported treatment schedule [6]. Briefly, lambs were stimulated with 4×40 mg FSH (Folltropin-V, Bioniche, Belleville, ON, Canada) that was administered in 12 h intervals. Sixty hours after the first FSH injection, the lambs underwent ovariectomy. To minimize genetic noise, the six lambs’ mother ewes formed the control group for gene expression analysis after superstimulation. The adult ewes were implanted with controlled internal drug release devices (CIDR, 0.3 g progesterone, Pharmacia and Upjohn Limited, Auckland, New Zealand) for oestrus synchronization on the first day. On the twelfth day, the ewes were given 4 progressively decreasing doses of FSH (120, 100, 80, and 60 mg) in 12 h intervals (Figure 1). Sixty hours after the first injection, the ewes were ovariectomized and the CIDR was removed.

All experiments were performed in accordance with relevant guidelines and regulations set by the Ministry of Agriculture of the People’s Republic of China.

Ovariectomies were performed using an abdominal surgery approach. Briefly, anesthesia was induced by intramuscular injection of xylidinothiazoline (2 mg/kg body weight). The lower abdomen was shaved, washed and disinfected with an iodine-based detergent solution. A small incision was made in the midline of the abdomen, close to the mammary gland. A hole was punctured in the peritoneum with a surgical knife and then enlarged manually, allowing direct access to and palpation of the reproductive tract. Local anesthesia was applied to the ovarian pedicle using gauze soaked with procaine and a plastic clip was then applied on the pedicle to minimize hemorrhage. Intestinal suture was used to ligate the ovarian pedicle, the ligature was slowly tightened and then the ovary was removed. The ovaries were placed in polyethylene bags, kept on ice, and were transported to the laboratory within 5 min after collection. In order to harvest antral and mural GCs, the growing antral follicles (3 to 4 mm in diameter) were stripped away from the ovaries using ophthalmic tweezers. The follicular fluid was collected by aspiration of the follicular antrum using a 4.5-gauge needle and a syringe. Follicles were flushed three times with Dulbecco’s phosphate buffered saline (D-PBS). The collapsed follicles were cut using a scalpel blade and small tweezers in D-PBS to expose the inner follicular wall; the inner follicular wall was subsequently scraped with a microbiology culture loop to harvest the mural layer of GCs. The mural GCs were collected by centrifugation at 1,500 rpm for 5 min. The cell pellet was resuspended in Sample Protector (TaKaRa, Biotechnology (Dalian) Co., LTD, Dalian, China) and stored at −80°C until RNA extraction.

Each individual sample total RNA was extracted using QIAzol (miRNeasy Mini Kit, QIAGEN, Hilden, Germany) according to the manufacturer’s instructions. Quality of the extracted RNA was assessed by gel electrophoresis and quantity was measured by NanoDrop 2000 spectrophotometer. Purity of extracted RNA was considered good when the ratio of optical density between 260 nm and 280 nm fell between 1.8 and 2.0, while the ratio of values measured at 260 nm and 230 nm was ≥2.0. The RNA integrity value required for RNA sequencing was >8.0. After total RNA extraction (three lambs or three ewes pooled in each group) and DNase treatment, construction of cDNA libraries and sequencing were performed as described previously [17].

After Illumina HiSeq2500 sequencing, raw data were processed using Perl and Python scripts provided by Novogene (Novogene Bioinformatics Technology Co. Ltd, Beijing, China). In this step, the clean data (clean reads) were obtained by removing reads containing adapter, reads containing poly-N and low quality reads from raw data. At the same time, quality parameters of clean data including Q20, Q30, GC-content and sequence duplication level were used for data filtering. All the succeeding analyses were carried out using high quality clean data.

Reference genome and gene model annotation files were downloaded from the sheep genome website at ftp://ftp.ensembl.org/pub/release-76/fasta/ovis_aries/dna/. Paired-end clean reads were aligned to the reference genome using TopHat2 [18]. TopHat was chosen as the mapping tool because it can generate a database of splice junctions based on the gene model annotation file, and thus give a better mapping result than other non-splice mapping tools.

The parameter FPKM (fragments per kilobase of exon per million fragments mapped) was used to evaluate the gene expression levels. HTSeq v0.5.3 (http://www-huber.embl.de/users/anders/HTSeq) was used to count the reads numbers mapped to each transcript. FPKM was calculated based on the mapped transcript fragments, transcript length and sequencing depth. Currently, this is the most commonly used method for estimating transcript expression [19].

Differential expression analysis of two groups was performed using the DESeq R package (1.12.0). DESeq provides statistical routines to determine differential expression in digital gene expression data using a model based on the negative binomial distribution. Absolute value of log2≥1 and q-value<0.005 were considered to reflect significant differences in gene expression [20].

Gene ontology (GO) enrichment analysis of differentially expressed transcripts was implemented using GOseq software package. GO terms with corrected P-value less than 0.05 were considered significantly enriched by differentially expressed transcripts.

A differentially expressed gene (DEG) data set was imported into ingenuity pathway analysis (IPA, http://www.ingenuity.com) to examine functional and molecular pathway enrichment for DEGs based on Fisher’s exact test. Analisis of canonical pathways and regulator effects as well as network analysis were performed as previously described [21]. The activation score (z-score) was used to infer the activation state of transcriptional regulators (z-score<−2 and >2). The expression and/or differential expression of transcriptional regulators in the data set provided evidence for biological mechanisms. Combining lamb and ewe GCs regulator analyses provided clues to the putative regulation involved in the lamb GCs on the follicle development stage.

In order to validate our RNA sequencing data, SYBR green real-time quantitative polymerase chain reaction were performed with the TransStart Green qPCR SuperMix (Transgen, Beijing, China) on a Roche HOLD CYCLE LightCycler 480 II. Three lambs and three ewes that have not been used for the RNA sequencing experiments were selected. Total RNA (1 μg) was used to synthesize first-strand cDNA with the TransScript-Uni One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen, China) with oligo (dT) primer. All qPCR experiments were run in duplicates and included no-RT and no-template controls. The qPCR reactions comprised a denaturation step of 3 min at 95°C followed by 40 cycles of product amplification (5 s at 95°C, 15 s at 60°C, and 15 s at 72°C). Quantities of mRNA corresponding to individual genes were normalized to glyceraldehyde-3phosphate dehydrogenase (GAPDH) expression.

Relative gene expression was calculated with the comparative Ct (ΔCt) method as follows: Ct target gene −Ct GAPDH. ΔCt-values for each target gene were then compared across samples (ΔΔCt) in the following way: ΔCt (sample 1) minus ΔCt (sample 2). Results of real-time quantitative polymerase chain reaction analysis were expressed as fold changes of RNA transcripts and were calculated by the 2^−ΔΔCt method [22].

Unless otherwise stated, all chemicals were obtained from Sigma Chemical Company (St. Louis, MO, USA).

We sequenced the RNA from two groups, superstimulated lamb follicles vs superstimulated (adult) ewe follicles, and performed sequencing quality assessment. In this study, 29627505 and 32774030 raw reads were generated for each sample. After quality control, 3.23 G and 3.56 G clean bases were obtained from the two groups (Supplementary Table S1). Approximately 85% of the total reads were mapped to sheep chromosomes and about 83% of the reads were uniquely mapped to the sheep genome (Supplementary Table S2). Regarding coding characteristics of the mapped reads, 54.1% (ewe) and 59.0% (lamb) corresponded to exon sequences, 8% (ewe) and 15.4% (lamb) to introns, and 26.0% (ewe) and 25.6% (lamb) were matched to intergenic genome regions (Supplementary Figure S1). Read density in chromosomal distrbution and number of mapped reads in chromosomes were shown in Supplementary Figure S2. The correlation of transcripts expression between samples is the most important indicator for reliability of experimental results and rationality of sampling. Generally, the correlation value should be up to 0.92 (r²≥0.92). In our study, the two samples correlation value was up to 0.966. This analysis validates the quality of the library and its utility for further analysis.

To explore the biological mechanism of follicle development following ovarian stimulation in lamb, we sought to identify the DEGs between lamb and ewe GCs. Comparing the two samples, we identified 311 DEGs with the fold change ≥2 and q-value<0.005 used as cutoff values (Supplementary Table S3). Of these, 175 genes (42 novel) were down-regulated and 136 genes (5 novel) were up-regulated in lamb versus ewe GCs (Figure 2). Among them we identified 33 DEGs with more than four-fold change in expression, including 20 downregulated and 13 upregulated in lamb versus ewe GCs (Table 1).

We hypothesized that the DEGs we identified were important for follicle development. The 311 identified DEGs were classified into three major functional categories: cellular component, molecular function, and biological process. The top five functional groups included binding, catalytic activity, cellular process, metabolic process, and organic substance metabolic process (Figure 3). Further enrichment analysis was related to cellular functions and subcellular locations, and the detected DEGs were enriched in different terms related to follicle development. For example, the DEGs cytochrome P450, family 19 (CYP19A1), steroidogenic acute regulatory protein (STAR), neuropilin 1 (NRP1), EGF containing fibulin-like extracellular matrix protein 1 (EFEMP1) and proliferating cell nuclear antigen (PCNA) were enriched in binding, lipid binding, protein binding, cellular process, catalatic activity and metabolic process, respectively (Supplementary Table S4).

The DEG list was subjected to IPA and pathway, network and upstream regulator analyses of lamb versus ewe were performed. Among 311 DEGs, the 234 genes (119 up-regulated and 115 down-regulated) have been annotated (Supplementary Table S3) and subjected to pathway network and upstream regulator analyses in IPA.

Forty-two IPA canonical pathways were identified in lamb compared to ewe with threshold −log (p-value) ≥ 1.3 (Supplementary Table S5). The fourteen pathways were with values greater than 2 (Figure 4). Among these pathways, the top three IPA canonical pathways were base excision repair (BER) pathway (ligase I [LIG1, DNA, ATP-dependent], PCNA, polymerase [DNA directed], epsilon, catalytic subunit [POLE], [flap structure-specific endonuclease 1]), cell cycle control of chromosomal replication (MCM3, MCM4, MCM5, MCM6 (minichromosome maintenance complex component 3, 4, 5, 6), CDC6 [cell division cycle 6]) and GADD45 (growth arrest and DNA-damage-inducible) signaling (PCNA, GADD45G, BRCA1 [breast cancer 1], CDK1 [cyclin-dependent kinase 1]) (Figure 3). The top three pathways by the number of included DEGs were axonal guidance signaling SEMA3A (semaphorin 3A), SEMA6D (semaphorin 6D), PAPPA (pregnancy-associated plasma protein A), NFAT5 (nuclear factor of activated T-cells 5), PTCH1 (patched 1), NTRK1 (neurotrophic tyrosine kinase, receptor, type 1), GNB2 (guanine nucleotide binding protein [G protein], beta polypeptide 2), TUBA1C (tubulin, alpha 1c), ADAMTS9 (ADAM metallopeptidase with thrombospondin type 1 motif, 9), ADAMTS4 (ADAM metallopeptidase with thrombospondin type 1 motif, 4), NRP1, SEMA7A (semaphorin 7A), protein ubiquitination pathway CDC20 (cell division cycle 20), USP27X (ubiquitin specific peptidase 27, X-linked), USP53 (ubiquitin specific peptidase 53), DNAJB11 (DnaJ [Hsp40] homolog, subfamily B, member 11), DNAJC27 (DnaJ [Hsp40] homolog, subfamily C, member 27), UBE2S (ubiquitin-conjugating enzyme E2S), NEDD4L (neural precursor cell expressed, developmentally down-regulated 4-like, E3 ubiquitin protein ligase), BRCA1, HSPA5 (heat shock 70 kDa protein 5), UBE2C (ubiquitin-conjugating enzyme E2C), and LPS (lipopolysaccharide)/IL-1 (Interleukin 1) mediated inhibition of retinoid × receptor (RXR) function MGST1 (microsomal glutathione S-transferase 1), ALDH1L2 (aldehyde dehydrogenase 1 family, member L2), IL1RL1 (interleukin 1 receptor-like 1), SREBF1 (sterol regulatory element binding transcription factor 1), CHST11 (carbohydrate [chondroitin 4] sulfotransferase 11), ALAS1 (aminolevulinate, delta-, synthase 1), CHST15 (carbohydrate [N-acetylgalactosamine 4-sulfate 6-O] sulfotransferase 15. In contrast, there were six pathways with only one DEG each, namely andamide degradation (FAAH), taurine biosynthesis (CDO1 [cysteine dioxygenase type 1]), glutathione biosynthesis (GCLC), arginine degradation I (arginase pathway) (OAT, ornithine aminotransferase), L-cysteine degradation I (CDO1), and myo-inositol biosynthesis (ISYNA1) (Supplementary Table S5).

On the IPA function analysis, these DEGs have been involved in 500 function cases (Supplementary Table S6). There were 17 functional groups with activation z-score greater than 2 and predicted activation state activated (Table 2). The most significantly affected functions in lamb follicles were found to be proliferation of cells, DNA replication, cellular movement, organ development and tissue development (Table 2). IPA upstream regulator analysis was used to identify upstream transcriptional regulators. We had identified 32 activated upstream regulators, including transcription regulators (MYC [v-myc avian myelocytomatosis viral oncogene homolog], E2F1 [E2F transcription factor 1], E2F2 [E2F transcription factor 2], TBX2 [T-box protein 2], CCND1 [cyclin D1], SREBF1, GATA4 [GATA binding protein 4], STAT6 [signal transducer and activator of transcription 6], MYOD1 [myoblast determination protein 1], XBP1 [X-box binding protein 1]), growth factors (IGF1 [insulin-like growth factor 1], HGF [hepatocyte growth factor], EGF, VEGFA [vascular endothelial growth factor A], INHA [inhibin, alpha]), enzymes and kinases (HRAS [Harvey rat Sarcoma viral oncogene homolog], CAT [catalase], KRAS [Kirsten rat Sarcoma viral oncogene homolog], ERBB2 [erb-b2 receptor tyrosine kinase 2]), and others (Table 3). The top activated upstream regulator was the growth factor EGF. The EGF receptor upstream regulator target molecules in the dataset included ACTB (actin, beta), ANGPT2 (angiopoietin 2), BIRC5 (baculoviral IAP repeat containing 5), CAV1 (caveolin 1), CDCP1 (CUB domain containing protein 1), CDK1, CLK1 (CDC-like kinase 1), CTGF (connective tissue growth factor), CYP19A1, FAP (fibroblast activation protein), FOSL2 (FOS-like antigen 2), HK2 (hexokinase 2), MYBL2 (v-myb avian myeloblastosis viral oncogene homolog-like 2), NOTCH1 (notch 1), PCNA, POSTN (periostin), and THBS1 (thrombospondin 1). They interacted within a mechanistic network which included Akt (serine/threonine-specific protein kinase), CTNNB1 (catenin [cadherin-associated protein], beta 1), EGFR (epidermal growth factor receptor), ERK1/2 (extracellular regulated protein kinases), ESR1 (estrogen receptor 1), FOXO1 (forkhead box O1), HIF1A (hypoxia inducible factor 1, alpha subunit), HRAS, JUN (Jun Proto-Oncogene), MYC, NFkB (complex), P38 MAPK (mitogen-activated protein kinase), PI3K (phosphatidylinositol-4,5-bisphosphate 3-kinase, complex), PPARG (peroxisome proliferator-activated receptor gamma), RELA (v-rel avian reticuloendotheliosis viral oncogene homolog A), RUNX2 (runt-related transcription factor 2), SMAD3 (SMAD family member 3), SP1 (specificity protein 1), STAT3 (signal transducer and activator of transcription 3), and TP53 (tumor protein P53). In addition, 23 upstream regulators were predicted to be inhibited, including transcription regulators (RB1 [retinoblastoma 1], RBL1 [retinoblastoma-like 1], TP53, CDKN2A [cyclin-dependent kinase inhibitor 2A], E2F6 [E2F transcription factor 6], SMARCB1 [SWI/SNF related, matrix associated, actin dependent regulator of chromatin, subfamily B, member 1], HDAC2 [histone deacetylase 2], SIM1 [single-minded family BHLH transcription factor 1], ARNT2 [aryl-hydrocarbon receptor nuclear translocator 2]), kinases (CDKN1A [cyclin-dependent kinase inhibitor 1A], DYRK1A [dual-specificity tyrosine-(Y)-phosphorylation regulated kinase 1A], TGFBR2 [transforming growth factor, beta receptor II]), a phosphatase (PTEN), micro-RNAs (miR-24-3p, let-7), and others (Table 4). Upstream regulator let-7 predicted target molecules including cell division cycle family members (CDC20, CDC6, CDCA3 [cell division cycle associated 3]), minichromosome maintenance complex components (MCM3, MCM4, MCM5, MCM6), ribonucleotide reductase (RRM1 [ribonucleotide reductase M1], RRM2 [ribonucleotide reductase M2]) and others transcripts (BRCA1, CCNA2 [cyclin A2], CDK1, THBS1). The top inhibited upstream regulator was TP53.

The networks generated by IPA based on the dataset above are shown in Figure 4. The network in Figure 4A shows a role for cellular movement, cell death and survival. There is considerable connectivity associated with Akt, which exerts direct effects on the extracellular matrix molecules and indirectly affects collagen. The other network in Figure 4B shows a role for cardiovascular system development and function indicating significant interaction between the extracellular matrix genes and ADM (adrenomedullin), which appears to mainly signal through the intracellular PI3K pathway. The extracellular matrix genes PTX3 (pentraxin 3), NOV (nephroblastoma overexpressed), EFEMP1, OXT (oxytocin-neurophysin 1 precursor), and ADM and membrane protein gene NRP1 were all downregulated in these networks in lamb.

To confirm the reliability of the high-throughput sequencing data, we used qPCR to compare the expression levels. We selected 8 differently expressed genes, two upregulated (CYP19A1, STAR) and 6 downregulated (EFEMP1, SPP1 [osteopontin, secreted phosphoprotein 1], PTX3, NRP1, ADM, OXT), for validating. The results of the qPCR analysis for these 8 genes were in a good agreement with the RNA-seq data. The primers used for each gene and the fold changes are summarized in Table 5.