Skip to main content

Advertisement

Comparative analysis of differentially expressed genes between the ovaries from pregnant and nonpregnant goats using RNA-Seq

Abstract

Background

A multitude of genes tightly regulate ovarian follicular development and hormone secretion. These complex and coordinated biological processes are altered during pregnancy. In order to further understand the regulatory role of these genes during pregnancy, it is important to screen the differentially expressed genes (DEGs) in the ovaries of pregnant and nonpregnant mammals. To detect the genes associated with the development of pregnancy in goats, DEGs from the ovaries from pregnant and nonpregnant Anhui white goats (pAWGs and nAWGs, respectively) were analyzed using RNA sequencing technology (RNA-Seq).

Results

In this study, 13,676,394 and 13,549,560 clean reads were generated from pAWGs and nAWGs, respectively, and 1724 DEGs were identified between the two libraries. Compared with nAWGs, 1033 genes were upregulated and 691 genes were downregulated in pAWGs, including PGR, PRLR, STAR and CYP19A1, which play important roles in goat reproduction. Gene Ontology analysis showed that the DEGs were enriched for 49 functional GO terms. Kyoto Encyclopedia of Genes and Genomes analysis revealed that 397 DEGs were significantly enriched in 13 pathways, including “cell cycle”, “cytokine–cytokine receptor interaction” and “steroid biosynthesis”, suggesting that the genes may be associated with cell cycle regulation, follicular development and hormone secretion to regulate the reproduction process. Additionally, quantitative real-time PCR was used to verify the reliability of the RNA-Seq data.

Conclusions

The data obtained in this work enrich the genetic resources of goat and provide a further understanding of the complex molecular regulatory mechanisms occurring during the development of pregnancy and reproduction in goats. The DEGs screened in this study may play an important role in follicular development and hormone secretion and they would provide scientific basis for related research in the future.

Background

The goat (Capra hircus) can provide high-quality wool, meat, and other products and it is thus an important domestic and commercial animal. However, compared with other commercial animals, goat fecundity is relatively low, which has hindered the development of the goat industry. Anhui white goats (AWGs) have been used as a model to study goat fecundity in recent years [1] because of its precocious puberty, higher fertility, and higher leather quality than other types of goats [2].

Pregnancy is one of the most complex reproductive activities of mammals and is tightly regulated by the external environment, various endocrine factors, and the expression of a large number of genes [3, 4]. The ovary plays a vital role during pregnancy. There are significant differences in the activity and endocrine characteristics of the ovary during pregnancy and nonpregnancy [5]. In the nonpregnant phase, ovulation is normal and estrogen secretion dominates to maintain normal estrous function; in contrast, with the formation of the corpus luteum, ovulation is suspended during pregnancy and progesterone secretion is gradually increased to maintain pregnancy [6]. These differences between pregnancy and nonpregnancy may be associated with the distinct regulation of gene expression, which is implicated in various reproductive activities, including follicular development [7, 8], hormone secretion [9], ovulation [10], luteinization [11], and pregnancy maintenance [12].

Many studies have indicated that the differential expression of ovarian mRNA hormone receptor genes (follicle-stimulating hormone receptor [FSHR] and luteinizing hormone/choriogonadotropin receptor [LHR]) and growth factor genes (insulin-like growth factor 1 [IGF-I] and bone morphogenetic protein 6 [BMP6]) during follicular development significantly impact the development of preantral follicles [13,14,15]. Pramod et al. showed that the differential expression of bone morphogenetic protein 15 (BMP-15) and growth/differentiation factor (GDF-9) might be involved in higher prolificacy in goats [16]. Suppression subtractive hybridization has been used to screen differentially expressed genes (DEGs) in ovarian tissues between polytocous and monotocous goats, obtaining 29 differentially expressed sequence tags [17]. Microarray analysis has been used to analyze gene expression during early folliculogenesis in goats, and 2466 genes and three cell death-related pathways were explored [18]. In 2014, the different expression of genes in the ovaries from uniparous and multiparous AWGs were analyzed by high-throughput sequencing. In total, 2201 DEGs were identified, and 12 genes that may be associated with AWG high prolificacy were explored [19]. However, the information on the functions of specific genes in regulating goat pregnancy and reproduction is insufficient.

RNA sequencing technology (RNA-Seq) is mainly used for quantitative gene expression analysis of biological processes in a particular tissue or cells in specific species. With the development of high-throughput technologies, RNA-Seq has been extensively used to study specific genes in numerous species, including mice, pig, cattle, human, sheep, and goat [19,20,21,22,23,24]. It was reported that 338 genes were upregulated in the Jining grey goats and 404 were upregulated in Laiwu blank goats, and these genes might be involved in the regulation of goat fecundity and prolificacy [25]. Ling et al. used RNA-Seq to study the DEGs between the ovaries of uniparous and multiparous goats, and identified candidate genes for goat prolificacy [19].

In the present study, we analyzed the gene expression levels in ovaries from pregnant and nonpregnant AWGs (pAWGs and nAWGs, respectively) using an Ion Proton™ RNA-Seq platform. The results shed further light on the role of specific genes in reproductive biological processes, including hormone secretion, ovulation, luteinization formation, and pregnancy maintenance, which will help to identify genes that could potentially be used to regulate hircine reproduction and breeding practices in the future.

Methods

Animal materials

