Introduction

Puberty is the stage of transition from infancy to sexual maturity. Numerous studies have shown that early-maturing livestock often have shorter reproductive cycles1,2, which can shorten the generation interval and improve the reproductive capacity. The age of pubertal onset is an important indicator reflecting the reproductive performance of goats. However, the basic molecular mechanisms that regulate the onset of puberty in goats are still unclear.

Research on the initiation of puberty has primarily focused on endocrine aspects3,4,5. Generally, the hypothalamic-pituitary-gonadal axis (HPGA) controls and drives pubertal initiation. The release of gonadotropin-releasing hormone (GnRH) from the hypothalamus leads to the secretion of follicle-stimulating hormone (FSH) and luteinizing hormone (LH) from the pituitary gland, which promotes the first ovulation and vaginal opening of female animals6. However, endocrine regulation is influenced by a complex gene network7. Therefore, in order to elucidate the mechanism of pubertal initiation, it is necessary to investigate the expression regulation of genes involved in regulating the onset of puberty.

In recent years, a large amount of genetic data related to pubertal onset has been obtained. For example, studies have shown that some rare mutations in genes such as GnRH receptor (GNRHR)8, KiSS-1 metastasis suppressor (KiSS1)9, Tachykinin Precursor 3(TAC3) and Tachykinin Receptor 3 (TACR3)10 can affect the normal onset of puberty in humans. Researchers found that methylation levels affect the expression of Calcium Voltage-Gated Channel Subunit Alpha1 C (CACNA1C) and Corticotropin Releasing Hormone (CRH) in Jining gray goats11. A study on the transcriptome of cow pituitary and ovaries suggested that Diencephalon/Mesencephalon Homeobox 1 (DMBX1), Complement Component 4 Binding Protein Alpha (C4BPA), and Myelin Associated Glycoprotein (MAG) may involve in the initiation of puberty12. These genes deserve further study to elucidate their roles in the onset of puberty in female mammals.

In addition to the HPGA, the pineal gland also plays an important role in mammalian pubertal onset. The mammalian pineal gland can convert photoperiod information into a night hormone signal, thereby controlling the synthesis and secretion of melatonin13. Before puberty, the pineal gland secretes melatonin and other bioactive substances that inhibit the activity of hypothalamic GnRH neurons. With the development of the body, the inhibitory activity substance secretion gradually reduces, and in turn relieves the inhibition of GnRH neurons in the hypothalamus near puberty14. Studies have shown that pineal gland dysfunction caused by pineal gland tumors, cysts, calcification, and other diseases, is associated with precocious puberty in children13.

The objective of this study is to profile gene expression and genome-scale DNA methylation of the hypothalamus in prepubertal and pubertal goats using RNA sequencing (RNA-seq) and whole genome bisulfite sequencing (WGBS). The results will provide useful insights into the epigenetic mechanisms involved in the onset of puberty in goats.

Materials and methods

Ethics statement

The animal experimental protocols and procedures for this study were approved by the Animal Care and Use Committee of Anhui Agricultural University and complied with the EU Directive 2010/63/EU for animal experiments.

Animals and sampling

We raised twelve healthy prepubertal female Wanlin White goats (WWGs) in a semi-enclosed shed, and fed them three times a day (7:00, 11:30 and 18:00). We fed them with a 1:1 mixture of fattening goat concentrate (Jiangsu Mudi Feeds Co, Ltd, Q/321324JSMD002) and corn in the morning and evening, and give them dried peanuts seedlings at midday. On average, each goat consumed about 0.15 kg of feed per meal. The amount of feed for each goat was gradually increased as body weight. All the goats had free access to water and natural light.

We conducted estrus testing on these goats at 8:00 a.m. and 4:00 p.m daily. According to reference15, we determined the onset of puberty as follows: (1) first acceptance of mounting behavior; (2) reddening and swelling of vulvas with mucus discharge (Fig. 1a2); (3) presence of mature follicles on the surface of their ovaries (Fig. 1a4, b4).

Fig. 1
figure 1

Vulvas, ovaries, and ovarian histological sections of prepuberta; and pubertal goats. (a, a1) is the vulva of the prepubertal goat; (a2) is the vulva of the adolescent goat; (a3) is the ovary of the prepubertal goat; (a4) is the ovary of the pubertal goat. (b) histological sections of the prepubertal ovaries had a large number of primordial follicles (PrF), primary follicles (PF) and secondary follicles (SF) (b1, b2). The ovaries of pubertal goats had a large number of tertiary follicles (TF) and mature follicles (MF) (b3, b4).

