A correction has been applied to this article in:
Corrigendum: Identification of mRNAs and lncRNAs Involved in the Regulation of Follicle Development in Goat
Zhifeng Zhao1,3†
Xian Zou2†
Tingting Lu1Ming Deng1
Yaokun Li1Yongqing Guo1
Baoli Sun1
Guangbin Liu1*
Dewu Liu1*- 1College of Animal Science, South China Agricultural University, Guangzhou, China
- 2State Key Laboratory of Livestock and Poultry Breeding, Guangdong Key Laboratory of Animal Breeding and Nutrition, Institute of Animal Science, Guangdong Academy of Agricultural Sciences, Guangzhou, China
- 3National Engineering Laboratory for Animal Breeding, College of Animal Science and Technology, China Agricultural University, Beijing, China
Follicular development and maturation has a significant impact on goat reproductive performance, and it is therefore important to understand the molecular basis of this process. The importance of long non-coding RNAs (lncRNAs) in mammalian reproduction has been established, but little is known about the roles of lncRNAs in different follicular stages, especially in goats. In this study, RNA sequencing (RNA-seq) of large follicles (>10 mm) and small follicles (<3 mm) of Chuanzhong black goats was performed to investigate the regulatory mechanisms of lncRNAs and mRNAs in follicular development and maturation. A total of 8 differentially expressed lncRNAs (DElncRNAs) and 1,799 DEmRNAs were identified, and the majority of these were upregulated in small follicles. MRO, TC2N, CDO1, and NTRK1 were potentially associated with follicular maturation. KEGG pathway analysis showed that the DEmRNAs involved in ovarian steroidogenesis (BMP6, CYP11A1, CYP19A1, 3BHSD, STAR, LHCGR, and CYP51A1) and cAMP signaling play roles in regulating follicular maturation and developmental inhibition respectively. Five target pairs of DElncRNA-DEmRNA, namely, ENSCHIT00000001255-OTX2, ENSCHIT00000006005-PEG3, ENSCHIT00000009455-PIWIL3, ENSCHIT00000007977-POMP, and ENSCHIT00000000834-ACTR3 in co-expression analysis provide a clue in follicular development and maturation of lncRNA-mRNA interaction. Our findings provide a valuable resource for lncRNA studies, and could potentially provide a deeper understanding of the genetic basis and molecular mechanisms of goat follicular development and maturation.
Introduction
In domestic animals such as cattle, sheep, and goats, follicles develop in the form of wave-like pattern during oestrous cycles (Evans, 2003). In each follicular wave, follicles in the early stage grow slowly (Fortune, 1994). However, with the antrum formation, the process speeded up for recruiting qualified follicles to ovulate. Most follicles in follicular wave undergo atresia and only a few dominant follicles develop to ovulate. Many studies have been done on the mechanisms of follicular maturation and atresia. While the question of how developmental follicles go into different fates, stop and turn to atretic or ovulate, still plague researchers.
Follicle development and maturation is regulated by various hormones secreted by the hypothalamus-pituitary-gonadal axis and several growth factors and cytokines expressed by the ovary and follicles. Over the last few decades, it has been reported that a great number of protein-coding genes are involved in the regulation of folliculogenesis, granulosa cell (GC) proliferation, and oocyte maturation (Miller and Auchus, 2011; Chang et al., 2016). For example, activation of epidermal growth factor receptor is required for phosphorylation of SMAD3 by the EGFR-SFKs-ERK1/2 pathway and the subsequent mediation of GC proliferation (Sasseville et al., 2010). Additionally, the estradiol and estrogen receptor system maintains meiotic arrest in oocytes by binding directly to the promoter regions of the NPPC and NPR2 genes (Liu et al., 2017). Resumption of meiosis leads to a decrease in HDAC3 in GCs, and acetylation of histone H3K14 and binding to the SP1 transcription factor to initiate AREG transcription and oocyte maturation (Wang et al., 2019). Apart from the protein-encoding genes, many lncRNAs have been found to play roles in follicular cell proliferation and oocyte maturation.
Recent studies have identified several lncRNAs in the ovary that play roles in multiple reproduction events across various species via regulation of protein-coding genes. LncRNA is defined as a non-protein coding transcript with a nucleotide length greater than 200 bp. In general, lncRNAs are characterized by traits such as transcription by RNA polymerase II, alternative splicing, polyadenylation, and capping (Ginger et al., 2006; Mehler and Mattick, 2007; Lagarde et al., 2017). For example, lncRNA-AMH2 was found to activate the expression of AMH2 in ovarian GCs by enhancing promoter activity (Kimura et al., 2017). Further, in rat polycystic ovarian syndrome, lncRNA-HOTAIR was found to upregulate the expression of IGF1 and eventually aggravate endocrine disorders and GC apoptosis through competitive binding to miR-130a (Jiang et al., 2020). Additionally, overexpression of another lncRNA, lncRNA-LET, could inhibit migration and promote apoptosis in GCs by upregulating the expression of TIMP2 and activating the Wnt/β-catenin and Notch signaling pathways (Han et al., 2018). Thus, lncRNAs, along with mRNAs, may be closely involved in the regulation of follicle growth and maturation.
The studies of lncRNAs in human and mice are in progress, which contributed to our understanding of non-coding regulatory mechanisms. However, because of the tissue-specificity of lncRNA expression and low conservation of sequences between species (Derrien et al., 2012; Hezroni et al., 2015; Muret et al., 2019), the roles of lncRNAs in ruminant reproduction, especially in goats, is still poorly understood. The present study aimed to shed light on this question by investigating lncRNA and mRNA expression in follicles of Chuanzhong black goat at different stages and attempting to elucidate the regulatory role of lncRNAs in goat folliculogenesis, providing a theoretical basis for improving ovulation in goat breeding.
Materials and Methods
Availability of Supporting Data
This study utilized previously published datasets of 17 samples of uniparous/multiple follicles from 9 Chuanzhong black goats that were classified as 8 large follicles (LF, d > 10 mm) and 9 small follicular pools (eight to ten small follicles per pool, SF, d < 3 mm) (Zou et al., 2020). The datasets generated and analyzed in the current study are available from the corresponding authors upon request. Sample collection, RNA extraction, preparation of libraries for sequencing, and data processing have been described in a previous study (Zou et al., 2020). After removing the low-quality sequences of adaptor sequences and sequences with quality score <Q20 using Cutadapt, the clean reads obtained were aligned using Bowtie2 (Langmead and Salzberg, 2012) and mapped to the goat reference genome (GCF_001704415.1_ARS1) Ensembl V96 using Tophat2 (Kim et al., 2013). The Q30 (the percentage of the number of bases with phred score >30 in the original data to the total number of bases) of each sample is more than 91%, and is relatively constant within each sample. Moreover, the proportion of clean reads in each library is more than 99%. The proportion of samples located in the reference genome of goat ranged from 83 to 91% (Zou et al., 2020).
lncRNA Identification
The transcripts were reconstructed using the Stringtie software, and the cuffcompare software was used to identify transcripts in the database that matched the reconstructed transcripts. Candidate lncRNAs were screened from the reconstructed transcripts based on the structural characteristics of the lncRNAs. The inclusion criteria were (1) sequence length of >200 bp, (2) presence of two or more exons, and (3) minimum coverage greater than 3. The CPC, CNCI, and Pfamscan software programs were used to evaluate the protein-coding ability of the candidate lncRNAs (Bateman et al., 2004; Kong et al., 2007; Sun et al., 2013), and transcripts without protein-coding potential in the three programs were considered as novel lncRNAs.
Expression Analysis
LncRNA and mRNA expression in each sample was evaluated based on the fragments per kilobase per million mapped reads (FPKM). Principal component analysis (PCA), implemented in prcomp function in R, was used to extract the main features of the data as indicators of the overall state of the data and the result was presented by ggplot2 package in R. LncRNAs and mRNAs for which |log2Fold-Change| was >1 and FDR was <0.05 were considered to be differentially expressed between small follicles and large follicles, and the prefix “DE” was added so that they are referred to as DElncRNAs and DEmRNAs henceforth. Hierarchical clustering analysis was performed on DElncRNAs and DEmRNAs using the pheatmap package in R1.
Functional Enrichment Analysis
In this study, mRNAs within the nearest upstream or downstream window of the obtained DElncRNAs were included in the cis-target mRNA datasets of the DElncRNAs. Gene Ontology (GO)2 and Kyoto Encyclopedia of Genes and Genomes (KEGG)3 analysis of DEmRNAs and DElncRNAs were performed using the clusterProfiler package in R. Both GO terms and KEGG pathways with FDR < 0.05 were considered to be significantly enriched.
Construction of DElncRNA–DEmRNA Co-expression Network
Pearson correlation calculations were performed for DEmRNA and DElncRNA (raw P value < 0.05) FPKM of 17 samples (8 large follicles and 9 small follicular pools). The correlation pairs for which the absolute value of Pearson correlation coefficient was lower than 0.60 were excluded from the co-localization pairs. The co-expression relationship between DElncRNAs and DEmRNAs was visualized using Cytoscape (V3.5.1) software (Smoot et al., 2011).
qRT-PCR Verification
The accuracy of RNA-Seq was validated by qRT-PCR by randomly selecting five DElncRNAs and seven DEmRNAs. cDNA for qRT-PCR was synthesized using the PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China), and qRT-PCR was performed using the SYBR® Green PCR Supermix (Bio-Rad, CA, United States). Specific primers for the lncRNAs and mRNAs were designed using Primer Premier 5.0 (shown in Table 1). Each 20-μL reaction mixture included 10 μL of SYBR® Green PCR Supermix, 1 μL of each primer, 1 μL of cDNA, and 7 μL of RNase-free H2O. The PCR reaction conditions were as follows: an initial single cycle at 95°C for 1 min, 34 cycles at three different settings (95°C for 30 s, 58°C for 30 s, and 72°C for 1 min), and a final step at 72°C for 1 min. The dissociation curve of amplified products was used to evaluate product specificity. Relative expression levels were normalized to GAPDH with the 2–ΔΔCt method.
TABLE 1
Table 1. Primer sequences used for qRT-PCR.
Results
Genomic Features of lncRNAs and mRNAs
In this study, a total of 20,985 mRNAs was detected by RNA-seq. The software programs CPC, CNCI, and Pfam were used to identify novel lncRNAs after removal of the transcripts of protein-coding RNA. Eventually, a total of 4,384 lncRNAs including 197 novel lncRNAs (Figure 1A and Supplementary File 1) were retained. The annotated and novel lncRNAs had fewer exons than mRNAs (Figure 1B), and the expression levels of the annotated and novel lncRNAs were lower than that of mRNAs (Figure 1C). Moreover, the average length of the lncRNAs was lesser than that of the mRNAs (Figure 1D).
FIGURE 1
Figure 1. Genomic characteristics of mRNAs and lncRNAs in goat follicles. (A) Screening of candidate lncRNAs by CPC, CNCI, and PFAM. (B) The density curve of exon numbers. (C) Boxplot of the expression level of mRNAs and lncRNAs [expressed as log10 (FPKM+1)]. (D) Boxplot of sequence length [expressed as log10 (length)].
Differentially Expressed lncRNAs and mRNAs
PCA of the mRNA cluster indicated reasonable grouping (Figure 2A), but this was not observed for the lncRNA cluster in small follicles. This is probably because the small follicles in the pools contains different small follicles, which results in the irregular PCA of lncRNAs in small follicular group (Figure 2B). A total of 1,799 mRNAs (1,375upregulated and 424 downregulated mRNAs in small follicles) were differentially expressed between the small and large follicles, based on the criteria applied (|log2FoldChange| > 1 and FDR < 0.05 criteria) (Figure 2C and Supplementary File 2). In addition, 8 lncRNAs were found to be up-regulated in small follicles (Figure 2D and Supplementary File 3). The 30 DEmRNAs and 8 DElncRNAs with the most significant differential expression are shown in Table 2.
FIGURE 2
Figure 2. Expression analysis of mRNAs and lncRNAs. (A) PCA analysis of mRNA expression with FPKM. (B) PCA analysis of lncRNA expression with FPKM. (C) Clustering analysis showing the expression profiles of DEmRNAs between SF and LF groups. (D) Clustering analysis showing the expression profiles of DElncRNAs between SF and LF groups.
TABLE 2
Table 2. The 30 DEmRNAs and 8 DElncRNAs with the most significant difference between large follicles and small follicles.
Functional Analysis of DEmRNAs and DElncRNAs
A total of 1,799 DEmRNAs were significantly enriched in 77 GO terms (FDR < 0.05, Supplementary File 4): biological processes, cellular components, and molecular functions account for 22, 33, and 22 of the terms respectively. The mRNAs targeted by DElncRNAs were significantly enriched in three GO terms in cellular components, including basement membrane, proteinaceous extracellular matrix and extracellular matrix component (FDR < 0.05, Supplementary File 4). KEGG pathway analysis of the DEmRNAs revealed 27 enriched pathways among which involved in steroid hormone biosynthesis, steroid biosynthesis, ovarian steroidogenesis, cAMP signaling, and cell adhesion (FDR < 0.05, Supplementary File 5). The pathway, glycosphingolipid biosynthesis ganglio series were the only pathway in KEGG pathway analysis of DElncRNAs (FDR < 0.05, Supplementary File 5). The top 10 GO terms and 27 KEGG pathways for the DEmRNAs are shown in Table 3.
TABLE 3
Table 3. The top 10 GO terms and 27 KEGG pathways for the DEmRNAs.
DElncRNA–DEmRNA Co-expression Network
The correlation between co-expressed DEmRNAs and lncRNAs was analyzed using Pearson correlation analysis. Because there may be differences among individuals and the FDR were too large after multiple testing control, we used DElncRNA and mRNA with raw P value < 0.05 to perform the co-expression. A total of 32 DEmRNA–lncRNA pairs were found to play a role in follicular development (Figure 3 and Supplementary File 6). Among them, five target pairs, ENSCHIT00000001255-OTX2, ENSCHIT00000006005-PEG3, ENSCHIT00000009455-PIWIL3, ENSCHIT00000007977-POMP and ENSCHIT00000000834-ACTR3 provide a clue in follicular development and maturation of lncRNA-mRNA interaction.
FIGURE 3
Figure 3. Co-expression network of DElncRNAs and DEmRNAs. The red label represents the up-regulation gene in small follicle, and the green label represents the down-regulation gene in small follicle. The circles represent mRNAs and the triangles represent lncRNAs. The sizes of the labels were arranged according to the absolute value of Pearson correlation coefficient; the larger the value, the larger is the size.
Validation of RNA-Seq by qRT-PCR
To validate the reliability of the RNA-seq results, five DElncRNAs (ENSCHIT00000001721, ENSCHIT00000005825, ENSCHIT00000009619, ENSCHIT00000005875, and ENSCHIT00000005031) and seven DEmRNAs (IGFBP5, PTCH1, TSHR, APCDD1, INHA, CYP19A1, and MRO) were randomly selected for qRT-PCR verification. The expression pattern of these DElncRNAs and DEmRNAs was found to be consistent with the patterns determined by RNAs-Seq analysis; therefore, the reliability of the sequencing results was confirmed (Figure 4).
FIGURE 4
Figure 4. Validation of the DEmRNAs and DElncRNAs by qPCR. (A) The qRT-PCR results of the DEmRNAs were compared with the RNA-Seq. (B) The qRT-PCR results of the DElncRNAs were compared with the RNA-Seq. Black: qRT-PCR; Gray: RNA-Seq. In each group, 4 samples were analyzed.
Discussion
In this study, goat follicle samples deposited in a database were analyzed to identify the mRNAs and lncRNAs that are potentially involved in the regulatory mechanisms underlying follicular development and maturation.
A total of 8 DElncRNAs and 1,799 DEmRNAs were identified by RNA-seq. KITLG, INHA, FGF10, and PAPPA were among the 30 most significant DEmRNAs that were found to participate in follicular maturation. These mRNAs have been previously found to play roles in follicle maturation in other species (Driancourt et al., 2000; Rivera and Fortune, 2003; Buratini et al., 2007; Zhang et al., 2010; Namwanje and Brown, 2016). Specifically, 43 DEmRNAs were identified based on the set criteria (FPKM > 10 and |log2FoldChange| > 2, Supplementary File 7). In total, we found that 6 of the 43 highly expressed genes (LHCGR, TMEM200A, STAR, 3BHSD, CDO1, and NTRK1) were also differentially expressed in the follicles between uniparous and multiple in the earlier paper. They were up-regulated in large follicles compared with small one, and were up-regulated in multiple small follicles compared with uniparous small follicles, suggesting that these genes are not only involved in follicular maturation regulation but also improving the litter size through acting on in small follicles of multiple goats (Zou et al., 2020). Among them, LHCGR, STAR, and 3BHSD are known regulatory genes responsible for LH signal acceptance, conversion of cholesterol into pregnenolone, and conversion of progesterone and testosterone in the process of ovarian steroidogenesis which is one of the most important hormonal pathways known for follicular maturation (Miller and Auchus, 2011). In addition, CDO1, NTRK1, MRO, and TC2N (the other two genes in the list of 43 genes) were all upregulated in large follicles and might therefore potentially participate in the regulation of follicular maturation in goats, which have not been identified in the context of follicular development and maturation before (as far as we know). The protein encoded by CDO1 affects mitochondrial function by converting cystine to cysteine sulfamic acid and inhibits the production of glutathione to promote apoptosis, and it has been used as a diagnostic marker in a variety of cancers (Kang et al., 2019; Nishizawa et al., 2019). The NTRK1 protein was found to facilitate follicle assembly and early follicular development in the mouse ovary (Kerr et al., 2009). MRO was first found to be expressed in the testis and was associated with testis differentiation, and it is considered as a candidate testis-determining gene (Smith et al., 2003). A recent study had found that MRO was also expressed abnormally in ovarian cumulus cells of polycystic ovaries, but its function in follicular maturation is unknown (Kenigsberg et al., 2017). TC2N is a tandem C2 domain-containing protein that was found to attenuate p53 signaling via Cdk5 degradation or disruption of the interaction between Cdk5 and p53, thus inhibiting cell proliferation and, subsequently, tumorigenesis (Hao et al., 2019). In summary, LHCGR, TMEM200A, STAR, 3BHSD, CDO1, and NTRK1 are not only involved in follicular maturation regulation but also improving the litter size through acting on in small follicles of multiple goats. The underlying mechanism of CDO1, NTRK1, MRO, and TC2N in follicular maturation needs to be further investigated.
In the present study, enrichment analysis of the identified DEmRNAs was used to understand the important metabolic and molecular differences between small follicles and large ones in goats. KEGG pathway analysis indicated that ovarian steroidogenesis, steroid biosynthesis, cAMP signaling, and cytokine–cytokine receptor interaction were closely related to follicular development. In particular, ovarian steroidogenesis and steroids biosynthesis were upregulated in large follicles. Ovarian steroidogenesis is one of the most important processes in follicular maturation; therefore, the related mRNAs could be used as important markers for goat follicular maturation. Many mRNAs, such as BMP6, CYP11A1, CYP19A1, 3BHSD, STAR, LHCGR, and CYP51A1, have been shown to promote follicular development by regulating steroid production in the ovary (Segaloff et al., 1990; Miller and Auchus, 2011; Rajesh et al., 2018). However, 1,375 of the 1,799 DEmRNAs identified in this study were found to be more actively expressed in small follicles. For example, a total of 29 mRNAs were involved in cAMP signaling. cAMP signaling has previously been implicated that it was expressed in mature follicles to maintain oocyte meiosis inhibition in ovarian (Gulappa et al., 2015; Lyga et al., 2016). In particular, the downstream genes of PTCH1, GLI1, HHIP, and AMH were found to be upregulated in small follicles, suggesting that these genes are implicated in Hedgehog and AMH transduction in small follicles. These findings might suggest that unselected small follicles have more complicated signal transduction or functions in follicular cells before they undergo apoptosis.
8 different expression lncRNAs, ENSCHIT00000002664, ENSCHIT00000009148, ENSCHIT00000007167, ENSCHIT 00000002197, ENSCHIT00000008243, ENSCHIT00000002077, ENSCHIT00000002148, and ENSCHIT00000010661 were all upregulated in the small follicles. Among them, three DEmRNAs, NTN1, AMIGO2, and ENSCHIG00000016389, were located near ENSCHIT00000007167, ENSCHIT00000002077, and ENSCHIT00000010661, respectively, which involved in the follicular process potentially. Study had shown that NTN1 could inhibit granulosa cell viability and estradiol 17β levels while it stimulates progesterone production (Basini et al., 2011). The potential core regulating factors, 32 target pairs of DElncRNA–DEmRNA with raw P value, were predicted using bioinformatics analysis in our study. 2 of target mRNA with negative correlation, ENSCHIT00000007977-POMP and ENSCHIT00000000834- ACTR3, were upregulated in large follicles. The rest of 30 target pairs with positive correlation were upregulated in small follicles. Among them, three target pairs, namely, ENSCHIT00000001255-OTX2, ENSCHIT00000006005-PEG3, and ENSCHIT00000009455-PIWIL3, might be involved in follicular development inhibition by transcriptional regulation of genes. In particular, OTX2 targeted by ENSCHIT00000001255 could affect animal reproductive ability by promoting GnRH promoter activity in GnRH neuron cells, as OTX2 mutation was associated with hypogonadotropism (Diaczok et al., 2011). OTX2 has also been found to be differentially expressed in follicles and oocytes, but its regulatory mechanisms remain unclear (Assou et al., 2006; Wu et al., 2016). Thus, the function of the lncRNA ENSCHIT00000001255 might be associated with the function of its target gene OTX2. Similarly, ENSCHIT00000006005 might play a role in follicular development inhibition through its target gene PEG3. Loss of PEG3-imprinted methylation was observed in individual blastocyst stage embryos after in vitro oocyte culture and hyperovulation in mice (Song et al., 2009; Market-Velker et al., 2010). Additionally, loss of PEG3 in a way of promoter methylation might stimulate clonogenic growth and contribute to the pathogenesis of a majority of ovarian cancers, suggesting the potential effect of PEG3 in female reproduction (Feng et al., 2008). The PIWIL3 protein belongs to the P-element induced wimpy testis-like family, and it plays a role in several stages of genome-wide DNA methylation and meiosis of germ cells (Juliano et al., 2011; Iwasaki et al., 2015). PIWIL3 has been found to affects oocyte development and embryogenesis (Roovers et al., 2015); therefore, it might play a similar role in goat too. Based on all these findings, it is possible that the lncRNAs ENSCHIT00000001255, ENSCHIT00000006005, and ENSCHIT00000009455 are involved in the regulation of promoter activity or methylation, potentially exerting co-regulatory effects on gene transcription in small follicles.
In the present study, the findings indicate that small follicles in goat might involve a complex mechanism of gene expression in follicular inhibition. Several mRNAs and lncRNAs that potentially play important roles in follicular maturation were identified, based on the expression of the classical protein-coding genes in large follicles. However, we were unable to identify a fairly large number of lncRNAs in the goat gene pool. Accordingly, our future work will focus on verifying the functions of the already identified lncRNAs, so as to provide a richer understanding of the molecular mechanism of goat follicular development and maturation.
In conclusion, the present study provides systematic information about the expression of mRNAs and lncRNAs at different stages of antral follicle in goats. Specifically, DEmRNAs and DElncRNAs involved in important biological processes and pathways associated with follicular development inhibition and maturation were identified, and this provides a valuable resource for lncRNA and mRNA studies of goat follicles.
Data Availability Statement
The datasets analyzed for this study are publicly available. This data can be found in NCBI SRA (accession code: PRJNA579007).
Ethics Statement
The animal study was reviewed and approved by Institutional Animal Care and Use Committee of South China Agricultural University (Approval No. 2018-P002).
Author Contributions
GL and DL designed and managed the project. ZZ, XZ, TL, and MD performed the experiment. YL, YG, and BS carried out the data processing and analysis. ZZ and XZ were major contributors in writing the manuscript. All authors have read and approved the manuscript.
Funding
This work was supported by the Guangdong Provincial Promotion Project on Preservation and Utilization of Local Breed of Livestock and Poultry, Modern Agricultural Industrial Technology System of Guangdong Province (2019KJ127), and the Science and Technology Program of Guangdong (2019A050505007).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.589076/full#supplementary-material
Supplementary File 1 | Novel lncRNAs identified in the goat follicles.
Supplementary File 2 | RNA-seq data of DEmRNA expression in large and small follicles.
Supplementary File 3 | RNA-seq data of DElncRNA expression in large and small follicles.
Supplementary File 4 | Gene ontology analysis of DEmRNAs and DElncRNAs.
Supplementary File 5 | Kyoto encyclopedia of genes and genomes analysis of DEmRNAs and DElncRNAs.
Supplementary File 6 | The DElncRNA–DEmRNA co-expression network basing on the co-localization (raw P value < 0.05).
Supplementary File 7 | The subset of 43 DEmRNAs whose FPKM is greater than 10 and absolute value of log2FoldChange is larger than 2.
Abbreviations
lncRNAs, long non-coding RNAs; RNA-seq, RNA sequencing; DElncRNAs, differentially expressed lncRNAs; DEmRNAs, differentially expressed mRNAs; GC, granulosa cell; FPKM, fragments per kilobase per million mapped reads; PCA, principal component analysis; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.
Footnotes
- ^ https://cran.r-project.org/web/packages/pheatmap/index.html
- ^ http://geneontology.org/
- ^ http://www.kegg.jp/
References
Assou, S., Anahory, T., Pantesco, V., Le Carrour, T., Pellestor, F., Klein, B., et al. (2006). The human cumulus–oocyte complex gene-expression profile. Hum. Reprod. 21, 1705–1719. doi: 10.1093/humrep/del065
Basini, G., Cortimiglia, C., Baioni, L., Bussolati, S., Grolli, S., Ramoni, R., et al. (2011). The axonal guidance factor netrin-1 as a potential modulator of swine follicular function. Mol. Cell. Endocrinol. 331, 41–48. doi: 10.1016/j.mce.2010.08.001
Bateman, A., Coin, L., Durbin, R., Finn, R. D., Hollich, V., Griffiths-Jones, S., et al. (2004). The Pfam protein families database. Nucleic Acids Res. 32, D138–D141. doi: 10.1093/nar/gkh121
Buratini, J. J., Pinto, M. G., Castilho, A. C., Amorim, R. L., Giometti, I. C., Portela, V. M., et al. (2007). Expression and function of fibroblast growth factor 10 and its receptor, fibroblast growth 5. Biol. Reprod. 77, 743–750. doi: 10.1095/biolreprod.107.062273
Chang, H. M., Qiao, J., and Leung, P. C. (2016). Oocyte-somatic cell interactions in the human ovary-novel role of bone morphogenetic proteins and growth differentiation factors. Hum. Reprod. Update 23, 1–18. doi: 10.1093/humupd/dmw039
Derrien, T., Johnson, R., Bussotti, G., Tanzer, A., Djebali, S., Tilgner, H., et al. (2012). The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 22, 1775–1789. doi: 10.1101/gr.132159.111
Diaczok, D., DiVall, S., Matsuo, I., Wondisford, F. E., Wolfe, A. M., and Radovick, S. (2011). Deletion of Otx2 in GnRH neurons results in a mouse model of hypogonadotropic hypogonadism. Mol. Endocrinol. 25, 833–846. doi: 10.1210/me.2010-0271
Driancourt, M. A., Reynaud, K., Cortvrindt, R., and Smitz, J. (2000). Roles of KIT and KIT LIGAND in ovarian function. Rev. Reprod. 5, 143–152. doi: 10.1530/ror.0.0050143
Evans, A. C. (2003). Characteristics of ovarian follicle development in domestic animals. Reprod. Domest. Anim. 38, 240–246. doi: 10.1046/j.1439-0531.2003.00439.x
Feng, W., Marquez, R. T., Lu, Z., Liu, J., Lu, K. H., Issa, J. P., et al. (2008). Imprinted tumor suppressor genes ARHI and PEG3 are the most frequently down-regulated in human ovarian cancers by loss of heterozygosity and promoter methylation. Cancer Am. Cancer Soc. 112, 1489–1502. doi: 10.1002/cncr.23323
Fortune, J. E. (1994). Ovarian follicular growth and development in mammals. Biol. Reprod. 50, 225–232. doi: 10.1095/biolreprod50.2.225
Ginger, M. R., Shore, A. N., Contreras, A., Rijnkels, M., Miller, J., Gonzalez-Rimbau, M. F., et al. (2006). A noncoding RNA is a potential marker of cell fate during mammary gland development. Proc. Natl. Acad. Sci. U.S.A. 103, 5781–5786. doi: 10.1073/pnas.0600745103
Gulappa, T., Menon, B., and Menon, K. M. (2015). Hypusination of eukaryotic initiation factor 5A via cAMP-PKA-ERK1/2 pathway is required for ligand-induced downregulation of LH receptor mRNA expression in the ovary. Mol. Cell. Endocrinol. 413, 90–95. doi: 10.1016/j.mce.2015.06.014
Han, Q., Zhang, W., Meng, J., Ma, L., and Li, A. (2018). LncRNA-LET inhibits cell viability, migration and EMT while induces apoptosis by up-regulation of TIMP2 in human granulosa-like tumor cell line KGN. Biomed. Pharmacother. 100, 250–256. doi: 10.1016/j.biopha.2018.01.162
Hao, X. L., Han, F., Zhang, N., Chen, H. Q., Jiang, X., Yin, L., et al. (2019). TC2N, a novel oncogene, accelerates tumor progression by suppressing p53 signaling pathway in lung cancer. Cell Death Differ. 26, 1235–1250. doi: 10.1038/s41418-018-0202-8
Hezroni, H., Koppstein, D., Schwartz, M. G., Avrutin, A., Bartel, D. P., and Ulitsky, I. (2015). Principles of long noncoding RNA evolution derived from direct comparison of transcriptomes in 17 species. Cell Rep. 11, 1110–1122. doi: 10.1016/j.celrep.2015.04.023
Iwasaki, Y. W., Siomi, M. C., and Siomi, H. (2015). PIWI-Interacting RNA: Its biogenesis and functions. Annu. Rev. Biochem. 84, 405–433. doi: 10.1146/annurev-biochem-060614-034258
Jiang, B., Xue, M., Xu, D., Song, J., and Zhu, S. (2020). Down-regulated lncRNA HOTAIR alleviates polycystic ovaries syndrome in rats by reducing expression of insulin-like growth factor 1 via microRNA-130a. J. Cell. Mol. Med. 24, 451–464. doi: 10.1111/jcmm.14753
Juliano, C., Wang, J., and Lin, H. (2011). Uniting germline and stem cells: The function of Piwi proteins and the piRNA pathway in diverse organisms. Annu. Rev. Genet. 45, 447–469. doi: 10.1146/annurev-genet-110410-132541
Kang, Y. P., Torrente, L., Falzone, A., Elkins, C. M., Liu, M., Asara, J. M., et al. (2019). Cysteine dioxygenase 1 is a metabolic liability for non-small cell lung cancer. eLife 8:e45572. doi: 10.7554/eLife.45572
Kenigsberg, S., Lima, P. D., Maghen, L., Wyse, B. A., Lackan, C., Cheung, A. N., et al. (2017). The elusive MAESTRO gene: Its human reproductive tissue-specific expression pattern. PLoS One 12:e174873. doi: 10.1371/journal.pone.0174873
Kerr, B., Garcia-Rudaz, C., Dorfman, M., Paredes, A., and Ojeda, S. R. (2009). NTRK1 and NTRK2 receptors facilitate follicle assembly and early follicular development in the mouse ovary. Reproduction 138, 131–140. doi: 10.1530/REP-08-0474
Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14:R36. doi: 10.1186/gb-2013-14-4-r36
Kimura, A. P., Yoneda, R., Kurihara, M., Mayama, S., and Matsubara, S. (2017). A long noncoding RNA, lncRNA-Amhr2, plays a role in amhr2 gene activation in mouse ovarian granulosa cells. Endocrinology 158, 4105–4121. doi: 10.1210/en.2017-00619
Kong, L., Zhang, Y., Ye, Z. Q., Liu, X. Q., Zhao, S. Q., Wei, L., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. doi: 10.1093/nar/gkm391
Lagarde, J., Uszczynska-Ratajczak, B., Carbonell, S., Perez-Lluch, S., Abad, A., Davis, C., et al. (2017). High-throughput annotation of full-length long noncoding RNAs with capture long-read sequencing. Nat. Genet. 49, 1731–1740. doi: 10.1038/ng.3988
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Liu, W., Xin, Q., Wang, X., Wang, S., Wang, H., Zhang, W., et al. (2017). Estrogen receptors in granulosa cells govern meiotic resumption of pre-ovulatory oocytes in mammals. Cell Death Dis. 8:e2662. doi: 10.1038/cddis.2017.82
Lyga, S., Volpe, S., Werthmann, R. C., Gotz, K., Sungkaworn, T., Lohse, M. J., et al. (2016). Persistent cAMP signaling by internalized LH receptors in ovarian follicles. Endocrinology 157, 1613–1621. doi: 10.1210/en.2015-1945
Market-Velker, B. A., Zhang, L., Magri, L. S., Bonvissuto, A. C., and Mann, M. R. (2010). Dual effects of superovulation: loss of maternal and paternal imprinted methylation in a dose-dependent manner. Hum. Mol. Genet. 19, 36–51. doi: 10.1093/hmg/ddp465
Mehler, M. F., and Mattick, J. S. (2007). Noncoding RNAs and RNA editing in brain development, functional diversification, and neurological disease. Physiol. Rev. 87, 799–823. doi: 10.1152/physrev.00036.2006
Miller, W. L., and Auchus, R. J. (2011). The molecular biology, biochemistry, and physiology of human steroidogenesis and its disorders. Endocr. Rev. 32, 81–151. doi: 10.1210/er.2010-0013
Muret, K., Desert, C., Lagoutte, L., Boutin, M., Gondret, F., Zerjal, T., et al. (2019). Long noncoding RNAs in lipid metabolism: Literature review and conservation analysis across species. BMC Genomics 20:882. doi: 10.1186/s12864-019-6093-3
Namwanje, M., and Brown, C. W. (2016). Activins and inhibins: roles in development, physiology, and disease. Cold Spring Harb. Perspect. Biol. 8:a021881. doi: 10.1101/cshperspect.a021881
Nishizawa, N., Harada, H., Kumamoto, Y., Kaizu, T., Katoh, H., Tajima, H., et al. (2019). Diagnostic potential of hypermethylation of the cysteine dioxygenase 1 gene (CDO1) promoter DNA in pancreatic cancer. Cancer Sci. 110, 2846–2855. doi: 10.1111/cas.14134
Rajesh, G., Mishra, S. R., Paul, A., Punetha, M., Vidyalakshmi, G. M., Narayanan, K., et al. (2018). Transcriptional and translational abundance of Bone morphogenetic protein (BMP) 2, 4, 6, 7 and their receptors BMPR1A, 1B and BMPR2 in buffalo ovarian follicle and the role of BMP4 and BMP7 on estrogen production and survival of cultured granulosa cells. Res. Vet. Sci. 118, 371–388. doi: 10.1016/j.rvsc.2018.04.002
Rivera, G. M., and Fortune, J. E. (2003). Selection of the dominant follicle and insulin-like growth factor (IGF)-binding proteins: Evidence that pregnancy-associated plasma protein a contributes to proteolysis of IGF-binding protein 5 in bovine follicular fluid. Endocrinology 144, 437–446. doi: 10.1210/en.2002-220657
Roovers, E. F., Rosenkranz, D., Mahdipour, M., Han, C. T., He, N., Chuva, D. S. L. S., et al. (2015). Piwi proteins and piRNAs in mammalian oocytes and early embryos. Cell Rep. 10, 2069–2082. doi: 10.1016/j.celrep.2015.02.062
Sasseville, M., Ritter, L. J., Nguyen, T. M., Liu, F., Mottershead, D. G., Russell, D. L., et al. (2010). Growth differentiation factor 9 signaling requires ERK1/2 activity in mouse granulosa and cumulus cells. J. Cell Sci. 123, 3166–3176. doi: 10.1242/jcs.063834
Segaloff, D. L., Wang, H. Y., and Richards, J. S. (1990). Hormonal regulation of luteinizing hormone/chorionic gonadotropin receptor mRNA in rat ovarian cells during follicular development and luteinization. Mol. Endocrinol. 4, 1856–1865. doi: 10.1210/mend-4-12-1856
Smith, L., Van Hateren, N., Willan, J., Romero, R., Blanco, G., Siggers, P., et al. (2003). Candidate testis-determining gene, Maestro (Mro), encodes a novel HEAT repeat protein. Dev. Dyn. 227, 600–607. doi: 10.1002/dvdy.10342
Smoot, M. E., Ono, K., Ruscheinski, J., Wang, P. L., and Ideker, T. (2011). Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27, 431–432. doi: 10.1093/bioinformatics/btq675
Song, Z., Min, L., Pan, Q., Shi, Q., and Shen, W. (2009). Maternal imprinting during mouse oocyte growth in vivo and in vitro. Biochem. Biophys. Res. Commun. 387, 800–805. doi: 10.1016/j.bbrc.2009.07.131
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41:e166. doi: 10.1093/nar/gkt646
Wang, H., Cai, H., Wang, X., Zhang, M., Liu, B., Chen, Z., et al. (2019). HDAC3 maintains oocyte meiosis arrest by repressing amphiregulin expression before the LH surge. Nat Commun. 10:5719. doi: 10.1038/s41467-019-13671-8
Wu, S. J., Cheng, Y. S., Liu, H. L., Wang, H. H., and Huang, H. L. (2016). Global transcriptional expression in ovarian follicles from Tsaiya ducks (Anas platyrhynchos) with a high-fertilization rate. Theriogenology 85, 1439–1445. doi: 10.1016/j.theriogenology.2016.01.005
Zhang, K., Hansen, P. J., and Ealy, A. D. (2010). Fibroblast growth factor 10 enhances bovine oocyte maturation and developmental competence in vitro. Reproduction 140, 815–826. doi: 10.1530/REP-10-0190
Keywords: RNA-seq, follicle, goat, lncRNA, mRNA
Citation: Zhao Z, Zou X, Lu T, Deng M, Li Y, Guo Y, Sun B, Liu G and Liu D (2020) Identification of mRNAs and lncRNAs Involved in the Regulation of Follicle Development in Goat. Front. Genet. 11:589076. doi: 10.3389/fgene.2020.589076
Received: 30 July 2020; Accepted: 23 November 2020;
Published: 17 December 2020.
Edited by:
Yadong Zheng, Lanzhou Institute of Veterinary Research (CAAS), ChinaReviewed by:
Elton J. R. Vasconcelos, University of Leeds, United KingdomYuichiro Watanabe, The University of Tokyo, Japan
Copyright © 2020 Zhao, Zou, Lu, Deng, Li, Guo, Sun, Liu and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guangbin Liu, gbliu@scau.edu.cn; Dewu Liu, dwliu@scau.edu.cn
†These authors have contributed equally to this work