In the present study, AWGs raised under identical feeding conditions at the College of Animal Science and Technology (Anhui Agricultural University, Hefei, China) were used as animal samples. Six ewes, which were 2 years old and had a similar appearance, were selected from two groups for this study. Three of them had been pregnant for 3 months, and the rest were not pregnant and were undergoing anestrus, which was determined by hormone levels. After slaughter, both whole ovaries from each ewe were immediately collected and frozen in liquid nitrogen, followed by storage at − 80 °C. Half of the ovaries were used for RNA-Seq, and the other half were used for real-time PCR. All experimental procedures involving AWGs had been approved by the ethics committee of Anhui Agricultural University, Anhui, China (permit no. AHAU20101025).

Total RNA extraction and analysis for RNA-Seq

To minimize differences among individuals, one ovary from each pair of ovaries from the three nAWGs and three pAWGs were grouped for RNA extraction using RNAiso Plus (Takara, Dalian, China) following the manufacturer’s protocol. Briefly, the correct amount of frozen tissue sample was dissolved in 1 ml RNAiso Plus. After incubation at room temperature for 5 min, the homogenate was centrifuged at 12,000g for 5 min at 4 °C. The supernatant was then transferred to a new centrifuge tube. Chloroform (1/5 volume of RNAiso Plus) was then added to the supernatant. The solution was fully emulsified and incubated for 5 min at room temperature before centrifugation at 12,000g for 15 min at 4 °C. The supernatant was then carefully transferred to another centrifuge tube. An equal amount of isopropanol was added to the supernatant. The solution was fully mixed and incubated for 10 min at room temperature before an additional centrifugation at 12,000g for 10 min at 4 °C. Then, 1 ml of 75% ethanol was slowly added to the white precipitate, which was centrifuged at 12,000g for 5 min at 4 °C. Finally, the supernatant was discarded, and the precipitate was dried and dissolved in RNase-free ddH2O. The quantity and quality of the total RNA were measured using the Agilent 2100 Bioanalyzer system (Santa Clara, CA), and the samples were stored at − 80 °C until analysis.

Library preparation and sequencing

Two groups of total RNA were used for library preparation and sequencing by pooling equal quantities (10 μg) of total RNA isolated from 3 individual pregnant or 3 individual nonpregnant goat ovaries. Briefly, the total RNA samples were treated with RNase-free DNase I (Ambion Inc., Austin, TX) to degrade any possible genomic DNA contamination. The digested products were then purified by magnetic beads. Next, the mRNA was enriched using oligo (dT) magnetic beads (for eukaryotes) and fragmented into short fragments (approximately 200 bp) using fragmentation buffer. First-strand cDNA was then synthesized with random hexamer primers using the mRNA fragments as templates. Buffer, dNTPs, RNase H, and DNA polymerase I were added to synthesize the second strand cDNA. The double-stranded cDNA was purified with the QIAQuick PCR extraction kit (Qiagen, Germany) and eluted with elution buffer for end repair and poly (A) addition. Sequencing adapters were then ligated to the 5′ and 3′ ends of the fragments. Next, the ligation products were size selected and purified by agarose gel electrophoresis. Finally, the fragments were enriched by PCR amplification (95 °C for 10 min, followed by 40 cycles at 95 °C for 15 s and 60 °C for 45 s), purified by magnetic beads, and dissolved in the appropriate amount of elution buffer. During the quality control steps, an Agilent 2100 Bioanalyzer (Agilent, USA) was used for qualification and quantification of the sample library. Finally, the cDNA libraries were sequenced using an Ion Proton platform (Life Technologies, USA) at the Beijing Genomics Institute (BGI) in Shenzhen, China.

Sequencing data analysis

