Identification of Differentially Expressed Long Noncoding RNAs and mRNAs Involved with Dominant Follicle Selection in Goats using RNA-seq

Hide Journal Menu

Current Issue

PJZ_50_1_47-56

Identification of Differentially Expressed Long Noncoding RNAs and mRNAs Involved with Dominant Follicle Selection in Goats using RNA-seq

Guang-Xin E1, Yong-Ju Zhao1, Yue-Hui Ma2, Ming-Xing Chu2, Jia-Hua Zhang1, Zhong-Quan Zhao1, Hui-Jiang Gao2, Huai-Zhi Jiang3, Di Liu4, Li Liu5, Yan-Bin Zhu6, Wang-Dui Basang6, Luo-Bu Danjiu7, Tian-Wu An8, Xiao-Lin Luo8, Shi-Cheng He7 and Yong-Fu Huang1,*

1Chongqing Key Laboratory of Forage & Herbivore, Chongqing Engineering Research Centre for Herbivores Resource Protection and Utilization, College of Animal Science and Technology, Southwest University, Chongqing 400716, China

2Institute of Animal Science, Chinese Academy of Agricultural Sciences (CAAS), Beijing 100193, China

3College of Animal Science and Technology, Jilin Agricultural University, Changchun, Jilin, China

4Institute of Animal Husbandry, Heilongjiang Academy of Agricultural Science, Harbin 150086, China

5College of Animal Science and Technology, China Agricultural University, Beijing 100083, China

6Tibet Academy of Agriculture and Animal Husbandry Science, Lasa 850001, China

7Nagqu grassland station, Naqu 852000, China

8Sichuan Academy of Grassland Sciences, Chengdu, Sichuan 611731, China

Guang-Xin E and Yong-Ju Zhao contributed equally in this article.

ABSTRACT

In this study, we used high-throughput technology to provide the first transcriptome dataset for differentially expressed genes in mixed pools of dominant and nondominant follicles of goats. These data will contribute to research on the molecular mechanisms of dominant follicle selection in goats. In this study, 90276370 and 115579236 clear reads in dominant and nondominant follicles of goat were generated through Illumina paired-end sequencing, and their mapping rate was 84.99% and 84.47%, respectively. A total of 12577 differentially expressed genes (DEGs) were identified, including 6009 upregulated and 6568 downregulated genes in dominant follicles compared with nondominant follicles of goat. Of the 1026 significantly differentially expressed, long noncoding RNAs (lncRNAs) found, 419 were upregulated and 607 were downregulated. The DEGs related to 56 GO categories, and pathway analysis revealed that these DEGs were significantly enriched in 41 of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, including signaling pathways regulating the pluripotency of stem cells, the p53 signaling pathway, and oxidative phosphorylation. The results of the present study confirmed that the selection of the dominant follicle involved the regulation of various physiological systems. These results provided helpful data to understand on the molecular mechanisms of dominant follicle selection in goats.

Article Information

Received 25 July 2017

Revised 02 September 2017

Accepted 07 October 2017

Available online 19 December 2017

Authors’ Contribution

Y-FH, G-XE, Z-QZ, DL, LL and H-ZJ participated in the experimental design and wrote the manuscript. Y-JZ, M-XC, Y-HM, J-HZ, H-JG, W-DB, Y-BZ, L-BD, T-WA, Z-LL and S-CH performed the RNAseq experiment, analysed the bioinformatic data and approved the final manuscript.

Key words

Long noncoding RNA, RNA-seq, Dominant follicle, Goat.

DOI: https://dx.doi.org/10.17582/journal.pjz/2018.50.1.47.56

* Corresponding author: [email protected]

0030-9923/2018/0001-0047 $ 9.00/0

Abbreviations

