Integrated Analysis of mRNA and miRNA Expression Profiles in the Ovary of Oryctolagus cuniculus in Response to Gonadotrophic Stimulation

Molecular mechanisms responsible for gonadotrophic control of ovarian follicle development and ovulation have not been fully delineated. In this study, prepubertal female rabbits were subjected to a combined PMSG/hCG treatment for the induction of follicle maturation and ovulation. Ovaries of 6 does at different time points during gonadotrophic stimulation were collected for histomorphological examination and genome-wide analysis of miRNA and mRNA transcriptomes, and the plasma were separated for detecting melatonin (MT), prostaglandin E2 (PGE2), estradiol (E2), and progesterone (P4) levels. The results suggested that PMSG promoted the development of the reproductive tract by decreasing plasma levels of E2 and slightly increasing those of MT and PGE2 and that hCG induced ovulation and corpus luteum formation by significantly increasing MT, PGE2, and P4 levels. At the transcriptomic level, a total of 1,122 differentially expressed genes (DEGs) and 12 DE miRNAs were identified using three-group comparisons. Meanwhile, pairwise comparisons revealed that 279 and 103 genes as well as 36 and 20 miRNAs were up- and down-regulated during PMSG-stimulated follicle development while 11 and 5 genes as well as 33 and 16 miRNAs were up- and down-regulated during hCG-induced luteinization. KEGG enrichment analysis of the DEGs derived from both three-group- and two-group comparisons as well as the predicted target genes of DE miRNAs highlighted the crucial roles of pathways involving tissue remodeling, energy metabolism, and regulation of cellular functions in mediating gonadotrophin-induced follicle maturation. Specifically, 3 genes including the matrix metallopeptidase 13 (MMP13), protein phosphatase 1 regulatory subunit 3C (PPP1R3C), and solute carrier family 2 member 12 (SLC2A12), together with 2 miRNAs including the miR-205-1 and miR-34c, were predicted to be the promising downstream targets of both PMSG and hCG. Significantly, the miRNA-mRNA interaction pairs containing top 10 up- and down-regulated mRNAs/miRNAs upon PMSG/hCG stimulation were established, and so were those involved in the PI3K-Akt, ECM-receptor interaction, and focal adhesion pathways during PMSG-induced follicle maturation. Finally, qRT-PCR analysis confirmed the results from RNA-Seq and Small RNA-Seq. Our work may contribute to a better understanding of the regulatory mechanisms of gonadotrophins on ovarian follicle development and ovulation.

Molecular mechanisms responsible for gonadotrophic control of ovarian follicle development and ovulation have not been fully delineated. In this study, prepubertal female rabbits were subjected to a combined PMSG/hCG treatment for the induction of follicle maturation and ovulation. Ovaries of 6 does at different time points during gonadotrophic stimulation were collected for histomorphological examination and genome-wide analysis of miRNA and mRNA transcriptomes, and the plasma were separated for detecting melatonin (MT), prostaglandin E 2 (PGE 2 ), estradiol (E 2 ), and progesterone (P 4 ) levels. The results suggested that PMSG promoted the development of the reproductive tract by decreasing plasma levels of E 2 and slightly increasing those of MT and PGE 2 and that hCG induced ovulation and corpus luteum formation by significantly increasing MT, PGE 2 , and P 4 levels. At the transcriptomic level, a total of 1,122 differentially expressed genes (DEGs) and 12 DE miRNAs were identified using three-group comparisons. Meanwhile, pairwise comparisons revealed that 279 and 103 genes as well as 36 and 20 miRNAs were up-and down-regulated during PMSG-stimulated follicle development while 11 and 5 genes as well as 33 and 16 miRNAs were up-and down-regulated during hCG-induced luteinization. KEGG enrichment analysis of the DEGs derived from both three-group-and two-group comparisons as well as the predicted target genes of DE miRNAs highlighted the crucial roles of pathways involving tissue remodeling, energy metabolism, and regulation of cellular functions in mediating gonadotrophin-induced follicle maturation. Specifically, 3 genes including the matrix metallopeptidase 13 (MMP13), protein phosphatase 1 regulatory subunit 3C (PPP1R3C), and solute carrier family 2 member 12 (SLC2A12), together with 2 miRNAs including the miR-205-1 and miR-34c, were predicted to be the promising downstream targets of both PMSG and hCG. Significantly, the miRNA-mRNA interaction pairs containing top 10 up-and down-regulated mRNAs/miRNAs upon PMSG/hCG stimulation were established, and so were those involved in the PI3K-Akt, ECM-receptor interaction,

