Identification of Differentially Expressed Genes Associated With Precocious Puberty by Suppression Subtractive Hybridization in Goat Pituitary Tissues

Background

The cost of reserve ewes occupies a large part of the cost of raising sheep. If the age at puberty could be shortened, the production cost of ewes will be greatly reduced. How puberty is initiated is an enigma that still captivates scientists. The onset of puberty is a complex biological process involving numerous factors under the control of the neuroendocrine pathways that are regulated as part of the hypothalamus-pituitary-gonadal (HPG) axis [1]. The HPG axis is a system comprised of endocrine glands whose function is vital to the regulation of reproduction and associated behaviors [2]. Central neurotransmitters, neurohormones, and environmental cues integrate on the HPG-axis and regulate reproduction and puberty onset [3-5]. The pituitary is an endocrine gland that dynamically regulates peripheral tissues to coordinate fundamental physiological functions such as growth, metabolism, sexual maturity and reproduction. The pituitary regulates these homeostatic processes by interpreting hypothalamic signals and, in response, releases hormones from specialized cells in the anterior pituitary. At the onset of puberty, GnRH stimulates the secretion of the pituitary gonadotropins luteinizing hormone (LH) and follicle-stimulating hormone (FSH). LH and FSH in turn stimulate the ovaries to initiate follicular growth and luteal formation that secretes the sex steroid hormones estrogen (E2) and progesterone (P4) [6].

The timing of puberty in humans and other mammals is strongly influenced by genetic regulation [7]. Currently, mutations in the kisspeptin system, MKRN3, and Neurokinin B (NKB) have been identified in sporadic and familial cases of precocious puberty [8]. The loss of function of MKRN3 results in early puberty, implying an inhibitory role of MKRN3 on GnRH secretion [9]. The discovery of the kisspeptin system as a crucial component for pubertal activation of the hypothalamic-pituitary-gonadal axis occurred in 2003, when loss-of-function mutations of the KISS1R (previously known as GPR54) gene were identified in individuals with isolated hypogonadotropic hypogonadism, establishing KISS1R inactivation as a cause of this disorder [10,11]. The KISS1 positive regulation is believed to be important for the LH surge in females of many mammalian species [12]. NKB belongs to a family of closely related peptides called tachykinins. NKB system is necessary for the activation of the hypothalamic-pituitary-gonadal axis in puberty [13].

Suppression subtractive hybridization (SSH) is an effective method for isolation of specific DNA fragments that can be used to differentiate two closely related species [14]. A key feature of this method is simultaneous normalization and subtraction steps that respectively equalize the abundance of DNA fragments within the target population and exclude sequences common to the two populations being compared (10). In this study, the domestics goats breeds Jining Grey goats (JG) and Liaoning Cashmere goats (LC) at different growth stages (30 day, 90 day and 180 day) were selected as the experiment samples, pituitaries were collected and total RNA was extracted. The differential expression genes library was constructed by SSH, and differential expression sexual precocity related genes were screened. The important candidate genes interaction networks were predicted by STING online database. Several genes were identified to further molecular study to reveal the biology mechanism of goat reproduction.

Results

Construction of cDNA libraries

The results of 1.2% agarose gel electrophoresis of double-stranded cDNA products under different cycles is shown in Fig. 1. Treated group have the brightest bands and the widest range in 23 cycles. Control group have the brightest bands in 33 cycles, but the widest range in 23 cycles.

Construction of SSH libraries

324, 295 and 288 white clones were randomly selected from the API, BPI and EPI libraries respectively. Partial result of bacteria liquid PCR is presented in Fig. 2, most of them are positive clones, and the size of the bands is concentrated at 200 bp to 1000 bp, which conforms to the enzyme digestion effect and meets the requirements of the suppression subtractive hybridization libraries.

Sequencing and sequence analysis

