Open Access

Screening and evaluating of long noncoding RNAs in the puberty of goats

  • Xiaoxiao Gao1,
  • Jing Ye1, 2, 3,
  • Chen Yang1,
  • Kaifa Zhang1,
  • Xiumei Li1, 2, 3,
  • Lei Luo1,
  • Jianping Ding1, 2, 3,
  • Yunsheng Li1, 2, 3,
  • Hongguo Cao1, 2, 3,
  • Yinghui Ling1, 2, 3,
  • Xiaorong Zhang1, 2, 3,
  • Ya Liu1, 2, 3,
  • Fugui Fang1, 2, 3Email author and
  • Yunhai Zhang1, 2, 3Email author
Contributed equally
BMC Genomics201718:164

https://doi.org/10.1186/s12864-017-3578-9

Received: 13 July 2016

Accepted: 9 February 2017

Published: 14 February 2017

Abstract

Background

Long noncoding RNAs (lncRNAs) are involved in regulating animal development, however, their function in the onset of puberty in goats remain largely unexplored. To identify the genes controlling the regulation of puberty in goats, we measured lncRNA and mRNA expression levels from the hypothalamus.

Results

We applied RNA sequencing analysis to examine the hypothalamus of pubertal (case; n = 3) and prepubertal (control; n = 3) goats. Our results showed 2943 predicted lncRNAs, including 2012 differentially expressed lncRNAs, which corresponded to 5412 target genes. We also investigated the role of lncRNAs that act cis and trans to the target genes and found a number of lncRNAs involved in the regulation of puberty and reproduction, as well as several pathways related to these processes. For example, oxytocin signaling pathway, sterol biosynthetic process, and pheromone receptor activity signaling pathway were enriched as Kyoto Encyclopedia of Genes and Genomes (KEGG) or gene ontology (GO) analyses showed.

Conclusion

Our results clearly demonstrate that lncRNAs play an important role in regulating puberty in goats. However, further research is needed to explore the functions of lncRNAs and their predicted targets to provide a detailed expression profile of lncRNAs on goat puberty.

Keywords

LncRNA Puberty Goat Hypothalamus Transcriptome

Background

Puberty is a pivotal stage in female goat development. It marks the first occurrence of ovulation and the onset of reproductive capability [1]. The mechanism of puberty onset is complex and thought to be associated with environmental factors, neuroendocrine factors, genetic factors and their interactions. In general, the secretion of gonadotropin-releasing hormone (GnRH) is considered a crucial factor in puberty onset for goats [2]. A popular view is that during the prepubertal period, secreting neurons suffer persistent trans-synaptic inhibition. This means that GnRH secretions increase as long as this inhibition is eliminated, which leads to puberty [2]. However, these influences are based on substantial genetic control [3].

It was previously reported that the initiation of puberty in female rats is regulated by epigenetic mechanism of transcriptional repression [4], whereby epigenetic control was composed of several mechanisms. Two well established mechanisms include: modification of chromatin and chemical modification of the DNA (including DNA methylation and hydroxymethylation). Non-coding RNA is the most recently unveiled mechanism of epigenetic control, which affords epigenetic information by lncRNAs or microRNAs [5]. Broadly, lncRNAs are known as transcripts greater than 200 nt in length that do not appear to code proteins [6].

During the past decades of transcriptome studies, multiple lncRNAs have been discovered, such as Xist and H19. The advent of RNA-seq has been a powerful tool in exploring and quantifying lncRNAs [7], which has led to the identification of many more lncRNAs that await functional validation. Most identified lncRNAs have primarily originated from human and mouse studies [8, 9]. Recent studies in bovine [1012] and porcine species [13, 14] have enriched the mammalian lncRNAs databases, providing a promising future for further lncRNAs studies.

LncRNAs have been shown to participate in the regulation of transcriptional and post-transcriptional control [15]. In recent years, lncRNAs have proven to play roles in lactation, ovary development, and embryo and sperm maturation. Therefore, we inferred that the onset of goat puberty is also regulated by lncRNAs. In this study, we applied RNA-seq and investigated the expression profiles of mRNA and lncRNAs in pubertal and prepubertal goats to explore the association of lncRNAs with the onset of puberty.

Methods

Preparation of animals and tissues

This study was authorized and endorsed by the Animal Care and Use Committee of Anhui Agricultural University. We housed three, prepubertal (aged 2.5-3.0 months) and three, pubertal (aged 4.5-5.0 months) female Anhuai goats under the same conditions on a farm in Anhui Province, China. We determined puberty goats in studied femal by male goats detecting estrus and the appearance changes of vulva [16]. The average weight of pubertal goats was 16.17 kg compared with the prepubertal goats 8.30 kg, and the average weight of the pubertal goats’ ovary was 0.76 g compared with the prepubertal goats 0.30 g. The animals were deeply anesthetized by intravenous administration of 3% pentobarbital sodium (30 mg/kg; Solarbio, P8410, China) and sacrificed by exsanguination in a healthy physiological stage when pubertal goats were in the late follicular phase. Hypothalamic tissues were surgically removed, and frozen in liquid nitrogen immediately. These tissues were stored at −80 °C until the RNA extraction [17].

