RNA-seq Transcriptome Analysis in Ovarian Tissue of Pelibuey Breed to Explore the Regulation of Prolificacy

The Pelibuey sheep (Ovis aries) is an indigenous breed distributed in the tropical regions of Mexico. The prolificacy of this sheep is on average from 1 to 1.5 lambs, being an important breeding characteristic that owners seek to increase with the purpose of economic improvements. New-generation RNA sequencing technology has been used to identify the genes that are expressed in the ovarian tissue of sheep that have two or more lambs per parturition, as well as to elucidate the metabolic pathways that are affected by the expression of these genes, with the purpose of better understanding the prolificacy in the sheep. In the present study, the transcriptional expression of multiparous and uniparous sheep was compared using RNA sequencing. Multiparous (M group) and uniparous (U group) sheep that had a genealogical record for three generations (M, n = 5 and U, n = 5) were selected. RNA was extracted from ovarian tissue and subsequently used to prepare the libraries that were sequenced using the Illumina NextSeq500 platform. A total of 31,575 genes were detected from the transcriptomic analysis of which 4908 were significantly expressed (p-value ≤ 0.001) in the ovary of sheep. Subsequently, a second filter was carried out to evaluate the false discovery rate (FDR) and select those genes with p-values ≤ 0.05 and values of expression ≥ 1 (log2), obtaining 354 differential expressed genes (DEG): 120 genes up-regulated and 234 genes down-regulated in the group M with respect to the group U. Through Gene Ontology (GO) and metabolic analysis, we obtained information on the function of differentially expressed genes, and its importance in the reproduction of multiparous sheep. This result suggest that genes identified in the present study participate in the development of the final stages of follicles.


Introduction
Reproductive efficiency is considered one of the most important factors affecting the sheep meat industry. It is directly related to traits such as ovulation rate and litter size, which are genetically affected by causative mutations in some minor and major genes [1] and their interactions with the environment. The litter size improvement is one of the most important objectives for sheep breeders as it greatly increases production efficiency and economic gains [2], and therefore its analysis with

Materials and Methods
Two groups of Pelibuey (n = 10) ewes were selected based on the criteria of number of lambs per partum, and with records of three consecutive births. All sheep in the study were from a commercial sheep farm "El Cortijo", located in Campeche, Mexico (19 • 43 48.92" N and 90 • 05 11.58" W). The animals were transported and kept in the facilities of CESyRO (Center of Selection and Ovine Reproduction) of Technological Institute of Conkal, Yucatan, México. The conditioning of the animals was supervised by the ethics committee of the institute following the official Mexican standard (NOM-033-SAG/ZOO-2014).

Synchronize Estrus and Sample Preparation
The 10 experimental sheep were divided into two groups of five animals, referred to as Group M and Group U. Group M corresponds to multiparous sheep with two lambs per birth (M, n = 5), and group U corresponds to uniparous sheep with a single lamb at birth (U, n = 5). The animals of both groups were fed the same diet of Taiwan grass (Pennisetum purpureum) and 250 g/ewe/d of a commercial supplement with 16% crude protein, mineral lick and water ad libitum. The estrous cycle of the sheep was synchronized using intravaginal sponges of fluorogestone acetate (40 mg; Chonogest ® , Intervet, Santiago Tianguistenco, México) for 14 days, followed by intramuscular application of 1000 UI of gonadotropina (Folligon ® ). Ovarian tissue samples were obtained by biopsy using a left unilateral ovariectomy, immediately stored at −80 • C until RNA extraction.

RNA Extraction
For RNA extraction, 50 mg of tissue was taken from each sample using the Trizol Reagent according to the manufacturer's protocols (Invitrogen™, Carlsbad, CA, USA). The RNA was treated with Turbo DNA-free (Ambion, Austin, TX, USA). The quantity and quality of RNA were measured using the Agilent 2100 BioAnalyzer system (Agilent Technologies, Santa Clara, CA, USA) with the RNA 6000 Nanochip. RNA Integrity values (RIN) were >7.0 for all samples used in the analysis.