Full size image

We collected pineal glands from six prepubescent WWGs (aged 2.5–3.0 months, with a weight of 8.10 ± 0.30 kg), and six pubescent WWGs (aged 4.0–4.5 months, with a weight of 20.16 ± 0.35 kg). To ensure that each sample has sufficient quantity for DNA methylation and transcriptome sequencing analysis, we pooled pineal gland samples from the same period in pairs, ultimately obtained 3 prepubertal samples and 3 pubertal samples. We took out half of each sample from these three prepubertal and three pubertal samples for RNA sequencing. Then, we pooled the remaining prepubertal and pubertal samples together, respectively, for methylome analysis.

RNA sequencing

We extracted the total RNA using TRIzol® reagent following the manufacturer’s instructions (ThermoFisher, A33250), determined the quality of RNA using1% agarose gel, Nanodrop One (ThermoFisher, USA) and Agilent 2100 bioanalyzer as described before16. After enriching mRNAs with Oligo dT magnetic beads and fragmenting them with divalent cations, we constructed mRNA library based on previous studies17 and performed high-throughput sequencing on the Illumina HiSeqTM2500.

Transcriptome analysis

We processed the raw data (raw reads) with in-house Perl scripts to obtain clean data by removing adapter-containing reads, poly-N-containing reads, and low-quality reads (the number of bases with a Qphred score ≤ 20 accounts for more than 50% of the entire read length). We downloaded the goat (Capra hircus) reference genome ARS1 and the corresponding annotation files from the following link: (https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/704/415/GCF_001704415.1_ARS1/), and built the index of the reference genome using Hisat2 v2.0.5 with default parameters18,19, then aligned the clean reads with the reference genome.

We used StringTie for new gene prediction and rMATS software for alternative splice (AS) analysis of RNA-Seq data. We classified AS events and evaluated the differences in AS events between different samples. The significance threshold for filtering differentially spliced events was set at FDR less than 0.05. We used GATK software for variant calling and utilized SnpEff for annotating information on each variant, including variant function, variant region, and variant impact.

Individual mapped reads were quantified to calculate fragments per kilobase of exon per million mapped fragments (FPKM) values. Genes with FPKM > 1 in at least one sample were kept for downstream analysis, including Pearson correlation, principal component analysis (PCA), and cluster analysis. Differential gene expression analysis was performed using the DESeq2 R package (1.20.0)20. Genes were considered differentially expressed if they had a false discovery rate (FDR) of P < 0.05. Ingenuity pathway analysis (IPA) was used to reveal the Gene Ontology (GO) and pathways21. Stage-specific genes were identified as previously described, and expression pattern clusters were generated by unsupervised hierarchical clustering analysis using R. The cluster Profiler R package22 was used to implement GO enrichment analysis of differentially expressed genes, correcting for gene length bias. GO terms with a corrected P < 0.05 were considered significantly enriched by differentially expressed genes. The cluster Profiler R package was also used to analyze the statistical enrichment of differentially expressed genes in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Genome-wide sulfite sequencing

We extracted genomic DNA with the Magnetic Universal Genomic DNA Kit (Tiangen, Beijing, China) according to manufacturer’s protocols and checked its quality using agarose gel electrophoresis and Nanophotometer® spectrophotometer (IMPLEN, CA, USA). Then, we performed library construction and sequencing as described before23. In brief, after mixing 5.2 μg of genomic DNA with 26 ng of λDNA, we broke it into 200–300 bp fragments by ultrasonication. Then, we performed end-repair on the sheared DNA fragments, added an A-tail, and ligated sequencing adapters that had been methylated at all cytosines. After that, we treated the DNA fragments twice with the EZ DNA Methylation-GoldTM Kit (Zymo Research) to convert them to disulfite and amplified the products by PCR. After passing the quality control, we then performed paired-end sequencing of the DNA library on Illumina HiSeq 2500 sequencing instrument at Novogene (Beijing, China).

Methylation analysis