RNA sequencing and quality control

We isolated total RNA from goat hypothalamus using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA), according to the standard extraction protocol. The contamination and degradation of RNA was detected by 1% agarose gels. The purity of RNA was monitored using the NanoPhotometer® spectrophotometer (Implen, Los Angeles, CA, USA). We measured the concentration of RNA using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA). The integrity of RNA was monitored as previous reported [18]. We used 3 μg RNA per sample for the RNA sample preparations. Firstly, ribosomal RNA in total RNA was removed [19], and then the residue was cleaned up by using ethanol precipitation. Then the libraries whith high strand-specificity for sequencing was generated [19], following manufacturer’s recommendations. Then the process was followed as previously described [20]. Illumina Hiseq 4000 platform was adopted on sequencing and 150 bp paired-end reads were generated. Raw reads were dealt with in-house perl scripts. The reads with more than 10% unknown bases, reads containing adapter and reads with more than 50% of low-quality bases (whose Phred scores were < 5%) were removed, yielding only the clean reads. Meanwhile, the quality of clean reads (Q20, Q30, and GC content) were detected. All the following analyses were based on high quality clean reads.

Transcriptome assembly

We used a GTF file (ftp://ftp.ncbi.nlm.nih.gov/genomes/Capra_hircus/GFF/) with the annotation of the goat genome. Index of the reference genome was created by Bowtie v2.0.6 [21, 22] and then we aligned paired-end clean reads to the reference genome using TopHat v2.0.9 [23]. The mapped reads of each sample were assembled by both Scripture (beta2) [24] and Cufflinks (v2.1.1) [25, 26] in a reference-based approach. Both methods determined exons connectivity by spliced reads. Scripture ran using default parameters, while Cufflinks ran with min-frags-per-transfrag = 0’ and–library-type fr-firststrand’. Other parameters were set as default.

Expression and coding potential analysis of transcripts

Gene expression was calculated using FPKMs of transcripts in each sample [27]. We confirmed differential expression in gene expression data using Cuffdiff as it based on the negative binomial distribution provides statistical routines [25]. Transcripts with a P < 0.05 were assigned as significantly differentially expressed between two groups. We used three analytic tools, including Coding-Non-Coding-Index (CNCI; v2) [28], Coding Potential Calculator (CPC; 0.9-r2) [29], Pfam Scan (v1.3) [30] to screen out candidate lncRNAs. CNCI (v2) profiles distinguished protein-coding and non-coding sequences effectively by adjoining nucleotide triplets, which was independent of known annotations. CPC (0.9-r2) was mainly used to detect the extent and quality of the Open Reading Frames (ORF) in a transcript and discover the sequences in known protein database, clarifying the coding and non-coding transcripts. Each transcript was translated in all three possible frames, then any of the known protein family was identified by Pfam Scan (v1.3) in the Pfam database (release 27; adopted both Pfam A and Pfam B). The coding potential of transcripts predicted by any of the three tools above were filtered out (non-annotated transcriptional activity by identifying novel transcripts), and those without coding potential were our candidate lncRNAs for further analysis.

Target gene prediction and functional enrichment analysis

The cis role refers to the lncRNA acting on neighboring target genes [31, 32]. To predict the cis-regulated target genes of lncRNAs, we screened protein-coding genes as potential targets 10 K/100 K upstream and downstream of lncRNAs andanalyzed their function. The trans role refers to the coexpression relationship between lncRNAs and mRNA. Expression levels of lncRNAs and mRNAs were calculated for Pearson’s correlation coefficients by custom scripts (r > 0.95 or r < −0.95). The target genes of lncRNAs were performed functional enrichment analysis by clustering the genes from various samples using the DAVID platform [33]. The significance was described as a P-value, measured by the EASE score (P < 0.05 was considered significant).

Quantitative real-time PCR

We have validated the RNA-seq data by selected eight lncRNAs, two novel transcripts, and two target genes to investigate the expression patterns in the samples using qRT-PCR. Werepeated the qRT-PCR experiments three times per sample on expression from the same hypothalamic tissues of three pre vs. three pubertal goats. We designed primers online using Primer5 software and evaluated using BLAST at NCBI. A list of the primer sequences is shown in Table 1. We performed qRT-PCR using SYBR green (Vazyme, China) method. Expression levels of genes were quantified through the cycle threshold (Ct) values and evaluated as 2-ΔΔCT. The data of expression was normalized to β-action.
Table 1

qRT-PCR primer and size of the amplification products of the target and housekeeping genes

Gene

Forward primer, 5’-3’

Reverse primer, 5’-3’

Product size, (bp)

XLOC_957527

ACACGACCAGAACATCAG

AATCACAGGAGAAGAGTAGG

162

XLOC_910648

TGCCATCCAGCCATCTCATC

CACACACCACAGTTCCTTTACC

175

XLOC_1101518

CTCCTGGGCTACCGAATGT

GCGGCTGTGAACTAAATGG

172

XLOC_1276445

TCGCTCCGTCTTCACCTAC

CGTTGCTCCATCACCCTTG

100

XLOC_2056339

CCTGTTGTTGGAATCACTC

CTCTTATGCCTCGGATGG

184

XLOC_960044

CAAGGAGTCGCACGCTAC

CTCTTACGCCTCTGAATCGG

128

XLOC_032671

TGATGCCAAGAGGTAGCC

TTATACAGACAGTGAAAGAGAGG

180

XLOC_688924

ATCTCCACTCTACAAACCTATACC

ACTCTCAAAGGGAAGCCAATG

142

Novel_000476

TTACCTACTACCATCCTTCAC

ATCAGCAGAACCAGAACC

178

Novel_000453

GAAGTGTCGTCTGGAGATTACC

TGTTGAGTGAGTCCTGTATTACC

181

CD38

CTACTGCCTTCTTCTGTG

TTCTGCTTCTGGAATACG

185

GnRH1

CTAATCCTGCTGACTTTCTGTG

ACCTCTTTGGCTATCTCTTGG

127

β-actin

CGTGACATCAAGGAGAAG

GAAGGAAGGCTGGAAGAG

171

GO and pathway analysis

In this study, Gene Ontology (GO) enrichment [34] analysis of targets was performed by the GOseq R package and corrected by P (P < 0.05 were considered significantly enriched). Pathway analysis is a functional analysis in KEGG (http://www.genome.jp/kegg) pathways [35]. We evaluated the statistical enrichment of lncRNAs target genes or differential expression genes in KEGG pathways using KOBAS [36] software.

Statistical analysis

We performed further analysis of RNA-seq data and graphical representations using the statistical R package (R, Auckland, NZL), adopting multiple testing and P corrections. We applied SPSS 17.0 software package (SPSS, Chicago, IL, USA) to analyze the qRT-PCR data. Differential expression levels of genes were calculated by independent-samples t-test between prepubertal and pubertal goats. Significance of data was defined as P < 0.05.

Results

Identification of lncRNAs

We used pubertal and prepubertal goats to perform RNA-seq analysis from the hypothalamus of six female Anhuai goats. In total, 774,998,560 raw reads were produced under the Illumina HiSeq 4000 platform. We obtained 636,544,196 reads maped to goat reference genome after discarding low-quality sequences and adaptor sequences. The percentage of mapped reads among clean reads in each library ranged from 80.45% - 84.32% (Additional file 1). After the analysis of coding potential using the CNCI, CPC and Pfam-scan software, we identified 2943 lncRNAs (Fig. 1), including 2426 large intergenic noncoding RNAs, 217 anti-sense_lncRNAs, and 300 intronic_lncRNAs.
Fig. 1

Screening of candidate lncRNAs in hypothalamus transcriptome. The coding potential of lncRNAs were analyzed by three tools (CPC, CNCI and PFAM)

Genomic features of lncRNAs

Overall, we observed a lower expression of lncRNAs compared with mRNA (Fig. 2a). The mean length of lncRNAs in our dataset was 1180 nt, and the mean mRNA length was 2869 nt (Fig. 2b). Furthermore, we detected an ORF mean length of 105 nt for our lncRNAs, which tended to be shorter than protein coding genes (Fig. 2c). We also found lncRNAs contained fewer exons than mRNA (Fig. 2d).
Fig. 2

The comparison of features between predicted lncRNAs and mRNA. a Expression of lncRNAs and mRNA. b Length distribution of 2943 predicted lncRNAs and 30162 coding transcripts. c ORF length distribution of lncRNAs and coding transcripts. d Exon number distribution of lncRNAs and coding transcripts

Differential expression cluster analysis

Further analysis identified 1165 significant differential expression transcripts (including lncRNAs, mRNAs, and novel transcripts) (P < 0.05), 770 up-regulated and 395 down-regulated transcripts (Fig. 3), including 57 novel transcripts. Furthermore, we detected 59 lncRNAs transcripts from 58 lncRNAs gene loci significant differentially expressed lncRNAs (P < 0.05), including 29 up-regulated and 30 down-regulated lncRNAs transcripts in pubertal samples compared with prepubertal samples (Additional file 2). We validated sequencing results using qRT-PCR analysis (Fig. 4).
Fig. 3

Volcano plots of differential expression transcripts. X-axis is fold change (log 2) and Y-axis is P (−log 10). Red points indicate up-regulated (X axis > 0) transcripts; green points indicate down-regulated (X axis < 0) transcripts

Fig. 4

Validation of RNA-seq results by using quantitative qRT-PCR. Some lncRNAs and target genes were examined using quantitative qRT-PCR. The data are expressed as the mean ± 1 SD (n = 3). *p < 0.05, **p < 0.01

Prediction of target genes of lncRNAs in cis and trans

LncRNAs can act on target genes, either in cis (neighbor the site of lncRNA production) or in trans to coexpression whith target genes [37]. To explore whether differences in lncRNAs affects functional regulation of goat puberty, we predicted the target genes of lncRNA using the cis and trans model. To analyze the cis role of lncRNA, we screened protein-coding genes as potential targets 10 K/100 K upstream and downstream of the lncRNAs. The results indicated that there were 2012 lncRNAs that corresponded to 5412 target genes (Additional file 3). Interestingly, we observed several genes related puberty such as PRLHR, EMC3, IGF2BP1, ZNF175, and ZNF444 [3841], which were respectively a target of XLOC_1486935, XLOC_1284300, XLOC_957527, XLOC_080674, and XLOC_910648, indicating that the onset of puberty is probably regulated by the lncRNA-tatget genes. Regarding the trans role of lncRNA, our results showed that the lncRNA, XLOC_957527, acted on GnRH1 (Table 2; Additional file 4).
Table 2

LncRNAs and its potential target genes associated with puberty

Target genes

Cis-lncRNA

Trans-lncRNA

ZNF175

XLOC_080674

XLOC_1594657, XLOC_720006, XLOC_1831092,

DNAJB2

XLOC_1101518

XLOC_1891931, XLOC_1516990, XLOC_1554449,

XLOC_1021724

EMC3

XLOC_1284300

XLOC_1429620, XLOC_1269384, XLOC_047864,

XLOC_1409381, XLOC_1988116, XLOC_1430952

PRLHR

XLOC_1486935

XLOC_559241

MEF2C

XLOC_2123676

XLOC_674601, XLOC_2253794, XLOC_2334337,

XLOC_1598694, XLOC_1334265, XLOC_2423839

XLOC_1178955, XLOC_1959074, XLOC_263620,

XLOC_1375793

ZNF444

XLOC_910648

XLOC_1529384, XLOC_083680, XLOC_1561175, XLOC_1876674, XLOC_2092234, XLOC_1371300, XLOC_1341798, XLOC_471078, XLOC_904488, XLOC_1248122, XLOC_070767, XLOC_1753539

XLOC_1985452, XLOC_1680696, XLOC_912734, XLOC_2285546, XLOC_1011371

IGF2BP1

XLOC_957527

XLOC_1787362, XLOC_1068007, XLOC_428323,

XLOC_1405012, XLOC_395262, XLOC_1970958, XLOC_228837

CD38

 

XLOC_523219

GnRH1

 

XLOC_957527

GO and KEGG analysis

Our GO analysis of predicted targets demonstrated 73 significantly enriched terms (P < 0.05). The top eight terms were as follows: pheromone receptor activity, hyalurononglucosaminidase activity, hexosaminidase activity, sensory perception of taste, viral genome packaging, helicase activity, sensory perception of chemical stimulus, and sensory perception (Additional file 5). Interestingly, the signaling pathway of the pheromone receptor activity was significantly enriched, which relates to goat estrus. In addition, the sterol biosynthetic process signaling pathway was significantly enriched. DNAJB2 was the differentially expressed target gene on the pathway, which suggests that it may be a new gene involved in the regulation of puberty onset via the sterol biosynthetic process signaling pathway.

KEGG pathway analysis of lncRNAs targets showed 90 terms were enriched (Additional file 6), in which oxytocin signaling pathway was related to puberty [42]. These results suggested that lncRNAs may be cis-acting on its target genes to regulate onset of puberty (Fig. 5).
Fig. 5

KEGG annotation for target gene functions of predicated lncRNAs. Red indicates higher expression and green indicates lower expression. The number of differentially expressed genes is shown in parentheses

We also evaluated the trans role of 2943 lncRNAs in protein-coding genes by its correlation coefficient of gene expression (Pearson correlation ≥ 0.95 or ≤ −0.95). The results showed that 2551 lncRNAs had interactions with target genes in trans of the goat genome. Functional analysis illustrated that target genes in trans were enriched (P < 0.05) in 158 GO terms including a variety of processes (Additional file 7), such as G-protein coupled receptor activity, transmembrane signaling receptor activity, receptor activity, and so on.

We identified 273 KEGG pathways (Additional file 8), several of which were associated with puberty, including ovarian steroidogenesis, GnRH signaling pathway, steroid biosynthesis, steroid hormone biosynthesis, oxytocin signaling pathway, GABAergic synapse, estrogen signaling pathway, oocyte meiosis, glutamatergic synapse, and others. These findings indicate that lncRNAs may act on the target genes associated with puberty of goat in trans.

Specific expression of lncRNAs

There were 187 specific expressions of lncRNAs in the pubertal samples, especially, XLOC_2409732, which has a lower P than other specific expression of lncRNAs. The targets of XLOC_2409732 were detected as ASB5, WDR17, SPATA4 and SPCS3 according to RNA-seq analysis. We found 243 specific expressions of lncRNAs in prepubertal samples, particularly XLOC_1498149, which has a lower P than other specific expression of lncRNAs, and CDR1 was targeted to XLOC_1498149 (Additional file 9). These two specifically expressed lncRNAs may play a pivotal role in goat puberty and, with further studies, provide crucial information regarding the regulation of puberty.

Discussion

We initially performed RNA-seq to analyze lncRNAs of hypothalamus from pubertal and prepubertal female Anhuai goats. Through sequencing, we acquired 2943 predicted lncRNAs and 30162 coding transcripts. Many studies have indicated that lncRNAs have unique features compared with mRNA; for example, lncRNAs are shorter in length than protein-coding transcripts [27]. Furthermore, we found that lncRNAs in hypothalamus were shorter than in skin (1809 bp on average); however, the number of exons is similar [20]. Interestingly, the length of lncRNAs in goat hypothalamus are longer than that in human (1 kb on average) and mouse (550 nt on average), containing fewer exons than human (2.9 exons on average) and mouse (3.7 exons on average) [43].

In this study, we screened out significant differentially expressed transcripts, including 59 lncRNA transcripts from 58 lncRNA gene loci. As previous research, the functions of lncRNAs were reflected by acting on the protein-coding genes. For example, in a recent study, a muscle-specific lncRNA, linc-MD1, influenced muscle development by targeting to MAML1 [44]. Moreover, the lncRNA, Neat1, could make a difference in pregnancy by acting on corpus luteum formation in mice [45]. Therefore, we could predict the role of mammalian lncRNAs by the relevant protein-coding genes.

Here, we predicted the potential functions of lncRNAs through the protein-coding genes in cis and trans. Several genes have been confirmed to be associated with puberty onset, including Kiss1/GPR54 [4649], IGFs [50], GABA [51] and FSHR. We discovered several differentially expressed targets in cis and trans for lncRNAs in pubertal and prepubertal hypothalamus. PRLHR, EMC3, IGF2BP1, ZNF175, ZNF444 have been reported involved in the regulation of puberty and reproduction of animal [40, 52]. For example, puberty of female rats is significantly advanced by GnRH release under the stimulation of IGF-1; IGF-1 can affect the puberty-related events by hypothalamic GnRH release [53].

Moreover, previous research has demonstrated that puberty is delayed after over expression of ZNFs in the arcuate nucleus (ARC) of female rats; subsequent oestrous cyclicity is also disrupted [40]. Our results also showed the lncRNA, XLOC_957527, acted on GnRH1 through trans interactions. Consequently, we confirmed that relevant lncRNAs might play a crucial role on regulation of puberty via the above targets. However, these predicted functions of lncRNAs need further experimental verification.

In the present study, oxytocin signaling pathways were enriched in KEGG pathways. MEF2C, as the target gene of lncRNA XLOC_2123676, is an important gene in oxytocin signaling pathway associated with puberty regulation. The age that vaginal opening occurs in female rats is delayed by treatment with an oxytocin antagonist, indicating that oxytocin enhances sexual maturation [42]. We also observed that DNAJB2, the target of lncRNA XLOC_1101518, has a crucial role in sterol biosynthetic process, which is involved in regulation puberty. ESR1 is essential for multiple estrogen feedback loops and required for puberty onset in female mouse [54].

Our GO analysis of the predicted targets indicates that pheromone receptor activity signaling pathway, which relates to goat estrus, is significantly enriched (Fig. 6) [55].
Fig. 6

GO enrichment analysis for target gene functions of predicated lncRNAs. (MF: molecular function)

In our study, the enriched KEGG pathways and GO pathways associated with reproduction and puberty clearly suggest that these lncRNAs play a vital role in regulation of puberty in goats. However, the functions of lncRNAs and their predicted targets analyses should be carefully evaluated by further experiments.

Conclusion

We performed RNA-seq analysis, and screened out differentially expressed lncRNAs of pubertal and prepubertal goats. We elucidated genomic differences between lncRNAs compared with mRNA. Then, we observed several target genes of lncRNAs related to puberty. Our results clearly demonstrate that lncRNAs play an important role in regulating puberty in goats.

Abbreviations

ASB5: 

ankyrin repeat and SOCS box containing 5

DNAJB2: 

DnaJ (Hsp40) homolog, subfamily B, member 2

EMC3: 

ER membrane protein complex subunit 3

FSHR: 

follicle stimulating hormone receptor

GnRH: 

gonadotropin-releasing hormone

IGF2BP1: 

insulin-like growth factor 2 mRNA binding protein 1

LncRNA: 

long noncoding RNA

MEF2C: 

myocyte enhancer factor 2C

PRLHR: 

prolactin releasing hormone receptor

WDR17: 

WD repeat domain 17

ZNF175: 

zinc finger protein 175

ZNF444: 

zinc finger protein 444

Declarations

Acknowledgments

We thank all members of the Anhui Provincial Laboratory of Animal Genetic Resources Protection and Breeding for their stimulating discussions.

Funding

This work was supported by the National Natural Science Foundation of China (Grant number 31472096), the National Transgenenic New Species Breeding Program of China (No. 2014ZX08008-005-004), the National Natural Science Foundation of China (Grant number 31301934), and the Anhui Provincial Natural Science Foundation (Grant number 1408085MKL40).

Availability of data and materials

The sequencing data were submitted to the Genome Expression Omnibus (Accession Numbers GSE84301) in NCBI. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=wjydswcgxrstfib&acc=GSE84301.

Authors’ contributions

GXX and YJ conceived of the study, participated in its design and coordination and drafted the manuscript; YC, ZKF, LXM and LL conducted qRT-PCR validation and statistical analysis; DJP and LYS performed the statistical analysis; CHG, LYH, ZXR and LY carried out the analysis of data, animal experiments, and surgical processes; FFG and ZYH participated in the design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval

The study was approved by the Animal Care and Use Committee of Anhui Agricultural University. The methods were carried out in accordance with the approved guidelines. All experimental procedures involving goats were performed according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China; revised in June 2004).

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Anhui Provincial Laboratory of Animal Genetic Resources Protection and Breeding, College of Animal Science and Technology, Anhui Agricultural University
(2)
Anhui Provincial Laboratory for Local Livestock and Poultry Genetic Resource Conservation and Bio-Breeding
(3)
Department of Animal Veterinary Science, College of Animal Science and Technology, Anhui Agricultural University

References

  1. Cao GL, Feng T, Chu MX, Di R, Zhang YL, Huang DW, Liu QY, Hu WP, Wang XY. Subtraction suppressive hybridisation analysis of differentially expressed genes associated with puberty in the goat hypothalamus. Reprod Fertil Dev. 2015;28(11):1781–1787.Google Scholar
  2. Smith MJ, Jennes L. Neural signals that regulate GnRH neurones directly during the oestrous cycle. Reproduction. 2001;122:1–10.View ArticlePubMedGoogle Scholar
  3. Abreu AP, Macedo DB, Brito VN, Kaiser UB, Latronico AC. A new pathway in the control of the initiation of puberty: the MKRN3 gene. J Mol Endocrinol. 2015;54(3):R131–9.View ArticlePubMedPubMed CentralGoogle Scholar
  4. Lomniczi A, Loche A, Castellano JM, Ronnekleiv OK, Bosch M, Kaidar G, Knoll JG, Wright H, Pfeifer GP, Ojeda SR. Epigenetic control of female puberty. Nat Neurosci. 2013;16(3):281–9.View ArticlePubMedPubMed CentralGoogle Scholar
  5. Lomniczi A, Wright H, Ojeda SR. Epigenetic regulation of female puberty. Front Neuroendocrinol. 2015;36:90–107.View ArticlePubMedGoogle Scholar
  6. Mattick JS, Rinn JL. Discovery and annotation of long noncoding RNAs. Nat Struct Mol Biol. 2015;22(1):5–7.View ArticlePubMedGoogle Scholar
  7. Atkinson SR, Marguerata S, Bähler J. Exploring long non-coding RNAs through sequencing. Semin Cell Dev Biol. 2012;23(2):200–5.View ArticlePubMedGoogle Scholar
  8. Volders PJ, Verheggen K, Menschaert G, Vandepoele K, Martens L, Vandesompele J, Mestdagh P. An update on LNCipedia: a database for annotated human lncRNA sequences. Nucleic Acids Res. 2015;43(Database issue):D174–180.View ArticlePubMedGoogle Scholar
  9. Quek XC, Thomson DW, Maag JL, Bartonicek N, Signal B, Clark MB, Gloss BS, Dinger ME. lncRNAdb v2.0: expanding the reference database for functional long noncoding RNAs. Nucleic Acids Res. 2015;43(Database issue):D168–173.View ArticlePubMedGoogle Scholar
  10. Huang W, Long N, Khatib H. Genome-wide identification and initial characterization of bovine long non-coding RNAs from EST data. Anim Genet. 2012;43(6):674–82.View ArticlePubMedGoogle Scholar
  11. Billerey C, Boussaha M, Esquerré D, Rebours E, Djari A, Meersseman C, Klopp C, Gautheret D, Rocha D. Identification of large intergenic non-coding RNAs in bovine muscle using next-generation transcriptomic sequencing. BMC Genomics. 2014;15:499.View ArticlePubMedPubMed CentralGoogle Scholar
  12. Weikard R, Hadlich F, Kuehn C. Identification of novel transcripts and noncoding RNAs in bovine skin by deep next generation sequencing. BMC Genomics. 2013;14:789.View ArticlePubMedPubMed CentralGoogle Scholar
  13. Zhou ZY, Li AM, Adeola AC, Liu YH, Irwin DM, Xie HB, Zhang YP. Genome-wide identification of long intergenic noncoding RNA genes and their potential association with domestication in pigs. Genome Biol Evol. 2014;6(6):1387–92.View ArticlePubMedPubMed CentralGoogle Scholar
  14. Zhao W, Mu Y, Ma L, Wang C, Tang Z, Yang S, Zhou R, Hu X, Li MH, Li K. Systematic identification and characterization of long intergenic non-coding RNAs in fetal porcine skeletal muscle development. Sci Rep. 2015;5:8957.View ArticlePubMedPubMed CentralGoogle Scholar
  15. Yang G, Lu X, Yuan L. LncRNA: A link between RNA and cancer. Biochim Biophys Acta. Biochim Biophys Acta. 2014;1839(11):1097–109.View ArticlePubMedGoogle Scholar
  16. Dantas A, Siqueira ER, Fernandes S, Oba E, Castilho AM, Meirelles PRL, Sartori MMP, Santos PTR. Influence of feeding differentiation on the age at onset of puberty in Brazilian Bergamasca dairy ewe lambs. Arq Bras Med Vet Zootec. 2016;68:22–8.View ArticleGoogle Scholar
  17. Pan L, Ma J, Pan F, Zhao D. Gao J: long non-coding RNA expression profiling in aging rats with erectile dysfunction. Cell Physiol Biochem. 2015;37(4):1513–26.View ArticlePubMedGoogle Scholar
  18. Kiewe P, Gueller S, Komor M, Stroux A, Thiel E, Hofmann WK. Prediction of qualitative outcome of oligonucleotide microarray hybridization by measurement of RNA integrity using the 2100 Bioanalyzer capillary electrophoresis system. Ann Hematol. 2009;88(12):1177–83.View ArticlePubMedGoogle Scholar
  19. Li P, Conley A, Zhang H, Kim HL. Whole-Transcriptome profiling of formalin-fixed, paraffin-embedded renal cell carcinoma by RNA-seq. BMC Genomics. 2014;15:1087.View ArticlePubMedPubMed CentralGoogle Scholar
  20. Ren H, Wang G, Chen L, Jiang J, Liu L, Li N, Zhao J, Sun X, Zhou P. Genome-wide analysis of long non-coding RNAs at early stage of skin pigmentation in goats (Capra hircus). BMC Genomics. 2016;17(1):67.View ArticlePubMedPubMed CentralGoogle Scholar
  21. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.View ArticlePubMedPubMed CentralGoogle Scholar
  22. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–U354.View ArticlePubMedPubMed CentralGoogle Scholar
  23. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.View ArticlePubMedPubMed CentralGoogle Scholar
  24. Guttman M, Garber M, Levin JZ, Donaghey J, Robinson J, Adiconis X, Fan L, Koziol MJ, Gnirke A, Nusbaum C, et al. Ab initio reconstruction of cell type-specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs (vol 28, pg 503, 2010). Nat Biotechnol. 2010;28(7):756.View ArticleGoogle Scholar
  25. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–515.Google Scholar
  26. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.View ArticlePubMedPubMed CentralGoogle Scholar
  27. Flegel C, Manteniotis S, Osthold S, Hatt H, Gisselmann G. Expression profile of ectopic olfactory receptors determined by deep sequencing. PLoS One. 2013;8(2):e55368.Google Scholar
  28. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.View ArticlePubMedPubMed CentralGoogle Scholar
  29. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35:W345–9.View ArticlePubMedPubMed CentralGoogle Scholar
  30. Bateman A, Birney E, Durbin R, Eddy SR, Howe KL, Sonnhammer ELL. The Pfam Protein Families Database. Nucleic Acids Res. 2000;28(1):263–6.View ArticlePubMedPubMed CentralGoogle Scholar
  31. Gomez JA, Wapinski OL, Yang YW, Bureau JF, Gopinath S, Monack DM, Chang HY, Brahic M, Kirkegaard K. The NeST Long ncRNA Controls Microbial Susceptibility and Epigenetic Activation of the Interferon-gamma Locus. Cell. 2013;152(4):743–54.View ArticlePubMedPubMed CentralGoogle Scholar
  32. Lai F, Orom UA, Cesaroni M, Beringer M, Taatjes DJ, Blobel GA, Shiekhattar R. Activating RNAs associate with Mediator to enhance chromatin architecture and transcription. Nature. 2013;494(7438):497–501.View ArticlePubMedPubMed CentralGoogle Scholar
  33. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13.View ArticleGoogle Scholar
  34. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.View ArticlePubMedPubMed CentralGoogle Scholar
  35. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–484.PubMedGoogle Scholar
  36. Mao XZ, Cai T, Olyarchuk JG, Wei LP. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.View ArticlePubMedGoogle Scholar
  37. Wang Kevin C, Chang Howard Y. Molecular Mechanisms of Long Noncoding RNAs. Mol Cell. 2011;43(6):904–14.View ArticlePubMedPubMed CentralGoogle Scholar
  38. Bachelot A, Carre N, Mialon O, Matelot M, Servel N, Monget P, Ahtiainen P, Huhtaniemi I, Binart N. The permissive role of prolactin as a regulator of luteinizing hormone action in the female mouse ovary and extragonadal tumorigenesis. Am J Physiol Endocrinol Metab. 2013;305(7):E845–52.View ArticlePubMedGoogle Scholar
  39. Greenwald-Yarnell ML, Marsh C, Allison MB, Patterson CM, Kasper C, MacKenzie A, Cravo R, Elias CF, Moenter SM, Myers Jr MG. ERalpha in Tac2 Neurons Regulates Puberty Onset in Female Mice. Endocrinology. 2016;157(4):1555–65.View ArticlePubMedGoogle Scholar
  40. Lomniczi A, Wright H, Castellano JM, Matagne V, Toro CA, Ramaswamy S, Plant TM, Ojeda SR. Epigenetic regulation of puberty via Zinc finger protein-mediated transcriptional repression. Nat Commun. 2015;6:10195.View ArticlePubMedPubMed CentralGoogle Scholar
  41. Zhu H, Shah S, Shyh-Chang N, Shinoda G, Einhorn WS, Viswanathan SR, Takeuchi A, Grasemann C, Rinn JL, Lopez MF, et al. Lin28a transgenic mice manifest size and puberty phenotypes identified in human genetic association studies. Nat Genet. 2010;42(7):626–U106.View ArticlePubMedPubMed CentralGoogle Scholar
  42. Parent A-S, Rasier G, Matagne V, Lomniczi A, Lebrethon M-C, Gérard A, Ojeda SR, Bourguignon J-P. Oxytocin Facilitates Female Sexual Maturation through a Glia-to-Neuron Signaling Pathway. Endocrinology. 2008;149(3):1358–65.View ArticlePubMedGoogle Scholar
  43. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn JL. Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 2015;25(18):1915–27.View ArticleGoogle Scholar
  44. Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, Tramontano A, Bozzoni I. A Long Noncoding RNA Controls Muscle Differentiation by Functioning as a Competing Endogenous RNA. Cell. 2011;147(2):358–69.View ArticlePubMedPubMed CentralGoogle Scholar
  45. Nakagawa S, Shimada M, Yanaka K, Mito M, Arai T, Takahashi E, Fujita Y, Fujimori T, Standaert L, Marine J-C, et al. The lncRNA Neat1 is required for corpus luteum formation and the establishment of pregnancy in a subpopulation of mice. Development. 2014;141(23):4618–27.View ArticlePubMedPubMed CentralGoogle Scholar
  46. Funes S, Hedrick JA, Vassileva G, Markowitz L, Abbondanzo S, Golovko A, Yang S, Monsma FJ, Gustafson EL. The KiSS-1 receptor GPR54 is essential for the development of the murine reproductive system. Biochem Biophys Res Commun. 2003;312(4):1357–63.View ArticlePubMedGoogle Scholar
  47. Kaiser UB, Kuohung W. KiSS-1 and GPR54 as new play-ers in gonadotropin regulation and puberty. Endocrine. 2005;26(3):277–84.View ArticlePubMedGoogle Scholar
  48. d’Anglemont de Tassigny XD, Fagg LA, Dixon JP, Day K, Leitch HG, Hendrick AG, Zahn D, Franceschini I, Caraty A, Carlton MB, et al. Hypogonadotropic hypogonadism in mice lacking a functional Kiss1 gene. Proc Natl Acad Sci U S A. 2007;104(25):10714–9.View ArticlePubMedPubMed CentralGoogle Scholar
  49. Lapatto R, Pallais JC, Zhang D, Chan YM, Mahan A, Cerrato F, Le WW, Hoffman GE, Seminara SB. Kiss1−/− mice exhibit more variable hypogonadism than Gpr54−/− mice. Endocrinology. 2007;148(10):4927–36.View ArticlePubMedGoogle Scholar
  50. Wolfe A, Divall S, Wu S. The regulation of reproductive neuroendocrine function by insulin and insulin-like growth factor-1 (IGF-1). Front Neuroendocrinol. 2014;35(4):558–72.View ArticlePubMedPubMed CentralGoogle Scholar
  51. Camille Melon L, Maguire J. GABAergic regulation of the HPA and HPG axes and the impact of stress on reproductive function. J Steroid Biochem Mol Biol. 2016;160:196–203.View ArticlePubMedGoogle Scholar
  52. Allen MP, Xu M, Zeng C, Tobet SA, Wierman ME. Myocyte enhancer factors-2B and -2C are required for adhesion related kinase repression of neuronal gonadotropin releasing hormone gene expression. J Biol Chem. 2000;275(50):39662–70.View ArticlePubMedGoogle Scholar
  53. Hiney JK, Srivastava VK, Volz CE, Dees WL. Alcohol alters insulin-like growth factor-1-induced transforming growth factor beta1 synthesis in the medial basal hypothalamus of the prepubertal female rat. Alcohol Clin Exp Res. 2014;38(10):2572–8.View ArticlePubMedPubMed CentralGoogle Scholar
  54. Cheong RY, Czieselsky K, Porteous R, Herbison AE. Expression of ESR1 in Glutamatergic and GABAergic Neurons Is Essential for Normal Puberty Onset, Estrogen Feedback, and Fertility in Female Mice. J Neurosci. 2015;35(43):14533–43.View ArticlePubMedGoogle Scholar
  55. Sakamoto K, Wakabayashi Y, Yamamura T, Tanaka T, Takeuchi Y, Mori Y, Okamura H. A population of kisspeptin/neurokinin B neurons in the arcuate nucleus may be the central target of the male effect phenomenon in goats. PLoS One. 2013;8(11):e81017.View ArticlePubMedPubMed CentralGoogle Scholar

Copyright

© The Author(s). 2017