Library Preparation and Sequencing
Samples were processed following the manufacturer's protocol for NEBNext Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs (NEB), Ipswich, MA, USA) in conjunction with the NEBNext Poly (A) mRNA Magnetic Isolation Module and the NEBNext Multiplex Oligos for Illumina (Index Primers Set 1) NEB. This step was followed by first-strand cDNA synthesis using reverse transcriptase and oligo-dT primers. Synthesis of ds-cDNA was done using the 2nd strand master mix provided in the kit, followed by end-repair and dA-tailing. This was done to add universal Illumina adaptors to the samples. Finally, the library was enriched (and barcoded) by 6-9 cycles of amplification, and purefied by Agencourt AMPure beads (Beckman Coulter, Brea, CA, USA). The library size and mass were assessed by analysis in the Agilent 2100 BioAnalyzer (Tape Station). Briefly, 10 µL of library (4 nM) was mixed with 10 µL 0.1 N NaOH for 5 mins, then the library was diluted to 20 pM in HT1 buffer and sequenced using the reagents provided in the Illumina 150-cycle kit and following the manufacturer's protocol for high output NextSeq500 sequencing.

Data Analysis and Differential Expression Analysis
The raw data obtained from the sequencing by RNA-seq were first filtered to remove adaptors, as well as low-quality reads, using the NGS QC Toolkit v2.3.3 [15] software, and its program IlluQC.pl for Ilumina data using default parameters. Subsequently, to carry out the transcriptomic study of the filtered sequences, an analysis of the RNA reads was carried out with the strategy of using a reference genome. The TopHat v2.1.1 program was used [16] to assemble the reads, and the genome of the sheep O. aries v4.1 was used as reference genome to map the filtered reads and downloaded from NCBI genomes database (https://www.ncbi.nlm.nih.gov/genome?term=Ovis%20aries). The expression level of the genes in each sample were estimated and normalized as RPKM (Reads Per Kilobase Million Mapped Reads) [17], which is a normalized measure of read density that allows transcript levels to be compared both within and between samples. Cuffdiff v2.2.1 software [16], was used to identify the DEGs between the two groups, and a p-value was assigned to each gene to evaluate its statistical significance. Then we determine the false discovery rate (FDR) of the test to account for Type I errors. Multiple-testing corrections were performed using the Benjamini and Hochberg step-up false-discovery rate (FDR)-controlling procedure to calculated adjusted p-values. Genes with an adjusted p-value ≤ 0.05 and an expression |log2 ratio| ≥ 1 were identified as DEGs.

Analysis of Gene Ontology Category and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway
To obtain information on the function of differentially expressed genes, and its importance in the reproduction of multiparous sheep, an enrichment analysis was performed for up-regulated and down-regulated genes, using the online platform DAVID [18]. The list of selected DEGs obtained from the transcriptomic analysis made with TopHa was sent to DAVID database to obtained their biological functions using Gene Ontology (GO) terms [19], which analyzes three functional categories: biological process (BP), cellular component (CC), and molecular function (MF). To obtain an overview of the DEGs and the different metabolic pathways in which they participate, the KEGG database was used [20]. To graph the GO enrichment analysis obtained from the DAVID database, the WEGO v2.0 program was used [21].

Validation of RNA-Seq Data by Quantitative Real-Time Polymerase Chain Reaction
The same tissues used for RNA-seq were used for quantitative real-time polymerase chain reaction (qRT-PCR) analysis. The RNA was treated with the DNase I enzyme (TURBO DNA-free™, Ambion ® ). For cDNA synthesis, DNase-treated RNA was used with the RevertAid First Strand ® kit (Thermo Scientific ® , Lithuania). The real-time PCR were performed in a Thermal Cycler (iCycler IQ5; Bio-Rad, Hercules, CA, USA) with specific primers designed (Table 1). Each analysis was performed in triplicate with the Maximum SYBR ® Green/Fluorescein qPCR kit (Thermo Scientific ® ) according to the manufacturer's instructions. The reaction of each sample resulted in a final volume of 25 µL, constituting of 12.5 µL of SYBR Green/Fluorescein (2X), 0.3 µM of each oligonucleotide and 2000 ng of RNA. The thermal protocol was 95 • C for 30 s and 35 cycles of 64 • C for 1 min and 72 • C for 30 s. The data were analyzed with iQ5 optical system software v.2 (Bio-Rad). Finally, the method of 2 −∆∆CT by Livak et al. [22], was used to calculate fold changes. The significance of the expression levels between the control group (β-actin) and the treatments (CA5A and FEM1B) were determined for analysis using the Mann-Whitney U-test and Kruskal-Wallis test. Statistical analyses were performed using Sigma Plot (SYSTAT, Chicago, IL, USA) and figures were plotted using Prism 5.0 (GraphPad Software Inc., San Diego, CA, USA). Table 1. Primer sequences used to amplify target genes using by quantitative real-time polymerase chain reaction (RT-qPCR). β-actin for the housekeeping as an internal control.

Name
Size (pb) Primer Sequence Tm (C • ) Efficiency

Results
Using the Illumina NextSeq500 sequencing platform, a total of 64,998,983 paired-end (PE) reads were obtained for group M, while 61,941,256 PE reads were obtained for group U. The quality control resulted in 60,528,631 (93.11%) reads for group M and 57,442,547 (92.73%) reads for group U (Table 2) and their results indicate that they were of high quality ( Figure 1).

Identification of Differentially Expressed Genes
The transcriptomic analysis of the two sheep groups showed a total of 31,575 genes detected using the reference genome of O. aries v4.1. After selecting those genes with a p-value ≤ 0.001, were obtained 4908 genes for group M, of these 2706 genes were up-regulated, and 2202 genes were down-regulated. In a subsequent step, a Benjamini and Hochberg test to control false discovery rate (FDR) was done. The genes with p-value ≤ 0.05 and fold change ratio (log2) ≥ |1| were identified as DEGs, obtaining 120 genes up-regulated and 234 genes down-regulated in the group M (Table S1). The top 10 differentially expressed genes involved in the ovary development process in ewes with two lambs at birth were: FGF18, CA5A, SVIP, CXCL14 and FST up-regulated; while the down-regulated were RENT, FEM1B, MMP9, CYP2S1 and HOXA9 (Table 3). The set of most relevant genes identified in the present analysis and that has been reported in other studies were IL1A, CYP2S1, BMP15, and INHIB [13,22], involved in reproductive processes as regulatory factors in ovarian development and in the ovulation process in female mammals [23,24].

Gene Ontology and KEGG Analysis
The 354 selected genes were categorized into 48 functional groups of the GO, distributed as follows: 232 genes were categorized in 24 functional groups for biological process (BP); 267 genes in 16 groups for cellular component (CC) category; and 191 genes in 8 groups for molecular function (MF) category. The highest abundance of genes was represented in the CC category ( Figure 2). The enrichment of DEGs in GO terms was tested to gain insights about the biological implications. The GO terms GO: 0009987 (Cellular process) for BP, GO: 0005623 (Cell) for CC, and GO: 0005488 (Binding) for MF, showed high abundance of DEGs into these categories (Figure 2). Regarding to functions associated with reproduction processes, we find that the GO: 0022414 (Reproductive process) and GO: 0000003 (Reproduction) are present in the BP category. The metabolic analysis of the DEGs using KEGG database, classified them into 116 metabolic pathways (Table S2). The top 10 metabolic pathways were identified (Table 4) and only six are related to reproductive processes: calcium-signaling pathway, TGF-β-signal pathway, insulin-signaling pathway, estrogen-signaling pathway, PI3K-Akt-signaling pathway, MAPK-signaling pathway (Table 4). Among these KEGG pathways, were the TGF-β-signaling pathway, estrogen-signaling pathway and PI3K-Akt-signaling pathway. These metabolic pathways are associated with gonadal development, ovarian steroidogenesis, oocyte maturation and steroid hormone biosynthesis [25]. These pathways may be playing an important role in the regulation of follicle development in sheep.

Validation of RNA-seq Data by qRT-PCR
The CA5A up-regulated and FEM1B down-regulated genes (Table 3), involved in the nitrogen pathway and the regulation of ubiquitin-protein transferase activity, respectively, were selected to confirm their expression profiles obtained by RNA-seq analysis using qRT-PCR analysis. The β-actin gene was used as housekeeping control. The experimental evaluation of the CA5A gene using qRT-PCR showed a positive expression, agreed with the results obtained from the RNA-seq analyzes. While in the case of the relative expression of the gene FEM1B it was not a negative expression, as was predicted in the RNA-seq analyzes, observing an increase in its expression in the results of qRT-PCR analysis, but its expression being less than in the tissue of the uniparous sheep ( Figure 3).

Discussion
The ovulation rate is determined by the regulation of intraovarian factors and their increase is associated with the expression of genes involved in the development of antral follicles in the preovulatory stage, resulting in an increase in prolificacy traits [8]. The transcriptomic analysis performed with the RNA-seq technique in this study yielded a number of reads similar to other studies (53,024,205) [14] and (86,966,348) [26]. Similarly, with our transcriptomic analysis in the Pelibuey sheep, we have identified genes that are associated with multiple calving sheep and have been previously reported in other sheep breeds such as the IL1A (Interleukin 1 α) and CYP2S1 (Metabolism of xenobiotics by Cytochrome P450) genes in Hu ewes, the function of these genes in the process of folliculogenesis is unknown [14]; the BMP15 and INHBA (Inhibin β A subunit) genes in Finnsheep ewes [26]. The BMP15 gene has activity in the oocytes (both immature and mature) within the ovary, stimulating the proliferation of GC in a follicle-stimulating hormone receptor (FSHR)-independent manner [2] and INHBA gene shows a considered level of mRNA expression in the GC of antral follicles during the stages of the estrous cycle in sheep [27]. While the gene IL1A has activity in the mechanism of apoptosis, which is responsible for the elimination of oocytes and a deficiency or blockade of increased expression of FSHR in GC in rodents [26]. The CYP2S1 gene acts in the synthesis of steroid hormones for the development of ovarian follicles from the preantral stage to the Graafian follicle stage in sheep [24]. These results suggest that genes identified in the present study participate in the development of the final stages of follicles, as well as in the low expression of IL1A, which decreases cellular apoptosis in the samples of oocytes.
The analysis of the DEGs, revealed that some of the main genes that are positively regulated are FGF18, CA5A, SVIP, CXCL14 and FST; while, that RENT, FEM1B, MMP9, ANKRD37 and CYP2S are down-regulating as shown in Table 2. The Fibroblast Growth Factor 18 (FGF18) is involved in the signaling of paracrine in the ovarian follicle, differentially activating the receptors depending on the properties of the extracellular matrix in the follicle, suggesting that it activates an apoptotic pathway through FGFR3c (fibroblast growth factor receptor 3) in GC as well as in other cell types [28]. Atresia of the subordinate follicle is also associated with an increase in the expression levels of FGF18 in theca cells and in the follicular fluid, which in turn acts on GC to inhibit steroidogenesis [29]. Carbonic anhydrase 5A (CA5A), is a zinc metalloprotein that catalyze the reversible conversion of CO 2 to bicarbonate ion and a proton [30]. In previous studies, CA5A activity has been reported in cellular functions in mammals, including processes in the ovary [30][31][32]. In these studies, its expression is reported on the surface of the epithelium of the ovary; and is absent from the corpora lutea in humans, pigs, rats, rabbits, mice, cats and dogs, and in the rabbit oocyte [32]. The activity of the CA5A gene in granulosa cells has been reported in rats and guinea pigs and is especially outstanding in small developing follicles before the stage of a recognizable antrum [33]. In our study the CA5A gene was identified in the zinc ion binding function in the MF of the GO categories. In humans, this gene is involved in the process of ureagenesis, where it supplies HCO 3− for carbamoyl phosphate synthetase and in gluconeogenesis for pyruvate carboxylase [33], and the low-expression or deficiency of CA5A results in impaired of bicarbonate which predisposes to hyperammonemia and hypoglycemia [34]. Follistatin (FST) is identified as member of the TGF-superfamily; the FST plays an important role in female physiology by regulating of follicle stimulating hormone (FSH) levels through blocking activin actions and conditional knockout females develop fertility defects, which include reduced litter sizes and, in more severe cases infertility [35]. In this study, the FST gene showed a significant expression, which is possibly related to reproductive processes in this breed.
Within the group of down-regulated genes, the Resistin (RETN) gene, it is associated with lipid metabolism, cause insulin resistance and decreased adipocyte differentiation [36]. The recombinant resistance (0.1, 1, 10 and 100 ng/mL) in the steroid hormone (i.e., progesterone, androstendione, testosterone, and estradiol) observed in prepubertal sows stimulates the steroidogenesis in the ovarian follicles, increasing the expression of the genes of the cytochrome family (CYP11A1 and CYP17A1) [37]. The Fem-1-Like Death Receptor-Binding Protein α (FEM1B) gene intervenes in the insulin secretion and acts on the degradation of Stem-Loop Binding protein (SLBP) [38]. In studies that used mice, an interaction was demonstrated between the products of the ANKRD37 (Ankyrin Repeat Domain 37) and FEM1B genes, which were highly expressed in male germ cells; the binding of the products of the FEM1B and ANKRD37 genes and their interaction resulted in degradation through the ubiquitin-mediated proteolysis pathway [39]. A functional variation (SNP rs10152450 and haplotype GGAAT) in the FEM1B gene increases insulin secretion and insulin resistance in mice, and the result is ovarian hyperandrogenism that occurs at the time of puberty in polycystic ovary syndrome (PCOS) [40]. In this study, the FEM1B gene showed a negative expression. In the GO functional annotation, the FEM1B gene is identified as a transcriptional regulator of ubiquitin-protein transferase activity. Reports by Lim et al. [41] suggest that chronic androgenization is a negative regulator of the growth of the antral follicle and the activity of the androgen receptor mediated by ubiquitin, resulting in antral follicle growth arrest in a chronically androgenized PCOS in rat. The Matrix Metalloproteinase 9 (MMP9) gene functions for homeostasis of the extracellular matrix and has a role in LH-induced ovarian steroidogenesis in mouse preovulatory follicles [42]. In Romanov ewes, the activity MMP9 is regulate by the events associated with follicular development, plays a role in the disappearance of P450 aromatase in CG and a progressive decrease in levels of P450 17alpha-hydroxylase in the theca internal, suggesting that MMPs can have a role in atresia of the ovary and is associated with distinct changes in levels of extracellular matrix following hypophysectomy [43]. We identified that the MMP-9 gene is down-regulated in prolific Pelibuey sheep, and metabolic pathway analyses showed that the MMP-9 gene have activity in the tumor necrosis factor (TNF)-and estrogen-signaling pathway. The Cytochrome P450, Family 2, Subfamily S, Polypeptide 1 (CYP2S1) has a role in the expressed during at all stages of development embryonic in mouse, endogenous substrates have as yet not been determined [44]. In the GO analysis, the CYP2S1 gene was not annotated, however Luo et al. [14] shown that the CYP2S1 gene was enriched in the metabolism of xenobiotics by cytochrome P450. The Homeobox protein gene HOX-A9 (HOXA9), is expressed in the reproductive tract of adult females, allowing further development and may play a role in the establishment and maintenance of pregnancy [45]. In recent human studies, up-expression of HOXA9 can stimulate the progression of ovarian cancer by promoting an immunosuppressive microenvironment through paracrine effects in peritoneal macrophages [46]. In our study, we observed the HOXA9 gene being expressed negatively, which may also be associated with the estrus stage.

Gene Ontology and KEGG Analysis: Metabolic Functions in Folliculogenesis
To investigate the functions of DEGs involved in the reproductive processes of prolific Pelibuey sheep, we performed the GO and metabolic pathway analysis, however only the six most important metabolic pathways are discussed here. Our analysis revealed that six genes identified with negative regulation are found in the calcium (Ca 2+ )-signaling pathway such as ATP2A3, ERBB3, MYLK4, PLCB2, PTGFR and PRKCB. The Ca 2+ regulates aspects of cell functions including cell cycle progression, maturation and fertilization of mammalian oocytes [47]. In the TGF-β-signaling pathway the up-regulated genes INHBA, FST and BMP15 were identified; and only the up-regulated gene BMP15 was identified in the ovarian steroidogenesis pathway, which is the metabolism of steroid hormones important for reproductive regulation [14]. The TGF-β-signaling pathway is one of the main metabolic pathways of folliculogenesis in mammals [13]. Previous research showed that INHBA functions in follicular maturation and steroidogenesis in Small Tail Han ewes, is associated with increased litter size [13], and BMP15 is expressed in oocytes regulating the proliferation and differentiation of GC in the final development in sheep [8]. In this study, the INHBA and BMP15 genes were up-regulated. The INHBA can act in a negative way by inhibiting the production of FSH, whereas the inhibitory effect of BMP15 on the action of FSH regulates the development of oocytes. However, a change of nucleotide modifies the protein, increasing the sensitivity to FSH and increasing the rate of ovulation and therefore prolificacy. Insulin signaling pathway, regulates blood glucose levels, maintenance of energy and glucose metabolism, also regulates ovarian androgens producction in humans, which suggests that it has an indirect role in early folliculogenesis, ovarian androgen is important to maintain the development of the follicle [48].
The NR4A1 gene was enriched in the pathway of the PI3K-Akt signaling pathway. The NR4A1 gene plays an important role in the development of the decomposition of the germinal vesicle (GVBD) which is regulated by FSH in cattle [49]. The IL1A gene was enriched in the MAPK signaling pathway, where has been reported as an important regulator of the maturation of oocytes in different species [50]. One study by Yu et al. [51] in ovaries of immature rats, suggests that FSH-activated p38 MAPK signal pathway regulates progesterone and estrogen production in GCs differentially.
Six enriched pathways are related to reproductive processes. However, four signaling pathways, the cAMP-signaling pathway, prolactin-signaling pathway, Rap1-signaling pathway and the insulin secretion pathway, have no association with reproductive function in the ovary. However, these pathways were derived from their corresponding physiological processes, we speculate that some of the genes involved in these pathways, may be related to some function of the reproductive processes of the ovary. This could agree with that not all the prolific sheep have point mutations reported as associated to litter size, and there are probably other molecular mechanisms that influence it.
The identifying genes in the Pelibuey sheep, previously reported in other sheep breeds, such as those that regulate the TGF-β-signaling pathway, are a sign that our analyzes are consistent with what was reported by other groups. We assume that the expression of the new genes identified in our study may also be due to other physiological factors, where a marked seasonality is not observed in this sheep and to nutrition, where they are an important factor that regulates reproductive performance of sheep and affects follicle development. Possibly these factors, may be stimulated some of the six enriched pathways found in the ovarian tissue, controlling the litter size. This information serves as a basis to design new experiments to evaluate prolificacy and control these factors in the prolificacy of Pelibuey sheep.

Conclusions
The analysis of ovarian tissue of sheep that gave birth to twin lambs and sheep that gave birth to a single lamb was performed using the new-generation sequencing technology RNA-seq. The most relevant set of genes identified such as BMP15, INHIB, SERPINE5, IL1A, CD5 and CA5A are associated with reproductive processes. In this study we report new genes, not previously documented in sheep such as ATP2A3, MMP9, CBL, SOCS3, FGF18, NR4A1, CHC4 and FST that are involved in the reproduction process. These genes can be considered candidates to regulate prolificacy, because they enriched six pathways with related to function in reproductive processes.
The information generated in this study provides the basis for understanding the reproductive processes in hair sheep such as Pelibuey, as well as providing additional data to aid in future research and breeding programs.