DF, dominant follicle; NF, nondominant follicles; RNAseq, transcriptome shotgun sequencing; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; lncRNA, long non-coding RNA; mRNA, messenger RNA; DEGs, differentially expressed genes; FDR, false discovery rate; FRKM, kilobase of exon per million mapped reads; qPCR, quantitative PCR; OXPHOS, mitochondrial oxidative phosphorylation; SDHC, succinate dehydrogenase complex subunit C; Cox11, cytochrome c oxidase assembly protein 11; Cyt1, cystatin; NDUFS1, NADH-ubiquinone oxidoreductase Fe–S protein 1; BMP15, bone morphogenetic protein 15; GDF9, growth differentiation factor 9; IGFBP5, insulin-like growth factor binding protein 5; BMPR2, bone morphogenetic protein receptor type 2; BAX, BCL2 associated X protein; BCL-2, B-cell CLL/lymphoma 2; BMP2, bone morphogenetic protein 2; INHA, inhibiting the alpha subunit of inhibin; PTGS1, prostaglandin-endoperoxide synthase 1; PTGS-2, prostaglandin-endoperoxide synthase 2.

Introduction

Recent advances in high-throughput technology have led to the demonstration that lncRNA is actively transcribed in the eukaryotic genome and plays an important role in all aspects of normal biological function, particularly in model animals (Yue et al., 2016).

Follicle selection, the process giving rise to the dominant follicle (DF), and the physiological state of the DF are important processes to study, especially the selection of an ovarian follicle for further differentiation and finally ovulation in an ovulating animal. However, the transcriptional regulation of ovarian follicles is very complex, and large differences have been observed during different developmental stages and in different breeds (Zhao et al., 2015). Therefore, the precise mechanisms regulating the growth of follicles and ovulation are still poorly understood.

In this study, we investigated the differential expression of lncRNAs and mRNAs in different development periods of follicles using RNA-seq technology. The data provide a large amount of useful information about RNAs that are related to mammal reproductive processes, and they help to understand the importance of different transcriptional factors in follicular development.

Materials and Methods

Animals

The experimental conditions of this study were approved by the Committee on the Ethics of Animal Experiments of the Southwest University (No. [2007] 3) and the Animal Protection Law of China. Ovaries of three samples of Dazu black goats were collected from domestic animal conservation field of the Southwest University. The follicles of ovaries were separated into DF (diameter > 5 mm) and nondominant follicles (NF, 3mm < diameter < 5mm), the detail standard protocol, animal description as Supplementary Material.

RNA extraction, library preparation and genome-wide resequencing

Total RNA was extracted using TRIzol® Reagent (Invitrogen, USA). Equal amounts of RNA from three different goats and the same follicle developmental phases (dominant and nondominant) were pooled. Libraries were generated using the rRNA-depleted RNA of the NEBNext Ultra Directional RNA Library Prep Kit (NEB, Ipswich, MA, USA). Libraries were sequenced on the Illumina HiSeq 2500 platform by Gene Denovo Technologies (Guangzhou, China). The detail libraries optiation and informatic analysis was displyed as Supplementary Material.

Quantitative PCR of candidate RNA factors

The samples used in the q-PCR analyses were the same as those used in the sequencing experiments. The cDNA was synthesized using the First Strand cDNA Synthesis Kit (GE Healthcare) and 1 mg total RNA. The primers are shown in Table I. Quantitative RT-PCR was performed in triplicate using 10 ml SYBR® Select Master Mix (LifeTechnologies), 6.4 μL H2O, 0.8 ml primer (10 pmol/ml) and 2 μL cDNA (~16 ng). The reaction was performed on an Applied Biosystems StepOnePlusTM Real-Time PCR System (LifeTechnologies). The cycle threshold (Ct) values were normalized to the control gene (GAPDH). The relative quantification was calculated using the 2-ΔΔCt method (Schmittgen and Livak, 2008). Each experiment was independently repeated three times, and each sample was evaluated in triplicate. Using the t test, a P value less than 0.05 was considered significant.

Results

Illumina sequencing and gene annotation

