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

Abstract The aim of this study was to identify genes related to precocious puberty expressed in the pituitary of goats at different growth stages by suppression subtractive hybridization (SSH). The pituitary glands from Jining Gray (JG) goats (early puberty) and Liaoning Cashmere (LC) goats (late puberty) at 30, 90, and 180 days were used in this study. To identify differentially expressed genes (DEGs) in the pituitary glands, mRNA was extracted from these tissues, and SSH libraries were constructed and divided into the following groups: juvenile group (30-JG vs. 30-LC, API), puberty group (90-JG vs. 180-LC, BPI), and control group (90-JG vs. 90-LC, EPI). A total of 60, 49, and 58 DEGs were annotated by 222 Gene Ontology (GO) terms and 75 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Most of the DEGs were significantly enriched in GO terms related to ‘structural constituent of ribosome’, ‘translation’ and ‘GTP binding’, and numerous DEGs were also significantly enriched in KEGG terms related to the Jak-STAT signaling and oocyte meiosis pathways. Candidate genes associated with precocious puberty and sexual development were screened from the SSH libraries. These genes were analyzed to determine if they were expressed in the pituitary tissues of the goats at different growth stages and to identify genes that may influence the hypothalamic-pituitary-gonadal (HPG) axis. In this study, we found precocious puberty-related genes (such as PRLP0, EIF5A, and YWHAH) that may be interesting from an evolutionary perspective and that could be investigated for use in future goat breeding programs. Our results provide a valuable dataset that will facilitate further research into the reproductive biology of goats.


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 hypothalamuspituitary-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][4][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 in uenced by genetic regulation [7].
Currently, mutations in the kisspeptin system, MKRN3, and Neurokinin B (NKB) have been identi ed 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 identi ed 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 speci c 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 identi ed to further molecular study to reveal the biology mechanism of goat reproduction.

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.

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 identi ed 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 in uence 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. classi ed 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 pro ling by qRT-PCR [18]. Nikishin et al. con rmed the conclusion of RPLP0 gene as the most stable reference genes in the case of vitri ed/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 774 in the case of human PKN1 (Thr 778 in mouse/rat) which is necessary for the catalytic activation of PKN1/2 14 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 ne-tuning the balance between axonal outgrowth and presynaptic differentiation in the parallel ber-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][27][28][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, identi ed 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 in uence 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-speci c 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.
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][44][45], nucleocytoplasmic transport of RNA [46], cell proliferation [47], in ammation, 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 identi ed 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][53][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-pituitarygonadal (HPG) axis to regulate the goats precocious puberty. [56,57]. The results indicate that oocyte-speci c 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-speci c 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.

Conclusion
In summary, the expression of many genes in pituitary tissues from JG and LC goats at different growth stages signi cantly affected precocious puberty, including 18 related genes. Many genes related with oocyte meiosis and steroid biosynthesis may have signi cant roles in the puberty onset, of which, PRLP0, EIF5A and YWHAH play a regulatory role in the precocious puberty of goats, thus will be our candidates for further research.  Table 1 in this study were housed in open sheepfolds and under the same nutrition condition. The pituitary tissues were collected from each goat, sacri ced 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 puri cation 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 e ciency 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 rst 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 2 O were then mixed with the two samples from the rst hybridization simultaneously and incubated at 68°C overnight.
The nal 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. Ampli ed 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 speci c ampli cation 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 puri ed 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.   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.
[ 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] been associated with early-onset schizophrenia and psychotic bipolar disorder.