According to the requirements of this experiment and the principles of standard bioinformatics analysis, some contaminant reads, such as low-quality reads, adaptor sequences and reads with a length less than the threshold (30 bp), should be removed from the raw reads [26]. After filtering, the valid sequences, also known as the clean reads, were retained for further analysis. To identify the gene expression patterns in each genotype from the two libraries, the clean reads were mapped to goat reference gene sequences (http://goat.kiz.ac.cn/GGD/download.htm) and goat reference genome sequences (http://goat.kiz.ac.cn/GGD/download.htm) using SOAPaligner/SOAP2 [27]. No more than two mismatches were allowed in the alignment.

Identification of DEGs

To identify the DEGs between pAWGs and nAWGs, the gene expression levels were calculated using the reads per kilobase per million reads (RPKM) as follows [20]:

$$RPKM(A) = \frac{{10^{6} C}}{{{\raise0.7ex\hbox{${NL}$} \!\mathord{\left/ {\vphantom {{NL} {10^{3} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${10^{3} }$}}}}$$

where RPKM (A) is the expression level of gene A, C is the number of reads uniquely aligning to gene A, N is the total number of reads uniquely aligning to all genes, and L is the number of bases in gene A. DEGs were identified using a strict algorithm developed by BGI based on the method described by Audic and Claverie [28]. To compare the differential expression of the two samples, fold changes and p values were used, which calculated from the normalized expression using the following formulas:

$${\text{Fold-change}} = \log_{2} \left( {{\text{RPKM}}\;{\text{of}}\;{\text{Pergnant}}\;{\text{Ovary/RPKM}}\;{\text{of}}\;{\text{Nonpregnant}}\;{\text{Ovary}}} \right)$$

Formula for the p value:

$$2\sum\limits_{i = 0}^{{i = {\text{y}}}} {p\left( {i\left| x \right.} \right)} \;or\;2 \times \left( {1 - \sum\limits_{i = o}^{i = y} {p\left( {i\left| x \right.} \right)} } \right)\;\left( {if\sum\limits_{i = 0}^{i = y} {p\left( {i\left| x \right.} \right) > 0.5} } \right)$$
$$p\left( {y\left| x \right.} \right) = \left( {\frac{{N_{2} }}{{N_{1} }}} \right)^{y} \frac{{\left( {x + y} \right)!}}{{x!y!\left( {1 + \frac{{N_{2} }}{{N_{1} }}} \right)^{{\left( {x + y + 1} \right)}} }}$$

N1 and x represent the total clean reads count and the normalized expression level of a given gene in the nonpregnant ovary library, respectively. N2 and y represent the total clean reads count and the normalized expression level of a given gene in the pregnant ovary library, respectively. The false discovery rate (FDR) was used to determine the threshold p values in the multiple tests and perform analyses through the manipulation of the FDR values. In the present study, FDR ≤ 0.001 and |log2Ratio| ≥ 1 were used as the thresholds to judge the significance of the gene expression differences [29].

Real-time quantitative PCR

Eight DEGs were randomly selected to validate the RNA-Seq (quantification) data through real-time quantitative PCR (qPCR). The primers were designed using the Primer Express 5.0 software (ABI, China) (Table 1). GAPDH was used as a reference gene for data standardization. The left three ovaries from each nAWG or pAWG were pooled for total RNA extraction as described above. First-strand cDNA was synthesized using 1 μg of total RNA and the First-Strand cDNA Synthesis Kit (Toyobo, Japan). The qPCR reactions were performed on a Step One Plus™ Real-Time PCR System (Life Technologies, USA) using THUNDERBIRD SYBR qPCR Mix (Toyobo, Japan). The reaction contained 2.0 μl cDNA, 12.5 μl qPCR Mix, 2.0 μl of each primer, and 8.5 μl RNase-free H2O. The thermal cycling conditions were 95 °C for 1 min, followed by 40 cycles of 95 °C for 15 s, 58 °C for 20 s, and 72 °C for 20 s. All reactions were performed in triplicate. Relative quantification analyses were performed using the comparative CT method, and the relative gene expression levels were calculated using the 2−ΔΔCT method [30]. All data are expressed as the mean ± standard deviation and were analyzed by Student’s t-test with SPSS 17.0 for Windows (SPSS Inc. Released 2008. SPSS Statistics for Windows, Version 17.0. SPSS Inc, Chicago, USA). Data were considered statistically significant at p < 0.05.

Table 1 The primers used for real-time quantitative PCR analysis

GO function enrichment and KEGG pathway analysis of the DEGs

Gene Ontology (GO) enrichment analysis provided all the GO terms that were significantly enriched for the DEGs compared with the genome background and filtered the DEGs that correspond to biological functions. First, the DEGs were mapped to the GO terms in the database (http://www.geneontology.org/), and the gene numbers for each term were calculated. A hypergeometric test was then used to find significantly enriched GO terms among the DEGs compared with the genome background. The formula for the calculation is as follows:

$$p = 1 - \sum\limits_{i = 0}^{m - 1} {\frac{{\left( {\begin{array}{*{20}c} M \\ i \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {N - M} \\ {n - i} \\ \end{array} } \right)}}{{\left( {\begin{array}{*{20}c} N \\ n \\ \end{array} } \right)}}}$$

where N is the number of genes with GO annotations, n is the number of DEGs in N, M is the number of genes that are annotated to specific GO terms, and m is the number of DEGs in M. The calculated p value underwent Bonferroni correction using a corrected p value ≤ 0.05 as a threshold. GO terms fulfilling this condition were defined as significantly enriched GO terms for the DEGs. This analysis is able to recognize the main biological functions performed by the DEGs.

Using nr annotation, the gene annotation of the DEGs was conducted using Blast2GO [31]. After obtaining the GO annotation for the DEGs, Web Gene Ontology Annotation Plot (WEGO) software was used to perform GO functional classification of the DEGs and understand the distribution of the gene functions from the species at the macro level [32].

The Kyoto Encyclopedia of Genes and Genomes (KEGG) is the major public pathway-related database [33]. Pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways for the DEGs compared with the whole genome background. The formula for calculating the pathways is the same as for the GO analysis, with N as the number of all genes with KEGG annotations, n as the number of DEGs in N, M as the number of genes annotated to specific pathways, and m as the number of DEGs in M. Pathways with Q values ≤ 0.05 are significantly enriched for the DEGs.

Results

Overview of sequencing data

In this study, we generated cDNA libraries of ovaries from pregnant and nonpregnant goat (Capra hircus) (Fig. 1). Approximately 14.39 Gb reads, including 7.12 Gb from pAWGs and 7.27 Gb from nAWGs, were obtained through the sequencing of two ovary cDNA libraries from pAWGs and nAWGs using the Ion Proton sequencing platform at BGI in Shenzhen, China. After discarding sequences shorter than 30 bp, filtering the adaptor sequences, and eliminating low-quality sequences from the raw data, 13,676,394 and 13,549,560 clean reads were ultimately obtained from the ovary libraries for pAWGs and nAWGs, respectively, and were retained for further analysis. The clean reads from each library were sufficient for the quantitative analysis of gene expression.

Fig. 1
figure1

Summary of genetic analysis procedures for pregnant and unpregnant goats

According to the mapping analysis, approximately 98% of the reads (98.01% for nAWGs and 97.95% for pAWGs) could be mapped to the goat genome sequence. When the unique matches were regarded, slightly more than 92% (92.50% for nAWGs and 92.72% for pAWGs) of the reads corresponded to the genome, and approximately 45% of the reads (45.03% for nAWGs and 43.47% for pAWGs) could be perfectly matched to the genome.

When the reference genes were used to analyze the data, only 75.91% of the reads in nAWGs and 75.64% of the reads in pAWGs could be mapped to reference genes. Regarding the unique matches, 54.14% of the reads in nAWGs and the 54.54% in pAWGs corresponded to the reference genes. Additionally, 37.83% of the reads in nAWGs and 36.43% of the reads in pAWGs could be perfectly matched to the reference genes. The major characteristics of the two libraries are described in Additional file 1.

Identification and analysis of the DEGs in the ovaries from pregnant and nonpregnant goats

In total, 20,651 genes were detected in the study, among which 17,181 genes were co-expressed in the two libraries and 2002 and 1468 genes were specifically expressed in the ovary libraries from pAWGs and nAWGs, respectively (Additional file 2).

Analysis of the gene coverage distribution indicated that 16% of genes had 70–90% coverage, while 26% of genes in the nAWG and 27% of genes in the pAWG libraries had 90–100% coverage (Additional file 3: Fig. S1), which means the read distributions were similar between the two libraries.

To detect pregnancy-related genes, the expression levels of all the detected genes from pAWGs and nAWGs were calculated using the RPKM method. An FDR ≤ 0.001 and the absolute value of the log2 ratio ≥ 1 were used as criteria to judge the significant differences in the expressed genes. A total of 1724 DEGs were detected between the nAWG and pAWG libraries, with 1033 upregulated genes and 691 downregulated genes in pAWGs compared to nAWGs (Additional file 4: Fig. S2 and Additional file 5). Some of the DEGs were related to reproduction and reproductive processes, such as progesterone receptor (PGR), prolactin receptor precursor (PRLR), steroidogenic acute regulatory protein (STAR), cytochrome P450 19A1 (CYP19A1) and so on (Table 2). Among the DEGs, 15 genes and 71 genes were found to be specifically expressed in the nAWG and pAWG libraries, respectively. Among the specifically expressed DEGs, we screened four DEGs related to reproduction, including two genes [beta-defensin 112 (DEFB112) and interferon epsilon (IFNE)] that were specifically expressed in nAWGs and two genes [somatostatin precursor (SST) and somatostatin receptor type 1 (SSTR1)] that were specifically expressed in pAWGs.

Table 2 Details of the differentially expressed genes (DEGs) associated with reproduction

QPCR validation

The expression profiles of eight randomly selected DEGs were validated using qPCR analysis. As shown in Fig. 2, the expression of SST, SLC1A2, TSPAN1, and OSTF1 was upregulated, whereas the expression of DEFB112, UGDH, APOA2, and FADS1 was downregulated in pAWGs compared with nAWGs. These expression patterns were largely consistent with the RNA-Seq (Quantification) results, which means that the RNA-Seq results were reliable and that RNA-Seq could be used to reliably and accurately perform mRNA differential expression analysis.

Fig. 2
figure2

Comparison of RNA-Seq (Quantification) and qPCR results. Results represent the mean (± SD) of three experiments. *p ≤ 0.05

GO functional enrichment and KEGG pathway analysis of the DEGs

GO is an international standardized classification system for gene function that supplies a set of controlled terms to comprehensively describe the properties of genes and gene products. In our present study, GO function enrichment was used to screen the DEGs that were significantly associated with different biological functions. The well-annotated gene sequences were assigned to 23, 15, and 11 functional groups in the “biological process”, “cellular component”, and “molecular function” categories, respectively (Additional file 6). Many of the enriched GO terms are reproduction-related, such as “signaling”, “reproduction”, “reproduction process”, “developmental process”, “growth”, “organelle”, “receptor activity”, “biological regulation”, “regulation of biological process”, “molecular transducer activity” and “enzyme regulator activity” (Fig. 3). These GO terms were mostly associated with cell proliferation, apoptosis, follicular development, hormone secretion, and pregnancy maintenance. Among the reproduction-related GO terms, some important genes related to reproduction were detected, including PGR, CYP19A1, CCNB1, STAR, SST, SSTR1, estrogen receptor beta 2 (ESR2), BMP receptor 1B (BMPR1B) and erythropoietin receptor (EPOR).

Fig. 3
figure3

Histogram of the Gene Ontology (GO) classifications. The results are summarized for the following three main categories: “biological process”, “cellular component”, and “molecular function”. Among these groups, the GO terms “cellular process”, “binding”, “cell”, and “cell part” predominated in each of the three main categories. The reproduction-related GO terms were underlined

Analysis of the KEGG pathway enrichment revealed that the 1413 DEGs participated in 242 pathways (Additional file 7), among which the largest category was “metabolic pathways” (ko01100, 11.18%), which had 158 annotated genes. Furthermore, 397 DEGs were significantly enriched (p ≤ 0.05) in 13 pathways, several of which were related to reproduction, such as “cytokine–cytokine receptor interaction”, “cell cycle”, and “steroid biosynthesis” (Figs. 4 and 5, Additional file 8: Figure S3, Additional file 9: Figure S4, Additional file 10: Figure S5 and Table 3). And these DEGs up- and down-regulated also shown in Fig. 5. The significantly enriched DEGs, including PRLR, BMPR1B, EPOR, kit Hardy (KIT), interleukin 1 beta (IL1B), cyclin A2 (CCNA2), cyclin B1 (CCNB1), and proliferating cell nuclear antigen (PCNA), have been reported to be involved in follicular development, hormone secretion, luteinization, and pregnancy maintenance.

Fig. 4
figure4

Top 20 pathway enrichment assignments for nonpregnant and pregnant ovaries. The Q value is the corrected p value ranging from 0 to 1. The lower the Q value, the greater the intensity. Red lines are significantly enriched signal pathways

Fig. 5
figure5

“Cytokine–cytokine receptor interaction”, “steroid biosynthesis” and “cell cycle” pathway and the identified regulator

Table 3 The significantly enriched pathways among the DEGs between the pregnant and nonpregnant

Discussion

We successfully identified the gene expression profiles of ovarian tissues from pregnant and nonpregnant AWGs using RNA-Seq (Quantification) and analyzed the gene expression differences between the two libraries. Furthermore, the GO enrichment and KEGG pathways for the DEGs were analyzed.

In our present study, 1724 genes were found to be significantly differentially expressed between the two libraries. Among them, a large number of DEGs were associated with prolificacy processes, follicular development, ovulation control, hormone secretion, and luteinization or pregnancy maintenance. For example, genes previously associated with reproduction, including PRLR, LHR, STAR, SPP1, BMP6, 3β-hydroxysteroid dehydrogenase estradiol (3BHSD), 17-beta-dehydrogenase 12 (HSD17B12), and BMPR1B, were found to be upregulated in pAWGs compared to nAWGs. The upregulated expression of PRLR mRNA during pregnancy may be associated with corpus luteum function for maintaining pregnancy [34]. LHR is a specific membrane receptor on the granulosa and theca cells that binds to luteinizing hormone (LH), resulting in androgen and progesterone production. LHR may impact follicular development from the primary follicle stage onwards [35]. 3BHSD plays a crucial role in the biosynthesis of all classes of hormonal steroids [36]. HSD17B12 plays an important role in female fertility [37]. STAR, which is the key rate-limiting enzyme in the production of steroid hormones, is involved in the regulation of steroid hormone biosynthesis and follicular development in the mammalian ovary [9]. SPP1, which is also known as OPN (osteopontin), is localized to luteal cells during the luteal phase of the estrous cycle and is involved in the development and regression of the corpus luteum [11]. Previous research has indicated that BMP6 may regulate ovary granulosa cell steroidogenesis, ovarian hormone secretion, follicular growth, and luteinization [38, 39]. A BMPRIB mutation may decrease the ability of BMP to inhibit differentiation, increase INHBA (inhibin beta A) mRNA expression in granular cells, and enhance the responsiveness of granule cells to gonadotropin [40]. PGR, FSHR, CYP19A1, ESR2, INHA, NHBA, estradiol 17-beta-dehydrogenase 1 (HSD17B1), and inhibin beta B (INHBB) were found to be downregulated in pAWG. Recent studies showed that PGR, ESR2, FSHR, CYP19A1, HSD17B1, INHA, INHBA and INHBB played a very important role in hormone secretion, follicle growth, cell proliferation and apoptosis [41,42,43].

Our data showed some uniquely transcribed genes in nAWGs and pAWGs. Among these specifically expressed genes, two genes (DEFB112 and IFNE) from nAWGs and two genes (SST and SSTR1) from pAWGs may be involved in related events during pregnancy. DEFB112 is a type of B-defensin that is preferentially expressed in the male reproductive system and actively involved in sperm maturation and capacitation but is poorly documented. Further study is required to determine whether DEFB112 is involved in the reproductive process in female animals [44]. IFNE is a novel type I IFN that has been detected in the lung, reproductive tissue, intestine and so on. IFNE may protect the female reproductive tract from viral and bacterial infection. Its expression was found to be at the lowest levels during diestrus and the highest levels at estrus [45]. IFNE can activate the JAK-STAT signaling pathway, which plays a very important role in reproduction [46]. SST is a 14-amino acid polypeptide with a short half-life that may inhibit the secretion of various hormones including GH and LH. SSTR1 is somatostatin receptor subtype 1. SST and its receptors are present in the ovary in various species. It was reported that SSTR2 and 5 activation modulate ovarian steroidogenesis by upregulating endogenous BMP activity in growing follicles [47]. Previous studies have shown that SSTR1 inhibits cAMP accumulation by coupling to pertussis toxin-sensitive G proteins [48]. However, the specific functions of SSTR1 are unclear and require further study.

Regarding GO enrichment analysis, we found that some of the correlated genes were enriched into reproduction-related GO terms. For example, CYP19A1, ESR2, PGR INHA, INHBA and STAR were enriched into GO terms for reproduction, reproduction process, molecular transducer activity, and transporter activity. The results indicated that these genes may play a specific role in goat reproduction. For the uniquely transcribed genes, the two genes (ENC1 and RGS17) from nAWGs were enriched for signaling, enzyme regulator activity, and molecular transducer activity GO terms, while the six genes (DCDC2, TBC1D15, PTHRP, SLC1A2, RDH16 and SST) from pAWGs were enriched for biological regulation, enzyme regulator activity, growth, metabolic process, and response to stimulus GO terms. Our data indicated that these genes may be related to cell proliferation and hormone secretion.

Three pathways, including “cytokine–cytokine receptor interaction”, “steroid biosynthesis” and “cell cycle”, were selected from the 13 significantly enriched pathways and were found to be intimately associated with mammalian reproduction. Most of the DEGs participating in cytokine–cytokine receptor interaction and steroid biosynthesis pathways were upregulated, including EPOR, KIT, IL1B, PRLR and BMPR1B, indicating that these pathways are activated in the ovaries of pregnant goats. It was shown that EPOR and IL1B might play a role in establishing pregnancy [49, 50]. KIT plays a very important role in the maintenance of the primordial follicle reserve and in primary to secondary follicle transition [51]. Most DEGs were downregulated in the cell cycle pathway, such as CCNA2, CCNB1 and so on, and only one gene (cyclin D1, CCND1) was unregulated. CCNA2 regulates cell cycle progression by interacting with CDK kinases and is involved in the G2/M transition. CCNB1 is a regulatory protein involved in mitosis. Research indicates that the downregulated expression of CCNA2 may promote primordial follicle activation, while CCNB1 plays a significant role in promoting the maturation of large-follicle oocytes [52, 53]. CCND1 is a protein required for progression through G1 phase of the cell cycle and its main function is to promote cell proliferation [54]. It can be inferred that these genes may be involved in cell cycle regulation, granulosa cell proliferation and follicular development.

Furthermore, some other signaling pathways associated with reproduction were also identified in our current analysis. These pathways included the ECM-receptor interaction, focal adhesion, and Jak-STAT, Wnt, insulin, and p53 signaling pathways, suggesting that some molecules in these pathways might also be involved in goat reproduction or maintaining the physiological activity of the ovary. However, the exact relationship between these signaling networks is not fully understood and further research is needed. Generally, the current KEGG pathway analyses have provided a solid basis for future work.

Conclusion

In conclusion, this study presents the transcriptomic profiles of ovarian tissue from pregnant and nonpregnant AWGs using RNA-Seq technology. 1724 genes were detected to be significantly differentially expressed in pAWGs compared with nAWGs. Among them, 1033 genes were upregulated and 691 were downregulated in pAWGs, suggesting that DEGs may play an important role in the regulation of goat ovarian function. GO annotation and KEGG pathway analyses were conducted on the DEGs between the two libraries, and 13 pathways were pointed out for significantly enriched, in which “cytokine–cytokine receptor interaction”, “cell cycle”, and “steroid biosynthesis” were found to be associated with reproduction. The DEGs and pathways identified could facilitate further to elucidate the functions of DEGs in goat follicular development and hormone secretion, which will help us to understand the relationship between genes and mammalian reproduction and provide a solid foundation for future studies.

References

  1. 1.

    Ling YH, Ren CH, Guo XF, Xu LN, Huang YF, Luo JC, et al. Identification and characterization of microRNAs in the ovaries of multiple and uniparous goats (Capra hircus) during follicular phase. BMC Genomics. 2014;15:339.

  2. 2.

    Chen S, Cheng GL, Zhu DJ, Jiang XC, Zhao HL. Determination on the body properties and meat performance of Anhui white goat. Anim Husb Feed Sci. 2009;30:150–2.

  3. 3.

    Gram A, Boos A, Kowalewski MP. Uterine and placental expression of canine oxytocin receptor during pregnancy and normal and induced parturition. Reprod Domest Anim. 2014;49:41–9.

  4. 4.

    Agaoglu OK, Agaoglu AR, Guzeloglu A, Kurar E, Kayis SA, Ozmen O, et al. Expression of hypoxia-inducible factors and vascular endothelial growth factor during pregnancy in the feline uterus. Theriogenology. 2015;84:24–33.

  5. 5.

    Stephens SM, Moley KH. Follicular origins of modern reproductive endocrinology. Am J Physiol Endocrinol Metab. 2009;297:E1235–6.

  6. 6.

    Ling YH, Guo XF, Chen T, Ding JP, Ma YH, Chu MX, et al. Characterization and analysis of differentially expressed microRNAs in hircine ovaries during the follicular and luteal phases. Anim Reprod Sci. 2016;166:47–57.

  7. 7.

    Orisaka M, Jiang JY, Orisaka S, Kotsuji F, Tsang BK. Growth differentiation factor 9 promotes rat preantral follicle growth by up-regulating follicular androgen biosynthesis. Endocrinology. 2009;150:2740–8.

  8. 8.

    Sargent KM, Lu N, Clopton DT, Pohlmeier WE, Brauer VM, Ferrara N, et al. Loss of vascular endothelial growth factor A (VEGFA) isoforms in granulosa cells using pDmrt-1-Cre or Amhr2-Cre reduces fertility by arresting follicular development and by reducing litter size in female mice. PLoS ONE. 2015;10:e0116332.

  9. 9.

    Yamashita H, Murayama C, Takasugi R, Miyamoto A, Shimizu T. BMP-4 suppresses progesterone production by inhibiting histone H3 acetylation of StAR in bovine granulosa cells in vitro. Mol Cell Biochem. 2011;348:183–90.

  10. 10.

    Galloway SM, McNatty KP, Cambridge LM, Laitinen MP, Juengel JL, McLaren RJ, et al. Mutations in an oocyte-derived growth factor gene (BMP15) cause increased ovulation rate and infertility in a dosage-sensitive manner. Nat Genet. 2000;25:279–83.

  11. 11.

    Poole DH, Ndiaye K, Pate JL. Expression and regulation of secreted phosphoprotein 1 in the bovine corpus luteum and effects on T cell lymphocyte chemotaxis. Reproduction. 2013;146:527–37.

  12. 12.

    Piccinni MP, Maggi E, Romagnani S. Role of hormone-controlled T-cell cytokines in the maintenance of pregnancy. Biochem Soc Trans. 2000;28:212–5.

  13. 13.

    Saraiva MV, Celestino JJ, Araújo VR, Chaves RN, Almeida AP, Lima-Verde IB, et al. Expression of follicle-stimulating hormone receptor (FSHR) in goat ovarian follicles and the impact of sequential culture medium on in vitro development of caprine preantral follicles. Zygote. 2011;19:205–14.

  14. 14.

    Magalhães-Padilha DM, Duarte ABG, Araújo VR, Saraiva MVA, Almeida AP, Rodrigues GQ, et al. Steady-state level of insulin-like growth factor-I (IGF-I) receptor mRNA and the effect of IGF-I on the in vitro culture of caprine preantral follicles. Theriogenology. 2012;77:206–13.

  15. 15.

    Frota IM, Leitão CC, Costa JJ, van den Hurk R, Saraiva MV, Figueiredo JR, et al. Levels of BMP-6 mRNA in goat ovarian follicles and in vitro effects of BMP-6 on secondary follicle development. Zygote. 2013;21:270–8.

  16. 16.

    Pramod RK, Sharma SK, Singhi A, Pan S, Mitra A. Differential ovarian morphometry and follicular expression of BMP15, GDF9 and BMPR1B influence the prolificacy in goat. Reprod Domest Anim. 2013;48:803–9.

  17. 17.

    An XP, Hou JX, Li G, Peng JY, Liu XQ, Liu HY, et al. Analysis of differentially expressed genes in ovaries of polytocous versus monotocous dairy goats using suppressive subtractive hybridization. Reprod Domest Anim. 2012;47:498–503.

  18. 18.

    Magalhães-Padilha DM, Geisler-Lee J, Wischral A, Gastal MO, Fonseca GR, Eloy YR, et al. Gene expression during early folliculogenesis in goats using microarray analysis. Biol Reprod. 2013;89:19.

  19. 19.

    Ling YH, Xiang H, Li YS, Liu Y, Zhang YH, Zhang ZJ, et al. Exploring differentially expressed genes in the ovaries of uniparous and multiparous goats using the RNA-Seq (Quantification) method. Gene. 2014;550:148–53.

  20. 20.

    Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5:621–8.

  21. 21.

    Gunawan A, Sahadevan S, Neuhoff C, Große-Brinkhaus C, Gad A, Frieden L, et al. RNA deep sequencing reveals novel candidate genes and polymorphisms in boar testis and liver tissues with divergent androstenone levels. PLoS ONE. 2013;8:e63259.

  22. 22.

    Huang W, Khatib H. Comparison of transcriptomic landscapes of bovine embryos using RNA-Seq. BMC Genomics. 2010;11:711.

  23. 23.

    Hackett NR, Butler MW, Shaykhiev R, Salit J, Omberg L, Rodriguez-Flores JL, et al. RNA-Seq quantification of the human small airway epithelium transcriptome. BMC Genomics. 2012;13:82.

  24. 24.

    Chen HY, Shen H, Jia B, Zhang YS, Wang XH, Zeng XC. Differential gene expression in ovaries of Qira Black sheep and Hetian sheep using RNA-Seq technique. PLoS ONE. 2015;10:e0120170.

  25. 25.

    Miao X, Luo Q, Qin X. Genome-wide transcriptome analysis in the ovaries of two goats identifies differentially expressed genes related to fecundity. Gene. 2016;582:69–76.

  26. 26.

    Kerpedjiev P, Frellsen J, Lindgreen S, Krogh A. Adaptable probabilistic mapping of short reads using position specific scoring matrices. BMC Bioinform. 2014;15:100.

  27. 27.

    Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009;25:1966–7.

  28. 28.

    Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997;7:986–95.

  29. 29.

    Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-Seq data. Bioinformatics. 2010;26:136–8.

  30. 30.

    Ishida-Takagishi M, Enomoto A, Asai N, Ushida K, Watanabe T, Hashimoto T, et al. The Dishevelled-associating protein Daple controls thenon-canonical Wnt/Rac pathway and cell motility. Nat Commun. 2012;3:859.

  31. 31.

    Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.

  32. 32.

    Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, et al. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;34:W293–7.

  33. 33.

    Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:D480–4.

  34. 34.

    Zi XD, Chen DW, Wang HM. Molecular characterization, mRNA expression of prolactin receptor (PRLR) gene during pregnancy, nonpregnancy in the yak (Bos grunniens). Gen Comp Endocrinol. 2012;175:384–8.

  35. 35.

    Phoophitphong D, Srisuwatanasagul S, Tummaruk P. Immunohistochemical localization of luteinizing hormone receptor in the cyclic gilt ovary. Anat Histol Embryol. 2017;46:94–100.

  36. 36.

    Lachance Y, Luu-The V, Labrie C, Simard J, Dumont M, de Launoit Y, et al. Characterization of human 3 beta-hydroxysteroid dehydrogenase/delta 5-delta 4-isomerase gene and its expression in mammalian cell. J Biol Chem. 1990;265:20469–75.

  37. 37.

    Kemiläinen H, Adam M, Mäki-Jouppila J, Damdimopoulou P, Damdimopoulos AE, Kere J, et al. The hydroxysteroid (17β) dehydrogenase family Gene HSD17B12 is involved in the prostaglandin synthesis pathway, the ovarian function, and regulation of fertility. Endocrinology. 2016;157:3719–30.

  38. 38.

    Campbell BK, Kendall NR, Baird DT. Effect of direct ovarian infusion of bone morphogenetic protein 6 (BMP6) on ovarian function in sheep. Biol Reprod. 2009;81:1016–23.

  39. 39.

    Khalaf M, Morera J, Bourret A, Reznik Y, Denoual C, Herlicoviez M, et al. BMP system expression in GCs from polycystic ovary syndrome women and the in vitro effects of BMP4, BMP6, and BMP7 on GC steroidogenesis. Eur J Endocrinol. 2013;168:437–44.

  40. 40.

    Fabre S, Pierre A, Pisselet C, Mulsant P, Lecerf F, Pohl J, et al. The Booroola mutation in sheep is associated with an alteration of the bone morphogenetic protein receptor-IB functionality. J Endocrinol. 2003;177:435–44.

  41. 41.

    Komatsu K, Masubuchi S. The concentration-dependent effect of progesterone on follicle growth in the mouse ovary. J Reprod Dev. 2017;63:271–7.

  42. 42.

    Chakravarthi VP, Khristi V, Ghosh S, Yerrathota S, Dai E, Roby KF, et al. ESR2 Is essential for gonadotropin-induced Kiss1 expression in granulosa cells. Endocrinology. 2018;159:3860–73.

  43. 43.

    Cadoret V, Frapsauce C, Jarrier P, Maillard V, Bonnet A, Locatelli Y, et al. Molecular evidence that follicle development is accelerated in vitro compared to in vivo. Reproduction. 2017;153:493–508.

  44. 44.

    Patil AA, Cai Y, Sang Y, Blecha F, Zhang G. Cross-species analysis of the mammalian beta-defensin gene family: presence of syntenic gene clusters and preferential expression in the male reproductive tract. Physiol Genomics. 2005;23:5–17.

  45. 45.

    Fung KY, Mangan NE, Cumming H, Horvat JC, Mayall JR, Stifter SA, et al. Interferon-ε protects the female reproductive tract from viral and bacterial infection. Science. 2013;339:1088–92.

  46. 46.

    Yang L, Xu L, Li Y, Li J, Bi Y, Liu W. Molecular and functional characterization of canine interferon-epsilon. J Interferon Cytokine Res. 2013;33:760–8.

  47. 47.

    Nakamura E, Otsuka F, Inagaki K, Tsukamoto N, Ogura-Ochi K, Miyoshi T, et al. Involvement of bone morphogenetic protein activity in somatostatin actions on ovarian steroidogenesis. J Steroid Biochem Mol Biol. 2013;134:67–74.

  48. 48.

    Hadcock JR, Strnad J, Eppler CM. Rat somatostatin receptor type 1 couples to G proteins and inhibition of cyclic AMP accumulation. Mol Pharmacol. 1994;45:410–6.

  49. 49.

    Geisert R, Fazleabas A, Lucy M, Mathew D. Interaction of the conceptus and endometrium to establish pregnancy in mammals: role of interleukin 1β. Cell Tissue Res. 2012;349:825–38.

  50. 50.

    Ji YQ, Zhang YQ, Li MQ, Du MR, Wei DD, Li DJ. EPO improves the proliferation and inhibits apoptosis of trophoblast and decidual stromal cells through activating STAT-5 and inactivating p38 signal in human early pregnancy. Int J Clin Exp Pathol. 2011;4:765–74.

  51. 51.

    John GB, Shidler MJ, Besmer P, Castrillon DH. Kit signaling via PI3K promotes ovarian follicle maturation but is dispensable for primordial follicle activation. Dev Biol. 2009;331:292–9.

  52. 52.

    Tong Y, Li F, Lu Y, Cao Y, Gao J, Liu J. Rapamycin-sensitive mTORC1 signaling is involved in physiological primordial follicle activation in mouse ovary. Mol Reprod Dev. 2013;80:1018–34.

  53. 53.

    Liu HL, Gao Y, Zhai B, Jiang H, Ding Y, Zhang L, et al. The Effects of polyadenylation status on MPFs during in vitro porcine oocyte maturation. Cell Physiol Biochem. 2016;39:1735–45.

  54. 54.

    Fu M, Wang C, Li Z, Sakamaki T, Pestell RG. Minireview: cyclin D1: normal and abnormal functions. Endocrinology. 2004;145:5439–44.

Download references

Authors’ contributions

These studies were designed by QQ, QZ and YHL. QQ and QZ carried out all the experimental analyses and prepared all figures and tables. YHL and QQ analyzed the data and drafted the manuscript. FFG, CMX, ZXR LY, and LWY contributed to revisions of the manuscript. YHL assisted in explaining the results and revised the final version of the manuscript. All authors read and approved the final manuscript.

Acknowledgements

No application.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Consent for publication

All authors have consented the manuscript been published.

Ethics approval and consent to participate

All procedures involving animals were approved by the Ethics Committee of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAHCAAS201535). We comply with all the points claimed by the committee.

Funding

This research was supported by the National Natural Science Foundation of China (31772566), the Natural Science Foundation of Anhui Province (1708085MC61), the Key Research Projects of Natural Science in Anhui Colleges and Universities (KJ2017A334) and the Agricultural Science and Technology Innovation Program of China (ASTIP-IAS13).

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Correspondence to Yinghui Ling.

Additional files

40709_2019_95_MOESM1_ESM.xls

Additional file 1. Summary of the sequenced reads mapped to the Capra hircus reference genes and genome.

40709_2019_95_MOESM2_ESM.xls

Additional file 2. All detected genes in the two libraries.

40709_2019_95_MOESM3_ESM.tif

Additional file 3. Gene coverage distribution in the two ovary libraries.

40709_2019_95_MOESM4_ESM.tif

Additional file 4. plot indicating the log-transformed gene expression levels and DEG distributions between the pregnant and nonpregnant samples.

40709_2019_95_MOESM5_ESM.xls

Additional file 5. Differential gene detected in two libraries.

40709_2019_95_MOESM6_ESM.xls

Additional file 6. Gene Ontology (GO) classifications of two libraries.

40709_2019_95_MOESM7_ESM.xls

Additional file 7. KEGG pathway analysis of two libraries.

40709_2019_95_MOESM8_ESM.tif

Additional file 8. Cytokine–cytokine receptor interaction signal pathway.

40709_2019_95_MOESM9_ESM.tif

Additional file 9. Steroid biosynthesis signal pathway.

40709_2019_95_MOESM10_ESM.tif

Additional file 10. Cell cycle pathway.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Differentially expressed genes
  • Ovary
  • Pregnant and nonpregnant
  • Goat
  • RNA-Seq