Introduction
Litter size is influenced not only by nutrition levels and environment but also by inheritance (Cui et al., 2009). The ovary is the most important organ for the normal reproductive function of goats. It secretes estrogen to maintain sexual characteristics and cyclical reproductive activity; further, oocytes and ovulation have a major impact on the fertility of goats (Barnett et al., 2006; Zhao et al., 2015). Studies have shown that the ovulation rate of goats is linked to high productivity (Pramod et al., 2013). lncRNAs play a chief role in reproduction-related processes in animals, but very limited information is available on the functions of lncRNAs in goats. In particular, in the context of reproduction in goats, few studies have explored the regulation of lncRNAs (Xing et al., 2014). Long non-coding RNAs are non-coding RNA transcripts of >200 nucleotides in length; they have a complex structure and lack the ability to code proteins (Jarroux et al., 2017). They can regulate gene expression and protein function to perform biological functions. Studies have reported that lncRNAs can regulate reproductive processes, such as ovarian development and maturation in female animals (Li et al., 2015; Ling et al., 2017; Liu et al., 2020). Therefore, it is crucial to study their role by exploring the function of key target genes.
High-throughput RNA sequencing and functional analyses have been used to elucidate the reproductive function of lncRNAs that were identified to be differentially expressed between the ovaries of multiparous and uniparous Anhui white goats; TCONS_00136407, TCONS_00146968, and TCONS_00320849, for example, were suggested to participate in oocyte meiosis (Ling et al., 2017). Using the same method to study the function of differentially expressed lncRNAs (DElncRNAs) in Chuanzhong black goats, ENSCHIT00000005909 and ENSCHIT00000005910 were suggested to regulate the viability and proliferation of keratinocyte-derived cells by influencing IL1R2 (interleukin 1 receptor type II) thereby affecting ovarian function (Bouckenheimer et al., 2018). Leizhou black goat is a special local goat breed in southern China, which shows excellent adaptability to the living circumstance with high humidity and high temperature, and using high-throughput sequencing and bioinformatics analysis can help us to explore the novel functional DElncRNAs in the ovaries of goats.
Litter size is one of the most important economic traits in goat production, determining the benefit of farming enterprises. To provide a theoretical basis for goat breeding and improve the production efficiency of goat industry, it is vital to conduct in-depth research on the mechanisms regulating litter size. We herein screened DElncRNAs between the ovaries of high and low fecundity Leizhou black goats and predicted the target genes of DElncRNAs. In addition, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were used to analyze the function of target genes. Our results not only enrich the transcriptomic data of the goat ovary but also provide a theoretical basis for combining molecular breeding and conventional breeding technologies.
Materials and Methods
Ethics Statement
All study protocols were approved by the Ethics Committee for the Care and Use of Laboratory Animals at the South China Agricultural University (permit no.: SYXK-2014-0136). Further, all experiments were performed in accordance with the guidelines of the South China Agricultural University.
Animals and Sample Collection
Seven healthy female Leizhou black goats (age, 3.5–4.5 years) were divided into high and low fecundity groups. The litter size of high-fecundity group (n = 3) and low-fecundity group (n = 4) were more than one and only one, respectively. Meanwhile, all of the samples in this study were from goats with three parity records of litter size. The female goats were injected with 0.1 mg cloprostenol to induce estrus (Perera et al., 1978; Martemucci and D’Alessandro, 2011; Fierro et al., 2013). The goats were kept under observation to determine whether they were in heat (bleating, searching for the male goat, frequent urination, hyperemia, edema, contraction of the vulva, and vaginal mucus discharge). The basis of estrus was the female goat shaking their tail, standing, and accepting to mate with the male goat (Taylor, 1978; Mekuriaw et al., 2016). The ovaries were collected within 24 h of estrus. The selected goats were killed and dissected, and both whole ovaries from each goat were collected immediately. The intact ovaries were collected and washed with 75% alcohol thrice. Then they were soaked into phosphate buffered saline. After the collection of the ovary, the ligaments and attached tissues were trimmed off under surgical anatomy microscope, ovarian follicles were isolated from the ovary, and the isolated ovarian tissue was frozen in liquid nitrogen and stored at –80°C.
Total RNA Isolation, cDNA Library Construction, and Transcriptome Sequencing
After thoroughly grinding the ovarian tissue, total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, United States). NanoDrop ND-2000 was used to measure RNA concentration (Thermo Science, Wilmington, DE, United States). RNA integrity was assessed by denaturing agarose gel electrophoresis. Further, the cDNA library was constructed using 3 μl of total RNA from each sample, and double-terminal sequencing was performed on the HiSeq X-TEN sequencing platform by Shanghai Parsons Biotech Co., Ltd.
Quality Control of Raw Sequences
We used Cutadapt to remove reads with an average quality score below Q20. The Q20 value referred to the error probability of 1% for the identified bases in the process of base recognition. The reference genome index (GCF_001704415.1_ARS1_genomic1) was established by Bowtie 2, and the filtered reads were compared with the reference genome using TopHat 2. If the mismatch between the reads and the reference genome sequence was within 2, we considered the alignment to be successful (Kim et al., 2013).
Assembly and Novel lncRNA Prediction
According to the TopHat 2 results, StringTie was used for transcript assembly, and candidate lncRNAs were then selected based on the splicing results and structural features of lncRNAs. The screening conditions to identify lncRNAs were as follows: (1) transcripts with low expression levels, low credible single exon transcripts, and exon numbers < 2 were filtered out and (2) transcripts < 200 bp in length were excluded (Trapnell et al., 2010; Cabili et al., 2011). Moreover, Coding-Non-Coding-Index v2 (Barnett et al., 2006), Coding Potential Calculator (0.9-r2; Sun et al., 2013), Pfam Scan v1.3 (Punta et al., 2012), and phylogenetic codon substitution frequency (v20121028; Lin et al., 2011) were used for coding potential analyses. Transcripts without coding potential comprised the candidate set of lncRNAs. lncRNA expression at the transcription level was analyzed with StringTie. DESeq was used to analyze the expression of lncRNAs; the screening conditions were | Log2FoldChange| > 1 and P < 0.05 (Love et al., 2014). The ggplot 2 software package was used to construct a volcano map of DElncRNAs, and the pheatmap software package was used to perform clustering according to the expression level of same lncRNAs in different samples and that of different lncRNAs in the same sample. Distance was calculated with the Euclidean method and clustering was performed using hierarchical agglomerative clustering (Wang et al., 2010).
Target Gene Prediction
To explore the functions of lncRNAs, we predicted the target genes of DElncRNAs. Because the reliability of the analysis results is not high when the sample number is small, the function of trans-regulation can not be predicted. We searched the genes 100 kb downstream and upstream of lncRNAs and analyzed their functions.
GO and KEGG Pathway Analyses for Target Genes of DElncRNAs
GO analysis was performed with the predicted target genes using DAVID2. Furthermore, we used the KEGG database to analyze the potential functions of these genes in pathways3 (Dennis et al., 2003; Han et al., 2012). A hypergeometric test was applied to discover the significant enrichment of GO terms and KEGG pathways so as to determine the main biological functions of differentially expressed genes (Wang et al., 2020; Huang et al., 2007). P < 0.05 indicated statistical significance.
Quantitative Real-Time PCR (qRT-PCR) for DElncRNAs
Total RNA (1 μg) was first reverse-transcribed using an RT Reagent Kit with gDNA Eraser (Takara, Dalian, China), according to manufacturer instructions. qRT-PCR was performed on a StepOnePlus Real-Time PCR System (Life Technologies, United StatesA), as per the standard protocol, using TB Green Fast qPCR Mix (Takara, Dalian, China). Primer Premier 5 used in primer design. Capra hircus β-actin served as the endogenous control for mRNA and lncRNA expression analyses.
Discussion
The elucidation of mechanisms regulating litter size can provide a theoretical basis for breeding technologies in goats. Therefore, in this study, we used the ovaries of high and low fecundity Leizhou black goats to identify DElncRNAs by high-throughput RNA sequencing; moreover, we analyzed the target genes of lncRNAs to explore the role of DElncRNAs in ovarian development.
We herein identified enriched terms and signaling pathways; subsequently, we analyzed them as well as pertinent target genes involved in the regulation of reproduction. Fibroblasts are the main cellular component of loose connective tissue (Cai et al., 2012; Yeung et al., 2013). Carcinoma-associated fibroblasts evidently regulate the development of epithelial ovarian cancer by affecting the proliferation, apoptosis, migration, and invasive activity of ovarian cancer cells (Zhang et al., 2011). Fibroblast growth factor (FGFs) is involved in follicular development and follicular atresia (Costa et al., 2009; Miyoshi et al., 2010; Asgari et al., 2015; Coticchio et al., 2015). The expression of FGFs related to follicular development changes with the development of follicles. FGFs need to bind to distinct receptors to physiologically function (Basu et al., 2014). Apoptosis inhibitor 5, which is encoded by API5, is involved in regulating the cell cycle. Further, it promotes DNA synthesis and cell cycle G1/S transition, and regulates cell growth, proliferation, and apoptosis (Garcia-Jove Navarro et al., 2013). API5 plays an important role in the termination of diapause and early embryonic development of Artemia sinica (Zhang et al., 2017). Through GO analysis, we found that the target gene API5 of MSTRG.3782 was involved in the regulation of fibroblast apoptotic process and FGF binding. MSTRG.3782 was significantly upregulated in the high fecundity group. Accordingly, we hypothesized that lncRNA participates in cell growth, thereby affecting follicular development.
The nucleus is the main repository of genetic information in eukaryotic cells, and the site of DNA replication and transcription; it consequently controls genetic and metabolic activities (Lynch and Marinov, 2017). Chromosomes are the most important structures in the nucleus and carry hereditary information (Zetterström, 2008). A study found that OPA-interacting protein 5 (OIP5) was enriched in centrosomes during the G1 phase of the cell cycle and mediated the regulation of cell division (Naetar et al., 2007). In addition, OIP5 reportedly has a fundamental role in maintaining the structure and function of centrosomes/centromeres (Fujita et al., 2007). KRR1 encodes proteins present in early 90 S precursor particles of the small ribosomal subunit, and its locus has been implied to contribute to the development of polycystic ovary syndrome (Gromadka and Rytka, 2000; Zheng et al., 2014; Pau et al., 2017). In this study, we found that the target gene OIP5 of MSTRG.1201 and the target gene KRR1 of ENSCHIT00000001883 were significantly enriched in the cellular component of nuclear lumen, chromatin, organelle lumen, intracellular organelle lumen, among others. In the high fecundity group, MSTRG.1201 was significantly upregulated and ENSCHIT00000001883 was significantly downregulated. Therefore, we believe that MSTRG.1201 and ENSCHIT00000001883 affect follicular development by regulating cell division.
Mature snoRNP particles are composed of a series of small nucleolar RNA and core proteins. snoRNPs regulate the processing and modification of pre-rRNA and play an important role in ribosomal biogenesis (Richard and Kiss, 2006). Nucleolar protein 56 (encoded by Nop56) is involved in the synthesis of snoRNP as a core protein (Lykke-Andersen et al., 2018). Furthermore, an increase in the ribosome biosynthesis rate can promote the expression of the proto-oncogene C-myc and enhance the proliferative ability of cancer cells (Tomczak et al., 2015). C-myc encodes a transcription factor with a direct role in controlling translation (Ruggero, 2009). Nol5a/Nop56 may be a critical gene involved in Myc-mediated oncogenic transformation (Cowling et al., 2014). According to our GO and KEGG pathway analyses, MSTRG.2938 was significantly upregulated in the high fecundity group, and its target gene NOP56 was involved in ribosomal biogenesis in eukaryotes and the pre-snoRNP complex. We thus speculate that MSTRG.2938 regulates ribosomal biogenesis in the pre-snoRNP complex as well as cell transformation in eukaryotes. However, the specific mechanism of regulation of each lncRNA remains to be further investigated.
According to GO analyses, the target genes [Guanine Nucleotide Binding Protein, alpha 13(GNA13), Mothers against decapentaplegic homolog 2 (SMAD2), and Fibronectin Leucine Rich Transmembrane Protein 2 (FLRT2)] of all lncRNAs were mainly enriched in positive regulation of transcription from RNA polymerase II promoters, patterning of blood vessels, palate development, and positive regulation of synapse assembly. RNA polymerase II plays a pivotal role in the transcription of protein-encoding genes in all eukaryotic cells (Bernecky et al., 2016). GNA13 participates in regulating cell movement and developmental angiogenesis (Offermanns et al., 1997). Moreover, SMAD2 overexpression has been reported to repair secondary cleft palate by increasing apoptosis of medial edge epithelial cells in the TGF-β3 pathway (Miyazono et al., 2018). We thus report that these genes play a major role in maintaining the healthy growth of goats.
The Wnt signaling pathway and its downstream effectors not only regulate physiological processes such as cell growth and differentiation, cell migration, and genetic material stability but are also important for cancer progression, including for regulating tumor growth, cell senescence, and cell death. The Wnt/β-catenin signaling pathway is involved in various important processes, such as the regulation of embryo development, cell proliferation, and cell migration (Nusse and Clevers, 2017; Peng et al., 2017; Steinhart and Angers, 2018). β-Catenin is an essential structural component of cadherin-based adherens junctions and is a key component of Wnt/β-catenin signal transduction. Aberrant expression of CTNNB1 and WNT5A has been observed to affect cell proliferation and lead to cancer occurrence. A mutation in CTNNB1 is one of the many causes of β-catenin degradation (Kim et al., 2018). The Wnt/CTNNB1 pathway is a pivotal signaling pathway that regulates steroid production (Abedini et al., 2016). WNT5A is a highly evolved conservative non-classical Wnt ligand, which is required for normal ovarian follicle development (Abedini et al., 2016). WNT5A is differentially expressed during the development of mouse follicles, and it can significantly inhibit steroid production in atretic follicles (Lapointe and Boerboom, 2011; Abedini et al., 2016; Kumawat and Gosens, 2016). By blocking the function of FSH (follicle-stimulating hormone) and luteinizing protein, WNT5A can induce the down-regulation of CTNNB1 and cAMP-response element binding protein (CREB), thus affecting follicle development and gonadotropin reactivity (Abedini et al., 2015). We found that the target gene CTNNB1 of ENSCHIG00000000641 is one of the two signaling pathway members of Wnt and proteoglycans in cancer. ENSCHIG00000000641 was downregulated in the high fecundity group. Further, the target gene WNT5A of ENSCHIG00000000774 participated in the negative regulation of canonical Wnt signaling pathway terms and the Wnt signaling pathway. ENSCHIG00000000774 was also downregulated in the high fecundity group. Thus, we believe that ENSCHIG00000000641 and ENSCHIG00000000774 affect follicular development by regulating cell proliferation and steroid production, respectively.
The TGF-β superfamily participates in many physiological activities in mammals via autocrine and paracrine pathways, and TGF-β is mainly produced locally in the ovary. It has been reported that TGF-β1 can promote the growth of mice follicles. Both TGF-β and activin A have proliferative action and cytodifferentiative action on granulosa cells (Liu et al., 1999). As a member of the transforming growth factor β family, TGF-β2 also plays an important role in the growth and development of follicles. TGF-β2 is located in follicular membrane cells and luteal cells and regulates the production of inhibitors and activins in granulosa cells and luteal cells (Knight and Glister, 2003). Smad2/Smad3 are key molecules in the TGF-β/Smad signaling pathway that regulate ovarian growth and development and maintain ovarian function (Coutts et al., 2008; AlMegbel and Shuler, 2020). Smad2 and Smad3 can maintain normal fertility in women, and further support the Smad2/3 pathway in the ovary to participate in the regulation of signals produced by oocytes, which plays an important role in the coordination of ovulation (Li et al., 2008). In this study, we found that TGF-β2, TGF-βR2, and Smad2 participated in the TGF-β signaling pathway, and they were regulated by ENSCHIG00000000886, ENSCHIG00000000609, and ENSCHIG00000002761, respectively. We accordingly speculate that these lncRNAs regulate follicle development, but the specific mechanism needs to be further studied.
To conclude, we found that target genes of all lncRNAs were mainly involved in protein transcription and played a role in maintaining the healthy growth of animals. In addition, the TGF-β and Wnt signaling pathways were found to be related to reproduction in animals. Based on functional analyses of target genes of DElncRNAs, fibroblast apoptotic process, FGF binding, pre-snoRNP complex, and ribosomal biogenesis in eukaryotes were associated with reproduction in goats. Our data improves the current understanding of the transcriptome of goats and provides valuable information for functional genomics resources and biological studies; moreover, we believe that our results are of great significance for in-depth studies of candidate lncRNAs in breeding techniques.