INTRODUCTION
Paramount activities of the mammalian ovary are characterized by steroid hormone production as well as follicle maturation and ovulation, which are coordinately regulated by a number of factors emanating from extra-ovarian tissues and the ovary per se (1). With regard to the factors of extra-ovarian origin, adenohypophyseal hormones [e.g., follicle-stimulating hormone (FSH) and luteinization hormone (LH)] play pivotal roles in regulating ovarian activities (2). In support of this, the canonical two-cell/two-gonadotrophin mechanism controlling follicular steroidogenesis has been demonstrated in the mammalian ovary, where LH directs the androgen synthesis of theca cells while estrogens are aromatized from androgens in FSHregulated granulosa cells. Furthermore, the preovulatory LH surge promotes the luteinization of granulosa cells, which then acquire the ability to secrete progesterone (3). Thus, it is of great value to unravel the mechanisms by which gonadotrophins orchestrate mammalian follicle development and ovulation.
The past century has seen a great advance in understanding such mechanisms. Although primordial follicle activation and preantral follicle development are independent of FSH but rely on factors secreted by the oocyte and granulosa cells, FSH is indispensable for the development of immature follicles into a preovulatory phenotype via activation of several downstream signaling cascades in granulosa cells, including the classical Gα s /adenylyl cyclase (AC)/cyclic AMP (cAMP)/protein kinase A (PKA), inositol triphosphate (IP3)/diacylglycerol (DAG)dependent protein kinase C (PKC), and phosphoinositide 3-kinase (PI3K)/protein kinase B (Akt) pathways (4,5). Cooperation of these pathways induces developmental stagespecific expression patterns of target genes related to granulosa cell proliferation (e.g., CREB, CCND2, and CDKs), differentiation (e.g., EGFs, BMPs, and TGFs), and steroidogenesis (e.g., LHR, CYP11A1, and CYP19) as well as formation and expansion of a fluid-filled antrum (e.g., AQPs and HAS2) (6)(7)(8). By contrast, LH triggers ovulation and luteinization of preovulatory follicles by downregulating cell proliferation-related genes (e.g., FSHR, IGF1, and CCND2) and upregulating genes related to cumulus cell expansion (e.g., HAS2, PGS2, and COX2), follicular rupture (e.g., PR, ADAMTS1, CTSL, and MMPs), and granulosa cell luteinization (e.g., STAR, CYP11A1, LHR, and C/EBPβ), which are mainly mediated through the extracellular signal-regulated kinase 1/2 (ERK1/2) pathway in a cAMP/PKA-dependent manner (9,10). These findings suggest that gonadotrophins modulate follicle development and ovulation by inducing both developmental stage-and cell type-specific expression of numerous genes related to different cellular functions depending on the complex interactions of multiple signaling cascades. Nevertheless, it remains a huge challenge to depict the global gonadotrophin-stimulated gene expression profile changes during follicle development and ovulation and to decipher the signaling networks.
In addition to protein-coding genes, microRNAs (miRNAs), a class of endogenous ∼22 nt-long noncoding RNAs, also regulate follicle development and ovulation, since they are widely expressed in various ovarian cell types and regulate oogenesis and folliculogenesis by targeting a myriad of genes related to cellular proliferation, differentiation, luteinization, etc. (11)(12)(13). Rabbit (Oryctolagus cuniculus) not only serves as an economically important animal by providing superior meat, wool, and fur products but also an experimental animal model in biomedical researches. However, to date, information on global changes in the landscape of miRNA expression as well as on the miRNA-mRNA interaction networks during gonadotrophin-stimulated follicle development and ovulation remains scarce in rabbits. To achieve this, in the present study, the dynamic process of follicle maturation and ovulation was in vivo reproduced in sexually immature rabbits using validated hormonal treatments able to promote antral follicle development (PMSG injection) and ovulation (hCG treatment). Subsequently, gonadotrophin-stimulated changes in circulating hormones, ovarian morphology, and histology as well as mRNA and miRNA transcriptomes profiling were evaluated through the application of enzyme linked immunosorbent assay (ELISA), haematoxylin and eosin (H&E) staining, next-generation sequencing (NGS), and quantitative real-time PCR (qRT-PCR) techniques.

Ethics Statement
All experimental procedures involving the use of animals were conducted according to the "Guidelines for Experimental Animals" of the Ministry of Science and Technology (Beijing, China). This study was reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Sichuan Agricultural University, under the permit No. DKY-B20141401.