In the three libraries, 265, 217 and 225 positive clones were successfully sequenced in the API, EPI and BPI groups, respectively. Excluding the low-quality sequences, the unmatched sequences, the repeat sequences, and the unannotated protein sequences, 91, 58 and 45 valid ESTs were obtained in API, EPI and BPI groups, respectively (Table S1), respectively. These EST sequences were submitted to NR and UNIPROT database to obtain the annotation information of genes (Fig.3). After clustering analysis of Gene ontology on Level 2, data showed that in the API group cDNA library, within the biological process category, 128 ESTs were classified into 11 categories, comprising: metabolic process (25%), cellular process (25.78%), single-organism process (13.28%), biological regulation (9.38%), developmental process (7.03%), multicellular organismal process (5.47%), response to stimulus (5.47%), signaling (3.13%), cellular component organization or biogenesis (3.13%),localization (1.56%) and biological adhesion(0.78%) (Fig.4 A) . Within the cellular component category, they were classified into 11 categories, comprising: cells (20.77%), cell part (20.77%), organelle (19.13%), macromolecular complex (17.49%), organelle part (12.02%), membrane (4.37%), membrane part (2.19%) etc. (Fig.4 B). Within the molecular function category, they were classified into 6 categories, comprising: binding (38.71%), structural molecule activity (37.1%), catalytic activity (14.52%), enzyme regulator activity (4.84%) channel regulator activity (3.23%) and transporter activity (1.61%) (Fig.4 C). In the EPI group subtracted cDNA library, within the biological process category, 88 ESTs were classified into 10 categories, comprising: cellular process (27.27%), metabolic process (26.14%), single-organism process (14.77%), biological regulation (9.09%), cellular component organization or biogenesis (6.82%), response to stimulus (5.68%), developmental process (3.41%) and localization (3.41%) etc. (Fig.5 A). Regarding the cellular component category, they were classified into 10 categories, including cell (21.85%), cell part (21.85%), macromolecular complex (16.81%), organelle (16.81%), organelle part (11.76%), membrane (6.72%) and membrane part (1.68%) etc. (Fig.4 B). In the molecular function category, they were classified into 4 categories, including binding (39.02%), structural molecule activity (34.15%), catalytic activity (24.39%) and channel regulator activity (2.44%) (Fig.5 C). In the BPI group subtracted cDNA library, within the biological process category, 148 ESTs were classified into 12 categories, comprising: cellular process (24.32%), metabolic process (22.97%), single-organism process (13.51%), biological regulation (10.81%), response to stimulus (8.11%), cellular component organization or biogenesis (4.73%), developmental process (4.73%), multicellular organismal process (4.05%) and localization (2.7%) etc. (Fig.6 A). Regarding the cellular component category, they were classified into 8 categories, including cell part (21.02%), cell (20.38%), organelle (18.47%), macromolecular complex (16.56%), organelle part (11.46%), membrane (7.01%) and membrane part (4.46%) etc. (Fig.6 B). In the molecular function category, they were classified into 6 categories, including structural molecule activity (34.92%), binding (26.98%), catalytic activity (23.81%) and enzyme regulator activity (7.94%) etc. (Fig.6 C).

KEGG Pathway analysis

The differentially expressed annotation genes in the three groups were submitted to KEGG Pathway database (Fig.7). The results revealed that most pathways that differentially expressed known genes involved in were ribosome, PI3K-Akt signaling pathway, RNA transport, measles, legionellosis, insulin signaling pathway and Jak-STAT signaling pathway. The most three pathways that the DEGs in API, EPI, BPI groups participated are ribosome, PI3K-Akt signaling pathway, and RNA transport. Some KEGG pathways in the three groups are related with precocious puberty, for instance, estrogen signaling pathway, cAMP signaling pathway, steroid biosynthesis, oocyte meiosis, neuroactive ligand-receptor interaction. In addition, the DEGs in API groups participate in estrogen signaling pathway. Parts of DEGs in EPI are related with oocyte meiosis and some DEGs in BPI group are related with steroid biosynthesis. Some pathways are related with energy metabolism, including PI3K-Akt signaling pathway, insulin signaling pathway, and mTOR signaling pathway (Table 3).

Integrate analysis of Hypothalamic-pituitary-gonadal (HPG) axis genes

Due to the goat database were not provided in the STRING, the high homology species sheep as the alternatives to analysis the correlation of DEGs. Eighteen candidate genes were related to precocious puberty directly or indirectly (Table 4). Most of the DEGs enriched in the one pathway (shown in Fig.8). In our previous study, the DEGs from hypothalamic and ovary tissues were identified by SSH libraries. Integrate all of the data from the three tissues of hypothalamic-pituitary-gonadal (HPG) axis, the result showed that almost 70% DEGs were enriched in the one pathway (Fig.9).

Discussion

In this study, SSH technology was used to construct the differential gene libraries of pituitary tissues in JG and LC goats at different growth stages. Through sequencing, sequence alignment and statistics, the SSH library of pituitary contains 707 ESTs, 37.48% of which were expressed in API group and 30.69% were expressed in EPI group and 31.82% were expressed in BPI group. At the same time, the GO annotation of these libraries showed that most of the genes were mainly involved in biological process, cellular component or molecular function such as metabolic process, cellular process, catalytic activity, binding, etc., and some genes were also involved in the developmental process. This indicates that the molecular mechanism of different growth stages affecting the precocious puberty of goats is comprehensive and complex, which might regulate the growth, development, metabolism and development of goats through multiple signaling pathways.

According to the GO annotation of SSH library, 18 precocious puberty-related differential expression genes were screened from three groups, 44.44% of them in API group, 11.11% in EPI group and 44.44% in BPI group. The results indicated that most precocious puberty related genes were shown in the same growth stages of different goat breeding. The influence of growth stages is more important for the goat precocious puberty.

Some of these precocious puberty related differential genes were differentially expressed in the three groups. Among them, eighteen genes were appeared more than twice in three groups including RPS25, GNB2L1,RPS12, GH1,PET100,PLP1,RPL10A,RPL37,RPLP0,RPS19,PKN1,RPL23A,RPS10, RPS14,RPS20.S,RPS7, HSPA8 and SLC25A5. RPLP0 protein (ribosomal protein lateral stalk subunit P0) encodes one out of approximately 80 ribosomal proteins in human, which are involved in protein synthesis and apoptosis processes [15]. In a previous study, Ragni et al. classified the RPLP0 as the most stable reference gene in expression assays performed in mesenchymal stem cell differentiation, osteoblasts precursor cells [16]. In the study interactions of cytotoxic T lymphocytes with tumor cells, RPLP0 is stably expressed in melanoma cells [17]. RPLP0 gene as reference genes are stably expressed in GCs in both controls and polycystic ovarian syndrome patients and can be used for normalization in gene expression profiling by qRT-PCR [18]. Nikishin et al. confirmed the conclusion of RPLP0 gene as the most stable reference genes in the case of vitrified/thawed human ovarian cortical tissue [19].

