Expression profile analysis of sheep ovary after superovulation and estrus synchronisation treatment

Abstract Superovulation is a widely used reproductive technique in livestock production, but the mechanism of sheep's superovulation is not yet clear. Here, a method of superovulation and estrus synchronisation was used to treat female Duolang sheep. After treatment, there were significant differences in serum FSH and LH levels and the number of dominant follicles between the two groups of sheep. We identified a total of 5021 differentially expressed genes (11, 13 and 15 days after treatment) and performed RT‐qPCR analysis to identify several mRNA expression levels. GO and KEGG enrichment analysis revealed that differentially expressed genes were involved in the regulation of signalling pathways of follicular development, cell cycle, material synthesis, energy metabolism, such as COL3A1, RPS8, ACTA2, RPL7 RPS6 and TNFAIP6 may play a key role in regulating the development of follicles. Our results show a comprehensive expression profile after superovulation and estrus synchronisation treatment. We provide the basis for further research on breeding techniques to improve the ovulation rate and birth rate of livestock.


INTRODUCTION
The ovaries are female's animal reproductive organ that contains a large number of follicles at different developmental stages (Eppig, 2001;Fortune, 1994;Matzuk et al., 2002). The development of follicles presents a dynamic process of 'follicular waves' including recruitment, selection, dominance, turnover or ovulation (Ginther et al., 1996). This process is regulated by a variety of hormones, and the imbalance of hormone expression levels can disrupt the follicular development affecting animal reproduction (Wallach & Hodgen, 1982). Follicle stimulating hormone (FSH) is the main stimulator for mammalian follicles development (Bex and Goldman, 1975;Kauffman et al., 1998). The FSH hormone is necessary for the further development of follicles, although it is not required in the early stages of development (Norman et al., This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2022 The Authors. Veterinary Medicine and Science published by John Wiley & Sons Ltd. 1994). FSH mainly stimulates the growth and development of follicles, affects the number of growing follicles and stimulates the final maturation of follicles. In addition, luteinising hormone (LH) is another essential hormone that promotes oocyte maturation and ovulation (Cooke et al., 1996). FSH and LH are necessary for the cascade of gene expression necessary to initiate ovarian function (Richards & Hedin, 1988).
Ruminants, especially sheep, are widely bred as livestock in the world, providing a large number of meat, fat, milk, wool and other products for humans (Pierson & Ginther, 1988). Increasing the litter size of sheep and animal reproductive performance of sheep has always been the focus of breeding research. Many genes related to sheep reproductive capacity including Booroola gene (FecB), bone morphogenetic protein receptor type IB (BMPRIB), growth differentiation factor 9 (GDF9), bone morphogenetic protein 15 (BMP15) and Prostaglandin G/H synthase 2 (PTGS2) have been identified (Bonnet et al., 2011;Galloway et al., 2000;Hanrahan et al., 2004;Piper & Bindon, 1983;Souza et al., 2001).
Studying the mechanism of follicular development and differentiation is necessary to increase the litter size of livestock. Clinically, FSH can be used to treat infertility caused by abnormal gonadotropin secretion. In animal husbandry production, FSH can be used to perform estrus synchronisation and superovulation of livestock, which can greatly improve the litter size of livestock. In recent years, a lot of research has been conducted on sheep breeding , but the mechanism of causing different physiological states of sheep's ovary after superovulation and estrus synchronisation treatment has not been fully understood.
Duolang sheep is an excellent dual-purpose sheep breed for meat and fat. It has the characteristics of large size, rapid growth and development, large meat production, fresh and tender meat, early sexual maturity and high reproducibility. It is an ideal breeding breed for mutton production that is widely raised in Xinjiang, China. We used Illumina HiSeq 2500 to study Duolang sheep ovarian genes at different times after superovulation and estrus synchronisation treatment. Our research provides gene expression profiles of sheep's ovaries at various times after superovulation and estrus synchronisation treatment, which will lay the foundation for further understanding of the mechanism of sheep follicular differentiation and development. From the 10th to the 15th day of sheep treatment, we selected three sheep from each group (Gao et al., 2017;Ren et al., 2017;Ullah et al., 2020) to collect blood into EDTA tubes from jugular vein on an empty stomach and separate serum stored at -20 • C to measure hormone content. At the same time, three ewes from each group were immediately slaughtered in the slaughterhouse, which was performed by severing the carotid artery and jugular vein, and complete ovaries were collected within 5 min. We measured and recorded diameter and num- according to the manufacturer's protocol. Library quality was assessed on an Agilent Bioanalyzer 2100 system after addition of the adaptor (Parkhomchuk et al., 2009). The post-test library was subjected to 125 bp/150 bp paired-end sequencing using the Illumina HiSeq 2500 (Illumina, San Diego, USA) platform.

