- Email: [email protected]
Exploring differentially expressed genes in the ovaries of uniparous and multiparous goats using the RNA-Seq (Quantification) method
Gene 550 (2014) 148–153
Contents lists available at ScienceDirect
Gene journal homepage: www.elsevier.com/locate/gene
Short Communication
Exploring differentially expressed genes in the ovaries of uniparous and multiparous goats using the RNA-Seq (Quantification) method Ying-Hui Ling a,1, Hao Xiang a,b,1, Yun-Sheng Li a, Ya Liu a, Yun-Hai Zhang a, Zi-Juan Zhang a, Jian-Ping Ding a,⁎, Xiao-Rong Zhang a,⁎ a b
College of Animal Science and Technology, Anhui Agricultural University, Hefei, Anhui 230036, China Local Animal Genetic Resources Conservation and Biobreeding Laboratory of Anhui Province, Hefei, Anhui 230036, China
a r t i c l e
i n f o
Article history: Received 22 January 2014 Received in revised form 30 July 2014 Accepted 4 August 2014 Available online 5 August 2014 Keywords: Goat RNA-Seq Ovary Differentially expressed genes Prolificacy
a b s t r a c t Fecundity improvement is one of the most important objectives for goat breeders as it greatly increases production efficiency. To investigate the genes associated with litter sizes in the Anhui White goat (AWG), gene expression differences in the ovaries of uniparous and multiparous AWG were assessed using the RNA-Seq (Quantification) method. This analysis generated 6,027,714 and 5,884,062 clean reads in uniparous and multiparous libraries, respectively. A total of 2201 differentially expressed genes (DEGs) were thereby identified (FDR ≤ 0.001, |log2Ratio| ≥ 1). There were 1583 up-regulated and 618 down-regulated genes in the multiparous samples compared with the uniparous samples. A large number of these DEGs were related to the terms cellular process, cell & cell part and binding. Twelve genes which may be associated with the high prolificacy of AWG were identified using a bioinformatic screen. In addition, pathway analysis revealed that these DEGs were significantly enriched in 11 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, including ECM–receptor interactions, focal adhesion, and regulation of the actin cytoskeleton among others. This suggested a role for these pathways in the prolificacy of AWG. These results provide a list of new candidate genes for goat prolificacy. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Litter size is one of the most important economic traits in goat production and varies among individual animals. It is directly related to the production efficiency and results from the interactions between multiple genes and the environment. There has been some recent progress in characterizing the major genes involved in the prolificacy of sheep i.e. BMPR-IB, BMP15, and GDF9 (Davis et al., 2002; Galloway et al., 2000; Hanrahan et al., 2004). Such information however is still relatively limited in goat. The lambing rates in goats are low, litters of one or two in general, and heritability is also low in this species (Sun et al., 2010), which restricts breeding. There has been an increasing demand for mutton and consequent rapid development of the goat industry in recent years for which traditional breeding methods have not been adequate. The Anhui White goat (AWG) is known for its precocious puberty, higher fertility, and higher leather quality compared with other types Abbreviations: AWG, Anhui White goat; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology; RPKM, reads per kilo bases per million reads; FDR, false discovery rate. ⁎ Corresponding authors. E-mail addresses: [email protected] (J.-P. Ding), [email protected] (X.-R. Zhang). 1 These authors contributed equally to this work.
http://dx.doi.org/10.1016/j.gene.2014.08.008 0378-1119/© 2014 Elsevier B.V. All rights reserved.
of goat. AWG ewes can estrus all year round and can lamb twice in one year or three times in two years. The average AWG lambing rate is 227–239% which belonged to varieties of high lambing rates in goats (Chen et al., 2009). AWG is therefore an ideal model for the study of goat traits. RNA-Seq (Quantification) is mainly used for the quantitative gene expression analysis of biological processes in a particular tissue or cell in certain species. This approach can be used to study genome-wide differences in gene expression by using new generation high-throughput sequencing methods and information analysis platforms. It has many advantages including more accurate quantization, higher repeatability, a wider testing range and thus more reliable analysis. Currently, RNASeq (Quantification) has been used in the study of specific genes in mice (Mortazavi et al., 2008), pig (Gunawan et al., 2013), cattle (Huang and Khatib, 2010), and human (Hackett et al., 2012), but not in the study of prolificacy-related genes in goat. As the ovary directly mediates ovulation and affects litter size, it has a significant impact on the fecundity of mammals (Toosi et al., 2010). In the present study, estrus ovaries of multiparous and uniparous AWGs were sampled in order to analyze the differential expression of genes using RNA-Seq (Quantification). The aim of this analysis was to screen for major genes that control prolificacy in goat, and thus provide a molecular basis for new breeding approaches in this important commercial animal.
Y.-H. Ling et al. / Gene 550 (2014) 148–153
2. Materials and methods 2.1. Experimental animals and sample collection The experimental AWG group in this study was obtained from Boda Livestock Technology Development Co., Ltd of Hefei (Anhui, China). Twelve adult ewes (six multiparous and six uniparous) with at least three births were randomly selected based on breeding records from the previous five years. All animals were similar in age and appearance. These animals were fed under the same conditions of feeding. The ewes in the estrus period were slaughtered and ovaries from both sides were obtained in each case. The separated ovaries were snap frozen in liquid nitrogen in a freezing tube and stored at −80 °C until needed. The ovaries from three multiparous and three uniparous goats were selected for analysis. 2.2. RNA extraction and quality analysis Total RNA was extracted from pooled ovarian tissue samples using RNAiso Plus (Takara, Dalian, China) in accordance with the manufacturer's instructions. The RNA concentration and quality were tested using the SMA1000 UV spectrophotometer (Merinton, Beijing, China). The RNA preparations were stored at −80 °C until use. 2.3. Library preparation and sequencing Total RNA from all ovary samples was pooled prior to library preparation in the two experimental groups. Equimolar quantities of RNA from each ovary sample were combined into one pool. For library preparation, mRNA was first extracted from total RNA using oligo (dT) magnetic beads and sheared into short fragments of about 200 bases. These fragmented mRNAs were then used as templates for cDNA synthesis. The cDNAs were then PCR amplified to complete the library. The cDNA library was sequenced using an Illumina HiSeq™ 2000 platform. 2.4. Mapping reads to the reference genome As the raw reads which were transfered from base calling may contain low quality reads and/or adaptor sequences, preprocessing is necessary before starting further analysis (Lindgreen, 2012; Kerpedjiev et al., 2014). Clean reads were mapped to goat reference gene sequences (http://goat.kiz.ac.cn/GGD/download.htm) and to goat reference genome sequences (http://goat.kiz.ac.cn/GGD/download.htm) set using SOAPaligner/SOAP2 (Li et al., 2009). No more than two mismatches were allowed in the alignment.
149
et al., 2000; Zhou and Su, 2007). GO enrichment analysis provides all GO terms that are significantly enriched in DEGs compared with the genome background, and filters the DEGs that correspond to biological functions. With nr annotation, the Blast2GO program was used to obtain the GO annotation of the DEGs. WEGO software (Ye et al., 2006) was then used to perform GO functional classification for DEGs and determine the distribution of gene functions of the species at the macro level. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is the major public pathway-related database (Kanehisa et al., 2008). Pathway enrichment analysis identifies significantly enriched metabolic pathways or signal transduction pathways in DEGs compared with the whole genome background. A Q value ≤ 0.05 was defined as a significantly enriched pathway in terms of DEGs. 2.7. Validation of RNA-Seq data To validate the DEGs revealed by RNA-Seq, 10 genes (Table S1) identified to be differentially expressed were chosen randomly for qPCR validation. The primers used, designed using Primer Express 3.0 software for qPCR analysis, are listed in Table 1 with β-actin used as a reference control. The qPCR reactions were carried out on an ABI StepOnePlus™ Real-Time PCR System (ABI, USA) using a SYBR Green qPCR Mix (Roche, Indianapolis, USA) in accordance with the manufacturer's instructions. The thermal cycling conditions were 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. There were three replicates for each amplification. Relative quantification analyses were performed using the comparative CT method, and relative gene expression levels were calculated using the 2−ΔΔCT method (Ishida-Takagishi et al., 2012). 3. Results 3.1. Analysis of RNA-Seq (Quantification) libraries Two RNA-Seq libraries were constructed from multiparous and uniparous samples and generated 5.9 million raw reads which are sufficient for the quantitative analysis of gene expression. The rates of unique match or multi-position match reads were 82.14% and 84.60% for the uniparous and multiparous samples, respectively (Table 2). After filtering the adaptor sequences, regions containing N sequences and low quality sequences, the two RNA-Seq libraries still generated over 5.8 million clean reads in each library, and the percentages of
Table 1 List of primers used in the qPCR analysis.
2.5. Determination of gene expression levels and detection of differentially expressed genes
Gene
Primer sequences
FER1L4
The expression level for each gene was determined by the numbers of reads uniquely mapped to the specific gene and the total number of uniquely mapped reads in the sample. The gene expression level was then calculated using the reads per kilo bases per million reads (RPKM) method (Mortazavi et al., 2008). The statistical significance of the differential expression of each gene was determined using the P value (Benjamini et al., 2001). The FDR (false discovery rate) ≤ 0.001 and |log2Ratio| ≥ 1 were used to identify differentially expressed genes (DEGs).
Itga4
F: 5′-CCGGCCTGGATACCTCTCTAT-3′ R: 5′-CAGTCGGCCTTGGAAACAAA-3′ F: 5′-ACTGCCAACTGGCTCTCCAA-3′ R: 5′-CCACAAGGATCTCCGAAAGG-3′ F: 5′-GGCATAATTGCTGGCACCAT-3′ R: 5′-AGTCGCCATCTCCAGGCTTT-3′ F: 5′-TTTGGTTTTCCATTTGCAATGA-3′ R: 5′-TGTTCAAGCATGGGCAAAAAG-3′ F: 5′-CTCCCGGCTAGCCAGTTCTC-3′ R: 5′-GTTGGCGATGAAACCCATCT-3′ F: 5′-GCTCAGGAAGCCTGGAGAAAT-3′ R: 5′-GGCATAGCCGATCCATTCAA-3′ F: 5′-AAGCCCCAGGAAACAGATGA-3′ R: 5′-AGCGGATGCAGGACAGAAGA-3′ F: 5′-CAAGGCCAACAAGCAGGTTT-3′ R: 5′-CAGTCGGCCACGGTATTGAC-3′ F: 5′-GAGACCACGTACCGGCAAAC-3′ R: 5′-CGTCCAAGCGGGTAGTTGTC-3′ F: 5′-CAGCAAAATCCCGGGTTAAC-3′ R: 5′-TAGCTGCCCTGAGCATCTGA-3′ F: 5′-GTCATTGAGAGCAATGCCAG-3′ R: 5′-GTGTTCCTACCCCCAATGTG-3′
ADAM22 SYCP1 CDC14A SRD5A2 NMU
2.6. Gene ontology and pathway enrichment analysis of differentially expressed genes
COMP
Gene ontology (GO) is an international standardized gene functional classification system which offers a dynamic-updated controlled vocabulary and strictly defined concepts to comprehensively describe the properties of genes and their products in any organism (Ashburner
WWC1
TNC
β-Actin
150
Y.-H. Ling et al. / Gene 550 (2014) 148–153
clean reads among the raw reads reached 98.60% and 98.16% in the uniparous and multiparous libraries, respectively (Fig. 1). This demonstrated that the two libraries were of high-quality. All data were submitted to NCBI (http://www.ncbi.nlm.nih.gov/), including one study (ID: SRP041866), two samples (uniparous ID: SRS605132, multiparous ID: SRS605133), two experiments (uniparous ID: SRX540103, multiparous ID: SRX540104) and two runs (uniparous ID: SRR1283211, multiparous ID: SRR1283212).
3.2. Gene expression levels and differentially expressed genes To detect prolificacy-related genes, the expression levels of all detected genes were calculated (Fig. 2 and Table S2). One-thousandand-nine genes expressed only in uniparous sample and 858 genes expressed only in multiparous sample in all of the 15,890 detected genes. A total of 2201 genes were differentially expressed between the uniparous and multiparous libraries, in which 1583 genes were upregulated and 618 genes were downregulated in the multiparous library (Table S3).
3.3. Functional enrichment analysis GO contains cellular component, molecular function and biological process domains. Based on sequence homology, the 2201 DEGs identified in our analyses could be categorized into 50 functional groups. In the three main GO classification categories, 23, 15 and 12 functional groups were identified (Fig. 3). Among these groups, the terms cellular process, cell & cell part and binding were predominant in each of the three main categories. A high percentage of genes in the following categories was also observed among the DEGs: biological regulation, cellular component organization or biogenesis, developmental process, establishment of localization, localization, metabolic processes, multicellular organismal process, positive regulation of biological processes, regulation of biological processes, response to stimulus, signaling and single-organism process, macromolecular complex, membrane, membrane part, membrane-enclosed lumen, organelle, organelle part and catalytic activity. KEGG pathway significant enrichment analysis revealed that the DEGs participated in 244 pathways (Table S4) which are known to be enriched in Huntington's disease, endocytosis, lysine degradation, inositol phosphate metabolism, ECM–receptor interaction, focal adhesion, oxidative phosphorylation, Alzheimer's disease, pathways in cancer, the phosphatidylinositol signaling system and regulation of the actin cytoskeleton (Table 3).
3.4. Confirmation of differential gene expression by qPCR Ten DEGs were randomly selected for qPCR analysis to validate the expression profiles obtained by RNA-Seq. The results verified that FER1L4, Itga4, ADAM22, SYCP1, and CDC14A were upregulated and SRD5A2, NMU, COMP, TNC, and WWC1 were downregulated in multiparous samples (Fig. 4), which was basically consistent with the RNA-Seq (Quantification) findings. Thus, RNA-Seq can be used to reliably and accurately perform mRNA differential expression analysis.
Fig. 1. Composition of total raw reads from the uniparous and multiparous libraries. After filtering the only adaptor sequences, which contain N sequences and low quality sequences, the two libraries still generated over 5.8 million clean reads each. The percentage of clean reads among the raw reads reached 98.60% and 98.16% in the uniparous and multiparous libraries, respectively.
4. Discussion 4.1. The functions of the 10 most differentially expressed goat genes found in this study The 10 most differentially expressed goat genes identified in the present study (Table S5) are an important pool of candidate genes for future analysis. EVI2A may form homocomplexes or associate with other proteins within the membrane, as part of a cell-surface receptor complex (Cawthon et al., 1990). The Eyes Absent (EYA) proteins, firstly described in the context of fly eye development, are now implicated in processes as disparate as organ development, innate immunity, DNA damage repair, photoperiodism, angiogenesis, cancer metastasis, and early embryonic development (Ishihara et al., 2008; Tadjuidje and Hegde, 2013). 5 alpha-SR2 (3-oxo-5-alpha-steroid 4dehydrogenase 2, encoded by SRD5A2) converts testosterone (T) into 5-alpha-dihydrotestosterone (DHT) and progesterone or corticosterone into their corresponding 5-alpha-3-oxosteroids. This protein plays a central role in sexual differentiation and androgen physiology (Chang et al., 2003; Seo et al., 2009). Somatostatin (encoded by SST) inhibits the release of somatotropin, and some research has found that somatotropin has effects in animal reproduction (Anderson et al., 1991). Bone sialoprotein 2 (encoded by IBSP) binds tightly to hydroxyapatite and appears to form an integral part of the mineralized matrix. It is possibly important for cell–matrix interactions and likely promotes Arg-GlyAsp-dependent cell attachment (Ganss et al., 1999). BT (butyrophilin
Table 2 Read numbers and mapping results for the two constructed ovary libraries. Sample ID
Total reads
Total mapped reads
Unique match
Multi-position match
Total unmapped reads
Uniparous
6,027,714 (100.00%) 5,884,062 (100.00%)
4,951,157 (82.14%) 4,977,880 (84.60%)
4,071,224 (67.54%) 4,192,641 (71.25%)
879,933 (14.60%) 785,239 (13.35%)
1,076,557 (17.86%) 906,182 (15.40%)
Multiparous
Y.-H. Ling et al. / Gene 550 (2014) 148–153
151
GVINP1 and GVIN1 in the cytosol and nucleus, and the molecular functions of GVINP1 and GVIN1 involve GTP binding. Hence, the function of FER1L4 may be in the exchange of information between the inside and outside of cells, and the functions of GVINP1 and GVIN1 may be in cell energy metabolism, DNA replication and transcription. In view of the high log2Ratio values of these three genes (N10), they may also be associated with high fecundity, and are possibly the key candidate genes to emerge from the current analysis. 4.2. Some of the validated DEGs have previously been associated with high fecundity
Fig. 2. Scatter plot indicating the comparative results of log transformed gene expression levels and differentially expressed gene distributions between uniparous and multiparous samples.
subfamily 1 member A1, encoded by BTN1A1) may have some role in the secretion of milk-fat droplets and act as a specific membraneassociated receptor for the association of cytoplasmic droplets with the apical plasma membrane (Jack and Mather, 1990). CRP3 (cysteine and glycine-rich protein 3) encoded by CSRP3, and the CRP family (CRP1, CRP2, and CRP3) have been implicated in the regulation of muscle development (Xu et al., 2010). CRP3 may be a scaffold protein that promotes the assembly of interacting proteins at Z-line structures (Knoll et al., 2002). In addition, CSRP3 has been confirmed as a HCM disease gene (Geier et al., 2008). In summary, the functions of EVI2A, EYA2, SRD5A2, SST and IBSP involve the exchange of information between cells, DNA repair and hormone secretion which may be closely related to the development of follicular cells and oocytes. It was initially considered that these five genes may be associated with high fecundity in goats. The functions of FER1L4, GVINP1 and GVIN1 are unknown. However, within the cellular component of the GO, FER1L4 is positioned in the membrane and
Genes that emerged from the current screen and have been previously associated with high fecundity were further analyzed. These genes such as GDF9, GnRHR, INH, FSHR, BMPs, and OPNm mainly encode hormones and cytokines. The genes Bmp6 (bone morphogenetic protein 6), INHBB (inhibin beta B), INHA (inhibin, alpha) and OPN (osteopontin) (Koji et al., 2010; Steinheuer et al., 2003) were also identified here as DEGs in the goat, suggesting that they may also be relevant to goat prolificacy (Table S6). The absolute value of the log2Ratio of OPN is relatively high (greater than 3.7) and it may also therefore be a key candidate gene. Other than these four genes, differential expression was non-significant for all of the genes identified in this study. Possible reasons for this include differences between species and an insufficient number of samples. 4.3. Signaling pathways in reproduction At least three of the 11 significantly enriched pathways identified in this study, ECM–receptor interaction (Fig. S1), focal adhesion (Fig. S2) and regulation of actin cytoskeleton (Fig. S3), are intimately associated with mammalian reproduction. The extracellular matrix (ECM) consists of a complex mixture of structural and functional macromolecules and plays an important role in tissue and organ morphogenesis. Specific interactions between cells and the ECM lead to the direct or indirect control of cellular activities such as adhesion, migration, differentiation, proliferation, and apoptosis (Bosman and Stamenkovic, 2003). Focal adhesion plays essential roles in important biological processes including cell motility, cell proliferation, cell differentiation, regulation of gene expression and cell survival (Kanehisa et al., 2008). Regulation of the dynamic actin cytoskeleton includes signaling to the cytoskeleton which
Fig. 3. Histogram of gene ontology classifications. The results were summarized in three main categories: biological process, cellular component and molecular function. Among these groups, the terms cellular process, cell & cell part and binding were predominant in each of the three main categories, respectively.
152
Y.-H. Ling et al. / Gene 550 (2014) 148–153
Table 3 Significantly enriched pathways among the differentially expressed genes between the uniparous and multiparous samples. Pathway
DEGs with pathway annotation (1825)
All genes with pathway annotation (16,032)
P value
Q value
Pathway ID
Huntington's disease Endocytosis Lysine degradation Inositol phosphate metabolism ECM–receptor interaction Focal adhesion Oxidative phosphorylation Alzheimer's disease Pathways in cancer Phosphatidylinositol signaling system Regulation of actin cytoskeleton
55 (3.01%) 73 (4%) 26 (1.42%) 24 (1.32%) 38 (2.08%) 70 (3.84%) 32 (1.75%) 44 (2.41%) 84 (4.6%) 31 (1.7%) 80 (4.38%)
253 (1.58%) 385 (2.4%) 99 (0.62%) 97 (0.61%) 190 (1.19%) 412 (2.57%) 153 (0.95%) 233 (1.45%) 518 (3.23%) 149 (0.93%) 514 (3.21%)
1.511431e −06 7.280254e −06 3.141498e −05 0.0001715156 0.0003662485 0.0003790018 0.0004581824 0.0004921621 0.0004933844 0.0006124814 0.002192947
0.0004141321 0.0009973948 0.0028692348 0.0117488186 0.0168984157 0.0168984157 0.0133761993 0.0168984157 0.0168984157 0.0186466560 0.0486435516
ko05016 ko04144 ko00310 ko00562 ko04512 ko04510 ko00190 ko05010 ko05200 ko04070 ko04810
can lead to diverse effects on cellular activity such as changes in cell shape, migration and proliferation (Lee and Dominguez, 2010). The cell will migrate and move towards the adluminal compartment to prepare for meiosis if the actin cytoskeleton dynamics have changed appropriately (Zhang et al., 2013). These signaling pathways are closely associated with cell differentiation, proliferation, and meiosis during the development of follicles. Most of the DEGs (22/38, 50/70 and 59/80, respectively) identified in this study which participate in the abovementioned three pathways were up-regulated. This indicates that the three pathways are activated in the estrus ovaries of multiparous goats. As a whole therefore, these genes may promote the development of follicles and oocytes but the mechanism is unclear. Previous studies have indicated that pathways involving the DEGs in this study also associate with reproduction except the Notch signaling pathway (Larkin et al., 1996; Lopez-Schier and St Johnston, 2001), TGF-beta signaling pathway (Shi and Massagué, 2003), steroid biosynthesis (Issop et al., 2013), MAPK signaling pathway (Xiao et al., 2007), Wnt signaling pathway (Habas and He, 2006; Ishida-Takagishi et al., 2012) and Jak–STAT signaling pathway (Kisseleva et al., 2002). However, the exact relationship between these signaling networks is not fully understood.
5. Conclusion This study screened for DEGs in the ovarian tissues of the multiparous and uniparous Anhui White goats using RNA-Seq (Quantification), and 2201 DEGs were thereby obtained. Among the pathways in which these DEGs are mainly involved, at least three are associated with reproduction, and 12 of the DEGs may be relevant to the prolificacy of goats. This new information provides a solid foundation for further studies of the molecular mechanisms underlying goat prolificacy.
Fig. 4. Validation of the RNA-Seq (Quantification) results by qPCR. Results represent the mean (±SD) of three experiments. Asterisk, P ≤ 0.05.
Further biochemical and physiological analyses of these candidate genes will be conducted in the future. Acknowledgments This research was supported by the National Natural Science Foundation of China (31301934, 31372310), the Natural Science Foundation of Anhui Province (1308085QC54) and Anhui Science and Technology Program (11010302108) funds. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2014.08.008. References Anderson, L.L., Ford, J.J., Klindt, J., et al., 1991. Growth hormone and prolactin secretion in hypophysial stalk-transected pigs as affected by growth hormone and prolactinreleasing and inhibiting factors. Proc. Soc. Exp. Biol. Med. 196, 194–202. Ashburner, M., Ball, C.A., Blake, J.A., et al., 2000. Gene ontology: tool for the unification of biology. Nat. Genet. 25, 25–29. Benjamini, Y.,Drai, D., Elmer, G., et al., 2001. Controlling the false discovery rate in behavior genetics research. Behav. Brain Res. 125, 279–284. Bosman, F.T.,Stamenkovic, I., 2003. Functional structure and composition of the extracellular matrix. J. Pathol. 200 (4), 423–428. Cawthon, R.M., O'Connell, P., Buchberg, A.M., et al., 1990. Identification and characterization of transcripts from the neurofibromatosis 1 region: the sequence and genomic structure of EVI2 and mapping of other transcripts. Genomics 7 (4), 555–565. Chang, B.L.,Zheng, S.L.,Isaacs, S.D., et al., 2003. Evaluation of SRD5A2 sequence variants in susceptibility to hereditary and sporadic prostate cancer. Prostate 56 (1), 37–44. Chen, S.,Cheng, G.L.,Zhu, D.J., et al., 2009. Determination on the body properties and meat performance of Anhui white goat. Anim. Husb. Feed Sci. 30 (4), 150–152. Davis, G.H., Galloway, S.M., Ross, I.K., et al., 2002. DNA tests in prolific sheep from eight countries provide new evidence on origin of the Booroola (FecB) mutation. Biol. Reprod. 66 (6), 1869–1874. Galloway, S.M.,McNatty, K.P.,Cambridge, L.M., et al., 2000. Mutations in an oocyte-derived growth factor gene (BMP15) cause increased ovulation rate and infertility in a dosage-sensitive manner. Nat. Genet. 25 (3), 279–283. Ganss, B., Kim, R.H., Sodek, J., 1999. Bone sialoprotein. Crit. Rev. Oral Biol. Med. 10 (1), 79–98. Geier, C.,Gehmlich, K.,Ehler, E., et al., 2008. Beyond the sarcomere: CSRP3 mutations cause hypertrophic cardiomyopathy. Hum. Mol. Genet. 17 (18), 2753–2765. Gunawan, A., Sahadevan, S., Neuhoff, C., et al., 2013. RNA deep sequencing reveals novel candidate genes and polymorphisms in boar testis and liver tissues with divergent androstenone levels. PLos One 8 (5), e63259. Habas, R., He, X., 2006. Activation of Rho and Rac by Wnt/frizzled signaling. Methods Enzymol. 406, 500–511. Hackett, N.R., Butler, M.W., Shaykhiev, R., et al., 2012. RNA-Seq quantification of the human small airway epithelium transcriptome. BMC Genomics 13, 82. Hanrahan, J.P., Gregan, S.M., Mulsant, P., et al., 2004. Mutations in the genes for oocytederived growth factors GDF9 and BMP15 are associated with both increased ovulation rate and sterility in Cambridge and Belclare sheep (Ovis aries). Biol. Reprod. 70, 900–909. Huang, W.,Khatib, H., 2010. Comparison of transcriptomic landscapes of bovine embryos using RNA-Seq. BMC Genomics 11, 711. Ishida-Takagishi, M.,Enomoto, A.,Asai, N., et al., 2012. The Dishevelled-associating protein Daple controls the non-canonical Wnt/Rac pathway and cell motility. Nat. Commun. 3, 859. Ishihara, T.,Ikeda, K., Sato, S., 2008. Differential expression of Eya1 and Eya2 during chick. Gene Expr. Patterns 8 (5), 357–367.
Y.-H. Ling et al. / Gene 550 (2014) 148–153 Issop, L., Rone, M.B., Papadopoulos, V., 2013. Organelle plasticity and interactions in cholesterol transport and steroid biosynthesis. Mol. Cell. Endocrinol. 371 (1–2), 34–46. Jack, L.J., Mather, I.H., 1990. Cloning and analysis of cDNA encoding bovine butyrophilin, an apical glycoprotein expressed in mammary tissue and secreted in association with the milk-fat globule membrane during lactation. J. Biol. Chem. 265 (24), 14481–14486. Kanehisa, M.,Araki, M.,Goto, S., et al., 2008. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36 (Suppl. 1), D480–D484. Kerpedjiev, P., Frellsen, J., Lindgreen, S., et al., 2014. Adaptable probabilistic mapping of short reads using position specific scoring matrices. BMC Bioinforma. 15, 100. Kisseleva, T., Bhattacharya, S., Braunstein, J., et al., 2002. Signaling through the JAK/STAT pathway, recent advances and future challenges. Gene 285 (1–2), 1–24. Knoll, R.,Hoshijima, M.,Hoffman, H.M., et al., 2002. The cardiac mechanical stretch sensor machinery involves a Z disc complex that is defective in a subset of human dilated cardiomyopathy. Cell 111, 943–955. Koji, S., You, Q.S., John, J.E., 2010. Does bone morphogenetic protein 6 (BMP-6) affect female fertility in the mouse? Biol. Reprod. 83 (6), 997–1004. Larkin, M.K., Holder, K., Yost, C., et al., 1996. Expression of constitutively active Notch arrests follicle cells at a precursor stage during Drosophila oogenesis and disrupts the anterior–posterior axis of the oocyte. Development 122 (11), 3639–3650. Lee, S.H., Dominguez, R., 2010. Regulation of actin cytoskeleton dynamics in cells. Mol. Cells 29, 311–325. Li, R.Q., Yu, C., Li, Y.R., et al., 2009. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25 (15), 1966–1967. Lindgreen, S., 2012. AdapterRemoval: easy cleaning of next-generation sequencing reads. BMC Res. Notes 5, 337. Lopez-Schier, H.,St Johnston, D., 2001. Delta signaling from the germ line controls the proliferation and differentiation of the somatic follicle cells during Drosophila oogenesis. Genes Dev. 15 (11), 1393–1405. Mortazavi, A.,Williams, B.A.,McCue, K., et al., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5 (7), 621–628.
153
Seo, Y.K., Zhu, B., Jeon, T.I., et al., 2009. Regulation of steroid 5-alpha reductase type 2 (Srd5a2) by sterol regulatory element binding proteins and statin. Exp. Cell Res. 315 (18), 3133–3139. Shi, Y.G., Massagué, J., 2003. Mechanisms of TGF-β signaling from cell membrane to the nucleus. Cell 113 (6), 685–700. Steinheuer, R., Drögemüller, C., Hamann, H., et al., 2003. Possible uses of genetic markers for improving fertility and health in swine production. Dtsch. Tierarztl. Wochenschr. 110 (6), 255–265. Sun, W., Chang, H., Musa, H.H., et al., 2010. Study on relationship between microsatellite polymorphism and producing ability on litter size trait of Hu sheep in China. Afr. J. Biotechnol. 9 (50), 8704–8711. Tadjuidje, E.,Hegde, R.S., 2013. The Eyes Absent proteins in development and disease. Cell. Mol. Life Sci. 70 (11), 1897–1913. Toosi, B.M., Seekallu, S.V., Rawlings, N.C., 2010. Effects of the rate and duration of physiological increases in serum FSH concentrations on emergence of follicular waves in cyclic ewes. Biol. Reprod. 83 (4), 648–655. Xiao, Y.,Yang, F.Q.,Li, S.P., et al., 2007. Furanodiene induces G2/M cell cycle arrest and apoptosis through MAPK signaling and mitochondria–caspase pathway in human hepatocellular carcinoma cells. Cancer Biol. Ther. 6 (7), 1044–1050. Xu, X.W., Qiu, H.F., Du, Z.Q., et al., 2010. Porcine CSRP3: polymorphism and association analyses with meat quality traits and comparative analyses with CSRP1 and CSRP2. Mol. Biol. Rep. 37 (1), 451–459. Ye, J., Fang, L., Zheng, H.K., et al., 2006. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 34 (Suppl. 2), W293–W297. Zhang, X.J.,Hao, L.L.,Meng, L.J., et al., 2013. Digital gene expression tag profiling analysis of the gene expression patterns regulating the early stage of mouse spermatogenesis. PLos One 8 (3), e58680. Zhou, X., Su, Z., 2007. EasyGO: Gene Ontology-based annotation and functional enrichment analysis tool for agronomical species. BMC Genomics 8, 246.
Contents lists available at ScienceDirect
Gene journal homepage: www.elsevier.com/locate/gene
Short Communication
Exploring differentially expressed genes in the ovaries of uniparous and multiparous goats using the RNA-Seq (Quantification) method Ying-Hui Ling a,1, Hao Xiang a,b,1, Yun-Sheng Li a, Ya Liu a, Yun-Hai Zhang a, Zi-Juan Zhang a, Jian-Ping Ding a,⁎, Xiao-Rong Zhang a,⁎ a b
College of Animal Science and Technology, Anhui Agricultural University, Hefei, Anhui 230036, China Local Animal Genetic Resources Conservation and Biobreeding Laboratory of Anhui Province, Hefei, Anhui 230036, China
a r t i c l e
i n f o
Article history: Received 22 January 2014 Received in revised form 30 July 2014 Accepted 4 August 2014 Available online 5 August 2014 Keywords: Goat RNA-Seq Ovary Differentially expressed genes Prolificacy
a b s t r a c t Fecundity improvement is one of the most important objectives for goat breeders as it greatly increases production efficiency. To investigate the genes associated with litter sizes in the Anhui White goat (AWG), gene expression differences in the ovaries of uniparous and multiparous AWG were assessed using the RNA-Seq (Quantification) method. This analysis generated 6,027,714 and 5,884,062 clean reads in uniparous and multiparous libraries, respectively. A total of 2201 differentially expressed genes (DEGs) were thereby identified (FDR ≤ 0.001, |log2Ratio| ≥ 1). There were 1583 up-regulated and 618 down-regulated genes in the multiparous samples compared with the uniparous samples. A large number of these DEGs were related to the terms cellular process, cell & cell part and binding. Twelve genes which may be associated with the high prolificacy of AWG were identified using a bioinformatic screen. In addition, pathway analysis revealed that these DEGs were significantly enriched in 11 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, including ECM–receptor interactions, focal adhesion, and regulation of the actin cytoskeleton among others. This suggested a role for these pathways in the prolificacy of AWG. These results provide a list of new candidate genes for goat prolificacy. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Litter size is one of the most important economic traits in goat production and varies among individual animals. It is directly related to the production efficiency and results from the interactions between multiple genes and the environment. There has been some recent progress in characterizing the major genes involved in the prolificacy of sheep i.e. BMPR-IB, BMP15, and GDF9 (Davis et al., 2002; Galloway et al., 2000; Hanrahan et al., 2004). Such information however is still relatively limited in goat. The lambing rates in goats are low, litters of one or two in general, and heritability is also low in this species (Sun et al., 2010), which restricts breeding. There has been an increasing demand for mutton and consequent rapid development of the goat industry in recent years for which traditional breeding methods have not been adequate. The Anhui White goat (AWG) is known for its precocious puberty, higher fertility, and higher leather quality compared with other types Abbreviations: AWG, Anhui White goat; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology; RPKM, reads per kilo bases per million reads; FDR, false discovery rate. ⁎ Corresponding authors. E-mail addresses: [email protected] (J.-P. Ding), [email protected] (X.-R. Zhang). 1 These authors contributed equally to this work.
http://dx.doi.org/10.1016/j.gene.2014.08.008 0378-1119/© 2014 Elsevier B.V. All rights reserved.
of goat. AWG ewes can estrus all year round and can lamb twice in one year or three times in two years. The average AWG lambing rate is 227–239% which belonged to varieties of high lambing rates in goats (Chen et al., 2009). AWG is therefore an ideal model for the study of goat traits. RNA-Seq (Quantification) is mainly used for the quantitative gene expression analysis of biological processes in a particular tissue or cell in certain species. This approach can be used to study genome-wide differences in gene expression by using new generation high-throughput sequencing methods and information analysis platforms. It has many advantages including more accurate quantization, higher repeatability, a wider testing range and thus more reliable analysis. Currently, RNASeq (Quantification) has been used in the study of specific genes in mice (Mortazavi et al., 2008), pig (Gunawan et al., 2013), cattle (Huang and Khatib, 2010), and human (Hackett et al., 2012), but not in the study of prolificacy-related genes in goat. As the ovary directly mediates ovulation and affects litter size, it has a significant impact on the fecundity of mammals (Toosi et al., 2010). In the present study, estrus ovaries of multiparous and uniparous AWGs were sampled in order to analyze the differential expression of genes using RNA-Seq (Quantification). The aim of this analysis was to screen for major genes that control prolificacy in goat, and thus provide a molecular basis for new breeding approaches in this important commercial animal.
Y.-H. Ling et al. / Gene 550 (2014) 148–153
2. Materials and methods 2.1. Experimental animals and sample collection The experimental AWG group in this study was obtained from Boda Livestock Technology Development Co., Ltd of Hefei (Anhui, China). Twelve adult ewes (six multiparous and six uniparous) with at least three births were randomly selected based on breeding records from the previous five years. All animals were similar in age and appearance. These animals were fed under the same conditions of feeding. The ewes in the estrus period were slaughtered and ovaries from both sides were obtained in each case. The separated ovaries were snap frozen in liquid nitrogen in a freezing tube and stored at −80 °C until needed. The ovaries from three multiparous and three uniparous goats were selected for analysis. 2.2. RNA extraction and quality analysis Total RNA was extracted from pooled ovarian tissue samples using RNAiso Plus (Takara, Dalian, China) in accordance with the manufacturer's instructions. The RNA concentration and quality were tested using the SMA1000 UV spectrophotometer (Merinton, Beijing, China). The RNA preparations were stored at −80 °C until use. 2.3. Library preparation and sequencing Total RNA from all ovary samples was pooled prior to library preparation in the two experimental groups. Equimolar quantities of RNA from each ovary sample were combined into one pool. For library preparation, mRNA was first extracted from total RNA using oligo (dT) magnetic beads and sheared into short fragments of about 200 bases. These fragmented mRNAs were then used as templates for cDNA synthesis. The cDNAs were then PCR amplified to complete the library. The cDNA library was sequenced using an Illumina HiSeq™ 2000 platform. 2.4. Mapping reads to the reference genome As the raw reads which were transfered from base calling may contain low quality reads and/or adaptor sequences, preprocessing is necessary before starting further analysis (Lindgreen, 2012; Kerpedjiev et al., 2014). Clean reads were mapped to goat reference gene sequences (http://goat.kiz.ac.cn/GGD/download.htm) and to goat reference genome sequences (http://goat.kiz.ac.cn/GGD/download.htm) set using SOAPaligner/SOAP2 (Li et al., 2009). No more than two mismatches were allowed in the alignment.
149
et al., 2000; Zhou and Su, 2007). GO enrichment analysis provides all GO terms that are significantly enriched in DEGs compared with the genome background, and filters the DEGs that correspond to biological functions. With nr annotation, the Blast2GO program was used to obtain the GO annotation of the DEGs. WEGO software (Ye et al., 2006) was then used to perform GO functional classification for DEGs and determine the distribution of gene functions of the species at the macro level. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is the major public pathway-related database (Kanehisa et al., 2008). Pathway enrichment analysis identifies significantly enriched metabolic pathways or signal transduction pathways in DEGs compared with the whole genome background. A Q value ≤ 0.05 was defined as a significantly enriched pathway in terms of DEGs. 2.7. Validation of RNA-Seq data To validate the DEGs revealed by RNA-Seq, 10 genes (Table S1) identified to be differentially expressed were chosen randomly for qPCR validation. The primers used, designed using Primer Express 3.0 software for qPCR analysis, are listed in Table 1 with β-actin used as a reference control. The qPCR reactions were carried out on an ABI StepOnePlus™ Real-Time PCR System (ABI, USA) using a SYBR Green qPCR Mix (Roche, Indianapolis, USA) in accordance with the manufacturer's instructions. The thermal cycling conditions were 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. There were three replicates for each amplification. Relative quantification analyses were performed using the comparative CT method, and relative gene expression levels were calculated using the 2−ΔΔCT method (Ishida-Takagishi et al., 2012). 3. Results 3.1. Analysis of RNA-Seq (Quantification) libraries Two RNA-Seq libraries were constructed from multiparous and uniparous samples and generated 5.9 million raw reads which are sufficient for the quantitative analysis of gene expression. The rates of unique match or multi-position match reads were 82.14% and 84.60% for the uniparous and multiparous samples, respectively (Table 2). After filtering the adaptor sequences, regions containing N sequences and low quality sequences, the two RNA-Seq libraries still generated over 5.8 million clean reads in each library, and the percentages of
Table 1 List of primers used in the qPCR analysis.
2.5. Determination of gene expression levels and detection of differentially expressed genes
Gene
Primer sequences
FER1L4
The expression level for each gene was determined by the numbers of reads uniquely mapped to the specific gene and the total number of uniquely mapped reads in the sample. The gene expression level was then calculated using the reads per kilo bases per million reads (RPKM) method (Mortazavi et al., 2008). The statistical significance of the differential expression of each gene was determined using the P value (Benjamini et al., 2001). The FDR (false discovery rate) ≤ 0.001 and |log2Ratio| ≥ 1 were used to identify differentially expressed genes (DEGs).
Itga4
F: 5′-CCGGCCTGGATACCTCTCTAT-3′ R: 5′-CAGTCGGCCTTGGAAACAAA-3′ F: 5′-ACTGCCAACTGGCTCTCCAA-3′ R: 5′-CCACAAGGATCTCCGAAAGG-3′ F: 5′-GGCATAATTGCTGGCACCAT-3′ R: 5′-AGTCGCCATCTCCAGGCTTT-3′ F: 5′-TTTGGTTTTCCATTTGCAATGA-3′ R: 5′-TGTTCAAGCATGGGCAAAAAG-3′ F: 5′-CTCCCGGCTAGCCAGTTCTC-3′ R: 5′-GTTGGCGATGAAACCCATCT-3′ F: 5′-GCTCAGGAAGCCTGGAGAAAT-3′ R: 5′-GGCATAGCCGATCCATTCAA-3′ F: 5′-AAGCCCCAGGAAACAGATGA-3′ R: 5′-AGCGGATGCAGGACAGAAGA-3′ F: 5′-CAAGGCCAACAAGCAGGTTT-3′ R: 5′-CAGTCGGCCACGGTATTGAC-3′ F: 5′-GAGACCACGTACCGGCAAAC-3′ R: 5′-CGTCCAAGCGGGTAGTTGTC-3′ F: 5′-CAGCAAAATCCCGGGTTAAC-3′ R: 5′-TAGCTGCCCTGAGCATCTGA-3′ F: 5′-GTCATTGAGAGCAATGCCAG-3′ R: 5′-GTGTTCCTACCCCCAATGTG-3′
ADAM22 SYCP1 CDC14A SRD5A2 NMU
2.6. Gene ontology and pathway enrichment analysis of differentially expressed genes
COMP
Gene ontology (GO) is an international standardized gene functional classification system which offers a dynamic-updated controlled vocabulary and strictly defined concepts to comprehensively describe the properties of genes and their products in any organism (Ashburner
WWC1
TNC
β-Actin
150
Y.-H. Ling et al. / Gene 550 (2014) 148–153
clean reads among the raw reads reached 98.60% and 98.16% in the uniparous and multiparous libraries, respectively (Fig. 1). This demonstrated that the two libraries were of high-quality. All data were submitted to NCBI (http://www.ncbi.nlm.nih.gov/), including one study (ID: SRP041866), two samples (uniparous ID: SRS605132, multiparous ID: SRS605133), two experiments (uniparous ID: SRX540103, multiparous ID: SRX540104) and two runs (uniparous ID: SRR1283211, multiparous ID: SRR1283212).
3.2. Gene expression levels and differentially expressed genes To detect prolificacy-related genes, the expression levels of all detected genes were calculated (Fig. 2 and Table S2). One-thousandand-nine genes expressed only in uniparous sample and 858 genes expressed only in multiparous sample in all of the 15,890 detected genes. A total of 2201 genes were differentially expressed between the uniparous and multiparous libraries, in which 1583 genes were upregulated and 618 genes were downregulated in the multiparous library (Table S3).
3.3. Functional enrichment analysis GO contains cellular component, molecular function and biological process domains. Based on sequence homology, the 2201 DEGs identified in our analyses could be categorized into 50 functional groups. In the three main GO classification categories, 23, 15 and 12 functional groups were identified (Fig. 3). Among these groups, the terms cellular process, cell & cell part and binding were predominant in each of the three main categories. A high percentage of genes in the following categories was also observed among the DEGs: biological regulation, cellular component organization or biogenesis, developmental process, establishment of localization, localization, metabolic processes, multicellular organismal process, positive regulation of biological processes, regulation of biological processes, response to stimulus, signaling and single-organism process, macromolecular complex, membrane, membrane part, membrane-enclosed lumen, organelle, organelle part and catalytic activity. KEGG pathway significant enrichment analysis revealed that the DEGs participated in 244 pathways (Table S4) which are known to be enriched in Huntington's disease, endocytosis, lysine degradation, inositol phosphate metabolism, ECM–receptor interaction, focal adhesion, oxidative phosphorylation, Alzheimer's disease, pathways in cancer, the phosphatidylinositol signaling system and regulation of the actin cytoskeleton (Table 3).
3.4. Confirmation of differential gene expression by qPCR Ten DEGs were randomly selected for qPCR analysis to validate the expression profiles obtained by RNA-Seq. The results verified that FER1L4, Itga4, ADAM22, SYCP1, and CDC14A were upregulated and SRD5A2, NMU, COMP, TNC, and WWC1 were downregulated in multiparous samples (Fig. 4), which was basically consistent with the RNA-Seq (Quantification) findings. Thus, RNA-Seq can be used to reliably and accurately perform mRNA differential expression analysis.
Fig. 1. Composition of total raw reads from the uniparous and multiparous libraries. After filtering the only adaptor sequences, which contain N sequences and low quality sequences, the two libraries still generated over 5.8 million clean reads each. The percentage of clean reads among the raw reads reached 98.60% and 98.16% in the uniparous and multiparous libraries, respectively.
4. Discussion 4.1. The functions of the 10 most differentially expressed goat genes found in this study The 10 most differentially expressed goat genes identified in the present study (Table S5) are an important pool of candidate genes for future analysis. EVI2A may form homocomplexes or associate with other proteins within the membrane, as part of a cell-surface receptor complex (Cawthon et al., 1990). The Eyes Absent (EYA) proteins, firstly described in the context of fly eye development, are now implicated in processes as disparate as organ development, innate immunity, DNA damage repair, photoperiodism, angiogenesis, cancer metastasis, and early embryonic development (Ishihara et al., 2008; Tadjuidje and Hegde, 2013). 5 alpha-SR2 (3-oxo-5-alpha-steroid 4dehydrogenase 2, encoded by SRD5A2) converts testosterone (T) into 5-alpha-dihydrotestosterone (DHT) and progesterone or corticosterone into their corresponding 5-alpha-3-oxosteroids. This protein plays a central role in sexual differentiation and androgen physiology (Chang et al., 2003; Seo et al., 2009). Somatostatin (encoded by SST) inhibits the release of somatotropin, and some research has found that somatotropin has effects in animal reproduction (Anderson et al., 1991). Bone sialoprotein 2 (encoded by IBSP) binds tightly to hydroxyapatite and appears to form an integral part of the mineralized matrix. It is possibly important for cell–matrix interactions and likely promotes Arg-GlyAsp-dependent cell attachment (Ganss et al., 1999). BT (butyrophilin
Table 2 Read numbers and mapping results for the two constructed ovary libraries. Sample ID
Total reads
Total mapped reads
Unique match
Multi-position match
Total unmapped reads
Uniparous
6,027,714 (100.00%) 5,884,062 (100.00%)
4,951,157 (82.14%) 4,977,880 (84.60%)
4,071,224 (67.54%) 4,192,641 (71.25%)
879,933 (14.60%) 785,239 (13.35%)
1,076,557 (17.86%) 906,182 (15.40%)
Multiparous
Y.-H. Ling et al. / Gene 550 (2014) 148–153
151
GVINP1 and GVIN1 in the cytosol and nucleus, and the molecular functions of GVINP1 and GVIN1 involve GTP binding. Hence, the function of FER1L4 may be in the exchange of information between the inside and outside of cells, and the functions of GVINP1 and GVIN1 may be in cell energy metabolism, DNA replication and transcription. In view of the high log2Ratio values of these three genes (N10), they may also be associated with high fecundity, and are possibly the key candidate genes to emerge from the current analysis. 4.2. Some of the validated DEGs have previously been associated with high fecundity
Fig. 2. Scatter plot indicating the comparative results of log transformed gene expression levels and differentially expressed gene distributions between uniparous and multiparous samples.
subfamily 1 member A1, encoded by BTN1A1) may have some role in the secretion of milk-fat droplets and act as a specific membraneassociated receptor for the association of cytoplasmic droplets with the apical plasma membrane (Jack and Mather, 1990). CRP3 (cysteine and glycine-rich protein 3) encoded by CSRP3, and the CRP family (CRP1, CRP2, and CRP3) have been implicated in the regulation of muscle development (Xu et al., 2010). CRP3 may be a scaffold protein that promotes the assembly of interacting proteins at Z-line structures (Knoll et al., 2002). In addition, CSRP3 has been confirmed as a HCM disease gene (Geier et al., 2008). In summary, the functions of EVI2A, EYA2, SRD5A2, SST and IBSP involve the exchange of information between cells, DNA repair and hormone secretion which may be closely related to the development of follicular cells and oocytes. It was initially considered that these five genes may be associated with high fecundity in goats. The functions of FER1L4, GVINP1 and GVIN1 are unknown. However, within the cellular component of the GO, FER1L4 is positioned in the membrane and
Genes that emerged from the current screen and have been previously associated with high fecundity were further analyzed. These genes such as GDF9, GnRHR, INH, FSHR, BMPs, and OPNm mainly encode hormones and cytokines. The genes Bmp6 (bone morphogenetic protein 6), INHBB (inhibin beta B), INHA (inhibin, alpha) and OPN (osteopontin) (Koji et al., 2010; Steinheuer et al., 2003) were also identified here as DEGs in the goat, suggesting that they may also be relevant to goat prolificacy (Table S6). The absolute value of the log2Ratio of OPN is relatively high (greater than 3.7) and it may also therefore be a key candidate gene. Other than these four genes, differential expression was non-significant for all of the genes identified in this study. Possible reasons for this include differences between species and an insufficient number of samples. 4.3. Signaling pathways in reproduction At least three of the 11 significantly enriched pathways identified in this study, ECM–receptor interaction (Fig. S1), focal adhesion (Fig. S2) and regulation of actin cytoskeleton (Fig. S3), are intimately associated with mammalian reproduction. The extracellular matrix (ECM) consists of a complex mixture of structural and functional macromolecules and plays an important role in tissue and organ morphogenesis. Specific interactions between cells and the ECM lead to the direct or indirect control of cellular activities such as adhesion, migration, differentiation, proliferation, and apoptosis (Bosman and Stamenkovic, 2003). Focal adhesion plays essential roles in important biological processes including cell motility, cell proliferation, cell differentiation, regulation of gene expression and cell survival (Kanehisa et al., 2008). Regulation of the dynamic actin cytoskeleton includes signaling to the cytoskeleton which
Fig. 3. Histogram of gene ontology classifications. The results were summarized in three main categories: biological process, cellular component and molecular function. Among these groups, the terms cellular process, cell & cell part and binding were predominant in each of the three main categories, respectively.
152
Y.-H. Ling et al. / Gene 550 (2014) 148–153
Table 3 Significantly enriched pathways among the differentially expressed genes between the uniparous and multiparous samples. Pathway
DEGs with pathway annotation (1825)
All genes with pathway annotation (16,032)
P value
Q value
Pathway ID
Huntington's disease Endocytosis Lysine degradation Inositol phosphate metabolism ECM–receptor interaction Focal adhesion Oxidative phosphorylation Alzheimer's disease Pathways in cancer Phosphatidylinositol signaling system Regulation of actin cytoskeleton
55 (3.01%) 73 (4%) 26 (1.42%) 24 (1.32%) 38 (2.08%) 70 (3.84%) 32 (1.75%) 44 (2.41%) 84 (4.6%) 31 (1.7%) 80 (4.38%)
253 (1.58%) 385 (2.4%) 99 (0.62%) 97 (0.61%) 190 (1.19%) 412 (2.57%) 153 (0.95%) 233 (1.45%) 518 (3.23%) 149 (0.93%) 514 (3.21%)
1.511431e −06 7.280254e −06 3.141498e −05 0.0001715156 0.0003662485 0.0003790018 0.0004581824 0.0004921621 0.0004933844 0.0006124814 0.002192947
0.0004141321 0.0009973948 0.0028692348 0.0117488186 0.0168984157 0.0168984157 0.0133761993 0.0168984157 0.0168984157 0.0186466560 0.0486435516
ko05016 ko04144 ko00310 ko00562 ko04512 ko04510 ko00190 ko05010 ko05200 ko04070 ko04810
can lead to diverse effects on cellular activity such as changes in cell shape, migration and proliferation (Lee and Dominguez, 2010). The cell will migrate and move towards the adluminal compartment to prepare for meiosis if the actin cytoskeleton dynamics have changed appropriately (Zhang et al., 2013). These signaling pathways are closely associated with cell differentiation, proliferation, and meiosis during the development of follicles. Most of the DEGs (22/38, 50/70 and 59/80, respectively) identified in this study which participate in the abovementioned three pathways were up-regulated. This indicates that the three pathways are activated in the estrus ovaries of multiparous goats. As a whole therefore, these genes may promote the development of follicles and oocytes but the mechanism is unclear. Previous studies have indicated that pathways involving the DEGs in this study also associate with reproduction except the Notch signaling pathway (Larkin et al., 1996; Lopez-Schier and St Johnston, 2001), TGF-beta signaling pathway (Shi and Massagué, 2003), steroid biosynthesis (Issop et al., 2013), MAPK signaling pathway (Xiao et al., 2007), Wnt signaling pathway (Habas and He, 2006; Ishida-Takagishi et al., 2012) and Jak–STAT signaling pathway (Kisseleva et al., 2002). However, the exact relationship between these signaling networks is not fully understood.
5. Conclusion This study screened for DEGs in the ovarian tissues of the multiparous and uniparous Anhui White goats using RNA-Seq (Quantification), and 2201 DEGs were thereby obtained. Among the pathways in which these DEGs are mainly involved, at least three are associated with reproduction, and 12 of the DEGs may be relevant to the prolificacy of goats. This new information provides a solid foundation for further studies of the molecular mechanisms underlying goat prolificacy.
Fig. 4. Validation of the RNA-Seq (Quantification) results by qPCR. Results represent the mean (±SD) of three experiments. Asterisk, P ≤ 0.05.
Further biochemical and physiological analyses of these candidate genes will be conducted in the future. Acknowledgments This research was supported by the National Natural Science Foundation of China (31301934, 31372310), the Natural Science Foundation of Anhui Province (1308085QC54) and Anhui Science and Technology Program (11010302108) funds. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2014.08.008. References Anderson, L.L., Ford, J.J., Klindt, J., et al., 1991. Growth hormone and prolactin secretion in hypophysial stalk-transected pigs as affected by growth hormone and prolactinreleasing and inhibiting factors. Proc. Soc. Exp. Biol. Med. 196, 194–202. Ashburner, M., Ball, C.A., Blake, J.A., et al., 2000. Gene ontology: tool for the unification of biology. Nat. Genet. 25, 25–29. Benjamini, Y.,Drai, D., Elmer, G., et al., 2001. Controlling the false discovery rate in behavior genetics research. Behav. Brain Res. 125, 279–284. Bosman, F.T.,Stamenkovic, I., 2003. Functional structure and composition of the extracellular matrix. J. Pathol. 200 (4), 423–428. Cawthon, R.M., O'Connell, P., Buchberg, A.M., et al., 1990. Identification and characterization of transcripts from the neurofibromatosis 1 region: the sequence and genomic structure of EVI2 and mapping of other transcripts. Genomics 7 (4), 555–565. Chang, B.L.,Zheng, S.L.,Isaacs, S.D., et al., 2003. Evaluation of SRD5A2 sequence variants in susceptibility to hereditary and sporadic prostate cancer. Prostate 56 (1), 37–44. Chen, S.,Cheng, G.L.,Zhu, D.J., et al., 2009. Determination on the body properties and meat performance of Anhui white goat. Anim. Husb. Feed Sci. 30 (4), 150–152. Davis, G.H., Galloway, S.M., Ross, I.K., et al., 2002. DNA tests in prolific sheep from eight countries provide new evidence on origin of the Booroola (FecB) mutation. Biol. Reprod. 66 (6), 1869–1874. Galloway, S.M.,McNatty, K.P.,Cambridge, L.M., et al., 2000. Mutations in an oocyte-derived growth factor gene (BMP15) cause increased ovulation rate and infertility in a dosage-sensitive manner. Nat. Genet. 25 (3), 279–283. Ganss, B., Kim, R.H., Sodek, J., 1999. Bone sialoprotein. Crit. Rev. Oral Biol. Med. 10 (1), 79–98. Geier, C.,Gehmlich, K.,Ehler, E., et al., 2008. Beyond the sarcomere: CSRP3 mutations cause hypertrophic cardiomyopathy. Hum. Mol. Genet. 17 (18), 2753–2765. Gunawan, A., Sahadevan, S., Neuhoff, C., et al., 2013. RNA deep sequencing reveals novel candidate genes and polymorphisms in boar testis and liver tissues with divergent androstenone levels. PLos One 8 (5), e63259. Habas, R., He, X., 2006. Activation of Rho and Rac by Wnt/frizzled signaling. Methods Enzymol. 406, 500–511. Hackett, N.R., Butler, M.W., Shaykhiev, R., et al., 2012. RNA-Seq quantification of the human small airway epithelium transcriptome. BMC Genomics 13, 82. Hanrahan, J.P., Gregan, S.M., Mulsant, P., et al., 2004. Mutations in the genes for oocytederived growth factors GDF9 and BMP15 are associated with both increased ovulation rate and sterility in Cambridge and Belclare sheep (Ovis aries). Biol. Reprod. 70, 900–909. Huang, W.,Khatib, H., 2010. Comparison of transcriptomic landscapes of bovine embryos using RNA-Seq. BMC Genomics 11, 711. Ishida-Takagishi, M.,Enomoto, A.,Asai, N., et al., 2012. The Dishevelled-associating protein Daple controls the non-canonical Wnt/Rac pathway and cell motility. Nat. Commun. 3, 859. Ishihara, T.,Ikeda, K., Sato, S., 2008. Differential expression of Eya1 and Eya2 during chick. Gene Expr. Patterns 8 (5), 357–367.
Y.-H. Ling et al. / Gene 550 (2014) 148–153 Issop, L., Rone, M.B., Papadopoulos, V., 2013. Organelle plasticity and interactions in cholesterol transport and steroid biosynthesis. Mol. Cell. Endocrinol. 371 (1–2), 34–46. Jack, L.J., Mather, I.H., 1990. Cloning and analysis of cDNA encoding bovine butyrophilin, an apical glycoprotein expressed in mammary tissue and secreted in association with the milk-fat globule membrane during lactation. J. Biol. Chem. 265 (24), 14481–14486. Kanehisa, M.,Araki, M.,Goto, S., et al., 2008. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36 (Suppl. 1), D480–D484. Kerpedjiev, P., Frellsen, J., Lindgreen, S., et al., 2014. Adaptable probabilistic mapping of short reads using position specific scoring matrices. BMC Bioinforma. 15, 100. Kisseleva, T., Bhattacharya, S., Braunstein, J., et al., 2002. Signaling through the JAK/STAT pathway, recent advances and future challenges. Gene 285 (1–2), 1–24. Knoll, R.,Hoshijima, M.,Hoffman, H.M., et al., 2002. The cardiac mechanical stretch sensor machinery involves a Z disc complex that is defective in a subset of human dilated cardiomyopathy. Cell 111, 943–955. Koji, S., You, Q.S., John, J.E., 2010. Does bone morphogenetic protein 6 (BMP-6) affect female fertility in the mouse? Biol. Reprod. 83 (6), 997–1004. Larkin, M.K., Holder, K., Yost, C., et al., 1996. Expression of constitutively active Notch arrests follicle cells at a precursor stage during Drosophila oogenesis and disrupts the anterior–posterior axis of the oocyte. Development 122 (11), 3639–3650. Lee, S.H., Dominguez, R., 2010. Regulation of actin cytoskeleton dynamics in cells. Mol. Cells 29, 311–325. Li, R.Q., Yu, C., Li, Y.R., et al., 2009. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25 (15), 1966–1967. Lindgreen, S., 2012. AdapterRemoval: easy cleaning of next-generation sequencing reads. BMC Res. Notes 5, 337. Lopez-Schier, H.,St Johnston, D., 2001. Delta signaling from the germ line controls the proliferation and differentiation of the somatic follicle cells during Drosophila oogenesis. Genes Dev. 15 (11), 1393–1405. Mortazavi, A.,Williams, B.A.,McCue, K., et al., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5 (7), 621–628.
153
Seo, Y.K., Zhu, B., Jeon, T.I., et al., 2009. Regulation of steroid 5-alpha reductase type 2 (Srd5a2) by sterol regulatory element binding proteins and statin. Exp. Cell Res. 315 (18), 3133–3139. Shi, Y.G., Massagué, J., 2003. Mechanisms of TGF-β signaling from cell membrane to the nucleus. Cell 113 (6), 685–700. Steinheuer, R., Drögemüller, C., Hamann, H., et al., 2003. Possible uses of genetic markers for improving fertility and health in swine production. Dtsch. Tierarztl. Wochenschr. 110 (6), 255–265. Sun, W., Chang, H., Musa, H.H., et al., 2010. Study on relationship between microsatellite polymorphism and producing ability on litter size trait of Hu sheep in China. Afr. J. Biotechnol. 9 (50), 8704–8711. Tadjuidje, E.,Hegde, R.S., 2013. The Eyes Absent proteins in development and disease. Cell. Mol. Life Sci. 70 (11), 1897–1913. Toosi, B.M., Seekallu, S.V., Rawlings, N.C., 2010. Effects of the rate and duration of physiological increases in serum FSH concentrations on emergence of follicular waves in cyclic ewes. Biol. Reprod. 83 (4), 648–655. Xiao, Y.,Yang, F.Q.,Li, S.P., et al., 2007. Furanodiene induces G2/M cell cycle arrest and apoptosis through MAPK signaling and mitochondria–caspase pathway in human hepatocellular carcinoma cells. Cancer Biol. Ther. 6 (7), 1044–1050. Xu, X.W., Qiu, H.F., Du, Z.Q., et al., 2010. Porcine CSRP3: polymorphism and association analyses with meat quality traits and comparative analyses with CSRP1 and CSRP2. Mol. Biol. Rep. 37 (1), 451–459. Ye, J., Fang, L., Zheng, H.K., et al., 2006. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 34 (Suppl. 2), W293–W297. Zhang, X.J.,Hao, L.L.,Meng, L.J., et al., 2013. Digital gene expression tag profiling analysis of the gene expression patterns regulating the early stage of mouse spermatogenesis. PLos One 8 (3), e58680. Zhou, X., Su, Z., 2007. EasyGO: Gene Ontology-based annotation and functional enrichment analysis tool for agronomical species. BMC Genomics 8, 246.