We used FastQC software24 (fastqc-v0.11.5) to perform basic quality assessment on the original reads. Next, we used fastp software (fastp 0.20.0)25 to process the original reads by filtering out reads containing adapter sequences, reads with a percentage of N ≥ 10%, and low-quality reads (HRED score ≤ 5, low-quality bases ≥ 50%.). We aligned the clean reads to the goat genome (Capra hircus, ARS1 assembly) previously mentioned using Bismark software (version 0.12.5)26, predicted repeat sequences with RepeatMasker (v4.0.6)27, and tracked a genome Center for Genome Imaging (CGI) with cpgIslandExt. Then, we identified differentially methylated regions (DMRs) using swDM software (http://122.228.158.106/swDMR/), with a window length set at 1000 bp and a step size of 100 bp. Genes overlapping with DMRs in gene body regions (transcription start site (TSS) to transcription end site (TES) or promoter regions (2 kb upstream of TSS) were defined as DMR-related genes. We performed GO enrichment analysis of DMR-related genes using the GOseq R package28 with gene length bias correction and identified KEGG pathway enrichment using KOBAS 3.029. A corrected P < 0.05 was considered significant for GO term and pathway enrichment analysis.

Integrated analysis screening for differentially methylated and expressed genes

For the integrated analysis of gene expression and DNA methylation, we used a hierarchical testing approach30 to identify DNA methylation sites surrounding differentially expressed genes. First, we identified all CpG sites located 5 kb upstream and downstream of the transcription start site of methylation-regulated genes. Second, we combined the methylation positions within a 200 bp distance to form methylation regions. Regions that contained less than 20% CpGs (BH-corrected p value < 0.05) were excluded, and the remaining regions were considered DMRs. We performed FDR correction on all CpGs within the retained regions (BH-corrected p < 0.05). Primer sequences are given in Additional file 1.

Results

Transcriptional changes in the goat pineal gland during the onset of puberty

To reveal the molecular mechanism underlying the initiation of goats, we extracted RNA and DNA from the pineal glands of prepubescent and pubescent WWGs. We obtained approximately 72.96 million clean reads and 68.84 million total maps from each sample on average, with Q30 of clean, full-length reads ranging from 93.60% to 93.98% (Additional file 2). We found a total of 21,322 genes expressed in the pineal gland during prepuberty and puberty, which could be used for further downstream analysis. Details of the sequencing data quality are shown in Additional file 3. Pearson correlation coefficients of biological replicates showed high correlations in both stages (Fig. 2a), suggesting the high reproducibility of sample preparation and sequencing data quality. PCA showed significant differences between groups and high repeatability of samples (Fig. 2b).

Fig. 2
figure 2

Transcriptional changes of goats’ pineal gland during the onset of puberty. (a) To detect relatedness between sample individuals, Pearson’s correlation coefficient (r2) was used as theevaluation index. The closer the number is to 1, the stronger the correlation. (b) Principal component (PC) analysis of the transcriptomes for two stages of pineal glands. Pineal glands from the same stage are shown in the same color. PC1 (38.84%) and PC2 (19.64%) represent the top two dimensions of detected genes among the prepubertal and pubertal pineal glands.

Full size image

Differentially expressed genes during goats in puberty

By comparing RNA-seq data between prepuberty and puberty, we found a total of 2,703 differentially expressed genes (fold change > 0 and p-adjusted value < 0.05). Of these, 1303 genes were downregulated and 1400 genes were upregulated in puberty compared with prepuberty (Fig. 3a). Cluster analysis of differentially expressed gene (DEG) sets showed that the red part represented high mRNA expression and the green part represented low mRNA expression. We found that the prepuberty and puberty data were clustered together because of the similar gene expression patterns, and the repeated expression patterns in the group were consistent (Fig. 3b). Table 1 shows the 10 upregulated genes and 10 downregulated genes with the largest differences. According to KEGG analysis, the upregulated genes are involved in signaling pathways including oxidative phosphorylation (OXPHOS), Parkinson’s disease, ribosomes, Huntington’s disease and thermogenesis (Fig. 3c);meanwhile, the downregulated genes are involved in signaling pathways related to axon guidance, growth hormone synthesis, secretion and action, thyroid hormone, prolactin, and stem cell pluripotency (Fig. 3c). A total of 313 signaling pathways were involved in the prepuberty and puberty DEGs, and the 20 signaling pathways with the most significant differences between prepuberty and puberty are shown in Fig. 3d. Based on the analysis, we found KEGG entries significantly enriched in the following signaling pathways: OXPHOS, thermogenesis, Huntington’s disease, Parkinson’s disease, Alzheimer’s disease, ribosomes, nonalcoholic fatty liver disease (NAFLD) and cardiac muscle contraction. In summary, these data provide a comprehensive transcriptome profile of the pineal gland during goat puberty (Table 2).

Fig. 3
figure 3

DEGs during goat development. (a) Volcano plot of the DEGs between prepuberty and puberty. The red spots indicate genes that were upregulated in puberty, the green spots indicate downregulated genes, and the blue spots indicate genes that were not differentially expressed. (b) Hierarchical clustering of DEGs of prepubertal and pubertal goats’ pineal glands. The color spectrum, ranging from red to green, indicates the normalization level of gene expression. (c) Significant functional pathways revealed from up- or downregulated genes in comparisons between prepubertal and pubertal goats’ pineal glands. (d) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs between prepuberty and puberty.

Full size image
Table 1 List of the top 10 up- and downregulated genes in the pineal glands of pubertal goats.
Full size table
Table 2 List of the predicted new transcripts.
Full size table

Developmental transition is defined by stage-specific gene expression profiles

To investigate the differences in transcriptional activities during puberty, we identified genes that were transcribed only in a specific stage of development. We found that among them, 409 genes were specifically expressed in prepuberty, and 280 genes were specifically expressed in puberty (Fig. 4a), suggesting that these genes may play important roles in specific stages of development. Through GO functional enrichment analysis of 409 genes that are specifically expressed in prepuberty, axons, condensed chromosome, centromeric region, synapse, neuronal projection and other cellular components were significantly enriched (Fig. 4b). Additionally, KEGG analysis showed that these genes were significantly enriched in signaling pathways related to Huntington disease, glycosphingolipid biosynthesis-lacto and neolacto series, spinocerebellar ataxia, synaptic vesicle cycle, insulin secretion, and parathyroid hormone synthesis and secretion (Fig. 4c).

Fig. 4
figure 4

Enrichment analysis of DEGs during goat development. (a) Venn chart of differentially expressed genes between prepuberty and puberty. (b) Gene Ontology of prepubertal-specifically expressed genes. (c) Scatter plot of prepubertal-specifically expressed genes with significant enrichment. (d) Gene Ontology of pubertal-specifically expressed genes. (e) Scatter plot of pubertal-specifically expressed genes with significant enrichment.

Full size image

Furthermore, GO enrichment analysis of the 280 genes specifically expressed in puberty showed that the biological processes mainly involved the positive regulation of T-cell proliferation and activation, interleukin-2 biosynthesis, and granulocyte migration. The cell components were significantly enriched in cell surface, basal plasma membrane, and basal part of cells, while the molecular functions mainly included fibroblast growth factor binding, neuropeptide receptor activity, integrin binding, and hydrolase activity (Fig. 4d). Based on the results of the KEGG analysis, these genes were found to be significantly enriched in leishmaniasis, phagosome, sphingolipid metabolism and signaling pathways regulating stem cells (Fig. 4e). Further details can be found in Tables 3 and 4.

Table 3 List of gene body region overlap genes with significantly enriched.
Full size table
Table 4 List of promoter region overlap genes with significantly enriched.
Full size table

Genome-wide DNA methylation analysis

We performed global DNA methylation analysis of the pineal gland using WGBS with 30 × genome coverage and > 99% conversion efficiency. On average, 73.66 and 71.34 Giga raw bases were generated for prepuberty and puberty pineal glands, respectively. After filtering out low-quality data, we generated approximately 270 million clean reads for each group, with the Q30 of clean, full-length reads ranging from 88.88 to 89.33%. The mapped reads with rates ranging from 88.33 to 89.87% were used for subsequent analysis. The details of the quality of the sequencing data are shown in Table 5.

Table 5 Sequencing data by whole genome bisulfite sequencing for goat prepuberty and puberty stages.
Full size table

All methylated genomic C sites were approximately 1.76% in the prepuberty group and 1.84% in the puberty group. The methylation levels of CG, CHH, and CHG (where H is A, C, or T) were different. In prepuberty samples, we found genome-wide methylated cytosine (mC) levels of 23.35% for CG, 0.22% for CHG, and 0.31% for CHH, and in puberty samples, 25.84% for CG, 0.11% for CHG, and 0.23% for CHH. The genome-wide methylated cytosine levels of CHG and CHH in puberty were significantly different from those in prepuberty (Fig. 5a).

Fig. 5
figure 5

DNA methylation levels across genomic elements in prepubertal and pubertal goats. (a) The average ratio of DNA methylation types in prepubertal and pubertal genomes of goats. The green, purple, and orange colors represent methylated mCG, mCHG, and mCHH, respectively. (b) Violin plot for the overall distribution of methylation levels for different methylation types: CG, CHG and CHH. H = A, C or T. The abscissa represents the different samples, the ordinate represents the level of methylation of the samples, the width of each violin represents the density of the point at that methylation level, and the boxplot shows the methylation levels in each violin. (c) Circos map of genome-wide differences in methylation levels and overall differences between prepuberty and puberty. The three circles from the outside to the inside represent the methylation level of puberty and the methylation level of prepuberty. Internal scales: DNA methylation level represents methylation level, and DNA methylation difference among samples represents heatmap. (d) DNA methylation levels across genomic elements in prepubertal and pubertal goats. The orange and green colors represent prepuberty and puberty, respectively. The abscissa represents different genomic elements: CGI, CGI, shore, promoter, utr5, exon, intron, utr3, and repeat. The left ordinate represents the methylation levels of CG contexts, and the right ordinate represents the methylation levels of the CHG/CHH context.

Full size image

We drew a violin graph with points representing different levels of methylation. The most of the methylation found were CG and not CHG and CHH (Fig. 5b). We then plotted chromosome methylation maps for prepuberty and puberty samples, which showed that most chromosomal cytosine hypermethylation was found in the CG context and that the chromosomal mC sites were different between the prepuberty and puberty stages (Fig. 5c).

To further compare the genome-wide distribution and methylation levels of various functional genomic elements between the two developmental stages, we analyzed the methylation status of six different regions, including the promoters, 5′UTRs (untranslated regions), exons, introns, 3′UTRs, and distal intergenic regions. No significant differences were observed among the different genetic elements for the three mC contexts. Overall, the methylation levels in the CG context were higher than those in the CHG and CHH contexts, where the CHG and CHH contexts were hypomethylated and stable in all functional elements. The DNA methylation levels in the CG context were highest in introns, followed by exons (except the first exon) and downstream regions, and the lowest level of DNA methylation was at the CGI site. The methylation levels gradually decreased from promoters to TSSs and increased from TSSs to introns. More detailed information is listed in Fig. 5d.

Characterization of DMRs

We identified 206,532 methylated regions and 1294 differentially methylated regions (DMRs), including 570 differentially methylated CG regions, 371 differentially methylated CHG regions, and 353 differentially methylated CHH regions. Among the DMRs, 986 (CG:359+CHG:310+CHH:317) were hypermethylated, and 308 (CG:211+CHG:61+CHH:36) were hypomethylated. The DMRs were mostly located in intergenic regions, followed by promoters, but were rarely located at TSSs and TESs (Fig. 6a). Additionally, we analyzed genome-wide methylation in goats at the prepuberty and puberty stages using hierarchical clustering, which revealed a clear separation between the two developmental stages (Fig. 6b). More detailed DMR results are listed in Additional file 4.

Fig. 6
figure 6

Enrichment analysis of CG-type differentially methylated genes. (a) The number of differentially methylated regions between different methylation types. Histograms show the distribution numbers of differentially methylated regions (DMRs) in different genomic elements in the CG, CHG, and CHH contexts. (b) Heatmap cluster analysis of DMRs between prepuberty and puberty. In the heatmap, highly methylated loci are displayed in red, and sparsely methylated loci are displayed in blue. (c) Counts of differentially methylated genes (DMGs) enriched for the top 30 Gene Ontology (GO) terms. The abscissa represents the number of DMGs, and the ordinate shows the GO pathway terms. (d) The scatterplot of 20 KEGG pathways. The ordinate represents the enriched pathways, and the abscissa represents the Rich factor of the corresponding pathways; the size of the spots represents the number of genes related to DMRs enriched in each pathway, while the color of the spot represents the corrected P value for each pathway. The Rich factors indicate the ratio of the number of DMGs mapped to a certain pathway to the total number of genes mapped to this pathway. A greater Rich factor means greater enrichment.

Full size image

To explore changes in the methylation status of the pineal gland during puberty initiation, we used the GO and KEGG databases to annotate the 1294 differentially methylated genes (DMGs) detected in the DMRs. As most of the DMGs were related to CG (CG:44.0%, CHG:28.7%, CHH:27.3%), we focused on CG methylation for the DMG functional enrichment analysis. Based on the GO database, we identified terms that play a key role in the regulation of growth, regulation of cell communication, and regulation of signaling (corrected P < 0.05). We used enrichment terms to preliminarily determine the signal transduction pathways mainly involved in genes in DMRs. The 30 most significantly differentially enriched onset-of-puberty-related GO terms for DMGs between prepuberty and puberty pineal gland samples are shown in Fig. 6c. According to the KEGG pathway analysis, DMGs were significantly enriched in the Alzheimer’s disease, Huntington’s disease, OXPHOS, Parkinson’s disease, and spliceosome signaling pathways. The 20 puberty-development-related KEGG terms with a corrected p value < 0.05 for the DMGs are listed in Fig. 6d. The results suggest that these DMGs, which are influenced by DNA methylation, can affect puberty onset. More detailed results of the GO and KEGG analyses of CG, CHG, and CHH methylation are shown in Additional files 5–10.

Screening of potential functional differentially methylated genes involved in puberty onset

To identify key genes involved in regulating puberty onset, we first screened 2703 differentially expressed genes using RNA-seq data and then screened 1914 candidate genes known to be associated with puberty initiation through GO and KEGG functional enrichment analysis of significant DMGs. Second, using our previous WGBS data between prepuberty and puberty, we further screened 129 overlapping genes between DMGs and DEGs. By means of GO analysis, it was found that it was mainly involved in the negative regulation of neuronal projection regeneration, regulation of neuronal differentiation, phospholipid metabolic processes, membrane protein complexes and membrane parts (Fig. 7a). According to the results of the KEGG enrichment analysis, the majority of the genes were enriched in metabolic pathways, followed by the Parkinson’s disease, OXPHOS, Huntington’s disease and Alzheimer’s disease signaling pathways (Fig. 7b). An interaction network was generated using STRING software, and a total of 395 candidate genes that interacted with each other were identified. As shown in Fig. 7c, ATP Synthase F1 Subunit Delta (ATP5F1D), Calcium voltage-gated channel auxiliary subunit beta 2 (CACNB2), and Phosphatase and Tensin Homolog (PTEN) were identified as hub genes in the interaction network related to puberty initiation. More detailed results on these genes are listed in Additional file 11. We also used RT‒qPCR to detect the transcription levels of some transcriptome differential genes and overlapping genes and found that the variation trend of these genes was consistent with the sequencing results (see Additional file 12).

Fig. 7
figure 7

Analysis of overlapping genes between DMGs and DEGs. (a) GO enrichment histogram of overlapping DMRs and DEGS genes in the promoter region. (b) GO enrichment histogram of overlapping DMRs and DEGS genes in the promoter region. (c) Construction of the network of DMGs related to puberty initiation. Analysis of the interaction between DMGs related to puberty initiation using STRING software according to the interplay index (confidence > 0.9). (d) Scatter diagram of DMR differential methylation and differential expression gene levels.

Full size image

To investigate the effect of DNA methylation on gene expression levels, we compared the trend between gene expression and methylation levels using the FPKM value for the RNA-seq data and the difference in methylation levels between prepuberty and puberty WGBS data samples. The results showed that regardless of whether the transcription levels of these genes were increased or decreased or whether the methylation sequences of these genes were CG, CHG, or CHH, DEGs had increased methylation levels in puberty compared with prepuberty (Fig. 7d).

Discussion

DNA methylation is an epigenetic regulatory form that plays a crucial role in gene expression and tissue development31. Previous studies have analyzed DNA methylation in the hypothalamus, pituitary gland, and ovary during puberty. However, whole-genome methylation analysis of the pineal gland has not been reported before. This study is the first to systematically compare the genome-wide DNA methylation profile and transcriptome changes of the prepubertal and pubertal goat pineal gland. In our study, we analyzed the differential expression of genes and pathways related to puberty in the goat pineal gland between the pre and pub stages by transcriptome analysis. We also studied the whole-genome DNA methylation spectrum of goat pineal glands by WGBS to clarify the relationship between transcription and DNA methylation at different developmental stages.

DNA methylation is mainly caused by the DNA Methyltransferase (DNMT) family, including DNA Methyltransferase 1(DNMT1), DNA methyltransferase 3A (DNMT3A), and DNA Methyltransferase 3 Beta (DNMT3B). DNMT1 mainly maintains DNA methylation, while DNMT3A and DNMT3B catalyze de novo methylation. In our study, we found that the expression levels of DNMT1, DNMT3A, and DNMT3B were higher in the pre stage than in the pub stage, indicating that DNMTs may play an important role in regulating the transcription of genes related to puberty in goats. In the functional genomic region, the whole-genome methylation patterns of the pre and pub stages are similar. However, there are differences among the three mC contexts that may be related to the differences in the sequences of various genetic elements5. After filtering the WGBS raw data, we obtained 270 million clean reads per sample. Approximately 1.8% of the cytosine sites were methylated, and CG had the highest methylation rate, which was similar to the results found in other species and tissues. Among the gene function regions, the methylation level of TSSs was the lowest, which was consistent with the results found by WGBS in sheep ovaries. In our study, we identified 30,000 DMRs and 22,161 genes associated with these DMRs. The proportion of DMRs in the promoter region, 3′UTR and 5′UTR was small, with the majority being concentrated in the distal intergenic region.

We attempted to reveal the roles of DMGs through functional annotations. According to the GO analysis of the DMGs in the promoter region, these are mainly terms related to the negative regulation of neuronal projection regeneration, regulation of neuronal differentiation, phospholipid metabolism processes, membrane protein complexes, membrane parts and phosphatidylinositol. Furthermore, KEGG analysis confirmed that several pathways of methylation differences between the two developmental stages may be related to the onset of puberty, with the largest number of genes enriched in metabolic pathways, followed by the Parkinson’s disease, OXPHOS, Huntington’s disease, and Alzheimer’s disease signaling pathways. OXPHOS is a metabolic pathway that uses the energy released by the structure, enzymes, and nutrients of intracellular mitochondria to transform ATP32, a molecule that provides energy for metabolism. In reproductive tissues and gametes, OXPHOS is crucial for early gametogenesis and normal embryonic development. Moreover, OXPHOS is essential for brain health, and mitochondrial damage and oxidation are part of a common neuropathological mechanism affecting neurological diseases, including Alzheimer’s and Parkinson’s diseases33.

The genes related to puberty in the hypothalamus of Jining gray goats and Liaoning cashmere goats are mainly enriched in metabolic pathways, Parkinson’s disease, and OXPHOS, which is consistent with the results of this study34. Meanwhile, metabolic factors strictly control the synthesis and release of kisspeptin, so sufficient energy storage is critical to successful reproduction35. Therefore, the metabolic pathway is also crucial to the onset of puberty. Methylation-regulated DMGs enriched in several signaling pathways, such as ATP5F1D, CACNB2 and PTEN, may play an important role in the onset of puberty. The ATP5F1D gene encodes ATP synthase, H+ transport, mitochondrial F1 complex, and mitochondrial ATP synthase d subunit, which play an essential role in proton-coupled electron transfer and ATP production. Patients with ATP5F1Dhomozygous missense mutations exhibit intermittent sleepiness, metabolic acidosis, 3-methylglutamic aciduria, and hyperammonemia36. CACNB2 encodes a subunit of voltage-gated calcium channel proteins. Studies have found genetic associations between members of the voltage-gated calcium channel supergene family and autism spectrum disorders37. An analysis of genome-wide methylation profiles in patients with schizophrenia found that CACNB2 was one of the most significant genes38,39. PTEN is considered a tumor suppressor, encoding phosphatidylinositol-3,4,5-triphosphate 3-phosphatase. Unlike most protein tyrosine phosphatases, the protein preferentially dephosphorylates inositol phosphate substrate. The gene has a high frequency of mutations in a large number of cancers and plays a tumor suppressor role by negatively regulating intracellular phosphatidylinositol-3,4,5-triphosphate levels and the AKT/PKB signaling pathway40. Studies have shown that the increase in hypothalamic basal pAKT in mice and the decrease in leptin-induced Phosphorylated transducer and activator of transcription-3 (PSTAT3) and PTEN levels inhibit the increase in PI3K sensitivity to leptin in metabolic regulation41. The pineal gland is an important neuroendocrine organ, and our research enriches metabolic pathways such as metabolism, auxin secretion, and signal transduction. Bioinformatics results show that DMGs related to these regulatory processes are significantly different between prepuberty and puberty, suggesting that they may play an important role in the onset of puberty.

However, little is known about how DNA methylation affects gene expression and how these genes work together. We generated DMG interaction networks to determine whether DMGs play a determinative role in the onset of puberty in goats. Network analysis shows that CACNB2, PTEN, Lin-28 Homolog B (LIN28B), Gastric Inhibitory Polypeptide (GIP), Opsin 1, Short Wave Sensitive (OPN1SW), and DCC Netrin 1 Receptor (DCC) are key nodes. This suggests that these genes might play important roles in puberty initiation, and the differentiation and regulation of these genes through DNA methylation might be one of the mechanisms that determine the difference in puberty initiation time and reproductive capacity of goats.

Many DEGs between prepubescent and pubertal goats were found in this study to be related to reproductive regulation or puberty. Numerous studies have shown that the Lin28B gene is involved in the onset of mammalian puberty42,43. A genome-wide association study (GWAS) confirmed that Lin28B mutation affects adult height and menarche age44. The Lin28B gene is widely expressed in various tissues of mouse fetuses, but its expression in adult animals is limited to the HPG axis and other tissues, and the expression level decreases45. Long-term dark environments increased the expression of Lin28B in the hypothalamus of 10- to 15-day-old rats, resulting in the delay of the opening of the vulva in female rats and the separation of the foreskin in male rats46. This may be related to the different levels of melatonin secretion in the pineal gland under different photoperiods. Studies have found that the Lin28B allele frequency of goat precocious puberty is significantly higher than that of delayed puberty47, and compared with juveniles, the Lin28B mRNA of goat hypothalamus in puberty is significantly increased42. This study also showed that the Lin28B gene in the goat pineal gland was significantly upregulated in puberty compared with prepuberty. However, the reproductive regulation of Lin28B in mammals such as goats and mice are opposite, which may be due to the different effects of photoperiod on the reproduction of different animals. Therefore, it is still necessary to further study the mechanism by which the Lin28B gene regulates puberty initiation.

In 2018, several rare pathogenic mutations of the DCC gene encoding theFibronectin type III ___domain (FN3) ___domain protein and its ligand Netrin-1 (NTN1) in patients with idiopathic hypogonadotropic hypogonadism were first confirmed in research48. Mouse genetics studies have also shown that NTN1 and DCC are necessary for the migration of GnRH neurons49, and the netrin signal mediated by DCC is involved in the maturation of the cerebral cortex dopamine circuit during puberty50. It was also found that the intracerebroventricular injection of Nesfatin-1 could advance the opening time of the vulva in rats, promote the onset of puberty and increase the transcription levels of GnRH, GnRHR, LHβ, LHR, FSHβ, and ERα genes on the reproductive axis. Combined with the results of these studies, it is suggested that DCC may affect the secretion of GnRH and the expression of related reproductive genes on the hypothalamic-pituitary-ovary axis by promoting the development of GnRH neurons, thereby affecting the onset of puberty.

In summary, our research findings demonstrate that DNA methylation of the pineal gland can affect the onset of puberty in goats. The relationship between gene DNA methylation and gene expression level is complex, and although DNA methylation in the promoter region can inhibit gene transcription, gene silencing increases with the increase in CpG dinucleotide density at the promoter 51. However, it is still unclear how DNA methylation in the gene affects gene expression. We propose that OPN1SW, DCC, Lin28B, and GIP are key regulatory genes in the onset of puberty and that their DNA methylation status may play a functional role in puberty initiation. Methylation of these genes may lead to significant differences in puberty initiation and pineal gland neuron development between prepubescent and pubertal goats. However, the epigenetic mechanisms of gene regulation and the genetic regions involved in the onset of puberty and neuronal development require further study.

Conclusion

In summary, we have provided a systematic description of the whole-genome DNA methylation patterns and transcriptomes of goat pineal gland during prepuberty and puberty, which is the first of its kind. We have identified several novel and important DMRs, DMGs, and pathways related to goat puberty. These findings provide valuable insights into the genetic and epigenetic mechanisms of economic traits in goats and can be used in marker-assisted selection procedures to promote goat breeding.