Data quality control and sequencing analysis
We filtered out the sequencing joints or the lower-quality reads to obtain clean data (clean reads). At the same time, the Q20, Q30 and GC content of clean data were calculated and evaluated (Andrews, 2010).

Differential expression analysis
We used the FeatureCounts (v1.5.0-p3) software to calculate the readings mapped to each gene by the Hisat2 software (Schmid & Grossniklaus, 2014), and then the fragments per kilobase million (FPKM) of each gene were calculated to quantify the expression of non-genes in combination with the length information of the gene (Trapnell et al., 2012). Cuffdiff (v2.1.1) was used to calculate the FPKM of genes in each sample. The FPKM of genes were calculated by summing the FPKM of the genes in each genome. Next, differential expression in gene expression data was calculated using a negative binomial distribution model in the DESeq2 R package (1.16.1) (Love et al., 2014). Benjamini and Hochberg were used to adjust the p value, and genes with p values <0.05 were designated as differentially expressed genes (Thissen et al., 2002).

Real-time quantitative RCR (RT-qPCR) analysis
Total RNAs of sheep ovaries were extracted using TRIzol (Invitrogen, CA, USA) according to the method at the time of library establishment. Required amount of RNAs was taken, and cDNA was synthesized in strict accordance with the protocol of the RT-PCR kit (Takara, Dalian, China) manufacturer. Primers required for the experiments were designed using Primer 5 software (Table S1). We used RT-qPCR (three copies each and three biological replicates per group) to quantify some genes (10 genes are randomly selected in each group, but only 3 differential expressions in the third group) by monitoring the product formation of the fluorescent dye SYBR Green (Takara Biotech, Dalian).
Data were then analysed using the equation 2 ΔΔCt (Livek & Schmittgen, 2001 p < 0.05 is considered to be significantly enriched by differentially expressed genes. At the same time, we performed PPI analysis on the differential genes based on the STRING database.

Ovarian tissue morphology
Comparing the DC group (SOV group) and the DT group (ES group), we found that there were significant differences in the ovary morphology of the sheep between the two groups ( Figure 1). The ovaries of proestrus had several follicles of different sizes (diameter <3 mm;

Number of dominant follicles
To quantify the effects of both treatments on follicular development, we calculated the number of dominant follicles with a follicle diameter greater than 3 mm in both groups of sheep ovary (Gonzalez et al., 2004). We found that there were significant differences between the two groups of ovaries. By comparison, the number of dominant follicles in the DC group was significantly greater than that in the DT group from the 12th day after treatment (average increase 4.66, 15.34, 11.67 and 10 dominant follicles, respectively, p < 0.01) (

FSH and LH hormone levels
We detected the levels of FSH and LH in the serum of sheep and found that the level of FSH in the serum of sheep in the DC group increased, while the overall level of FSH in the DT group remained unchanged.
From the 13th to 15th day after treatment, the serum FSH hormone levels in the DC group were significantly higher than those in the DT group. The serum LH hormone levels in the two groups showed a tendency to increase gradually, but on the 11th-14th days after treatment, the serum LH hormone levels in the DC group increased rapidly, which was significantly higher than that in the DT group (p < 0.05).
Interesting, on the 15th day after treatment, there were no significant difference in serum LH hormone levels between the two groups of sheep (Table 3).

Quality and statistics of the sequencing data
According to the procedure in the pipeline, we analysed the expression patterns of genes in sheep's ovary at different time point (Figure 2a).
We obtained a total of 668,709,449 raw reads and 656,489,500 clean reads (over 196G data) (Table S2). After quality control, we found that the high-quality data (Q30) accounted for 94.46% of the total data, indicating that the quality of data produced by sequencing was good. We found that 81.12% of the RNA-seq data was aligned with the exons, 10.09% was aligned with intron and 8.80% was aligned with intergenic sequences by mapping to the reference genome ( Figure 2b). A total of 15,025 genes were identified by further analysis. Detailed information on all newly identified genes is listed in Table S3.  Table S4.

Verify differential expression genes
To determine the accuracy of the RNA-seq results, we confirmed the expression levels of these genes by further experiments. We randomly selected 10 differentially expressed genes from each sample to verify (except DCC vs DTC, only 3 differentially expressed genes). By gel electrophoresis analysis, the reverse transcription PCR (RT-PCR) products showed a single band of the expected size, which indicated that the primers for these genes are accurate and specific (Figure 4a). We further validated these amplification products using DNA sequencing.
Compared with RNA-seq data, the results of qPCR suggested that all 23 genes amplified showed the expected difference (Figure 4b). There is a strong agreement between the results of these experiments and the RNA-seq data, indicating that the RNA-seq data can truly reflect the true differential expression in sheep.

F I G U R E 3
The information of genes in sheep at different time points after estrus synchronisation and superovulation. (a) Under different experimental conditions, the horizontal box plot is expressed, where the abscissa is the sample (group) name, the ordinate is log2 (FPKM+1), and the box plot of each region is for the five statistics (top to bottom: value, upper quartile, median, lower quartile and minimum, respectively). (b) Sample correlation heat map. (c) Differential gene Venn diagram. The sum of all the numbers in the circle represents the total number of differential genes in the comparison combination, and the overlapping regions represent the differential genes shared between the combinations. Red representative DCA vs. DTA, green representative DCB vs. DTB and blue representative DCC vs. DTC (p < 0.05)

Functional enrichment analysis of differentially expressed genes
To further understand the biological functions of these differentially expressed genes in the ovary between DC and DT group, we used the GO and KEGG pathways to enrich these differentially expressed genes. The results of GO annotation showed 46 (DCA vs DCT group) and 240 (DCB vs DTB group) significantly different genes in biological process, molecular function and cellular component (Table S5) (p < 0.05). The differential genes in the DCA vs DTA group were mostly enriched in molecular function (endopeptidase complex, substrate-specific transporter activity, cation transmembrane transporter activity and ion transmembrane transporter activity), while DCB vs DTB, the most enriched in the group, were the cellular composition (protein complex, peptidase complex and catalytic complex) (Figure 5a and b).
Especially in the DCB vs DTB group, the biological process of sheep ovary was the most active during this period. KEGG enrichment analysis was used to analyse the biological functions involved in these differential genes (Table S6). At different stages after treatment, differentially expressed genes were significantly enriched in signalling pathways involved in follicular development and maturation, such as the FoxO signalling pathway, Apelin signalling pathway, and ovarian F I G U R E 4 Validation of genes and their differential expression in sheep ovary. (a) RT-PCR amplification of genes with specific primes. The order of electrophoresis are Marker (Takara DL500 500 bp, 400 bp, 300 bp, 200 bp, 150 bp, 100 bp and 50 bp), N is negative control, lane 1 is β-actin, lanes 2-22 represent SLC27A2, AMIGO2,ADAM9,LDLR,EIF4EBP1,MGARP,COBLL1,PDE1C,CYP11A1,CAPN6,RRM2,KPNA2,LMNB1,PPIA,TNNI3,BIRC5,GTSE1,AURKB,TNFAIP6,TNFAIP6,KLF2 and APLP2,respectively. (b) Expression of 23 differentially expressed genes by RT-qPCR. Blue is the result of qPCR, and red is the result of RNA-seq (p < 0.05). The red gene name represents the DCA vs. DTA group; the black gene name represents the DCB vs. DTB group and the blue gene name represents the DCB vs. DTB group. The ordinate is the ratio of the relative expression levels of the gene in the treatment group and the estrus synchronisation (DTN/DCN, N stands for a, b and c) steroidogenesis (Table 4). These results suggest that sheep ovary after FSH treatment has a strong biological response ( Figure 5C-E). We found a large number of protein-protein interactions on 11th and 13th days after treatment by PPI analysis of differentially expressed genes in sheep ovaries at different stages of treatment, some of the main genes were shown in Figure 6 and all the interaction networks were shown in Figure S2 (Table S7).

DISCUSSION
In recent years, with the rapid breakthrough and development of highthroughput technology, many studies have revealed important issues in many livestock production processes by RNA-seq (Clark et al., 2017;Jäger et al., 2011;Fan et al., 2013). Increasing the reproductive rate of livestock can bring huge economic benefits (Nienaber & Hahn, 2007), estrus synchronisation and superovulation are widely used in livestock production to increase the reproductive capacity of livestock (Mendoza-Tesarik & Tesarik, 2018). In this study, the follicle development rate and the number of ovulations between the two groups were significantly different. This is consistent with the known function of FSH, but it was possible to promote follicular development and maturation, and we used RNA-seq to further study its molecular level regulation.
After FSH treatment, we found that the number of differentially expressed genes first increased and then decreased. There were 546 differentially expressed genes between the two groups immediately after FSH treatment (DCA vs. DTA). However, KEGG pathway F I G U R E 5 Enrichment analysis of differentially expressed genes sheep at different time periods after estrus synchronisation and superovulation. GO enrichment analysis of differentially expressed genes at DCA vs. DTA (a) and DCB vs. DTB (b) groups (p < 0.05). The abscissa is the GO macromolecules and GO term at the next level, and the ordinate is the number and proportion of genes annotated to the term. KEGG enrichment analysis was performed on the hosting genes of differentially expressed genes in DCA vs. DTA (c), DCB vs. DTB (d) and DCC vs. DTC (e) groups. The vertical axis and the horizontal axis represent the name of the access and the access factor, respectively. The size of the dots represents the number of genes that are enriched in the access, and the colours correspond to different q-values enrichment analysis showed that only Cardiac muscle contraction, Carbon metabolism and Biosynthesis of amino acids were significantly enriched. Although some differential expressed genes were enriched in PPAR signalling pathway, synthesis and degradation of ketone bodies, ovarian steroidogenesis and cell cycle, these differences had not reached a significant degree. Biosynthesis of amino acids is an important part of protein synthesis. These results indicated that the protein biosynthesis of the ovary is enhanced after FSH treatment, which is consistent with the previous report that FSH can accelerate protein biosynthesis within 1 h (Means & Hall, 1967). FSH may promote ovar-ian protein biosynthesis and may be an important way to regulate follicular development.
The number of differentially expressed genes in the two groups of ovaries was then increased after receiving more FSH treatment (DCB vs. DTB). We found that these genes involved in follicular development and premature activation were significantly upregulated, including COL3A1, RPS8 and ACTA2, RPL7 and RPS6, indicating high concentrations of FSH promote the development of follicles by stimulating their early entry into the developmental stage, which is consistent with previous reports (Diao et al., 2004;Mork et al., 2012;Reddy et al., TA B L E 4 Differentially expressed genes are significantly enriched in pathways primarily related to follicular development and maturation  (Assidi et al., 2008;Manna et al., 2016).
Previous studies have shown that CYP11A1 will lead to an increase in aromatase expression and serum estradiol levels in female mice, which will lead to an increase in the number of mature follicles. At the same time, the expression level of CYP11A1 in the follicles of multiple goats is higher than that of uniparous goats (Bakhshalizadeh et al., 2018;Zou et al., 2020). This may be the reason why the expression level of CYP11A1 in sheep ovary increased significantly after treatment with FSH. In addition, studies have reported that GDF9 is hardly expressed when oocytes develop into primordial follicles. And the oocytes expressing GDF9 in the later development are obviously larger. GDF9 sends signals through SMAD 2 and 3 to regulate oocyte development (Coutts et al., 2008). This is consistent with our results that the expression levels of GDF9 and SMAD3 were not different between the two groups immediately (DCA vs DTA) after FSH treatment, but they increased significantly with the increase of FSH treatment (DCB vs DTB). At the same time, these differentially expressed genes extensively enrich in signalling pathways related to metabolism, such as energy metabolism, cell cycle, cell proliferation and biomass synthesis, such as phosphorylation, carbon metabolism, RNA transport, cell cycle and DNA replication. These results indicated that FSH regulates the development of follicles through a wide range of regulatory pathways.
But, there was a significant difference in the expression of only three genes (APLP2, TNFAIP6 and KLF2) between the two groups as ovulation occurred (DCC vs. DTC). This indicated the expression level of differentially expressed genes caused by FSH treatment decreased. The occurrence of ovulation after removal of FSH treatment led to the rapid recovery of differentially expressed genes between the two groups, indicating that FSH plays a key role in superovulation. This also verifies that the above-mentioned difference between DTB and DTC is mainly caused by FSH treatment, which further indicates that FSH treatment can promote the maturation of follicles by regulating cell proliferation, cell cycle and energy metabolism.
Interestingly, we found that TNFAIP6 was the only gene that was upregulated at all three treatment time points. Previous reports have shown that the expression of TNFAIP6 in follicles of multiple goats during estrus is eight times higher than uniparous goats (Zou et al., 2020), and the lack of TNFAIP6 will cause mice to be sterile (Yu et al., 2010).
These studies have shown that TNFAIP6 is deeply related to the ovulation number and fertility of female animals. LH treatment inhibits the expression of FSHr and Lhcgr and weakens the effect of FSH; even in the early stage of follicular development, it may cause follicular atresia and cysts (Farhat et al., 2011). We believe that TNFAIP6 plays a key role in superovulation, although this requires further experimental verification. Simultaneously, CYP11A1 and TNFAIP6 are the core of the interaction network and interact with many important genes that regulate follicle development. Once again, they show the key role in these candidate genes. The specific mechanism by which these genes regulate superovulation in sheep will be the focus of our next research.
In conclusion, our ovarian transcriptome studies of the ovary at different time points after superovulation and estrus synchronisation treatment indicate that TNFAIP6 and CYP11A1 may play an important role in promoting development of sheep follicle. At the same time, LDLR, SMAD3, GDF9, ACTA2, RPS6 and CDK2 also play an active role in this physiological activity. According to GO and KEGG analysis, FSH can also promote the development of follicles through a combination of signalling pathways that regulate cell proliferation, cell cycle and energy metabolism. Our research provides valuable resources for further understanding of the mechanism of superovulation during sheep follicular maturation.