- Email: [email protected]
Transcriptome studies of granulosa cells at different stages of ovarian follicular development in buffalo
Animal Reproduction Science 187 (2017) 181–192
Contents lists available at ScienceDirect
Animal Reproduction Science journal homepage: www.elsevier.com/locate/anireprosci
Transcriptome studies of granulosa cells at different stages of ovarian follicular development in buffalo J. Lia,b, Z. Lia,b, S. Liua,b, R. Ziaa, A. Lianga,b, L. Yanga,b, a b
MARK
⁎
College of Animal Science and Technology, Huazhong Agricultural University, Wuhan, Hubei, China Hubei Province’s Engineering Research Center in Buffalo Breeding & Products, China
AR TI CLE I NF O
AB S T R A CT
Keywords: Buffalo Differentially expressed genes Follicle growth RNA-seq Transcriptome
The normal maturation and ovulation from ovarian follicles is important in ensuring conception and improving fertility of buffalo. The molecular regulation mechanism of buffalo follicles growth, however, remains unknown. This study analyzed the gene expression profiles associated with buffalo ovarian follicle growth. According to the analysis of RNA sequencing, 17,700 unigenes and 13,672 differentially expressed genes (DEGs) were detected. A total of 30 common DEGs were identified during four stages of follicle growth, and the expression patterns are basically synchronized, suggesting the products as a result of expressions of these genes may cooperate to regulate follicular development. Furthermore, GO and KEGG enrichment analyses revealed that the majority of DEGs in early stage of follicular growth were enriched in ribosomal and oxidative phosphorylation signaling pathways, and the expression patterns of these DEGs are basically up-regulated at the beginning of follicular growth (< 8 mm, diameter), and then downregulated (8–12 mm) in the following stages of follicular development. The pathway of immune signaling, including allograft rejection, chemokine signaling pathway, natural killer cell mediated cytotoxicity, phagosome, and antigen processing and presentation, was significantly enriched in the last stage of follicular development (> 12 mm), which indicates that the immune system has an important role in the last stage of follicular maturation and ovulation. This study provided a gene expression profile of buffalo follicle growth, and provided an insight into biological processes associated with molecular regulation of ovarian follicle growth.
1. Introduction Buffalo, one of the most important economic livestock worldwide, is known for its resistance to stress and its milk product with high milk fat and protein contents. The number of buffalo stock worldwide is estimated to be approximately 172 million (FAO: http://faostat.fao.org/). Unlike swine, sheep, and goats, however, where multiple follicles undergo maturation and from which ovulation occurs, buffalo is a monovular species, and thus, it normally has only one follicle develop into mature follicle from which ovulation occurs (Ginther et al., 2001) resulting in a lesser fecundity in buffalo. This is an important factor restricting the development of the buffalo breeding industry. Ovaries are critical organs for follicular development and ovulation. The follicle is the structural and functional units of the ovaries for which normal development during the developmental stages is important for lifelong fertility of buffalo. Follicle development results with granulosa cell growth, development, differentiation, and apoptosis in follicles, which is closely related to growth of oocyte and follicular atresia (Jain et al., 2016; Yeung et al., 2017). Two or three waves of ovarian follicular development emerge ⁎
Corresponding author at: No.1 Shizishan Street, Wuhan, Hubei Province, 430070, China. E-mail address: [email protected] (L. Yang).
http://dx.doi.org/10.1016/j.anireprosci.2017.11.004 Received 7 August 2017; Received in revised form 17 October 2017; Accepted 3 November 2017 Available online 06 November 2017 0378-4320/ © 2017 Elsevier B.V. All rights reserved.
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Table 1 Results of RNA-seq, assembly, and annotion. Stage
Raw reads
Clean reads
Assembled reads
Unigenes
DEGs
SDEGs
GC1 GC2 GC3 GC4 Total
50,760,648 50,760,170 50,760,644 54,035,416 206,316,878
46,219,792 46,675,486 46,192,560 46,066,398 185,154,236(89.74%)
29,117,053 29,120,413 30,062,845 28,098,004 116,398,315(62.87%)
16,249 16,339 16,167 16,349 17,700
– 13160 13048 13228 13,672(77.24%)
– 299 303 691 1076(7.87%)
DEG, differentially expressed genes based on the previous stage; SDEG, significantly differentially expressed genes based on the previous stage with FDR ≤ 0.05 and |log2Ratio| ≥ 1.
during each estrous cycle (Fortune et al., 1991), and the emergence of each follicular wave is synchronous with the surge of a folliclestimulating hormone (FSH) (Evans et al., 1997). The FSH has a pivotal role during the early stage of follicular development, ensuring follicle recruitment and growth of these structures (Adams et al., 1992) until the dominant follicle reaches a certain size (i.e., 7.2 mm in buffalo) (Gimenes et al., 2011). Furthermore, a fertilizable oocyte is extruded from the follicle as a result of the pre-ovulatory surge of gonadotropins. Knowledge about follicular substances and the effects on follicle growth, however, as well as the diameter deviation mechanism of the largest from subordinate follicles during a wave of ovarian follicular development is limited. The roles of several key regulator genes involved in follicle growth have been identified (Evans et al., 2004; Tsuiko et al., 2016), however, the mechanism of the underlying global regulatory networks at the transcriptome level is still poorly understood. Transcriptome assembly knowledge has been applied to exploring transcriptional regulation mechanisms of many species, such as dairy cattle (Salleh et al., 2017), sheep (Jäger et al., 2011), and pigs (Du et al., 2014). In the present study, four independent cDNA libraries were constructed, respectively, representing four buffalo follicle growth stages by Illumina RNA-seq. The transcriptome changes during follicle growth were analyzed so as to reveal the molecular mechanism of follicular development. Furthermore, this is the initial attempt to report the dynamic development of buffalo follicles at transcriptome level. 2. Materials and methods 2.1. Material and sample collection In the present study, ovaries were collected at local abattoir in Guangxi, China, from non-pregnant hybrid buffalo. The ovaries with luteinized and large cystic follicles were detected by macroscopical examination and were not used in the study. Only the one largest follicle from each ovary was collected. The morphological characteristics were used to assess the developmental stages of follicles. Size is the most obvious morphological characteristics. In this study, size of follicles were designated into four categories: GC1 (< 5 mm in diameter), GC2 (5–8 mm), GC3 (8–12 mm) and GC4 (> 12 mm) with number in each category being 44 (Ginther et al., 2003; Pandey et al., 2011). Granulosa cells and follicular fluid were collected (Hatzirodos et al., 2014) from all follicles. Mixed granulosa cell samples representing each stage (n = 44) were frozen immediately in liquid nitrogen, and stored at −80 °C until RNAseq and RT-PCR. The follicular fluid samples of each stage were collected with three replicate samples for hormone measurement. 2.2. RNA extraction, cDNA synthesis and sequencing Total RNA was extracted from the granulosa cells of all follicles by using E.Z.N.A total RNA Kit I (R6834-02) according to the manufacturer’s instructions, and RNase-free DNase (Takara, China) was used to remove genomic DNA contamination. A total of 4 μg purified RNA per sample was sent to Hanzhou One Gene Co., where the samples were used to construct cDNA library. The library was constructed by using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) according to the manufacturer’s instructions. In addition, the Agilent Bioanalyzer 2100 system was used for assessing the quality of library. The, Illumina HiSeq™ 4000 was subsequently used for sequencing the amplified fragments and 150-bp paired-end reads were obtained by Hanzhou One Gene Co. (Hanzhou, China). 2.3. Assembly and annotation It is necessary to perform quality control on raw data, including removal of the raw data containing adapter, low quality reads with quality value Q < 20 base greater than 30%, and more than 5% unknown nucleotides. The software SOAPaligner/SOAP2 (Li et al., 2009) was used to compare the “clean” data of each sample with the reference gene of the species, allowing up to five base mismatches. The RPKM was used to obtain the relative gene expression amounts (Mortazavi et al., 2008). The analysis of the differentially expressed genes was performed by using the edgeR package in Bioconductor. The DEGs were screened among comparison groups according to the principle of FDR ≤ 0.05 and | log2Ratio | ≥1 so as to control the false discovery rate (Benjamini and Yekutieli 2001). All the DEGs were combined to perform K-means clustering, and the genes were grouped according to various amounts of expression. Venn diagrams and clustering heat maps in the present study were generated using Venn diagram and Pheatmap packages in R based on the extent of gene expression of DEGs for each comparison. Unigenes were compared to the NR library using the blast software, while GO annotation information for all genes was extracted from the Gene Ontology database (http://www.geneontology.org/). Then, GO and KEGG pathway analyses were performed (Moriya et al., 2007; Young et al., 2010). 182
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig 1. DEGs involved in follicle development process of buffalo. a. Venn diagram of differentially expressed genes during transition between stages. b. Number of DEGs during follicle growth. c, Expression profiles of the 30 common DEGs which differentially expressed in the four stages of follicular development. Different amounts of expression are indicated in different colors, red indicates expression up-regulation, and green indicates expression down-regulation; color from blue to red indicates expression from small to large, and the greater the color, the greater the amount of expression. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
183
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 2. K-means clustering map of DEGs.
Fig. 3. GO functions of the DEGs during stage GC1 to GC2. The red represents the up-regulated gene, and green represents the down-regulated genes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
2.4. Real-time PCR validation Single-stranded cDNA were obtained from every 1 μg RNA which was isolated respectively from each of four follicles growth stages. Primers for quantitative reverse transcription PCR (qPCR) were designed using Primer 5.0 software and synthesized by Sangon Biotech (Shanhai) Co. Ltd. Using a SYBR Green-based PCR assay, qPCR was performed with LightCycler 480 II (Roche). The qPCR 184
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 4. GO functions of the DEGs during stage GC2 to GC3. The red represents the up-regulated gene, and green represents the down-regulated gene. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. The GO functions of the DEGs during stage GC3 to GC4. The red represents the up-regulated gene, and green represents the down-regulated gene. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
mixture system and programs are similar to those described in the previous literature (Hatzirodos et al., 2014). The expression values of selected genes were determined by method 2−△△Ct after normalization against values of GADPH.
2.5. Measurements of various hormones Hormones are important for ovarian follicular development. The concentrations of FSH, LH, estrogen, and progesterone in follicular fluid were measured using a Bovine ELISA Kit (catalogue number: JYM0045B, JYM0055Bo, JYM0054Bo, JYM0056Bo, respectively) according to the manufacturer’s instructions. Each sample was subjected to three repeated measurements.
185
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 6. Top 20 data of KEGG pathway enrichment for GC1 compared with GC2.
2.6. Statistical analysis For determination of significant differences between values, the one-way ANOVA analysis was conducted with Tukey-Kramer’s post hoc test in SPSS (Version17.0; SPSS, Chicago, IL, USA). The values of hormones and the extent of gene expression of selected genes represent three biological replicates. 3. Results 3.1. Sequencing, assembly and annotation of buffalo transcriptome A total of 206,316,878 raw reads with the length of 150 bp were obtained from four different granulosa cell groups (Table 1). After eliminating low quality data by quality control, all clean reads (185,154,236) were used for further analysis (Table 1). Furthermore, “clean” data were assembled and then were aligned with reference genes, and more than 60% of a unique match in each sample was obtained. A total of 17,700 unigenes with average length of 2960 bp were annotated using NCBI Nr database (Table 1). 3.2. Identification of differentially expressed genes In summary, 13,672 DEGs were identified in the samples from the four stages of follicle development (Table 1). Based on the corrected P-value, 1076 DEGs (FDR < 0.05) were detected in the transitional process between stages (Table 1; Fig. 1b). Transition from stage GC3 to GC4 was when the greatest number of DEGs (691) was detected with 671 being up-regulated and 20 downregulated genes. The smallest difference in number of DEGs (299) with 215 up-regulated and 84 down-regulated genes was observed in the transition from stage GC1 to GC2. For the GC2 compared with GC3 comparison, a total of 303 DEGs were detected with 56 genes being up-regulated and 247 down-regulated. The calculation result of the number of overlapping DEGs in the whole transitional process from GC1 to GC4 was presented in a Venn diagram (Fig. 1a). In detail, there were 30 shared DEGs during the four 186
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 7. Top 20 data of KEGG pathway enrichment for GC2 compared GC3.
stages of follicular development. There were 161, 133 and 595 DEGs for comparisons of GC1 with GC2, GC2 with GC3 and GC3 with GC4, respectively. 3.3. Cluster analysis of DEGs To characterize the major trends and important genes during the follicular development process, the DEGs were distributed to nine clusters (Fig. 2) by the K-means analysis. The extent of expression of DEGs of cluster 2, 5, 7 and 8 fluctuated without marked changes during the first three stages of follicular development, whereas in the last stage there was a marked increase in expression. The amount of expression of DEGs of cluster 1 increased continuously during the entire follicular development process, while the amount of expression of DEG of cluster 4 tended to decrease during the four stages of follicular development. The amount of expression of DEGs of cluster 7 was stage-specific in stage GC4. Similarly, that of cluster 3 was also stage-specific with the minimum gene expression in stage GC1. A total of 30 DEGs were shared during four follicular development stages, and the expression patterns of these shared genes are depicted in Fig. 1c at different stages. All of these genes were up-regulated in stages GC2 and GC4, and down-regulated in dominant follicular formation stage GC3 with an exception of one gene. 3.4. GO enrichment of DEGs All DEGs in four stages were subjected to GO analysis. The most significantly enriched GO terms (including biological process, cellular component, and molecular function) and the number of the up- and down-regulated genes during the four transitions were depicted in Figs. 3 through 5 . During the transition from stage GC1 to GC2 (Fig. 3), the largest terms within the biological process were protein targeting to membrane, viral gene expression, multi-organism metabolic process, single-organism localization, nuclear-transcribed mRNA catabolic process, and RNA catabolic process. Within the cellular component category, ribosomal subunit, ribosome, and 187
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 8. Top 20 data of KEGG pathway enrichment for GC3 vs. GC4.
ribonucleoprotein complex were the most common terms. Within the molecular function category, the most highly represented term was structural molecule activity. The process of determining the dominant follicles mainly emerged in the transitional period from GC2 to GC3 (Fig. 4). The most significantly enriched biological processes were protein localization to membrane, protein targeting to membrane, single-organism localization, membrane organization, nuclear-transcribed mRNA catabolic process, and viral life cycle. The most significant cellular component involved ribosome and ribosomal subunit. Structural molecule activity and RNA binding were most significantly enriched within molecular function. DEGs in the transition from GC3 to GC4 may have a pivotal role in follicle maturation of buffalo (Fig. 5). The biological process was mainly enriched in the immune system process, response to external stimulus, defense response, leukocyte activation, cell activation, leukocyte migration, immune effector process, regulation of immune system process, and lymphocyte activation. The most highly represented terms were extracellular region, plasma membrane, cell periphery, MHC protein complex, vacuole, and membrane within the cellular component. Receptor binding, cytokine activity, and cytokine receptor binding were the most abundant terms within the molecular function.
3.5. KEGG pathway enrichment analysis of DEGs In vivo, the biological process of follicle growth is coordinated by multiple genes. The 174, 176 and 203 DEGs, respectively, from comparisons of GC1 with GC2, GC2 with GC3 and GC3 with GC4 were enriched in the KEGG pathway, suggesting a biological function of DEGs for follicle growth. The greatest 20 data outputs for significant KEGG pathway enrichment are depicted in Figs. 6 through 8 . The largest number of DEGs (36) in GC1 compared with GC2 occurred in ribosome pathway (ko03010). The significant KEGG pathways enrichment in GC2 compared with GC3 occurred in the ribosome pathway (ko03010) with 34 DEGs (Fig. 7), however, the largest number of DEGs (37) was in metabolic pathways (ko01100). The largest number of DEGs (42) in GC3 compared with GC4 was observed in the chemokine signaling pathway (ko04062), followed by osteoclast differentiation (ko04380), cytokine–cytokine receptor interaction (ko04060), phagosome (ko04145), tuberculosis (ko05152), natural killer cell mediated cytotoxicity 188
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 9. Measurement of hormone content in follicular fluid during follicle growth. a. Progesterone. b. Luteinizing hormone (LH). c. Follicle stimulating hormone (FSH). d. Estradiol.
(ko04650), and B cell receptor signaling pathway (ko04662). 3.6. Endogenous hormone measurements The endogenous hormones, including FSH, LH, estradiol and progesterone, during the four stages of follicle growth were analyzed (Fig. 9). The concentration of progesterone and FSH rapidly decreased from GC1 to GC2 and then increased. Meanwhile, the concentration of LH and estradiol increased from GC1 to GC3 and then decreased. Both the concentration of progesterone and estradiol were greatest at stage GC4. During the follicular development process, the concentration of FSH significantly decreased from 48.84 (GC1) to 37.19 ng/ml(GC2). The concentration of LH, however, significantly increased from 4.31 (GC2) to 5.24 ng/ml(GC3), and the concentration of progesterone significantly increased from 9.48 (GC2) to 12.74 ng/ml (GC4). 3.7. Confirmation of selected genes expression by using qPCR Six unigenes were selected randomly to validate the RNA-seq data. The expression values were normalized against values of GADPH. The expression values of these selected genes in four stages obtained by qPCR were basically consistent with the RPKM values obtained by RNA-seq (Fig. 10). 4. Discussion Follicle maturation and ovulation is the key to ensuring conception and improving the fertility of buffalo. Extensive information has been obtained about the regulation of buffalo’s follicular development by genomics (Purohit et al., 2003) and hormones (Presicce et al., 2004), whereas little is known about transcriptomics. In the present study, the general features of transcription were evaluated with high quality sequencing and assembly data. Based on the transcript profiles in four stages of ovarian follicular development, the candidate genes regulating the follicle growth in buffalo were investigated, providing the insight into the complex mechanisms underlying follicular development of buffalo. A total of 1076 significant DEGs (FDR < 0.05) were detected through comparisons. Interestingly, there were a large number of stage-specific DEGs in the GC4 stage, suggesting that these DEGs participated in numerous and complex biological processes at this stage of follicular development. It is well established that the last stage of follicle growth has an important role in follicular atresia or ovulation, and the timing of ovulation is related to the pre-ovulatory LH peak in yaks (Sarkar and Prakash, 2005). Moreover, the vascularity of the dominant follicle was greater than that of the subordinate follicle, possibly determining whether there was ovulation from or atresia of the follicle (Acosta, 2007). It is assumed that the DEGs identified in stage GC4 are directly related to the fertility of buffalo. In the early stage of follicular growth (< 8 mm, diameter), the molecular function of DEGs mainly includes RNA-binding and 189
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 10. Comparison of amounts of expression of six genes obtained by qPCR analysis and by RNA-seq (RPKM values).
structural molecule activity. In addition, the biological processes of DEGs during this period mainly involve membrane organization, nuclear-transcribed mRNA catabolic process, and RNA catabolic process. Findings in the present study are consistent with the observation of the large amount of proliferation and differentiation of granulosa cells in follicles during this period. Notably, KRT18 is one of the genes that is enriched in a number of biological functions, and the encoded protein is the trophoblast cell marker factor (Hou et al., 2015). The expression of the KRT18 gene was reduced by RNA interference techniques, resulting in a lesser developmental competence at the blastocyst stage in cattle (Goossens et al., 2010). By comparing the transcriptome profiles of buffalo embryos that had normal or retarded growth on Day 25 after mating, the expression of the KRT18 gene was significantly up-regulated with the normally growing embryos, suggesting KRT18 gene expression is associated with endometrial attachment and establishment of pregnancy (Strazzullo et al., 2014). In the present study, the expression of the KRT18 gene was down-regulated when comparing the GC2 GC3 stages of development, which indicates that KRT18 is a candidate gene affecting formation of dominant follicles. The immune system has an important role in the last stages of follicular maturation and ovulation. In the present study, there was significantly enriched signal pathways including allograft rejection, chemokine signaling, natural killer cell mediated cytotoxicity, phagosome, and antigen processing and presentation, indicating that immune-related molecules have an important role in dominant follicles developing the ovulation capacity. Cattle are the closest species to buffalo in evolutionary relationship. The comparison of transcriptome profiles of granulosa cells between small (< 5 mm) and large (> 10) follicles of cattle were performed using Affymetrix microarrays (Hatzirodos et al., 2014), and the result suggested that the processes of immune signaling were enriched in large follicles to the greatest extent. When ovulation occurs, this is indicative of the final stage of follicle development. This process was regarded as one type of localized physiological inflammatory reaction. The greatest concentrations of inflammatory substances have been found in the follicle during the peri-ovulatory phase of follicular development (Brannstrom and Enskog 2002). Previous studies demonstrated that inflammatory substances including leukocytes (Kh et al., 2000), macrophages (Takaya et al., 1997), neutrophils (Ujioka et al., 1998), and mast cells (Gaytan et al., 1991) have an important role in the ovulatory process to facilitate ovulation. The knowledge about the exact roles of these cells in ovulation process is, however, limited. Interestingly, there are many DEGs that are enriched in the last stage of the allograft rejection signaling process. A possible interpretation of results from the present study is during stage GC4 of follicle development there is that there is preparation occurring for the protecting of growth of oocytes after ovulation and movement of the oocytes into the fallopian tubes. Of the previously-
190
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
mentioned DEGs, only TNF, GZMB and IL6 have been reported to be associated with female fertility in previous studies. The same situation also exists in natural killer cell mediated cytotoxicity and the phagosome signaling process. It is noteworthy that PAQR7, is enriched in the natural killer cell mediated cytotoxicity signaling process, and has been reported to be related to the induction of oocyte maturation in fish (Shi et al., 2016). Another interesting that the products of the expression of the gene, CTSS, which is enriched in the phagosome signaling process, have been reported to participate in the process of cumulus cell apoptosis and oocyte development (Boruszewska et al., 2014). 5. Conclusion A total of four independent cDNA libraries, respectively, for four stages of ovarian follicle growth, namely, GC1, GC2, GC3 and GC4, were obtained by using RNA-seq. Numerous DEGs were identified during follicle growth. The stage-specifically expressed genes in follicle growth were characterized by GO and KEGG enrichment. The present study revealed overall changes in gene expression during follicle growth. It also indicated that follicle maturation and ovulation might be subject to immune control which involves allograft rejection, chemokine signaling pathway, natural killer cell mediated cytotoxicity, phagosome, and antigen processing and presentation. Competing interests The authors declare that they have no competing interests. Acknowledgements This study was financially supported by the National Natural Science Foundation of China (No. 31772602 and No. 31772604). References Acosta, T.J., 2007. Studies of follicular vascularity associated with follicle selection and ovulation in cattle. J. Reprod. Dev. 53, 39–44. Adams, G.P., Matteri, R.L., Kastelic, J.P., Ko, J.C., Ginther, O.J., 1992. Association between surges of follicle-stimulating hormone and the emergence of follicular waves in heifers. J. Reprod. Fertil. 94, 177–188. Benjamini, Y., Yekutieli, D., 2001. The control of the false discovery rate in multiple testing under dependency. Ann. Stat. 29, 1165–1188. Boruszewska, D., Torres, A.C., Kowalczyk-Zieba, I., Diniz, P., Batista, M., Lopes-da-Costa, L., Woclawek-Potocka, I., 2014. The effect of lysophosphatidic acid during in vitro maturation of bovine oocytes: embryonic development and mRNA abundances of genes involved in apoptosis and oocyte competence. Mediators Inflamm. 670670. Brannstrom, M., Enskog, A., 2002. Leukocyte networks and ovulation. J. Reprod. Immunol. 57, 47–60. Du, Z.Q., Eisley, C.J., Onteru, S.K., Madsen, O., Groenen, M.A.M., Ross, J.W., Rothschild, M.F., 2014. Identification of species-specific novel transcripts in pig reproductive tissues using RNA-seq. Anim. Genet. 45, 198–204. Evans, A.C., Komar, C.M., Wandji, S.A., Fortune, J.E., 1997. Changes in androgen secretion and luteinizing hormone pulse amplitude are associated with the recruitment and growth of ovarian follicles during the luteal phase of the bovine estrous cycle. Biol. Reprod. 57, 394–401. Evans, A.C.O., Ireland, J.L.H., Winn, M.E., Lonergan, P., Smith, G.W., Coussens, P.M., Ireland, J.J., 2004. Identification of genes involved in apoptosis and dominant follicle development during follicular waves in cattle1. Biol. Reprod. 70, 1475–1484. Fortune, J.E., Sirois, J., Turzillo, A.M., Lavoir, M., 1991. Follicle selection in domestic ruminants. J. Reprod. Fertil. Suppl. 43, 187–198. Gaytan, F., Aceitero, J., Bellido, C., Sánchez-Criado, J.E., Aguilar, E., 1991. Estrous cycle-related changes in mast cell numbers in several ovarian compartments in the rat. Biol. Reprod. 45, 27–33. Gimenes, L.U., Carvalho, N.A., Sa Filho, M.F., Vannucci, F.S., Torres-Junior, J.R., Ayres, H., Ferreira, R.M., Trinca, L.A., Sartorelli, E.S., Barros, C.M., Beltran, M.P., Nogueira, G.P., Mapletoft, R.J., Baruselli, P.S., 2011. Ultrasonographic and endocrine aspects of follicle deviation: and acquisition of ovulatory capacity in buffalo (Bubalus bubalis) heifers. Anim. Reprod. Sci. 123, 175–179. Ginther, O.J., Beg, M.A., Bergfelt, D.R., Donadeu, F.X., Kot, K., 2001. Follicle selection in monovular species. Biol. Reprod. 65 638-147. Ginther, O.J., Beg, M.A., Donadeu, F.X., Bergfelt, D.R., 2003. Mechanism of follicle deviation in monovular farm species. Anim. Reprod. Sci. 78 239-157. Goossens, K., Tesfaye, D., Rings, F., Schellander, K., Holker, M., Van Poucke, M., Van Zeveren, A., Lemahieu, I., Van Soom, A., Peelman, L.J., 2010. Suppression of keratin 18 gene expression in bovine blastocysts by RNA interference. Reprod. Fertil. Dev. 22, 395–404. Hatzirodos, N., Irving-Rodgers, H.F., Hummitzsch, K., Harland, M.L., Morris, S.E., Rodgers, R.J., 2014. Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes. BMC Genomics 15, 24. Hou, D., Su, M., Li, X., Li, Z., Yun, T., Zhao, Y., Zhang, M., Zhao, L., Li, R., Yu, H., Li, X., 2015. The efficient derivation of trophoblast cells from porcine in vitro fertilized and parthenogenetic blastocysts and culture with ROCK inhibitor Y-27632. PLoS One 10, e0142442. Jäger, M., Ott, C.E., Grünhagen, J., Hecht, J., Schell, H., Mundlos, S., Duda, G.N., Robinson, P.N., Lienau, Jasmin, 2011. Composite transcriptome assembly of RNA-seq data in a sheep model for delayed bone healing. BMC Genomics 12, 158. Jain, A., Jain, T., Kumar, P., Kumar, M., De, S., Gohain, M., Kumar, R., Datta, T.K., 2016. Follicle-stimulating hormone-induced rescue of cumulus cell apoptosis and enhanced development ability of buffalo oocytes. Domest. Anim. Endocrinol. 55, 74–82. Kh, V.D.H., Maddocks, S., Woodhouse, C.M., Van, R.N., Robertson, S.A., Norman, R.J., 2000. Intrabursal injection of clodronate liposomes causes macrophage depletion and inhibits ovulation in the mouse ovary. Biol. Reprod. 62, 1059–1066. Li, R., Yu, C., Li, Y., Lam, T.W., Yiu, S.M., Kristiansen, K., Wang, J., 2009. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25, 1966–1977. Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A.C., Kanehisa, M., 2007. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 35, W182–185. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. Pandey, A.K., Dhaliwal, G.S., Ghuman, S.P., Agarwal, S.K., 2011. Impact of pre-ovulatory follicle diameter on plasma estradiol, subsequent luteal profiles and conception rate in buffalo (Bubalus bubalis). Anim. Reprod. Sci. 123, 169–174. Presicce, G.A., Senatore, E.M., Bella, A., De Santis, G., Barile, V.L., De Mauro, G.J., Terzano, G.M., Stecco, R., Parmeggiani, A., 2004. Ovarian follicular dynamics and hormonal profiles in heifer and mixed-parity Mediterranean Italian buffaloes (Bubalus bubalis) following an estrus synchronization protocol. Theriogenology 61, 1343–1355. Purohit, G.N., Duggal, G.P., Dadarwal, D., Kumar, D., Yadav, R.C., Vyas, S., 2003. Reproductive biotechnologies for improvement of buffalo: the current status. AsianAustralas. J. Anim. Sci. 16, 1071–1086.
191
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Salleh, M.S., Mazzoni, G., Hoglund, J.K., Olijhoek, D.W., Lund, P., Lovendahl, P., Kadarmideen, H.N., 2017. RNA-seq transcriptomics and pathway analyses reveal potential regulatory genes and molecular mechanisms in high- and low-residual feed intake in Nordic dairy cattle. BMC Genomics 18, 258. Sarkar, M., Prakash, B.S., 2005. Timing of ovulation in relation to onset of estrus and LH peak in yak (Poephagus grunniens L.). Anim. Reprod. Sci. 86, 353–362. Shi, B., Liu, X., Thomas, P., Pang, Y., Xu, Y., Li, X., Li, X., 2016. Identification and characterization of a progestin and adipoQ receptor (PAQR) structurally related to Paqr7 in the ovary of Cynoglossus semilaevis and its potential role in regulating oocyte maturation. Gen. Comp. Endocrinol. 237, 109–120. Strazzullo, M., Gasparrini, B., Neglia, G., Balestrieri, M.L., Francioso, R., Rossetti, C., Nassa, G., De Filippo, M.R., Weisz, A., Di Francesco, S., Vecchio, D., D'Esposito, M., D'Occhio, M.J., Zicarelli, L., Campanile, G., 2014. Global transcriptome profiles of Italian Mediterranean buffalo embryos with normal and retarded growth. PLoS One 9, e90027. Takaya, R., Fukaya, T., Sasano, H., Suzuki, T., Tamura, M., Yajima, A., 1997. Macrophages in normal cycling human ovaries; immunohistochemical localization and characterization. Hum. Reprod. 12, 1508–1512. Tsuiko, O., Noukas, M., Zilina, O., Hensen, K., Tapanainen, J.S., Magi, R., Kals, M., Kivistik, P.A., Haller-Kikkatalo, K., Salumets, A., Kurg, A., 2016. Copy number variation analysis detects novel candidate genes involved in follicular growth and oocyte maturation in a cohort of premature ovarian failure cases. Hum. Reprod. 31, 1913–1925. Ujioka, T., Matsukawa, A., Tanaka, N., Matsuura, K., Yoshinaga, M., Okamura, H., 1998. Interleukin-8 as an essential factor in the human chorionic gonadotropininduced rabbit ovulatory process: interleukin-8 induces neutrophil accumulation and activation in ovulation. Biol. Reprod. 58, 526–530. Yeung, C.K., Wang, G., Yao, Y., Liang, J., Tenny Chung, C.Y., Chuai, M., Lee, K.K., Yang, X., 2017. BRE modulates granulosa cell death to affect ovarian follicle development and atresia in the mouse. Cell. Death. Dis. 8, e2697. Young, M.D., Wakefield, M.J., Smyth, G.K., Oshlack, A., 2010. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, R14.
192
Contents lists available at ScienceDirect
Animal Reproduction Science journal homepage: www.elsevier.com/locate/anireprosci
Transcriptome studies of granulosa cells at different stages of ovarian follicular development in buffalo J. Lia,b, Z. Lia,b, S. Liua,b, R. Ziaa, A. Lianga,b, L. Yanga,b, a b
MARK
⁎
College of Animal Science and Technology, Huazhong Agricultural University, Wuhan, Hubei, China Hubei Province’s Engineering Research Center in Buffalo Breeding & Products, China
AR TI CLE I NF O
AB S T R A CT
Keywords: Buffalo Differentially expressed genes Follicle growth RNA-seq Transcriptome
The normal maturation and ovulation from ovarian follicles is important in ensuring conception and improving fertility of buffalo. The molecular regulation mechanism of buffalo follicles growth, however, remains unknown. This study analyzed the gene expression profiles associated with buffalo ovarian follicle growth. According to the analysis of RNA sequencing, 17,700 unigenes and 13,672 differentially expressed genes (DEGs) were detected. A total of 30 common DEGs were identified during four stages of follicle growth, and the expression patterns are basically synchronized, suggesting the products as a result of expressions of these genes may cooperate to regulate follicular development. Furthermore, GO and KEGG enrichment analyses revealed that the majority of DEGs in early stage of follicular growth were enriched in ribosomal and oxidative phosphorylation signaling pathways, and the expression patterns of these DEGs are basically up-regulated at the beginning of follicular growth (< 8 mm, diameter), and then downregulated (8–12 mm) in the following stages of follicular development. The pathway of immune signaling, including allograft rejection, chemokine signaling pathway, natural killer cell mediated cytotoxicity, phagosome, and antigen processing and presentation, was significantly enriched in the last stage of follicular development (> 12 mm), which indicates that the immune system has an important role in the last stage of follicular maturation and ovulation. This study provided a gene expression profile of buffalo follicle growth, and provided an insight into biological processes associated with molecular regulation of ovarian follicle growth.
1. Introduction Buffalo, one of the most important economic livestock worldwide, is known for its resistance to stress and its milk product with high milk fat and protein contents. The number of buffalo stock worldwide is estimated to be approximately 172 million (FAO: http://faostat.fao.org/). Unlike swine, sheep, and goats, however, where multiple follicles undergo maturation and from which ovulation occurs, buffalo is a monovular species, and thus, it normally has only one follicle develop into mature follicle from which ovulation occurs (Ginther et al., 2001) resulting in a lesser fecundity in buffalo. This is an important factor restricting the development of the buffalo breeding industry. Ovaries are critical organs for follicular development and ovulation. The follicle is the structural and functional units of the ovaries for which normal development during the developmental stages is important for lifelong fertility of buffalo. Follicle development results with granulosa cell growth, development, differentiation, and apoptosis in follicles, which is closely related to growth of oocyte and follicular atresia (Jain et al., 2016; Yeung et al., 2017). Two or three waves of ovarian follicular development emerge ⁎
Corresponding author at: No.1 Shizishan Street, Wuhan, Hubei Province, 430070, China. E-mail address: [email protected] (L. Yang).
http://dx.doi.org/10.1016/j.anireprosci.2017.11.004 Received 7 August 2017; Received in revised form 17 October 2017; Accepted 3 November 2017 Available online 06 November 2017 0378-4320/ © 2017 Elsevier B.V. All rights reserved.
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Table 1 Results of RNA-seq, assembly, and annotion. Stage
Raw reads
Clean reads
Assembled reads
Unigenes
DEGs
SDEGs
GC1 GC2 GC3 GC4 Total
50,760,648 50,760,170 50,760,644 54,035,416 206,316,878
46,219,792 46,675,486 46,192,560 46,066,398 185,154,236(89.74%)
29,117,053 29,120,413 30,062,845 28,098,004 116,398,315(62.87%)
16,249 16,339 16,167 16,349 17,700
– 13160 13048 13228 13,672(77.24%)
– 299 303 691 1076(7.87%)
DEG, differentially expressed genes based on the previous stage; SDEG, significantly differentially expressed genes based on the previous stage with FDR ≤ 0.05 and |log2Ratio| ≥ 1.
during each estrous cycle (Fortune et al., 1991), and the emergence of each follicular wave is synchronous with the surge of a folliclestimulating hormone (FSH) (Evans et al., 1997). The FSH has a pivotal role during the early stage of follicular development, ensuring follicle recruitment and growth of these structures (Adams et al., 1992) until the dominant follicle reaches a certain size (i.e., 7.2 mm in buffalo) (Gimenes et al., 2011). Furthermore, a fertilizable oocyte is extruded from the follicle as a result of the pre-ovulatory surge of gonadotropins. Knowledge about follicular substances and the effects on follicle growth, however, as well as the diameter deviation mechanism of the largest from subordinate follicles during a wave of ovarian follicular development is limited. The roles of several key regulator genes involved in follicle growth have been identified (Evans et al., 2004; Tsuiko et al., 2016), however, the mechanism of the underlying global regulatory networks at the transcriptome level is still poorly understood. Transcriptome assembly knowledge has been applied to exploring transcriptional regulation mechanisms of many species, such as dairy cattle (Salleh et al., 2017), sheep (Jäger et al., 2011), and pigs (Du et al., 2014). In the present study, four independent cDNA libraries were constructed, respectively, representing four buffalo follicle growth stages by Illumina RNA-seq. The transcriptome changes during follicle growth were analyzed so as to reveal the molecular mechanism of follicular development. Furthermore, this is the initial attempt to report the dynamic development of buffalo follicles at transcriptome level. 2. Materials and methods 2.1. Material and sample collection In the present study, ovaries were collected at local abattoir in Guangxi, China, from non-pregnant hybrid buffalo. The ovaries with luteinized and large cystic follicles were detected by macroscopical examination and were not used in the study. Only the one largest follicle from each ovary was collected. The morphological characteristics were used to assess the developmental stages of follicles. Size is the most obvious morphological characteristics. In this study, size of follicles were designated into four categories: GC1 (< 5 mm in diameter), GC2 (5–8 mm), GC3 (8–12 mm) and GC4 (> 12 mm) with number in each category being 44 (Ginther et al., 2003; Pandey et al., 2011). Granulosa cells and follicular fluid were collected (Hatzirodos et al., 2014) from all follicles. Mixed granulosa cell samples representing each stage (n = 44) were frozen immediately in liquid nitrogen, and stored at −80 °C until RNAseq and RT-PCR. The follicular fluid samples of each stage were collected with three replicate samples for hormone measurement. 2.2. RNA extraction, cDNA synthesis and sequencing Total RNA was extracted from the granulosa cells of all follicles by using E.Z.N.A total RNA Kit I (R6834-02) according to the manufacturer’s instructions, and RNase-free DNase (Takara, China) was used to remove genomic DNA contamination. A total of 4 μg purified RNA per sample was sent to Hanzhou One Gene Co., where the samples were used to construct cDNA library. The library was constructed by using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) according to the manufacturer’s instructions. In addition, the Agilent Bioanalyzer 2100 system was used for assessing the quality of library. The, Illumina HiSeq™ 4000 was subsequently used for sequencing the amplified fragments and 150-bp paired-end reads were obtained by Hanzhou One Gene Co. (Hanzhou, China). 2.3. Assembly and annotation It is necessary to perform quality control on raw data, including removal of the raw data containing adapter, low quality reads with quality value Q < 20 base greater than 30%, and more than 5% unknown nucleotides. The software SOAPaligner/SOAP2 (Li et al., 2009) was used to compare the “clean” data of each sample with the reference gene of the species, allowing up to five base mismatches. The RPKM was used to obtain the relative gene expression amounts (Mortazavi et al., 2008). The analysis of the differentially expressed genes was performed by using the edgeR package in Bioconductor. The DEGs were screened among comparison groups according to the principle of FDR ≤ 0.05 and | log2Ratio | ≥1 so as to control the false discovery rate (Benjamini and Yekutieli 2001). All the DEGs were combined to perform K-means clustering, and the genes were grouped according to various amounts of expression. Venn diagrams and clustering heat maps in the present study were generated using Venn diagram and Pheatmap packages in R based on the extent of gene expression of DEGs for each comparison. Unigenes were compared to the NR library using the blast software, while GO annotation information for all genes was extracted from the Gene Ontology database (http://www.geneontology.org/). Then, GO and KEGG pathway analyses were performed (Moriya et al., 2007; Young et al., 2010). 182
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig 1. DEGs involved in follicle development process of buffalo. a. Venn diagram of differentially expressed genes during transition between stages. b. Number of DEGs during follicle growth. c, Expression profiles of the 30 common DEGs which differentially expressed in the four stages of follicular development. Different amounts of expression are indicated in different colors, red indicates expression up-regulation, and green indicates expression down-regulation; color from blue to red indicates expression from small to large, and the greater the color, the greater the amount of expression. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
183
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 2. K-means clustering map of DEGs.
Fig. 3. GO functions of the DEGs during stage GC1 to GC2. The red represents the up-regulated gene, and green represents the down-regulated genes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
2.4. Real-time PCR validation Single-stranded cDNA were obtained from every 1 μg RNA which was isolated respectively from each of four follicles growth stages. Primers for quantitative reverse transcription PCR (qPCR) were designed using Primer 5.0 software and synthesized by Sangon Biotech (Shanhai) Co. Ltd. Using a SYBR Green-based PCR assay, qPCR was performed with LightCycler 480 II (Roche). The qPCR 184
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 4. GO functions of the DEGs during stage GC2 to GC3. The red represents the up-regulated gene, and green represents the down-regulated gene. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. The GO functions of the DEGs during stage GC3 to GC4. The red represents the up-regulated gene, and green represents the down-regulated gene. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
mixture system and programs are similar to those described in the previous literature (Hatzirodos et al., 2014). The expression values of selected genes were determined by method 2−△△Ct after normalization against values of GADPH.
2.5. Measurements of various hormones Hormones are important for ovarian follicular development. The concentrations of FSH, LH, estrogen, and progesterone in follicular fluid were measured using a Bovine ELISA Kit (catalogue number: JYM0045B, JYM0055Bo, JYM0054Bo, JYM0056Bo, respectively) according to the manufacturer’s instructions. Each sample was subjected to three repeated measurements.
185
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 6. Top 20 data of KEGG pathway enrichment for GC1 compared with GC2.
2.6. Statistical analysis For determination of significant differences between values, the one-way ANOVA analysis was conducted with Tukey-Kramer’s post hoc test in SPSS (Version17.0; SPSS, Chicago, IL, USA). The values of hormones and the extent of gene expression of selected genes represent three biological replicates. 3. Results 3.1. Sequencing, assembly and annotation of buffalo transcriptome A total of 206,316,878 raw reads with the length of 150 bp were obtained from four different granulosa cell groups (Table 1). After eliminating low quality data by quality control, all clean reads (185,154,236) were used for further analysis (Table 1). Furthermore, “clean” data were assembled and then were aligned with reference genes, and more than 60% of a unique match in each sample was obtained. A total of 17,700 unigenes with average length of 2960 bp were annotated using NCBI Nr database (Table 1). 3.2. Identification of differentially expressed genes In summary, 13,672 DEGs were identified in the samples from the four stages of follicle development (Table 1). Based on the corrected P-value, 1076 DEGs (FDR < 0.05) were detected in the transitional process between stages (Table 1; Fig. 1b). Transition from stage GC3 to GC4 was when the greatest number of DEGs (691) was detected with 671 being up-regulated and 20 downregulated genes. The smallest difference in number of DEGs (299) with 215 up-regulated and 84 down-regulated genes was observed in the transition from stage GC1 to GC2. For the GC2 compared with GC3 comparison, a total of 303 DEGs were detected with 56 genes being up-regulated and 247 down-regulated. The calculation result of the number of overlapping DEGs in the whole transitional process from GC1 to GC4 was presented in a Venn diagram (Fig. 1a). In detail, there were 30 shared DEGs during the four 186
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 7. Top 20 data of KEGG pathway enrichment for GC2 compared GC3.
stages of follicular development. There were 161, 133 and 595 DEGs for comparisons of GC1 with GC2, GC2 with GC3 and GC3 with GC4, respectively. 3.3. Cluster analysis of DEGs To characterize the major trends and important genes during the follicular development process, the DEGs were distributed to nine clusters (Fig. 2) by the K-means analysis. The extent of expression of DEGs of cluster 2, 5, 7 and 8 fluctuated without marked changes during the first three stages of follicular development, whereas in the last stage there was a marked increase in expression. The amount of expression of DEGs of cluster 1 increased continuously during the entire follicular development process, while the amount of expression of DEG of cluster 4 tended to decrease during the four stages of follicular development. The amount of expression of DEGs of cluster 7 was stage-specific in stage GC4. Similarly, that of cluster 3 was also stage-specific with the minimum gene expression in stage GC1. A total of 30 DEGs were shared during four follicular development stages, and the expression patterns of these shared genes are depicted in Fig. 1c at different stages. All of these genes were up-regulated in stages GC2 and GC4, and down-regulated in dominant follicular formation stage GC3 with an exception of one gene. 3.4. GO enrichment of DEGs All DEGs in four stages were subjected to GO analysis. The most significantly enriched GO terms (including biological process, cellular component, and molecular function) and the number of the up- and down-regulated genes during the four transitions were depicted in Figs. 3 through 5 . During the transition from stage GC1 to GC2 (Fig. 3), the largest terms within the biological process were protein targeting to membrane, viral gene expression, multi-organism metabolic process, single-organism localization, nuclear-transcribed mRNA catabolic process, and RNA catabolic process. Within the cellular component category, ribosomal subunit, ribosome, and 187
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 8. Top 20 data of KEGG pathway enrichment for GC3 vs. GC4.
ribonucleoprotein complex were the most common terms. Within the molecular function category, the most highly represented term was structural molecule activity. The process of determining the dominant follicles mainly emerged in the transitional period from GC2 to GC3 (Fig. 4). The most significantly enriched biological processes were protein localization to membrane, protein targeting to membrane, single-organism localization, membrane organization, nuclear-transcribed mRNA catabolic process, and viral life cycle. The most significant cellular component involved ribosome and ribosomal subunit. Structural molecule activity and RNA binding were most significantly enriched within molecular function. DEGs in the transition from GC3 to GC4 may have a pivotal role in follicle maturation of buffalo (Fig. 5). The biological process was mainly enriched in the immune system process, response to external stimulus, defense response, leukocyte activation, cell activation, leukocyte migration, immune effector process, regulation of immune system process, and lymphocyte activation. The most highly represented terms were extracellular region, plasma membrane, cell periphery, MHC protein complex, vacuole, and membrane within the cellular component. Receptor binding, cytokine activity, and cytokine receptor binding were the most abundant terms within the molecular function.
3.5. KEGG pathway enrichment analysis of DEGs In vivo, the biological process of follicle growth is coordinated by multiple genes. The 174, 176 and 203 DEGs, respectively, from comparisons of GC1 with GC2, GC2 with GC3 and GC3 with GC4 were enriched in the KEGG pathway, suggesting a biological function of DEGs for follicle growth. The greatest 20 data outputs for significant KEGG pathway enrichment are depicted in Figs. 6 through 8 . The largest number of DEGs (36) in GC1 compared with GC2 occurred in ribosome pathway (ko03010). The significant KEGG pathways enrichment in GC2 compared with GC3 occurred in the ribosome pathway (ko03010) with 34 DEGs (Fig. 7), however, the largest number of DEGs (37) was in metabolic pathways (ko01100). The largest number of DEGs (42) in GC3 compared with GC4 was observed in the chemokine signaling pathway (ko04062), followed by osteoclast differentiation (ko04380), cytokine–cytokine receptor interaction (ko04060), phagosome (ko04145), tuberculosis (ko05152), natural killer cell mediated cytotoxicity 188
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 9. Measurement of hormone content in follicular fluid during follicle growth. a. Progesterone. b. Luteinizing hormone (LH). c. Follicle stimulating hormone (FSH). d. Estradiol.
(ko04650), and B cell receptor signaling pathway (ko04662). 3.6. Endogenous hormone measurements The endogenous hormones, including FSH, LH, estradiol and progesterone, during the four stages of follicle growth were analyzed (Fig. 9). The concentration of progesterone and FSH rapidly decreased from GC1 to GC2 and then increased. Meanwhile, the concentration of LH and estradiol increased from GC1 to GC3 and then decreased. Both the concentration of progesterone and estradiol were greatest at stage GC4. During the follicular development process, the concentration of FSH significantly decreased from 48.84 (GC1) to 37.19 ng/ml(GC2). The concentration of LH, however, significantly increased from 4.31 (GC2) to 5.24 ng/ml(GC3), and the concentration of progesterone significantly increased from 9.48 (GC2) to 12.74 ng/ml (GC4). 3.7. Confirmation of selected genes expression by using qPCR Six unigenes were selected randomly to validate the RNA-seq data. The expression values were normalized against values of GADPH. The expression values of these selected genes in four stages obtained by qPCR were basically consistent with the RPKM values obtained by RNA-seq (Fig. 10). 4. Discussion Follicle maturation and ovulation is the key to ensuring conception and improving the fertility of buffalo. Extensive information has been obtained about the regulation of buffalo’s follicular development by genomics (Purohit et al., 2003) and hormones (Presicce et al., 2004), whereas little is known about transcriptomics. In the present study, the general features of transcription were evaluated with high quality sequencing and assembly data. Based on the transcript profiles in four stages of ovarian follicular development, the candidate genes regulating the follicle growth in buffalo were investigated, providing the insight into the complex mechanisms underlying follicular development of buffalo. A total of 1076 significant DEGs (FDR < 0.05) were detected through comparisons. Interestingly, there were a large number of stage-specific DEGs in the GC4 stage, suggesting that these DEGs participated in numerous and complex biological processes at this stage of follicular development. It is well established that the last stage of follicle growth has an important role in follicular atresia or ovulation, and the timing of ovulation is related to the pre-ovulatory LH peak in yaks (Sarkar and Prakash, 2005). Moreover, the vascularity of the dominant follicle was greater than that of the subordinate follicle, possibly determining whether there was ovulation from or atresia of the follicle (Acosta, 2007). It is assumed that the DEGs identified in stage GC4 are directly related to the fertility of buffalo. In the early stage of follicular growth (< 8 mm, diameter), the molecular function of DEGs mainly includes RNA-binding and 189
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Fig. 10. Comparison of amounts of expression of six genes obtained by qPCR analysis and by RNA-seq (RPKM values).
structural molecule activity. In addition, the biological processes of DEGs during this period mainly involve membrane organization, nuclear-transcribed mRNA catabolic process, and RNA catabolic process. Findings in the present study are consistent with the observation of the large amount of proliferation and differentiation of granulosa cells in follicles during this period. Notably, KRT18 is one of the genes that is enriched in a number of biological functions, and the encoded protein is the trophoblast cell marker factor (Hou et al., 2015). The expression of the KRT18 gene was reduced by RNA interference techniques, resulting in a lesser developmental competence at the blastocyst stage in cattle (Goossens et al., 2010). By comparing the transcriptome profiles of buffalo embryos that had normal or retarded growth on Day 25 after mating, the expression of the KRT18 gene was significantly up-regulated with the normally growing embryos, suggesting KRT18 gene expression is associated with endometrial attachment and establishment of pregnancy (Strazzullo et al., 2014). In the present study, the expression of the KRT18 gene was down-regulated when comparing the GC2 GC3 stages of development, which indicates that KRT18 is a candidate gene affecting formation of dominant follicles. The immune system has an important role in the last stages of follicular maturation and ovulation. In the present study, there was significantly enriched signal pathways including allograft rejection, chemokine signaling, natural killer cell mediated cytotoxicity, phagosome, and antigen processing and presentation, indicating that immune-related molecules have an important role in dominant follicles developing the ovulation capacity. Cattle are the closest species to buffalo in evolutionary relationship. The comparison of transcriptome profiles of granulosa cells between small (< 5 mm) and large (> 10) follicles of cattle were performed using Affymetrix microarrays (Hatzirodos et al., 2014), and the result suggested that the processes of immune signaling were enriched in large follicles to the greatest extent. When ovulation occurs, this is indicative of the final stage of follicle development. This process was regarded as one type of localized physiological inflammatory reaction. The greatest concentrations of inflammatory substances have been found in the follicle during the peri-ovulatory phase of follicular development (Brannstrom and Enskog 2002). Previous studies demonstrated that inflammatory substances including leukocytes (Kh et al., 2000), macrophages (Takaya et al., 1997), neutrophils (Ujioka et al., 1998), and mast cells (Gaytan et al., 1991) have an important role in the ovulatory process to facilitate ovulation. The knowledge about the exact roles of these cells in ovulation process is, however, limited. Interestingly, there are many DEGs that are enriched in the last stage of the allograft rejection signaling process. A possible interpretation of results from the present study is during stage GC4 of follicle development there is that there is preparation occurring for the protecting of growth of oocytes after ovulation and movement of the oocytes into the fallopian tubes. Of the previously-
190
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
mentioned DEGs, only TNF, GZMB and IL6 have been reported to be associated with female fertility in previous studies. The same situation also exists in natural killer cell mediated cytotoxicity and the phagosome signaling process. It is noteworthy that PAQR7, is enriched in the natural killer cell mediated cytotoxicity signaling process, and has been reported to be related to the induction of oocyte maturation in fish (Shi et al., 2016). Another interesting that the products of the expression of the gene, CTSS, which is enriched in the phagosome signaling process, have been reported to participate in the process of cumulus cell apoptosis and oocyte development (Boruszewska et al., 2014). 5. Conclusion A total of four independent cDNA libraries, respectively, for four stages of ovarian follicle growth, namely, GC1, GC2, GC3 and GC4, were obtained by using RNA-seq. Numerous DEGs were identified during follicle growth. The stage-specifically expressed genes in follicle growth were characterized by GO and KEGG enrichment. The present study revealed overall changes in gene expression during follicle growth. It also indicated that follicle maturation and ovulation might be subject to immune control which involves allograft rejection, chemokine signaling pathway, natural killer cell mediated cytotoxicity, phagosome, and antigen processing and presentation. Competing interests The authors declare that they have no competing interests. Acknowledgements This study was financially supported by the National Natural Science Foundation of China (No. 31772602 and No. 31772604). References Acosta, T.J., 2007. Studies of follicular vascularity associated with follicle selection and ovulation in cattle. J. Reprod. Dev. 53, 39–44. Adams, G.P., Matteri, R.L., Kastelic, J.P., Ko, J.C., Ginther, O.J., 1992. Association between surges of follicle-stimulating hormone and the emergence of follicular waves in heifers. J. Reprod. Fertil. 94, 177–188. Benjamini, Y., Yekutieli, D., 2001. The control of the false discovery rate in multiple testing under dependency. Ann. Stat. 29, 1165–1188. Boruszewska, D., Torres, A.C., Kowalczyk-Zieba, I., Diniz, P., Batista, M., Lopes-da-Costa, L., Woclawek-Potocka, I., 2014. The effect of lysophosphatidic acid during in vitro maturation of bovine oocytes: embryonic development and mRNA abundances of genes involved in apoptosis and oocyte competence. Mediators Inflamm. 670670. Brannstrom, M., Enskog, A., 2002. Leukocyte networks and ovulation. J. Reprod. Immunol. 57, 47–60. Du, Z.Q., Eisley, C.J., Onteru, S.K., Madsen, O., Groenen, M.A.M., Ross, J.W., Rothschild, M.F., 2014. Identification of species-specific novel transcripts in pig reproductive tissues using RNA-seq. Anim. Genet. 45, 198–204. Evans, A.C., Komar, C.M., Wandji, S.A., Fortune, J.E., 1997. Changes in androgen secretion and luteinizing hormone pulse amplitude are associated with the recruitment and growth of ovarian follicles during the luteal phase of the bovine estrous cycle. Biol. Reprod. 57, 394–401. Evans, A.C.O., Ireland, J.L.H., Winn, M.E., Lonergan, P., Smith, G.W., Coussens, P.M., Ireland, J.J., 2004. Identification of genes involved in apoptosis and dominant follicle development during follicular waves in cattle1. Biol. Reprod. 70, 1475–1484. Fortune, J.E., Sirois, J., Turzillo, A.M., Lavoir, M., 1991. Follicle selection in domestic ruminants. J. Reprod. Fertil. Suppl. 43, 187–198. Gaytan, F., Aceitero, J., Bellido, C., Sánchez-Criado, J.E., Aguilar, E., 1991. Estrous cycle-related changes in mast cell numbers in several ovarian compartments in the rat. Biol. Reprod. 45, 27–33. Gimenes, L.U., Carvalho, N.A., Sa Filho, M.F., Vannucci, F.S., Torres-Junior, J.R., Ayres, H., Ferreira, R.M., Trinca, L.A., Sartorelli, E.S., Barros, C.M., Beltran, M.P., Nogueira, G.P., Mapletoft, R.J., Baruselli, P.S., 2011. Ultrasonographic and endocrine aspects of follicle deviation: and acquisition of ovulatory capacity in buffalo (Bubalus bubalis) heifers. Anim. Reprod. Sci. 123, 175–179. Ginther, O.J., Beg, M.A., Bergfelt, D.R., Donadeu, F.X., Kot, K., 2001. Follicle selection in monovular species. Biol. Reprod. 65 638-147. Ginther, O.J., Beg, M.A., Donadeu, F.X., Bergfelt, D.R., 2003. Mechanism of follicle deviation in monovular farm species. Anim. Reprod. Sci. 78 239-157. Goossens, K., Tesfaye, D., Rings, F., Schellander, K., Holker, M., Van Poucke, M., Van Zeveren, A., Lemahieu, I., Van Soom, A., Peelman, L.J., 2010. Suppression of keratin 18 gene expression in bovine blastocysts by RNA interference. Reprod. Fertil. Dev. 22, 395–404. Hatzirodos, N., Irving-Rodgers, H.F., Hummitzsch, K., Harland, M.L., Morris, S.E., Rodgers, R.J., 2014. Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes. BMC Genomics 15, 24. Hou, D., Su, M., Li, X., Li, Z., Yun, T., Zhao, Y., Zhang, M., Zhao, L., Li, R., Yu, H., Li, X., 2015. The efficient derivation of trophoblast cells from porcine in vitro fertilized and parthenogenetic blastocysts and culture with ROCK inhibitor Y-27632. PLoS One 10, e0142442. Jäger, M., Ott, C.E., Grünhagen, J., Hecht, J., Schell, H., Mundlos, S., Duda, G.N., Robinson, P.N., Lienau, Jasmin, 2011. Composite transcriptome assembly of RNA-seq data in a sheep model for delayed bone healing. BMC Genomics 12, 158. Jain, A., Jain, T., Kumar, P., Kumar, M., De, S., Gohain, M., Kumar, R., Datta, T.K., 2016. Follicle-stimulating hormone-induced rescue of cumulus cell apoptosis and enhanced development ability of buffalo oocytes. Domest. Anim. Endocrinol. 55, 74–82. Kh, V.D.H., Maddocks, S., Woodhouse, C.M., Van, R.N., Robertson, S.A., Norman, R.J., 2000. Intrabursal injection of clodronate liposomes causes macrophage depletion and inhibits ovulation in the mouse ovary. Biol. Reprod. 62, 1059–1066. Li, R., Yu, C., Li, Y., Lam, T.W., Yiu, S.M., Kristiansen, K., Wang, J., 2009. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25, 1966–1977. Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A.C., Kanehisa, M., 2007. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 35, W182–185. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. Pandey, A.K., Dhaliwal, G.S., Ghuman, S.P., Agarwal, S.K., 2011. Impact of pre-ovulatory follicle diameter on plasma estradiol, subsequent luteal profiles and conception rate in buffalo (Bubalus bubalis). Anim. Reprod. Sci. 123, 169–174. Presicce, G.A., Senatore, E.M., Bella, A., De Santis, G., Barile, V.L., De Mauro, G.J., Terzano, G.M., Stecco, R., Parmeggiani, A., 2004. Ovarian follicular dynamics and hormonal profiles in heifer and mixed-parity Mediterranean Italian buffaloes (Bubalus bubalis) following an estrus synchronization protocol. Theriogenology 61, 1343–1355. Purohit, G.N., Duggal, G.P., Dadarwal, D., Kumar, D., Yadav, R.C., Vyas, S., 2003. Reproductive biotechnologies for improvement of buffalo: the current status. AsianAustralas. J. Anim. Sci. 16, 1071–1086.
191
Animal Reproduction Science 187 (2017) 181–192
J. Li et al.
Salleh, M.S., Mazzoni, G., Hoglund, J.K., Olijhoek, D.W., Lund, P., Lovendahl, P., Kadarmideen, H.N., 2017. RNA-seq transcriptomics and pathway analyses reveal potential regulatory genes and molecular mechanisms in high- and low-residual feed intake in Nordic dairy cattle. BMC Genomics 18, 258. Sarkar, M., Prakash, B.S., 2005. Timing of ovulation in relation to onset of estrus and LH peak in yak (Poephagus grunniens L.). Anim. Reprod. Sci. 86, 353–362. Shi, B., Liu, X., Thomas, P., Pang, Y., Xu, Y., Li, X., Li, X., 2016. Identification and characterization of a progestin and adipoQ receptor (PAQR) structurally related to Paqr7 in the ovary of Cynoglossus semilaevis and its potential role in regulating oocyte maturation. Gen. Comp. Endocrinol. 237, 109–120. Strazzullo, M., Gasparrini, B., Neglia, G., Balestrieri, M.L., Francioso, R., Rossetti, C., Nassa, G., De Filippo, M.R., Weisz, A., Di Francesco, S., Vecchio, D., D'Esposito, M., D'Occhio, M.J., Zicarelli, L., Campanile, G., 2014. Global transcriptome profiles of Italian Mediterranean buffalo embryos with normal and retarded growth. PLoS One 9, e90027. Takaya, R., Fukaya, T., Sasano, H., Suzuki, T., Tamura, M., Yajima, A., 1997. Macrophages in normal cycling human ovaries; immunohistochemical localization and characterization. Hum. Reprod. 12, 1508–1512. Tsuiko, O., Noukas, M., Zilina, O., Hensen, K., Tapanainen, J.S., Magi, R., Kals, M., Kivistik, P.A., Haller-Kikkatalo, K., Salumets, A., Kurg, A., 2016. Copy number variation analysis detects novel candidate genes involved in follicular growth and oocyte maturation in a cohort of premature ovarian failure cases. Hum. Reprod. 31, 1913–1925. Ujioka, T., Matsukawa, A., Tanaka, N., Matsuura, K., Yoshinaga, M., Okamura, H., 1998. Interleukin-8 as an essential factor in the human chorionic gonadotropininduced rabbit ovulatory process: interleukin-8 induces neutrophil accumulation and activation in ovulation. Biol. Reprod. 58, 526–530. Yeung, C.K., Wang, G., Yao, Y., Liang, J., Tenny Chung, C.Y., Chuai, M., Lee, K.K., Yang, X., 2017. BRE modulates granulosa cell death to affect ovarian follicle development and atresia in the mouse. Cell. Death. Dis. 8, e2697. Young, M.D., Wakefield, M.J., Smyth, G.K., Oshlack, A., 2010. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, R14.
192