Illumina sequencing of C. hircus yielded a total of 11401202.75 kb and 14600735 kb raw data from the mixed pools of DF and NF, respectively. The transcriptome sequencing data from the mixed pools of goat follicles were deposited in the NCBI Sequence Read Archive database (accession number, SRA342894). A reference dataset, which included protein data from the goat genome (CHIR_1.0), was constructed for gene annotation. A total of 90276370 kb (DF) and 115579236 kb (NF) clear reads were generated by Illumina paired-end sequencing, and their mapping rate was 84.99% (DF) and 84.47% (NF), respectively. These results demonstrated that both libraries were high quality. Transcripts were blasted against this dataset. Of the total 20763 (DF) and 21821 (NF) transcripts, those with known isoform numbers were 15756 (51.34%, DF) and 16511 (53.80%, NF).

Differential gene expression (DGE) of LncRNA and mRNA between DF and NF

In the analysis of the digital coding gene expression of the follicle tissue in the DF and NF of goats, we found 12,577 significantly differentially expressed genes, including 6,009 genes from the DF that were upregulated compared to the NF and 6,568 that were downregulated. Differential lncRNA expression in the tissues of the DF and NF from goats, we found 1,026 significant differentially expressed lncRNA, including 419 that were upregulated and 607 that were down regulated.

Functional annotation of the DEGs

A GO analysis was performed to classify the functions of DEGs that had hits in the NCBI NR database. These were summarized under three main GO categories, cellular component, molecular function, and biological process (Fig. 1) and were associated with 56 GO terms in total, which could then be used to implicate the biological processes involved in DF selection.

Table I.- Primer information of candidate gene using q-PCR.

Target	Primer name	Sequences (5’-3’)	Ampicon size (bp)
BAX#	BAX-F	GCATCCACCAAGAAGCTGAG	130
	BAX-R	CCGCCACTCGGAAAAAGAC
BCL2#	BCL2-F	ATGTGTGTGGAGAGCGTCA	182
	BCL2-R	AGAGACAGCCAGGAGAAATC
BMPR1B#	BMPR1B-F	CCTGTGGTCACTTCTGGATGTC	107
	BMPR1B-R	TTCCGTTCTGTGCAGCATTC
BMP15@	BMP15-F	TGAGGCCGCTGGCTAGTG	147
	BMP15-R	GGGAATGAGTTAGGTGAAGCTGAT
GDF9@	GDF9-F	ACAGACGCCACCTCTACAACACT	136
	GDF9-R	TTCCACAACAGTAACACGATCCA
IGFBP5#	IGFBP5-F	AGAGAGACTCCCGTGAGCAT	159
	IGFBP5-R	ACGAACTTGGACTGGGTCAG
INHA#	INHA-F	CTTCCCTCTGCTGACCCATC	250
	INHA-R	ATTGAGGGCGGCTCTGTG
BMP2#	BMP2-F	A TCACCTGAACTC CACGAA	140
	BMP2-R	TACCACCTTCTCATTCTCATC
INSIG1#	INSIG1-F	GTGGGGAACATAGGACGAC	unknown
	INSIG1-R	ACGAGTCATTTGTACAGTCAGCCC
BMP6@	BMP6-F	CGCCTCAGACTACAACAGCA	165
	BMP6-R	TTCATGTGTGCGTTGAGAGG
BMPRII@	BMPRII-F	AGACCCAAGTTCCCAGAAGC	96
	BMPRII-R	AGCCTCTGCATCCTGGTC
TGFBR1@	TGFBR1-F	CGGAAAGCCGTCATCTGGCCTC	164
	TGFBR1-R	CTCGATGGTGAATGACAGTGCGGT
PTGS-2#	Ptgs2-F	AGGAGGTCTTTGGTCTGGTG	126
	Ptgs2-R	TCTGGAACAACTGCTCATCG
GAPDH	GAPDH-F	AGGCTGGGGCTCACTTGAAG	223
	GAPDH-R	ATGGCGTGGACAGTGGTCAT
LOC10218 3322#	LOC102183 322-F2	CGGGAGGTTACTGTAGGGC	189
	LOC102183 322-R2	TTGCTTAGTTGTCGTCGTGC