Animal Manipulation and Sample Collection
Sixty healthy and sexually immature females of Tianfu black rabbit (a Chinese indigenous rabbit breed with superior meat quality, Sichuan, China), aged 60 d and having similar body weights, were used in this study. The hormonal protocol used for induction of follicle maturation and ovulation was depicted in Supplementary Figure 1. Briefly, all does were kept at cages under a controlled light regimen of 14 L: 10 D, had free access to food and water, and were allowed to acclimate for 10 days before gonadotrophic administrations. At the age of 70 d, does were firstly injected intramuscularly with 100 IU of pregnant mare serum gonadotrophin (PMSG; Ningbo Second Hormone Factory, Zhejiang, China), and were then administered with 100 IU of human chorionic gonadotrophin (hCG; Ningbo Second Hormone Factory, Zhejiang, China) through marginal ear vein after 72 h of PMSG treatment. Six does were randomly selected and slaughtered at each time point, including just before PMSG injection (Control), 24 and 72 h after PMSG injection (P24 and P72), and 12 and 48 h after hCG administration (H12 and H48). Blood samples were collected with ethylenediaminetetraacetic acid (EDTA) vacutainers, and plasma was separated after centrifugation at 3000 × g for 15 min and then stored at −20 • C until further analysis. After slaughter, ovaries were promptly removed and washed with 0.9% physiological saline, with the left ovary being used for histology and the right one for RNA extraction. Body weight and reproductive tract weight were recorded at all the time points.

Histological Observation
All left-side ovaries were 4% formaldehyde-fixed for 72 h at room temperature, dehydrated through a graded ethanol series, transferred to xylene, and embedded in paraffin-wax. Paraffin sections of 5 µm thickness from each ovary were stained with H&E and photographed under a Nikon 90i microscope (Nikon, Japan). Ovarian follicles were generally classified into primordial, primary, secondary, tertiary, and mature follicles.

Total RNA Isolation
Total RNA was extracted from the right-side ovaries using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and treated with DNase I (Invitrogen, Carlsbad, CA, USA) following the manufacturers' protocol. The RNA concentration, purity, and integrity number (RIN) were determined using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA), with the RNA concentration >434 ng/µl in a final volume of 35 µl and a range of RIN between 9.0 and 10.0. Thereafter, two RNA samples from each of the Control, P72, and H48 groups were randomly pooled in equal quantities to generate three RNA preparations per group, each of which was then divided into two aliquots processed for construction of either RNA-Seq or Small RNA-Seq library.