PKN1 (protein kinase N1) is a stress-responsive kinase and a member of the protein kinase novel (PKN) family also known as protein kinase C-related kinases (PRKs) [20, 21]. The interaction of Rho with PKN1 induces a conformational change in PKN1 leading to activation loop phosphorylation by 3-Phosphoinositide-Dependent Kinase-1(PDK1) on Thr⁷⁷⁴ in the case of human PKN1 (Thr⁷⁷⁸ in mouse/rat) which is necessary for the catalytic activation of PKN1/2¹⁴ and critical for the stability of the protein [22,23]. PKN1 is shown to provide a basal cardioprotective function in the face of I/R injury [24]. In the currently study found that PKN1 as a novel key player in fine-tuning the balance between axonal outgrowth and presynaptic differentiation in the parallel fiber–forming (PF-forming) cerebellar granule cells (Cgcs) [25]. In neurons PKN1 is the most abundant isoform and has been implicated in a variety of functions including cytoskeletal organization and neuronal differentiation [26-29]. In the precocious puberty related studies found that PKN1 gene modulated TGFβ and EGF dependent regulation of the responses of HEC-A-1 endometrial cancer cells proliferation, migration, and invasiveness, and therefore is a component of the network signaling downstream of TGFβ and EGF [30]. Attarha et al. used proteomics, systems biology, and immunohistochemistry to explore protein expression in human endometrial tumours, identified PKN1 may be considered as predictive biomarkers of endometrial cancer [31].

SLC25A5 (solute carrier family 25 member 5) facilitated the transport of molecules involved in the urea and citric acid cycles, oxidative phosphorylation, DNA maintenance and iron metabolism processes [32, 33]. In the dysregulated expression of androgen metabolism genes and genetic analysis in hypospadias, the data provided direct evidence that SLC25A5 may contribute the genetic etiology of hypospadias [34]. In the bovine reproduction study, SLC25A5 gene were found high expressed in the SSH library of bovine blastocyst [35]. RPL37 (ribosomal protein L37), as a component of the 60S subunit of ribosomes and belongs to the L37E family of ribosomal proteins. In the Zebra Finch Song System study, RPL37 might all influence sexual differentiation, perhaps with the hormone and proteins interacting, such that an appropriate balance is required for normal development [36]. PRL37 protein as a histotype-specific prognostic biomarker was found in the Ovarian Carcinomas Using Immunohistochemistry [37]. These genes above mentioned were potential molecular markers for the study of goat precocious puberty.

In our experiment, 18 DEGs were related to precocious puberty directly or indirectly. Among them, RPL10A, RPL23A, RPL24, RPL34, RPL37, RPLP0, RPS3A , RPS6 and RPS7 are all common protein family ribosomal protein, they play an important role in the development of ovarian and breast cancers related to precocious puberty [38,39]. And EIF5A, HNRNPDL, HSPA8, PKN1, PRKAR2B, SLC40A1, SLC25A5, SNRPD2 and YWHAH, play key roles in oocyte meiosis, steroid biosynthesis and conduction of reproduction process [40–42].

Precocious puberty is regulated via hormones of the hypothalamic-pituitary-gonadal (HPG) axis. In our study, we analysis of the hypothalamic-pituitary-gonadal (HPG) axis related genes by the online tool STRING. We found that most of the DEGs from the three tissues enriched in one pathway, the EIF5A and YWHAH gene play an important role in the network (Fig.9).

eIF5A (eukaryotic translation initiation factor 5A) is a highly conserved 17-kDa protein that is expressed ubiquitously in all cells. eIF5A participates in diverse cellular processes, including protein translation [43-45], nucleocytoplasmic transport of RNA [46], cell proliferation [47], inflammation, and apoptosis[48,49]. Stimulation of follicular growth by follicle-stimulating hormone (FSH) is associated with increased expression of luteinizing hormone/choriogonadotropin receptors (LHCGRs) in granulosa cells [50, 51]. EIF5A was identified as one of the proteins that interact with LRBP and lead to the degradation of LHCGR mRNA by facilitating the transport of LHCGR mRNA-LRBP complex to P bodies for degradation [52-54]. In the currently, a series studies were carried out to determine the mechanism by which inhibition of hypusination of EIF5A causes an increase in LHCGR mRNA expression [55]. The EIF5A gene usually expressed in ovary tissue, in our study the expression of EIF5A gene were detected in the pituitary tissue. The result indicated that this gene might as an important gene expressed in the hypothalamic-pituitary-gonadal (HPG) axis to regulate the goats precocious puberty.

YWHAH (tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein eta) is an isoform of YWHA protein. Numerous studies implicate YWHA as a critical regulator of the cell cycle in meiotic and mitotic cells as well as other cellular processes [56, 57]. The results indicate that oocyte-specific or global elimination of YWHAH protein does not result in abnormal fertility, oocyte maturation or development. Breeding and development of pups was normal in the absence of YWHAH or YWHAE in females with oocyte-specific knockout of these genes. Global inactivation of Ywhah in female mice does not appear to alter oogenesis, oocyte maturation and early development [58]. The previously noted that YWHAH may play an important role in meiotic spindle formation by employing antisense morpholino knockdown approaches [59].