LOC10218 1730#	LOC10218 1730-F2	ACTGCGAGGACATAACACCA	207
	LOC10218 1730-R2	CCACTTTCTTCCCCGCTTTC
LOC10218 4274#	LOC10218 4274-F2	AGAAGTACAACCCCTCTGCC	230
	LOC10218 4274-R2	GTCTTCCTTGGGCCTTTTCG

Note: Gene name with the upper standard. #, identified significant differences in the expression of genes in this study using transcriptome re-sequence; @, the known candidate gene from previous study (Lima et al., 2012; Paradis et al., 2009).

The KEGG database was used to identify potentially important biological pathways in our dataset. A total of 4,318 genes were assigned to 292 KEGG pathways. Furthermore, most of the KEGG pathways involved immune function (e.g., the B-cell receptor signaling pathway, the autoimmune thyroid disease pathway, the Toll like receptor signaling pathway, and reproductive signaling pathways, including the PPAR and GnRH signaling pathways. Of the 292 KEGG pathways, 41 were enriched (P < 0.05). These annotations provide a valuable resource for investigating specific processes, functions, and pathways in DF selection. In addition, the 20 top pathway enrichments, including signaling pathways regulating the pluripotency of stem cells, the p53 signaling pathway, and the oxidative phosphorylation pathway, were identified (Fig. 2). Some reproductive related pathways were identified, including oocyte meiosis, ovarian steroidogenesis, progesterone-mediated oocyte maturation, and the estrogen signaling pathway.

Differential gene expression profiling and qRT-PCR validation

To identify the differentially expressed transcription factors (DETs) associated with DF and NF from the two mixed pools in the current study, we selected some DEGs for real time quantitative PCR (q-PCR) (Fig. 3). Using this procedure, we found similar gene expression profiles for INHA, IGFBP5, BMPR1B, BMP2, INSIG, PTGS2, LOC102181730 and LOC102184274, which indicated the high credibility of these genes in transcript abundance; however, some factors in the RNA-seq results, including BAX, LOC102183322, and BCL2, were inconsistent (Fig. 4).

Discussion

A total of 1026 lncRNAs (884 novel lncRNAs) and 12577 coding mRNAs with significant differential expression were detected in DF compared to NF. These transcripts were assigned to a total of 56 GO categories. The largest proportion was biological process (22), followed by cellular process (19) and molecular function (15). The results of the present study confirmed that the initial follicle was selected to be the DF by the regulation of various physiological systems

Differentially expressed genes were involved in 292 pathways, of which the 20 most enriched were predicted using KEGG analysis. Including the signaling pathways regulating the pluripotency of stem cells, and research has revealed that the stem cell factor promotes the development of ovarian follicles in vitro (Thuwanut et al., 2016). The p53 signaling pathway, which is a known factor in the regulation of ovarian function, such as the proliferation and the apoptosis of ovarian cells (Sirotkin et al., 2014). The HIF-1 signaling pathway, has a key role in the development of ovarian tissue, including luteal development (Yu et al., 2015) and ovarian follicular growth (Zhang et al., 2015). In addition, it has indicated that the Hippo signaling pathway has a spatio-temporal correlation with the size of the primordial follicle pool (Xiang et al., 2015).

Ovarian follicles enter a massive growth phase during which they become highly dependent on gonadotrophic factors and nutrients. At the same time, the FoxO-dependent process causes prompt and efficient adaptation to nutrient supplies, thereby contributing to maintaining the balance between reproduction and nutrient availability (Jouandin et al., 2014). In addition, there is evidence that indicates that the FoxO pathway factor plays a role in follicular development (Li et al., 2014) and in female reproductive disorders (Christian et al., 2011). Thus, the FoxO pathway may be important for the regulation of DF selection, too.

However, some pathways were not enriched in this study but may be related to the regulation of DF selection. For example, an interesting pathway, mitochondrial oxidative phosphorylation (OXPHOS), takes place inside mitochondria and is highly efficient in releasing energy for biological action. Many studies have found that mitochondrial function is related to reproductive biology (Ben-Meir et al., 2015), particularly the OXPHOS pathway is inseparably related to oocyte maturation, fertilization and embryo development (Harvey et al., 2002; Ge et al., 2013). In this study, some genes in the mitochondrial OXPHOS pathway, such as succinate dehydrogenase complex subunit C (SDHC), cytochrome c oxidase assembly protein 11 (Cox11), and cystatin (Cyt1) were differentially expressed.

In addition, the target gene prediction for lncRNAs showed that lncRNA LOC102181730 and its potential target gene NDUFS1 were significantly differentially expressed between dominant and non-dominant follicles. In addition, lncRNA LOC102183322 and its potential target gene NDUFB6 were also significantly expressed in the DF pool. The NADH-ubiquinone oxidoreductase Fe–S protein 1 (NDUFS1) and the NADH-ubiquinone oxidoreductase subunit B6 (NDUFB6) belong to NADH dehydrogenase (ubiquinone) is related to oxidative-phosphorylation (OXPHOS). Some studies have found that a variant in NDUFS1 resulted in a complex I deficiency (Björkman et al., 2015). NDUFB6 is required for complex I activity, and it defines conditions suitable for a systematic and stable exclusion of the different supernumerary subunits in human cells (Loublier et al., 2011). These results indicate that the interaction of LOC102183322 and LOC102181730 with NDUFS1 and NDUFB6 might regulate OXPHOS. Here, we believe that this biological process plays a role in DF selection.

Some candidate genes were identified from previous study were controversy, including growth differentiation factor 9 (GDF9) promotes growth of the oocyte at the primary (Cook-Andersen et al., 2016; Kona et al., 2016; Pramod et al., 2013). In addition, numerous studies have revealed that bone morphogenetic protein 15 (BMP15) is differentially expressed in small and large antral follicles during the development in vitro of cultured preantral follicles (Pramod et al., 2013; Lima et al., 2012). However, some reported these genes (BMP15 or GDF9) to be inactive in reproductive functions because no differences were identified that were associated with enhanced ovulation rate (Feary et al., 2007), it was consistent with this study.

In previous study, it has been shown that the bone morphogenetic protein receptor type 2 (BMPR2), the insulin-like growth factor binding protein 5 (IGFBP5), the BCL2 associated X protein (BAX) and the B-cell CLL/lymphoma 2 (BCL-2) protein have different transcription patterns in different stages of follicles cultured in vitro (Lima et al., 2012; Yang et al., 2012; DeBem et al., 2014). These patterns are consistent with the expression patterns observed in the current study. However, the BMP6, BMP15 and GDF9 mRNAs were most abundant in the oocyte, and their expression remained relatively constant during follicular development, whereas the BMPR1B and TGFBR1 were temporally regulated in different stages in pigs (Paradis et al., 2009). However, the expression of GDF9, BMP15, BMP6, and TGFBR1 exhibited some differential expression in our study, and BMPR1B showed significantly different expression in DF and NF. It indicated that the regulated transcription in the selection of DF was species-specific.

Bone morphogenetic protein 2 (BMP2), which showed significantly different expression in DF and NF in this study, has been shown to regulate primordial follicle formation by promoting the germ cell to oocyte transition and the somatic cell to pre-granulosa cell formation (Chakraborty and Roy , 2015). Inhibiting the alpha subunit of inhibin (INHA), which is related to infertility (Rah et al., 2014) and ovarian insufficiency (Li et al., 2015). Prostaglandin-endoperoxide synthase 2 (PTGS-2) and prostaglandin-endoperoxide synthase 1 (PTGS1), also known as the cyclooxygenases, are key enzymes in prostaglandin biosynthesis and act both as dioxygenases and as peroxidases. Both genes showed significant differential expression in this study. Therefore, BMP2, INHA and PTGS-2 would be novel candidate genes for DF selection in goats.

The results of q-PCR that were not complete consistent with the results of RNA-seq (Consistency rate = 72.72%) indicated that this study design was not optimal. The lack of replication was due to using two mixed pooled samples for the transcriptome sequencing analysis, resulting in a certain number of false positives. Despite such limitations, our results enhance the understanding of the regulation of the selection of the DF and lay the foundation for future studies.

CONCLUSION

In this study, the lncRNA and mRNA profiles of DF selection in goats were investigated using RNAseq. A total of 1026 lncRNAs and 12577 coding RNAs were detected that were differentially expressed in the two mixed pools of follicles. Further functional annotation analysis indicated that some significant biological categories and signaling pathways were involved in DF selection. In addition, some candidate RNAs that provide a new perspective on the mechanisms regulating this physiological process were identified by q-PCR. In short, the results of the current study provide a valuable basis for understanding the molecular mechanisms of DF selection in goats.

ACKNOWLEDGMENTS

This study was supported by the Fundamental Research Funds for the Central Universities (No. XDJK2016B003), the National Natural Science Foundation of China (No. 31172195), Fundamental Research Funds for the Central Universities (No. SWU114023) and the 2013 Innovation Team-Building Program in Chongqing Universities (KJTD201334).

Supplementary material

There is supplementary material associated with this article. Access the material online at: https://dx.doi.org/10.17582/journal.pjz/2018.50.1.47.56

Statement of conflict of interest

The authors declare no conflict of interest.

References

Ben-Meir, A., Burstein, E., Borrego-Alvarez, A., Chong, J., Wong, E., Yavorska, T., Naranian, T., Chi, M., Wang, Y., Bentov, Y., Alexis, J., Meriano, J., Sung, H.K., Gasser, D.L., Moley, K.H., Hekimi, S., Casper, R.F. and Jurisicova, A., 2015. Coenzyme Q10 restores oocytemitochondrial function and fertility during reproductive aging. Aging Cell, 14: 887-895. https://doi.org/10.1111/acel.12368

Björkman, K., Sofou, K., Darin, N., Holme, E., Kollberg, G., Asin-Cayuela, J., Holmberg-Dahle, K.M., Oldfors, A., Moslemi, A.R. and Tulinius, M., 2015. Broad phenotypic variability in patients with complex I deficiency due to mutations in NDUFS1 and NDUFV1. Mitochondrion, 21: 33-40. https://doi.org/10.1016/j.mito.2015.01.003

Chakraborty, P. and Roy, S.K., 2015. Bone morphogenetic protein 2 promotes primordial follicle formation in the ovary. Scient. Rep., 5: 12664. https://doi.org/10.1038/srep12664

Christian, M., Lam, E.W., Wilson, M.S. and Brosens, J.J., 2011. FOXO transcription factors and their role in disorders of the female reproductive tract. Curr. Drug Targets, 12: 1291-1302. https://doi.org/10.2174/138945011796150253

Cook-Andersen, H., Curnow, K.J., Su, H.I., Chang, R.J. and Shimasaki, S., 2016, Growth and differentiation factor 9 promotes oocyte growth at the primary but not the early secondary stage in three-dimensional follicle culture. J. Assist. Reprod. Genet, 33: 1067-1077. https://doi.org/10.1007/s10815-016-0719-z

De Bem, T.H., Adona, P.R., Bressan, F.F., Mesquita, L.G., Chiaratti, M.R., Meirelles, F.V. and Leal, C.L., 2014. The influence of morphology, follicle size and Bcl-2 and Bax transcripts on the developmental competence of bovine oocytes. Reprod. Domest Anim., 49: 576-583. https://doi.org/10.1111/rda.12325

Feary, E.S., Juengel, J.L., Smith, P., French, M.C., O’Connell, A.R., Lawrence, S.B., Galloway, S.M., Davis, G.H. and McNatty, K.P., 2007. Patterns of expression of messenger RNAs encoding GDF9, BMP15, TGFBR1, BMPR1B, and BMPR2 during follicular development and characterization of ovarian follicular populations in ewes carrying the Woodlands FecX2W mutation. Biol. Reprod., 77: 990-998. https://doi.org/10.1095/biolreprod.107.062752