RNA Sequencing and Analysis
The mRNA was isolated from total RNA using poly-T oligo attached magnetic beads (Invitrogen, Carlsbad, CA, USA), followed by purification and interruption into short fragments with the fragmentation buffer. Subsequently, 9 cDNA libraries for RNA sequencing were prepared with an Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacture's instruction. Then, all libraries were sequenced on an Illumina Hiseq x-ten platform that generated >6G 150-bp paired-end (PE) raw reads per sample by Beijing Genome Institute (BGI, Shenzhen, China). After sequencing, the raw reads were evaluated with FastQC software and filtered by removing adaptor sequences as well as contamination and low quality reads with >5% unknown bases and a quality score < 20%. Then, the clean reads were mapped to the rabbit reference genome OryCun2.0 (https://www. ncbi.nlm.nih.gov/assembly/GCF_000003625.3) using HISAT2 (v2.0.0), and the Bowtie2 (v2.2.9) software was further used to map the reads to the reference genes for better control of sequencing quality.
The abundance of each gene was expressed as the number of clean reads mapped to its sequence that was normalized to fragments per kilobase of exon per million fragments mapped (FPKM) using RSEM (v1.2.30). Differentially expressed genes (DEGs) in two-group comparisons were screened using DESeq2 (v1.16.1), and the cutoff for significant differential expression was |log2 fold change (FC)| ≥ 1 with a false discovery rate (FDR) < 0.05. DEGs in three-group comparisons were screened using the TCC package (v1.20.0) in R (v3.4.4) software, and a Bayesian method called BaySeq was applied to estimate significant differential expression with a FDR of 0.1 and an iteration of 1 (14). Time-series expression pattern analysis was performed by the Short Time-series Expression Miner (STEM) (v1.3.11) software using the STEM Clustering Method, with "Maximum Number of Model Profiles" and "Maximum Unit Change in Model Profiles between Time Points" set at 50 and 2, respectively (15). Hierarchical clustering and principle component analysis (PCA) were conducted in R (v3.4.4) software. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the selected DEGs were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.8 (http://david.abcc.ncifcrf.gov/), with an EASE score at 0.05.

Small RNA Sequencing and Analysis
Nine small RNA cDNA libraries were prepared with the TruSeq Small RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) following the manufacture's instruction. All library preparations were checked for quality using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and then were sequenced on an Illumina Hiseq 4000 platform that generated >10 M 50-bp single-end (SE) raw reads per sample by Beijing Genome Institute (BGI, Shenzhen, China). After removal of adaptors, junk, and low quality sequences, the clean reads with the length of 18-30 nucleotides were mapped to the rabbit reference genome using Bowtie2 (v2.2.9) combined with the Rfam database to identify different categories of small RNAs.
Identification of known miRNAs and prediction of novel miRNAs were performed with the miRBase v22 database (http:// www.mirbase.org/) and miRDeep2 software, respectively. The expression level of each miRNA was expressed as the tags per million reads (TPM) using DESeq2, and significantly differentially expressed miRNAs (DE miRNAs) in two-group comparisons were identified under the criteria of |log2FC| ≥ 1 and FDR < 0.05. DE miRNAs in three-group comparisons were screened using the TCC package (v1.20.0) in R (v3.4.4) software, and BaySeq was applied to estimate significant differential expression with a FDR of 0.1 and an iteration of 1 (14). The target genes of DE miRNAs were predicted using miRanda and RNAhybrid with respective default parameters, and their overlapping gene list was used for GO and KEGG analyses. Finally, the miRNA-mRNA pairs involving PMSG/hCG-induced follicle development and ovulation as well as specific signaling pathways were identified according to their anti-correlated expression profiles, and their interaction networks were depicted using Cytoscape (v3.7.1) software.

Validation Experiments
For qRT-PCR analysis and validation of our sequencing data, total RNA were isolated from 30 ovarian tissue samples containing the Control, P24, P72, H12, and H48 groups (n = 6/group). 100 ng of total RNA were reversed transcribed using the PrimeScript TM RT reagent Kit with gDNA Eraser (Takara Biotechnology Co., Ltd., Dalian, China) according to the manufacture's instruction. The PCR reactions were performed on the CFX96 TM Real-Time PCR Detection System (Bio-Rad, USA) using the SYBR Premix Ex Taq TM II (Takara Biotechnology Co., Ltd., Dalian, China). Reactions were conducted with the following conditions: pre-denaturation at 95 • C for 5 min, followed by 40 cycles of denaturation at 95 • C for 15 s and annealing/extension at corresponding temperature of each primer set for 30 s. The no-template controls and negative controls without reverse transcriptase were also included in all qPCR runs. Target specificity for each primer set was validated by melting curve analyses, and the identity of all amplicons was verified by sequencing. All samples were amplified in triplicate and relative expression levels of target genes were normalized to the reference gene GAPDH using the comparative Cq method ( Cq) (16). First-strand cDNA of miRNA was synthesized from total RNA using the Mir-X TM miRNA First-Strand Synthesis Kit (Takara Biotechnology Co., Ltd., Dalian, China) according to the manufacture's instruction, and quantification of miRNAs of interest was also performed on the CFX96 TM Real-Time PCR Detection System (Bio-Rad, USA) using the Mir-X TM miRNA qRT-PCR SYBR R Kit (Takara Biotechnology Co., Ltd., Dalian, China). The PCR reactions were conducted under a similar condition to target gene expression detection. Relative expression levels of target miRNAs were normalized to the reference gene U6 snRNA using the Cq method (16). The qRT-PCR primers of selected genes and miRNAs are listed in Supplementary Table 1.

Statistical Analyses
All data were expressed as mean ± SEM. Statistical comparisons between groups were analyzed by analysis of variance (ANOVA) followed by Tukey's test using SAS 9.4 (SAS Institute, Cary, USA). A p < 0.05 was considered statistically significant.

Histomorphological Responses to Gonadotrophic Stimulation in Prepubertal Female Rabbits
Significant morphological and histological changes were seen in the reproductive tracts of does after hormone administrations, as shown in Table 1 and Figure 1. In detail, gonadotrophic stimulation had no significant effect on body weight, fallopian tube weight, and the proportion of fallopian tube weight to body weight of prepubertal does (P > 0.05), but it is notable that both fallopian tube weight and its proportion to body weight increased gradually from Control to P72 and reached the maxima in H12. An around 3-to 5-fold increase in both ovary weight and its proportion to body weight occurred after 72 h of PMSG treatment, and their highest levels were observed at 12 h following hCG treatment (P < 0.05). For uterus weight and its proportion to body weight, both increased continuously from Control to H48, reaching the maxima in H48 (P < 0.05). H&E staining of ovarian sections revealed that both the composition and number of follicles of different size categories differed before and after hormone injection. Compared to the Control group whose ovaries were only composed of primordial, primary, secondary, and tertiary follicles, mature follicles appeared at 72 h after PMSG injection until 48 h following hCG treatment, ovulating follicles were seen in the H12 group, and corpora lutea were present in both the H12 and H48 groups. Additionally, injection of PMSG increased the number of tertiary follicles in P24 and that of mature follicles in P72, while hCG treatment induced follicular rupture in H12 and increased the number of corpora lutea in H48.

Plasma Hormonal Responses to Gonadotrophic Stimulation in Prepubertal Female Rabbits
Changes in circulating levels of MT, PGE 2 , E 2 , and P 4 were determined before and after gonadotrophic stimulation. As shown in Figure 2, plasma concentrations of MT increased significantly (P < 0.05) from Control (248.07 ± 44.83 pg/ml) to H12 (644.89 ± 107.44 pg/ml), but then decreased to 526.90 ± 64.48 pg/ml in H48. A continuous increase from Control to H48 was seen in plasma levels of PGE 2 , reaching the highest level (110.11 ± 6.61 pg/ml) in H48 (P < 0.05). However, levels of E 2 decreased gradually (P > 0.05) from Control (26.17 ± 7.59 pg/ml) to P72 (19.55 ± 2.97 pg/ml), but then increased dramatically to 41.14 ± 3.08 pg/ml in H12 and peaked in H48 (50.87 ± 9.11 pg/ml; P < 0.05). Similar to PGE 2 , levels of P 4 increased continuously from Control to H48 and reached the maxima (7.70 ± 0.63 pg/ml) in H48; and compared to the groups before P72, hCG administration significantly (P < 0.05) increased its levels in both H12 and H48.

Mapping and Annotation of RNA Sequencing Data
In order to reveal the mRNA transcriptomic response to gonadotrophic stimulation, nine cDNA libraries representing rabbit ovaries from the control group (just before PMSG:

Dynamics of Rabbit Ovarian mRNA Transcriptome in Response to Gonadotrophic Stimulation
As shown in Supplementary Figure 2, hierarchical clustering and PCA analyses showed that nine cDNA libraries were sorted into three distinct clusters (i.e., C, P, and H) corresponding to respective treatment and that relatively higher similarity in the expression profiles of all expressed genes was found among three cDNA libraries from the same group, which was indicative of good reproducibility of high-throughput sequencing. After standardization of raw data using the Trimmed Mean of Mvalues (TMM) method, 1122 DEGs were obtained among the C, P, and H groups using TCC in R. Hierarchical clustering analysis suggested that gonadotrophic treatment had a significant effect on the expression profiling of these DEGs and that two distinct profiles were present before and after PMSG stimulation ( Figure 3A). The STEM results showed that these 1122 DEGs were categorized into 16 profiles according to their temporal expression patterns during follicular maturation and ovulation. Among the top seven profiles ranked by decreasing number of genes, profiles 11, 15, and 12 represented a total of 250 genes, whose expression were significantly unregulated through time; while profiles 0 and 4 contained 103 genes, presenting a significant time-dependent decrease in the mRNA levels ( Figure 3B). Three hundred fifty-three genes up-or down-regulated through time were annotated to three main GO categories, including biological process, molecular function, and cellular component. For biological process, most genes were enriched in regulation of gene expression, cell proliferation, cell adhesion, apoptotic process, and inflammatory response as well as angiogenesis and morphogenesis. For molecular function, most enriched terms included extracellular exosome and extracellular space. For cellular component, most enriched terms included calcium ion binding and endopeptidase activity ( Figure 3C). KEGG enrichment analysis revealed that most of these 353 genes were mainly enriched in the PI3K-Akt signaling, focal adhesion, regulation of actin cytoskeleton, ECM-receptor interaction, cytokine-cytokine receptor interaction, phagosome, TNF signaling, and ovarian steroidogenesis pathways (Figure 3D).
In addition to three-group comparison, we also performed pairwise comparisons among the C, P, and H groups using DESeq2 in R. The results showed that a total of 531 DEGs were identified between C and H, 279 up-and 103 down-regulated genes were identified in C vs. P, and 11 and 5 genes were upand down-regulated in P vs. H, respectively ( Figure 4A). In particular, the top 10 up-and down-regulated DEGs upon PMSG stimulation were shown in Supplementary Table 3 A Venn diagram showed that no common DEGs were present in all three pairwise comparisons, whereas, FPKM of 3 genes (i.e., PPP1R3C, MMP13, and SLC2A12) were significantly altered by both PMSG and hCG treatments (Figures 4B,C). KEGG analysis suggested that the representative pathways enriched by these DEGs included the PI3K-Akt signaling, protein digestion and absorption, complement and coagulation cascade, ECM-receptor interaction, focal adhesion, TNF signaling, mineral absorption, and ovarian steroidogenesis ( Figure 4D). Particularly, there were 8 DEGs enriched in ovarian steroidogenesis, including STAR, 3βHSD, CYP19A1, adenylate cyclase 8 (ADCY8), cytochrome P450 1B1 (CYP1B1), bone morphogenetic protein 15 (BMP15), BMP6, and prostaglandin-endoperoxide synthase 2 (PTGS2).

Mapping and Annotation of Small RNA Sequencing Data
To further identify global changes in the miRNA transcriptome before and after gonadotrophic stimulation, nine small RNA libraries were constructed using total RNAs extracted from the FIGURE 2 | Changes in plasma concentration of reproduction-related hormones in prepubertal female rabbits upon gonadotrophic stimulation. Values were expressed as the mean ± SEM of six rabbits per group. Different lowercase letters indicated significant differences among groups at P < 0.05. same rabbit ovaries for mRNA-Seq and sequenced. As shown in Supplementary Table 5, the Q20 ratio of each library was more than 99.63%, and the GC content ranged from 44.43 to 45.27%. The length distribution of clean reads was similar among all libraries, presenting that most reads had 21-23 nt in length, with the 22 nt RNAs being the most abundant. It was observed that 89.54-91.11% clean reads from each library were mapped to the reference rabbit genome. An average number of 591, 592, and 595 known miRNAs together with 325, 326, and 334 novel miRNAs were identified in the C, P, and H group, respectively.

Dynamics of Rabbit Ovarian miRNA Transcriptome in Response to Gonadotrophic Stimulation
Similar to the mRNA transcriptome analysis, 12 DE miRNAs were identified among the C, P, and H groups using TCC in R. Hierarchical clustering analysis verified that expression of these DE miRNAs significantly changed before and after gonadotrophic treatment (Figure 5A). KEGG enrichment analysis of the predicted target genes of these DE miRNAs showed that their targets were mainly enriched in the PI3K-Akt signaling, HIF-1 signaling, prostate cancer, hepatitis B, and alpha-linolenic acid metabolism (Figure 5B).
In addition to three-group comparison, we also performed pairwise comparisons among C, P, and H using DESeq2 in R. The results showed that a total of 60 DE miRNAs were identified between C and H, 36 up-and 20 down-regulated miRNAs were found in C vs. P, and 33 and 16 miRNAs were up-and down-regulated in P vs. H, respectively (Figure 5C). In particular, the top 10 up-and down-regulated DE known miRNAs following either PMSG or hCG administration were listed in Supplementary Tables 6, 7. Of note, miR-205-1 and miR-34c were not only in the list of 12 DE miRNAs screened by three-group comparison, but they were also identified to be significantly differentially expressed in any of three pairwise comparisons (Figures 5A,D). Furthermore, TPM of these two miRNAs decreased following 72 h of PMSG treatment, but later on, increased by hCG (Figure 5E).

Correlation Between qRT-PCR and Illumina Sequencing
From DE mRNAs and DE miRNAs screened by RNA-Seq and small RNA-Seq, 9 genes (MMP13, TIMP1, SLC2A12, PPP1R3C, RASD1, SPP1, ADIPOQ, ACSL6, and STAR) and 6 miRNAs (miR-22-5p, miR-542-5p-1, miR-7b, miR-129a-3p, miR-30c-1-3p-1, and miR-34a-5p) were randomly selected for qRT-PCR validation (Supplementary Figure 5). Moreover, the main miRNA-mRNA pairs involved in the PI3K-Akt, ECM-receptor interaction, and focal adhesion pathways were also validated by qRT-PCR (Figure 7). The quantitative results showed that despite differences in the magnitude of fold-changes, expression of almost all these selected mRNAs and miRNAs determined by qRT-PCR displayed changes in the same direction with that observed by the sequencing data, indicating the true reliability of our RNA-Seq and small RNA-Seq methods. Additionally, expression of these selected mRNAs and miRNAs at other sampling points (i.e., P24 and H12) were also detected by qRT-PCR, providing a better understanding of their roles during follicle maturation and CL formation.

DISCUSSION
Ovarian follicle maturation and ovulation are achieved by gonadotrophic stimulation (2). However, the mechanisms of gonadotrophins regulating these two events have not been fully deciphered. Herein, the effects of combined PMSG/hCG administrations on plasma hormone levels as well as ovarian histomorphology and mRNA and miRNA transcriptomes were systematically evaluated in sexually immature does. Similar to the observations that in primiparous and multiparous does a combination of PMSG and hCG used for either estrus synchronization or superovulation markedly stimulated follicular growth and ovulation (17,18), we demonstrated that PMSG was effective in stimulating the reproductive tract development, as manifested by an increase in its weight and the number of mature follicles, and moreover, ovulation of PMSGtreated does was triggered with hCG. At the plasma level, E 2 reached the minima at 72 h after PMSG treatment, while hCG significantly increased MT and P 4 levels of PMSG-treated does. Furthermore, PMSG/hCG combinations significantly increased PGE 2 levels. MT is mainly secreted by the pineal gland and regulates follicle development and ovulation by affecting both adenohypophyseal hormone release and ovarian steroid production (19). PGE 2 , E 2 , and P 4 are primarily synthesized by ovarian follicles (20,21), and among them, PGE 2 is recognized as a key mediator in gonadotrophin-stimulated ovulation and fertilization (22) while E 2 and P 4 are indispensable for follicle and CL development (2). Thus, it was postulated that follicle maturation is mainly associated with altering levels of E 2 and PGE 2 while rising MT, PGE 2 , and P 4 levels induce ovulation and CL formation.
To reveal the underlying molecular mechanisms, we determined genome-wide profiling of ovarian mRNA and miRNA expression. Results from three-group-and pairwisecomparisons identified 1122 and 672 DEGs throughout follicle development and ovulation, respectively, and 353 of 1,122 DEGs were further screened by STEM analysis. Since GO classification and KEGG enrichment analysis are widely used in the annotation of gene functions (23,24), these 353 DEGs were further annotated to be mainly enriched into the BP terms associated with regulation of gene expression, cellular functions, inflammatory process, and face morphogenesis. Meanwhile, both these 353 DEGs and those 672 DEGs were enriched into the PI3K-Akt signaling, focal adhesion, protein digestion and absorption, complement and coagulation cascades, ECM-receptor interaction, TNF signaling, mineral absorption, and ovarian steroidogenesis pathways, most of which are known to be associated with regulation of cellular functions and tissue remodeling (25)(26)(27)(28). For instance, those 8 DEGs enriched in the "ovarian steroidogenesis" pathway have been demonstrated to regulate ovarian steroid production in rabbits or other mammals (29)(30)(31)(32). Similarly, three DEGs identified during hCG-induced luteinization, including HPGD, PTGFR, and AKR1C5, were reported to be important for CL formation by regulating steroid hormone and prostaglandin metabolism (3,22,29,30). Hence, we proposed that differential expression of these DEGs could be responsible for gonadotrophin-stimulated plasma hormonal fluctuations, thereby promoting rabbit ovarian follicle development and ovulation. As for those three genes regulated by both PMSG and hCG treatments, MMP13 is crucial for ECM remodeling and follicular-luteal transition (33), overexpression of PPP1R3C increased hepatic gluconeogenesis and glycogen storage by targeting protein phosphatase 1 (34), and SLC2A12 was reported to modulate glucose utilization in chicken skeletal muscles (35), supporting the roles of ECM remodeling and glucose metabolism during follicle maturation and ovulation. However, the physiological significance of PPP1R3C and SLC2A12 in the mammalian ovary awaits further investigations. These data suggest that gonadotrophin-stimulated follicle development and ovulation are finely tuned by a wide array of genes related to steroidogenesis, morphogenesis, and functional differentiation of ovarian cells.
Numerous miRNAs have been characterized to regulate follicle maturation and ovulation in a range of domestic animals (12,36,37). However, to the best of our knowledge, this study represents the first to depict global miRNAome changes in rabbit ovaries responding to gonadotrophic stimuli. The predicted target genes of the 12 DE miRNAs identified by three-group comparison were mainly enriched into the PI3K-Akt signaling, prostate cancer, HIF-1 signaling, hepatitis B, and alpha-linolenic acid metabolism pathways, which were evidenced to activate the transcription of genes related to angiogenesis, inflammation, FIGURE 7 | qRT-PCR validation of expression of the main mRNAs (A) and miRNAs (B) involving the PI3K-Akt, ECM-receptor interaction, and focal adhesion pathways in rabbit ovaries before and after gonadotrophic stimulation. RNA-seq was performed only in the Control, P72, and H48 groups, and the values were expressed as the mean ± SEM of 3 pooled ovaries per group. In contrast, qRT-PCR was performed in all five groups, and the values were expressed as the mean ± SEM of 6 individual ovaries per group. The RNA-seq and qRT-PCR results were depicted as the bar and line charts, respectively. energy metabolism, and cellular functions (25,38,39), indicating that these miRNAs play important roles in mediating the actions of gonadotrophins in the rabbit ovary. Among them, both miR-205-1 and miR-34c levels were reduced by PMSG but were enhanced by hCG. Since the miR-205 and miR-34 families were recognized as key regulators of cell survival, proliferation, differentiation, apoptosis, and cell cycle arrest in ovarian cancer cells (40,41) and more than 10 predicted target genes of both miR-205-1 and miR-34c matched our mRNA-seq data upon PMSG stimulation but only one matched target gene for either one upon hCG stimulation, we speculated that both miR-205-1 and miR-34c may function differently between PMSG-stimulated follicle development and hCG-induced ovulation by targeting different genes. Besides, as the two most downregulated miRNAs by PMSG, miR-34b and miR-205-5p were reported to regulate cell apoptosis and cell-cycle arrest in many types of cancers, and among their common predicted target genes, LGALS9, TNC, and LIF are known as regulators of cell apoptosis, proliferation, and adhesion (42)(43)(44), suggesting that reduced levels of both miRNAs may promote PMSG-stimulated follicle maturation by regulating ovarian cellular proliferation and apoptosis. By contrast, less is known about the physiological functions of the two most unregulated miRNAs by PMSG as well as the two most up-and down-regulated miRNAs by hCG in the ovary.
The key miRNA-mRNA pairs involved in the PI3K-Akt, ECM-receptor interaction, and focal adhesion pathways were identified, because these pathways have essential roles in regulating PMSG-stimulated follicle maturation (26)(27)(28). Furthermore, expression of eight DEGs commonly enriched in these pathways and their predicted regulatory miRNAs were validated by qRT-PCR, showing reduced miRNAs levels in parallel with rising mRNAs levels upon PMSG stimulation. Among them, collagens and tenascins are known as components of ECM that promote and counteract follicular cell adhesion, respectively (45). SPP1 is an adhesion molecule that plays important roles in regulating developmental processes, cell-cell interactions, inflammatory responses, and carcinogenesis (46), and moreover, it promoted the progress of human ovarian cancer through the Integrin β1/FAK/AKT signaling pathway (47) and its expression was regulated by estrogen in both the porcine uterus (48) and the chicken oviduct (49). Integrins are cell surface receptors that mediates interactions between cells and surrounding ECM (45), and they were able to modulate cell cycle progression and cell proliferation via chk1 and Rb/E2F pathways and were also regulated by E 2 in breast cancer cells (50,51). Hence, we postulated that these eight DEGs are important for PMSG-stimulated follicle maturation by not only constituting follicular cell microenvironment but also regulating cellular functions through multiple signaling pathways, and their upregulated mRNA levels may be associated with a reduction in plasma E 2 concentration. Nevertheless, the real roles of their regulatory miRNAs during mammalian follicle development remain largely unknown and need to be further studied. In addition, regarding that only 16 DEGs were identified during hCG-induced ovulation, the miRNAs predicted to regulate all these DEGs that matched DE miRNAs were fully identified. Among them, since PTGFR mediates PGF 2α -induced luteolysis and parturition (22) and MMP13 regulates follicular-luteal transition (33), it was inferred that miR-199a-5p-1 and miR-22-5p could target PTGFR while miR-542-5p, miR-107-1, and miR-30c-1-3p could target MMP13 to regulate ovulation and CL formation.
We also recognized the limitations of using the rabbit as our experimental model because rabbits have distinct reproductive characteristics from humans. Specifically, compared to humans, rabbits are induced ovulators, do not have a pronounced reproductive cycle, and normally ovulate 8-12 oocytes (52). Nevertheless, rabbits share a similar process of gonadotrophinstimulated follicle development and ovulation with humans and other domestic animals, which, together with the high consistency between the qRT-PCR results and transcriptome sequencing, highlight the importance and accuracy of our data in fully dissecting the regulatory mechanisms of gonadotrophins in the mammalian ovary. These identified DE mRNAs and miRNAs as well as their interaction networks also provide us a new avenue to reveal the molecular mechanisms controlling rabbit ovarian activities, which will be beneficial for improving reproductive efficiency in the rabbit industry.

DATA AVAILABILITY STATEMENT
The transcriptome sequencing data are available in the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) at NCBI, with the BioProject ID: PRJNA552172 and SRA Accession Number: SRR9639659-9639676.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee (IACUC), Sichuan Agricultural University.