The effect of differential expressed genes from pituitary tissues at different growth stages may be related to the process of precocious puberty. However, the exact pathway remains unknown, and there are few reports on how these differential genes affect the precocious puberty of goats. In the future, more function analysis will be needed to further explore the effects of these genes on precocious puberty in goats.

Methods

Ethical statement

All experimental procedures mentioned in the present study were approved by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China). Ethical approval was provided by the animal ethics committee of IASCAAS (No. IAS2019-63).

Animals and tissue collection

All the Jining Gray (JG) goats (Jining Grey goats conservation base, Jiaxiang County, Shandong Province, People’s Republic of China) and Liaoning Cashmere (LC) goats (Liaoning Cashmere Goat Breeding Center, Liaoyang City, Liaoning Province, People’s Republic of China)displayed in Table 1 in this study were housed in open sheepfolds and under the same nutrition condition. The pituitary tissues were collected from each goat, sacrificed after anesthesia with 3% pentobarbital sodium salt injection (20mg/kg body weight) (Merck, Darmstadt, Germany), and preserved in RNA later RNA Stabilization reagent (Qiagen, Hilden, Germany) and kept at -20°C until RNA isolation.

poly A⁺RNA purification and concentration

According to the instruction of Dynabeads mRNA DIRECT™ Kit (Invitrogen, Inc. USA), the poly A⁺ RNA was isolated from pituitary. After mixed in equal quantity from three goats in each group, the poly A⁺RNA concentrated using RNA clean & concentration-5 mRNA (Zymo Research, Orange, CA, USA) to a suitable concentration.

Suppression subtractive hybridization (SSH)

The concentrated poly A⁺ RNA from the pituitaries of JG goats was tester and that of LC goats was driver. As shown in Table 2, the poly A⁺ RNA from the pituitaries of 30-JG was hybridized with that of 30-LC goats and this hybridization was named as A group. The hybridizations between 90-JG and 180-LC, 90-JG and 90-LC were names as B group and E group, respectively.

The regents for SSH were from PCR-Select cDNA Subtraction Kit (Clontech, palo alto, CA) and the procedure were carried out according to the protocols. Simply, the double strand cDNA (dscDNA) were synthesized from 2 μg poly A⁺ RNA of tester and driver respectively and then were digested by Rsa Ⅰ restriction endonuclease to obtain shorter, blunt-ends dscDNA fragments, required for adaptor ligation and optimal for subtraction. After analysis of Rsa Ⅰ digestion, the digested blunt-ends of tester cDNA, but not adaptor cDNA, were divided into two parts and ligated with two different cDNA adaptors (adaptor 1 5′-CTAATACGACTCACTATAGGGCTCGAGCGGCCGCCCGGGCAGGT-3′ and adaptor 2R 5′-CTAATACGACTCACTATAGGGCAGCGTGGTCGCGGCCGAGGT-3′), respectively. After ligation, the ligation efficiency analysis was performed by PCR experiment with GAPDH primers (GAPDHF: 5′- AGGCTGGGGCTCACTTGAAG-3′, GAPDHR: ATGGCGTGGACAGTGGTCAT-3′) (goat GAPDH mRNA sequence-GenBank: AJ431207) and PCR primer I (5′-CTAATACGACTCACTATAGGGC-3′) using the Advantage cDNA PCR kit (Clontech, palo alto, CA).

Two cycles of hybridization were followed after the ligation. In the first hybridization, 1.5 μL Rsa I-Digested Driver cDNA and 1.5μl Adaptor-Ligated Tester were hybridized in 1.0 μL 4×Hybridization buffer at 68°C for 8 hours. For the second hybridization process, 1.0 μL fresh denatured driver cDNA, 1.0 μL 4×Hybridization Buffer and 2.0 μL ddH₂O were then mixed with the two samples from the first hybridization simultaneously and incubated at 68°C overnight.

The final hybridization solution (also called the subtracted library) was employed as a template to amplify the differentially expressed sequences in the tester population by using a set of PCR primer1 and was followed by nested PCR primers (Nested PCR primer 1 5′-TCGAGCGGCCGCCCGGGCAGGT-3′, Nested PCR primer2R 5′-AGCGTGGTCGCGGCCGAGGT-3′), which does not exponentially amplify the non-adaptor (derived from driver cDNA), cDNA with the one adaptor on either end (derived from tester cDNA hybridized with driver cDNA), or cDNA with the same adaptor on both ends (derived from relatively abundant tester cDNA). For the primary PCR, 1 μL sample was added to 24 μL PCR master mix prepared using the reagents supplied in the kit, and cycling conditions commenced as follows: 75 °C for 5 min to extend the adaptors; 94 °C for 25 sec; and 27 cycles at 94 °C for 10 sec, 66 °C for 30 sec, and 72 °C for 1.5 min. Amplified products were diluted 10-fold in sterile water and 1 μL of diluted primary PCR products were added to 24 μL of secondary PCR master mix containing nested primers, 1 and 2R, to ensure specific amplification of double-stranded templates containing both adaptors. Secondary PCR was performed 12 cycles at 94 °C for 10 sec, 68 °C for 30 sec, and 72 °C for 1.5 min. Primary and secondary PCR products were analyzed on a 2% agarose gel. And now the PCR mixture is enriched for differentially expressed cDNA in tester sample.

Cloning, sequencing and sequence analysis

3 μL of purified PCR products was ligated into pGEM-T Easy Vector (Promega, Madison, WI, USA) and transformed into the competent E. coli DH5α. Positive mono clones were detected by PCR and sequenced by Invitrogen Corporation.

The vector and adaptor sequence of sequences obtained by sequencing were remove using the UniVec database (http://www.ncbi.nlm.nih.gov/VecScreen/) to get clean sequence.

The differentially expressed ESTs were classified into known genes, known ESTs and unknown ESTs according to homologies analyzed by BLAST (www.ncbi.nlm.nih.gov/blast/blast.cgi). The differentially expressed genes (DEGs) function predicted in the Cluster of Orthologous Groups of proteins database (COG) (http://www.ncbi.nlm.nih.gov/COG/) and were classified according to the COG database. Then the DEGs were submitted to the KEGG pathway database (http://www.genome.jp/kegg/pathway.html) to analysis which pathway they participate in. The protein-protein interactions were predicted in the STRING (Search Tool for Retrieval of Interacting Genes/Proteins) net (http://string.embl.de/).

Declarations

Acknowledgement

Not Applicable.

Funding

This research was funded by the National Natural Science Foundation of China (30871773, 31201778), the Earmarked Fund for China Agriculture Research System (CARS-38), the Agricultural Science and Technology Innovation Program of China (ASTIP-IAS13), Shandong Provincial Natural Science Foundation (ZR2017BC054).

Author information

Affiliation

Key Laboratory of Animal Genetics, Breeding and Reproduction of Ministry of Agriculture and Rural Affairs, China Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing 100193, China

Yufang Liu&Mingxing Chu

College of Life Sciences and Food Engineering ,Hebei University of Engineering, Handan 056021, China

Yufang Liu

College of Agriculture, Liaocheng University, Liaocheng 252059, China

Guiling Cao& Yujing Xie

Contributions

Conceived and designed the experiments: GLC, MXC; Performed the experiments: GLC, YJX, YFL; Analyzed the data: GLC, YJX, YFL; Contributed reagents: MXC; Acquisition of data: GLC, YJX, YFL; YFL and MXC drafted the manuscript. YFL, GLC and MXC contributed in writing the paper; YFL, GLC and MXC substantively revised the paper; Study initiation: GLC, YJX, YFL and MXC; Read and approved the final version of the paper: all co-authors.

Corresponding author

Correspondence to Mingxing Chu.

Ethics declarations

Ethics approval and consent to participate

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

The datasets used and analyzed during the current study available from the corresponding author on reasonable request.

References

Tena-Sempere M. Deciphering puberty: novel partners, novel mechanisms. Eur J Endocrinol. 2012; 167(6):733-47.
Sower SA, Freamat M, Kavanaugh SI. The origins of the vertebrate hypothalamic-pituitary-gonadal (HPG) and hypothalamic-pituitary-thyroid (HPT) endocrine systems: new insights from lampreys. General and Comparative Endocrinology. 2009; 161(1): 20-29.
Parent AS, Teilmann G, Juul A, et al. The timing of normal puberty and the age limits of sexual precocity: variations around the world, secular trends, and changes after migration. Endocr Rev. 2003; 24(5):668-93.
Ullah R, Su Y, Shen Y, et al. Postnatal feeding with high-fat diet induces obesity and precocious puberty in C57BL/6J mouse pups: a novel model of obesity and puberty. Front Med. 2017; 11(2):266-276.
Ullah R, Shen Y, Zhou YD, et al. Expression and actions of GnIH and its orthologs in vertebrates: current status and advanced knowledge. Neuropeptides. 2016, 59:9-20.
Ojeda SR, Lomniczi A, Mastronardi C, et al. Minireview: the neuroendocrine regulation of puberty: is the time ripe for a systems biology approach? Endocrinology. 2006; 147(3):1166-74.
Howard S.R. The genetic basis of delayed puberty. Front Endocrinol (Lausanne). 2019; 10: 423.
Bulcao MD, Nahime BV, Latronico AC. New causes of central precocious puberty: the role of genetic factors. Neuroendocrinology. 2014; 100:1-8
Hughes IA. Releasing the brake on puberty. N Engl J Med. 2013; 368(26):2513-5.
De Roux N, Genin E, Carel JC, et al. Hypogonadotropic hypogonadism due to loss of function of the KiSS1-derived peptide receptor GPR54. Proc Natl Acad Sci U S A. 2003; 100(19):10972-6
Seminara SB, Messager S, Chatzidaki EE, et al. The GPR54 gene as a regulator of puberty. N Engl J Med. 2003; 349(17):1614-27.
Adachi S, Yamada S, Takatsu Y, et al. Involvement of anteroventral periventricular metastin/kisspeptin neurons in estrogen positive feedback action on luteinizing hormone release in female rats. J Reprod Dev. 2007; 53(2):367-78.
Gianetti E, Tusset C, Noel SD, et al. TAC3/TACR3 mutations reveal preferential activation of gonadotropin-releasing hormone release by neurokinin B in neonatal life followed by reversal in adulthood. J Clin Endocrinol Metab. 2010; 95(6):2857-67.
Rebrikov DV, Desai SM, Siebert PD, et al. Suppression subtractive hybridization. Methods Mol Biol. 2004; 258(258): 107-134.
Artero-Castro A, Castellvi J, García A, et al. Expression of the ribosomal proteins Rplp0, Rplp1, and Rplp2 in gynecologic tumors. Hum Pathol. 2011; 42(2):194-203.
Ragni E, Viganò M, Rebulla P, et al. What is beyond a qRT-PCR study on mesenchymal stem cell differentiation properties: how to choose the most reliable housekeeping genes. J Cell Mol Med. 2013; 17(1):168-80.
Lv Y, Zhao SG, Lu G, et al. Identification of reference genes for qRT-PCR in granulosa cells of healthy women and polycystic ovarian syndrome patients. Sci Rep. 2017; 7: 6961.
Neubert NJ, Soneson C, Barras D, et al. A well-controlled experimental system to study interactions of cytotoxic T lymphocytes with tumor cells. Front Immunol. 2016; 7: 326.
Nikishin DA, Filatov MA, Kiseleva MV, et al. Selection of stable expressed reference genes in native and vitrified/thawed human ovarian tissue for analysis by qRT-PCR and western blot. J Assist Reprod Genet. 2018; 35(10): 1851-60.
Palmer RH, Ridden J, Parker PJ. Identification of multiple, novel, protein kinase C-related gene products. FEBS Lett. 1994; 356(1):5-8.
Balendran A, Biondi RM, Cheung PC, et al. A 3-phosphoinositide-dependent protein kinase-1 (PDK1) docking site is required for the phosphorylation of protein kinase Czeta (PKCzeta ) and PKC-related kinase 2 by PDK1. J Biol Chem. 2000, 275(27):20806-13.
Francois AA, Obasanjo-Blackshire K, Clark JE, et al. Loss of Protein Kinase Novel 1 (PKN1) is associated with mild systolic and diastolic contractile dysfunction, increased phospholamban Thr17 phosphorylation, and exacerbated ischaemia-reperfusion injury. Cardiovasc Res. 2018; 114(1): 138-57.
Nedden S, Eith R, Schwarzer C, et al. Protein kinase N1 critically regulates cerebellar development and long-term function. J Clin Invest. 2018; 128(5): 2076-2088.
Mukai H. The structure and function of PKN, a protein kinase having a catalytic domain homologous to that of PKC. J Biochem. 2003; 133(1):17-27.
Gudi T, Chen JC, Casteel DE, et al. cGMP-dependent protein kinase inhibits serum-response element-dependent transcription by inhibiting rho activation and functions. J. Biol. Chem. 2002; 277(40):37382-93.
Sihag RK, Inagaki M, Yamaguchi T, et al. Role of phosphorylation on the structural dynamics and function of types III and IV intermediate filaments. Exp. Cell Res. 2007; 313(10):2098-2109.
Thauerer B, zur Nedden S, Baier-Bitterlich G. Vital role of protein kinase C-related kinase in the formation and stability of neurites during hypoxia. J. Neurochem. 2010; 113(2):432-46.
Sanaz Attarha, Ravi Kanth Rao Saini, Sonia Andersson, et al. PKN1 modulates TGFβ and EGF signaling in HEC-1-A endometrial cancer cell line. Onco Targets Ther. 2014; 7: 1397-1408.
Attarha S, Andersson S, Mints M, et al. Individualized proteome profiling of human endometrial tumours improves detection of new prognostic markers. Br J Cancer. 2013; 109(3): 704-713.
Palmieri F. The mitochondrial transporter family SLC25: identification, properties and physiopathology. Mol Aspects Med. 2013; 34(2-3):465-84.
Chen ZZ, Lin XL, Wang YP, et al. Dysregulated expression of androgen metabolism genes and genetic analysis in hypospadias. Mol Genet Genomic Med. 2020; 8(8): e1346.
Goossens K, Van Soom A, Van Poucke M, et al. Identification and expression analysis of genes associated with bovine blastocyst formation. BMC Dev Biol. 2007; 7: 64.
Tang YP, Wade JL. Effects of estradiol on incorporation of new cells in the developing zebra finch song system: potential relationship to expression of ribosomal proteins L17 and L37. Dev Neurobiol. 2009; 69(7): 462-475.
Engqvist H, Parris TZ, Kovács A, et al. Validation of novel prognostic biomarkers for early-stage clear-cell, endometrioid and mucinous ovarian carcinomas using immunohistochemistry. Front Oncol. 2020; 10: 162.
Gari HH, Gearheart CM, Fosmire S, et al. Genome-wide functional genetic screen with the anticancer agent AMPI-109 identifies PRL-3 as an oncogenic driver in triple-negative breast cancers. Oncotarget. 2016; 7(13): 15757-15771.
Shan NC, Zhou W, Zhang SF, et al. Identification of HSPA8 as a candidate biomarker for endometrial carcinoma by using iTRAQ-based proteomic analysis. Onco Targets Ther. 2016; 9: 2169-2179.
Saini P, Eyler DE, Green R, et al. Hypusine-containing protein eIF5A promotes translation elongation. Nature. 2009; 459(7243):118-121.
Rosorius O, Reichart B, Krätzer F, et al. Nuclear pore localization and nucleocytoplasmic transport of eIF-5A: evidence for direct interaction with the export receptor CRM1. J Cell Sci. 1999; 112(14):2369-2380.
Park MH, Joe YA, Kang KR, et al. The polyamine-derived amino acid hypusine: its post-translational formation in eIF-5A and its role in cell proliferation. Amino Acids. 1996; 10(2):109-121.
Taylor CA, Sun Z, Cliche DO, et al. Eukaryotic translation initiation factor 5A induces apoptosis in colon cancer cells and associates with the nucleus in response to tumors necrosis factor alpha signaling. Exp Cell Res. 2007; 313(3):437-449.
Lu DL, Peegel H, Mosier SM, et al. Loss of lutropin/human choriogonadotropin receptor messenger ribonucleic acid during ligand-induced down-regulation occurs post transcriptionally. Endocrinology. 1993; 132(1):235-240.
Gulappa T, Menon B, Menon KMJ. 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. 2015; 413:90-95.
Gulappa T, Menon B, Menon KMJ. LHCGR expression during follicle stimulating hormone–induced follicle growth is negatively regulated by eukaryotic initiation factor 5A. Endocrinology. 2017; 158(8): 2672-2679.
Mackintosh C. Dynamic interactions between 14-3-3 proteins and phosphoproteins regulate diverse cellular processes. Biochem J. 2004; 381(Pt 2):329-42.
Eisa AA, De Santanu, Detwiler A, et al. YWHA (14-3-3) protein isoforms and their interactions with CDC25B phosphatase in mouse oogenesis and oocyte maturation. BMC Dev Biol. 2019; 19: 20.
De S, Kline D. Evidence for the requirement of 14-3-3eta (YWHAH) in meiotic spindle assembly during mouse oocyte maturation. BMC Dev Biol. 2013; 13:10.
Mémin E, Hoque M, Jain MR, et al. Blocking eIF5A modification in cervical cancer cells alters the expression of cancer-related genes and suppresses cell proliferation. Cancer Res. 2014; 74(2): 552-562.
Wang CX, Zheng MJ, Wang S, et al. Whole genome analysis and prognostic model construction based on alternative splicing events in endometrial cancer. Biomed Res Int. 2019; 2019: 2686875.
James RG, Bosch KA, Kulikauskas RM, et al. Protein Kinase PKN1 Represses Wnt/β-Catenin signaling in human melanoma cells. J Biol Chem. 2013; 288(48): 34658-34670.
Bettina Thauerer, Stephanie zur Nedden, Gabriele Baier-Bitterlich. Protein Kinase C-Related Kinase (PKN/PRK). Potential key-role for PKN1 in protection of hypoxic neurons. Curr Neuropharmacol. 2014; 12(3): 213-218.
Manser C, Stevenson A, Banner S, et al. Deregulation of PKN1 activity disrupts neurofilament organization and axonal transport. FEBS Lett. 2008; 582(15): 2303-2308.
Zhang CX, Chen ZY, Yin QZ, et al. The chromatin remodeler Snf2h is essential for oocyte meiotic cell cycle progression. Genes & Development. 2020; 34(3-4): 166–178.
Monique Tomazele Rovani, Bernardo Garziera Gasperin, Gustavo Freitas Ilha, Rogério Ferreira, et al. Expression and molecular consequences of inhibition of estrogen receptors in granulosa cells of bovine follicles. Journal of Ovarian Research. 2014; 7: 96.
Kunwadee Palasin, Tamayo Uechi, Maki Yoshihama, et al. Abnormal development of zebrafish after knockout and knockdown of ribosomal protein L10a. Scientific Reports. 2019; 9: 18130.
Denisa Jansova, Anna Tetkova, Marketa Koncicka, et al. Localization of RNA and translation in the mammalian oocyte and embryo. PLoS ONE. 2018; 13(3): e0192544.
Nasser Ghanem, Michael Hölker, Franca Rings, et al. Alterations in transcript abundance of bovine oocytes recovered at growth and dominance phases of the first follicular wave. BMC Dev Biol. 2007; 7: 90.
Pokharel K, Peippo J, Weldenegodguad M, et al. Gene expression profiling of corpus luteum reveals important insights about early pregnancy in domestic sheep. Genes (Basel). 2020; 11(4): 415.
Lin PF, Lan XL, Chen FL, et al. Reference gene selection for real-time quantitative PCR analysis of the mouse uterus in the peri-implantation period. PLoS ONE. 2013; 8(4): e62462.
Cai J, Li T, Huang BX, et al. The use of laser microdissection in the identification of suitable reference genes for normalization of quantitative real-time PCR in human FFPE epithelial ovarian tissue samples. PLoS ONE. 2014; 9(4): e95974.
Vijay S Baddela, Anja Baufeld, Vengala R Yenuganti, et al. Suitable housekeeping genes for normalization of transcript abundance analysis by real-time RT-PCR in cultured bovine granulosa cells during hypoxia and differential cell plating density. Reprod Biol Endocrinol. 2014; 12: 118.
Li N, Qian SH, Li B, et al. Quantitative analysis of the human ovarian carcinoma mitochondrial phosphoproteome. Aging (Albany NY). 2019; 11(16): 6449-6468.
Kim M, Sim C, Denlinger DL. RNA interference directed against ribosomal protein S3a suggests a link between this gene and arrested ovarian development during adult diapause in Culex pipiens. Insect Mol Biol. 2010; 19(1): 27-33.
Chermuła B, Brązert M, Jeseta M, et al. The unique mechanisms of cellular proliferation, migration and apoptosis are regulated through oocyte maturational development—a complete transcriptomic and histochemical study. Int J Mol Sci. 2019; 20(1): 84.
Yang XQ, Xu LX, Yang Y, et al. Knockdown of ribosomal protein S6 suppresses proliferation, migration, and invasion in epithelial ovarian cancer. J Ovarian Res. 2020; 13: 100.
Wang ZL, Hou J, Lu LL, et al. Small ribosomal protein subunit S7 suppresses ovarian tumorigenesis through regulation of the PI3K/AKT and MAPK pathways. PLoS ONE. 2013; 8(11): e79117.
Manoel Francisco de Sá Filho, Angela Maria Gonella-Diaza, Mariana Sponchiado, et al. Impact of hormonal modulation at proestrus on ovarian responses and uterine gene expression of suckled anestrous beef cows. J Anim Sci Biotechnol. 2017; 8: 79.

Tables

Table 1 The detail information of experiment samples

Goat breed	Age	Abbreviation	Number	Development stage
Jining Grey goats	30-day-old	30-JG	3	Juvenile
Jining Grey goats	90-day-old	90-JG	3	Puberty
Liaoning Cashmere goats	30-day-old	30-LC	3	Juvenile
	90-day-old	90-LC	3	Juvenile
	180-day-old	180-LC	3	Puberty

Table 2 Suppression subtractive hybridization groups. JG: Jining Grey goats; LC: Liaoning Cashmere goats.

Hybridization Groups	Tester	Driver
A	30-JG	30-LC
B	90-JG	180-LC
E	90-JG	90-LC

Table 3 The genes enriched in different KEGG function in the three SSH groups

KOG function	AO	EO	BO
PI3K-Akt signaling pathway	GH1	LOC102322233 GH1	LOC102322233 bPKN1 RPS6
Jak-STAT signaling pathway	GH1	LOC102322233 GH1	LOC102322233
MAPK signaling pathway	HSPA8
Wnt signaling pathway	LOC102250913 LOC102324853
Estrogen signaling pathway	HSPA8
cGMP-PKG signaling pathway	LOC102250913	SLC25A5	Slc25a5
mTOR signaling pathway	LOC100073610		RPS6
TGF-beta signaling pathway	LOC102250913
Insulin signaling pathway			PRKAR2B EXOC7 RPS6
Steroid biosynthesis			SC5D
Oocyte meiosis		YWHAH

Table 4 Precocious puberty related important DEGs among three groups

Gene name	Description	Reference
EIF5A (eukaryotic translation initiation factor 5A)	mRNA-binding protein involved in translation elongation	[54-55,60]
HNRNPDL(heterogeneous nuclear ribonucleoprotein D like)	Involved in regulation of gene expression	[61]
HSPA8(heat shock protein family A (Hsp70) member 8)	A chaperone, and binds to nascent polypeptides to facilitate correct folding	[42]
PKN1 (protein kinase N1)	mediate insulin signals to the actin cytoskeleton	[30-31, 62-64]
PRKAR2B (protein kinase cAMP-dependent type II regulatory subunit beta)	important for a variety of cellular functions	[65,66]
RPL10A (ribosomal protein L10a)	the organelles that catalyze protein synthesis	[67]
RPL23A (ribosomal protein L23a)	one of the target molecules involved in mediating growth inhibition by interferon	[38]
RPL24 (ribosomal protein L24)	The protein belongs to the L24E family of ribosomal proteins.	[68,69]
RPL34 (ribosomal protein L34)	Overexpression of this gene has been observed in some cancer cells.	[70]
RPL37 (ribosomal protein L37)	The protein contains a C2C2-type zinc finger-like motif.	[37]
RPLP0 (ribosomal protein lateral stalk subunit P0)	It is a neutral phosphoprotein	[18-19,71-74]
RPS3A (ribosomal protein S3A)	The protein belongs to the S3AE family of ribosomal proteins.	[75,76]
RPS6 (ribosomal protein S6)	It is the major substrate of protein kinases in the ribosome	[77]
RPS7 (ribosomal protein S7)	Exhibits mRNA binding activity	[78]
Slc25a5 (solute carrier family 25 member 5)	Suppressed expression of this gene has been shown to induce apoptosis and inhibit tumor growth.	[35]
SLC40A1 (solute carrier family 40 member 1)	The product of this gene functions as a gated pore that translocates ADP from the cytoplasm into the mitochondrial matrix.	[41]
SNRPD2 (small nuclear ribonucleoprotein D2 polypeptide)	It is required for pre-mRNA splicing and small nuclear ribonucleoprotein biogenesis.	[79]
YWHAH (tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein eta)	This gene contains a 7 bp repeat sequence in its 5' UTR, and changes in the number of this repeat have been associated with early-onset schizophrenia and psychotic bipolar disorder.	[56-59]