Abstract

Background

   The hypothalamus is a critical organ that regulates sexual development
   in animals. However, current research on the hypothalamic regulation of
   sexual maturation in female goats remains limited. In this study, we
   conducted metabolomic and transcriptomic analyses on the hypothalamic
   tissues of female Jining grey goats at different stages of sexual
   development (1 day old (neonatal, D1, n = 5), 2 months old (prepuberty,
   M2, n = 5), 4 months old (sexual maturity, M4, n = 5), and 6 months old
   (breeding period, M6, n = 5)).

Results

   A total of 418 differential metabolites (DAMs) were identified in this
   study, among which the abundance of metabolites such as anserine,
   L-histidine, carnosine, taurine, and 4-aminobutyric gradually increased
   with the progression of sexual development. These metabolites may
   regulate neuronal development and hormone secretion processes by
   influencing the metabolism of histidine and phenylalanine. Through
   combined transcriptomic and metabolomic analyses, we identified that
   differentially expressed genes such as mitogen-activated protein kinase
   kinase kinase 9 (MAP3K9), prune homolog 2 with BCH domain (PRUNE2), and
   potassium voltage-gated channel interacting protein 4(KCNIP4) may
   jointly regulate the development and energy metabolism of hypothalamic
   Gonadotropin-releasing hormone neurons in conjunction with DAMs,
   including LPC22:5, 2-Arachidonyl Glycerol ether, LPE22:5, and
   Lysops22:5. Additionally, we elucidated the molecular mechanism through
   which glutathione metabolism regulates sexual maturation in goats.

Conclusions

   In summary, this study illustrates the dynamic changes in metabolites
   and mRNA within hypothalamic tissue during postnatal sexual maturation
   in female Jining grey goats. This research may provide significant
   scientific insights for future animal breeding.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12864-025-11492-2.

   Keywords: Precocious goat, Multi-omics analysis, Hormone secretion,
   Glutathione metabolism

Introduction

   Sexual maturation is a dynamic developmental process that is regulated
   by a variety of factors such as genetics, endocrine, and nutrition, and
   is primarily achieved through the hypothalamic-pituitary-gonadal (HPG)
   axis [[34]1, [35]2]. Gonadotropin-releasing hormone (GnRH) stimulates
   the pituitary gland to synthesize and secrete luteinizing hormone (LH)
   and follicle-stimulating hormone (FSH) [[36]3, [37]4], activates the
   hypothalamic-pituitary signaling cascade, promotes ovarian development,
   and promotes gonadotropin-induced gene expression [[38]5, [39]6,
   [40]7], while steroids secreted by the gonads provide feedback
   regulation of GnRH secretion [[41]8, [42]9]. This complex
   transcriptional network that regulates reproduction has been
   extensively studied, with the highest levels of network regulation
   controlled by key genes and metabolic pathways, such as genes and
   effectors associated with Kisspeptin, GABAergic, and glutamatergic
   neurons and glial cells, which act as upstream regulators of GnRH
   neurosecretion [[43]10]. In cases of negative energy balance, such as
   energy/nutrient deficiencies due to fasting, the onset of puberty and
   the acquisition of fertility are often delayed or inhibited [[44]11].
   In addition, Kisspeptin neurons play a major role in energy
   expenditure, which may be mediated by Kisspeptin and/or its
   cotransmitters such as glutamate [[45]12].

   Mutations in the kisspeptin (KISS1) or its receptor (KISS1R) have been
   shown to lead to decreased hypothalamic GnRH secretion and inability to
   enter puberty [[46]13, [47]14]. GABAA receptor activation, on the other
   hand, produces excitatory responses primarily to embryonic and
   prepubertal GnRH neurons [[48]15]. In addition to the regulation of
   gene expression, the process of sexual development is also tightly
   regulated by metabolic factors [[49]16], such as lipid metabolism
   [[50]17], amino acid metabolism [[51]18, [52]19], and others. Compared
   with young mice, there were significant reductions in multiple
   metabolites and neurotransmitters in the hypothalamus of older mice,
   and there was an increase in the levels of genes associated with the
   production and utilization of monocarboxylates [[53]20].
   3-Hydroxybenzoic acid and proline glutamine have been reported to be
   biomarkers of precocious puberty in women [[54]21], while ceramides
   affect ovarian steroid hormone production and the synthesis of central
   nervousamides promotes precocious puberty in women [[55]22].
   Metabolomics and transcriptomics techniques have been widely used to
   study the regulatory mechanisms of mammalian reproduction [[56]23].

   As the main source of meat, milk, and fur, goats play a vital role in
   animal husbandry [[57]24]. However, the low fecundity of goats poses a
   major challenge to the development of the goat industry. Precocious
   puberty can accelerate the intergenerational succession of livestock,
   reduce the cost of feeding, and effectively improve the fertility of
   goats [[58]25]. Understanding the physiological and molecular pathways
   involved in the sexual maturation of goat is crucial for genetic
   improvement and breeding practice. Compared with other goat breeds, The
   Jining grey goat is a nationally recognized indigenous breed in China,
   renowned for its exceptionally early sexual maturity and high
   reproductive efficiency. This breed initiates puberty as early as 2–3
   months postpartum and achieves full sexual maturity by 4 months of age
   [[59]26], which is an ideal model for studying the regulation mechanism
   of reproduction [[60]27, [61]28].

   Although the study of hypothalamic tissue has been abundant, the
   specific mechanism of its initiation of complex neuroendocrine
   signaling cascade, regulation of neuronal development and GnRH
   secretion to promote sexual development has not been fully elucidated.
   This study hypothesizes that during the development of goats after
   birth, the transcriptomic and metabolomic expression profiles of the
   hypothalamus undergo significant changes, which drive sexual maturation
   by regulating GnRH secretion and related neuroendocrine signals. In
   this study, metabolomics and transcriptomics were combined to analyze
   the dynamic changes of hypothalamic metabolites and mRNA at different
   stages after birth (1 day old (neonatal stage, D1, n = 5), 2 months old
   (prepuberty, M2, n = 5), 4 months old (sexual maturity stage, M4,
   n = 5) and 6 months old (reproductive stage, n = 5). This study
   revealed the transcriptional regulation mechanism and metabolic
   characteristics of the hypothalamus during goat sexual development, and
   provided important insights for genetic improvement and breed
   selection.

Materials and methods

Animals

   In this study, the hypothalamic tissues of 20 female Jining grey goats
   were collected from the original breeding farm of Jining grey goats
   (Jiaxiang County, Shandong Province, China). The goats were categorized
   into the following age groups:1 day old (neonatal stage, D1, n = 5), 2
   months old (prepuberty, M2, n = 5),4 months old (sexual maturity, M4,
   n = 5), and 6 months old (breeding period, M6, n = 5). All goats were
   raised under the same environmental conditions. 1-day-old and
   2-month-old kid goats lived with their mother goats; 4-month-old and
   6-month-old goats were separated from their mother goats and housed in
   the same shed but in different pens. The conditions of water source,
   lighting, temperature, and feed in the shed were completely consistent.
   All selected goats are healthy, disease-free, and are able to eat and
   drink freely. Starting from 4 weeks of age, the goats were weaned from
   milk to solid feed, and fully weaned at 10 weeks of age. All
   experimental goats were slaughtered at the same time, with the D1 group
   goats having already ingested colostrum. The selected goats were either
   in anestrus or in the diestrus and had similar weights in the same
   group. Experimental goats were slaughtered on the same day after
   electroshock anesthesia. According to the method of Sesti and Britt
   [[62]29], the hypothalamus was divided into the preoptic area (POA),
   mediobasal hypothalamus (MBH), and stalk median eminence (SME). The POA
   spanned ∼ 5 mm anterior to the optic chiasm, dorsally bounded by the
   MBH’s anterior edge. The MBH extended from the optic chiasm’s posterior
   edge to the mammillary bodies’ anterior edge, laterally limited by the
   hypothalamic sulcus and dorsally by a 5 mm posterior depth. The SME was
   ventrally dissected from the MBH along the third ventricle base
   midline. After sampling, the tissues were pooled, flash-frozen in
   liquid nitrogen, and stored at -80 °C for metabolomic and transcriptome
   analysis. To improve the accuracy of metabolome analysis, reduce
   pretreatment and instrument errors, and better perform subsequent
   combined analysis of transcriptome and metabolome data, we divided each
   biological replicate into two technical replicates [[63]30].

Extraction, detection and quantification of hypothalamic metabolites

   One hundred milligrams of hypothalamic tissue were ground in liquid
   nitrogen and transferred to an Eppendorf tube, where it was combined
   with 500 μL of an 80% methanol aqueous solution. The mixture was
   vortexed and then kept in an ice bath for 5 min before being
   centrifuged at 15,000 g for 20 min at 4 °C. The supernatant was
   collected and diluted with mass spectrometry-grade water to achieve a
   methanol concentration of 53%. The prepared supernatant was centrifuged
   again under the same conditions [[64]31].

   A UHPLC-MS/MS analysis was performed using a Vanquish UHPLC system
   coupled with an Orbitrap Q Exactive™ HF mass spectrometer (Thermo
   Fisher, Germany) [[65]32]. Sample separation occurred on a Hypersil
   Gold column (100 × 2.1 mm, 1.9 μm), applying a linear gradient for
   17 min at a flow rate of 0.2 ml/min. In positive polarity mode, solvent
   A was 0.1% formic acid in water and solvent B was methanol; in negative
   polarity mode, solvent A consisted of 5 mM ammonium acetate (pH 9.0)
   with methanol as B. The gradient protocol included holding at 2% B for
   1.5 min, increasing to 85% B over 3 min, ramping to 100% B in 10 min,
   then returning to 2% B at 10.1 min for an additional 12 min. The mass
   spectrometer was set at a spray voltage of 3.5 kV and a capillary
   temperature of 320 °C. For quality control, all samples were mixed in
   equal volumes to create mixed QC samples.

Data processing and metabolite quantification

   The UHPLC-MS/MS raw data files were processed using Compound Discoverer
   3.1 (CD3.1, ThermoFisher) for peak alignment, extraction, and
   metabolite quantification. Parameters included a retention time
   tolerance of 0.2 min, mass tolerance of 5 ppm, 30% signal strength
   tolerance, a signal-to-noise ratio of 3, and a minimum intensity
   threshold [[66]33]. Peak intensities were normalized to the total
   spectral intensity, enabling the prediction of molecular formulas from
   adducts, molecular ion peaks, and fragment ions. Finally, identified
   peaks were cross-referenced with mzCloud
   ([67]https://www.mzcloud.org/), mzVault, and MassList databases for
   accurate qualitative and relative quantitative analyses [[68]34].

Metabolome data analysis

   We conducted metabolite annotation for the samples and utilized the
   KEGG ([69]https://www.genome.jp/kegg/pathway.html), HMDB
   ([70]https://hmdb.ca/metabolites), and LIPIDMaps
   ([71]http://www.lipidmaps.org/) databases to gather further metabolic
   information. To examine the overall metabolic differences among the
   groups, we employed principal component analysis (PCA) and orthogonal
   partial least squares discriminant analysis (OPLS-DA) [[72]35].
   Statistical significance was determined using univariate analysis
   (t-test), and p-values were adjusted using the Benjamini-Hochberg
   correction method [[73]36]. The thresholds for differentially expressed
   metabolites (DAMs) were set as VIP > 1, false discovery rate
   (FDR) < 0.05, and fold change (FC) ≥ 1.5 (upregulated) or VIP > 1,
   FDR < 0.05, and FC ≤ 0.67 (downregulated). The identified differential
   metabolites were then clustered and visualized using a heatmap
   generated by the Pheatmap package in R (version 3.3.2). To investigate
   the expression patterns of metabolites, we conducted k-means clustering
   analysis using the mfuzz package [[74]37]. Metabolite annotations were
   performed using a dual matching strategy based on metabolite names and
   molecular formulas against the KEGG database (version 110.0, released
   on April 1, 2024). Data visualization was carried out using the R
   package ClusterGVis (V.0.1.1).

Transcriptome data analysis

   According to the manufacturer’s protocol, total RNA was extracted from
   hypothalamic tissues at four developmental stages using TRIzol reagent
   (Takara Bio, Dalian, China). The RNA integrity number (RIN) of the
   sequencing RNA was assessed using the Agilent Bioanalyzer 2100 system
   (Agilent Technologies, Santa Clara, CA, USA). Total RNA concentration
   was measured with the NanoDrop 2000 spectrophotometer (Thermo Fisher
   Scientific, Waltham, MA, USA), and the quality of total RNA was
   evaluated through agarose gel electrophoresis. RNA sequencing was
   conducted on the Illumina NovaSeq 6000 platform, with raw data filtered
   through FastQC (v0.11.8) to ensure high-quality reads. A reference
   genome index was created using HISAT2 (v2.0.5) [[75]38], and clean
   reads were aligned to the goat reference genome
   (GCF_001704415.2_ARS1.2). The data was assembled with StringTie
   (1.3.3b) [[76]39], and gene expression was normalized to Fragments Per
   Kilobase of transcript per Million mapped reads (FPKM) using
   featureCounts (v1.6.2). Differentially expressed genes (DEGs) were
   identified with DESeq2 (v1.22.1) using |log2(Fold Change)| ≥ 1 and
   p.adj < 0.05. K-means clustering analysis was performed with the mfuzz
   package [[77]37] in R, and KEGG enrichment analysis was conducted using
   clusterProfiler for DEGs exhibiting different expression patterns
   [[78]40].

Metabolome and transcriptome WGCNA analysis

   The co-expression network modules for all DEGs and DAMs were created
   using the WGCNA package (version 1.70-3) in R. The transcriptome and
   metabolome used soft thresholds of β = 7 and β = 6, respectively, with
   parameters minModuleSize = 30 and mergeCutHeight = 0.25 to build the
   topological network. Combining this data with the serum hormone indexes
   of the same batch of goats published in our laboratory in the past
   [[79]41], we evaluated the module-trait relationship. A significant
   correlation between modules and features was defined as p ≤ 0.05. KEGG
   enrichment analysis of DEGs and DAMs in key modules was performed using
   clusterProfiler. We filtered hub DEGs and hub DAMs in the module based
   on the top 20 of Module Eigenvector connectivity (KME). The correlation
   coefficient was then calculated using Spearman correlation analysis.

Integrative analysis of metabolomic and transcriptomic data

   Perform Spearman correlation analysis using R packages to analyze the
   hub DEGs and hub DAMs identified by WGCNA. Additionally, based on the
   expression pattern analysis and KEGG pathway analysis results in this
   study, compare those with similar or opposite expression patterns and
   their corresponding KEGG outcomes to identify shared KEGG pathways.
   This will determine the key signaling pathways jointly regulated by
   DAMs and DEGs.

Result

Metabolic profiling of hypothalamic tissue

   In order to better understand the dynamics of hypothalamic metabolites
   during the sexual maturation of goats, we performed a non-targeted
   metabolomic analysis of goat hypothalamic tissue using UPLC-MS/MS.
   Metabolomic data under the combination of positive and negative modes
   were pooled and analyzed. A total of 925 metabolites were identified
   from the hypothalamus (D1, M2, M4, M6) at different stages of sexual
   development after birth in goats (Figure [80]S1). Excluding unannotated
   metabolites, lipids and lipid-like molecules (292) and organic acids
   and derivatives (159) represented the most abundant metabolite
   categories. PCA analysis showed that the D1, M2, M4, and M6 groups had
   an obvious tendency for separation, and the M4 and M6 groups were close
   to each other. PC1 and PC2 together accounted for 41.32% of the total
   variation (Fig. [81]1a). The results of PLS-DA are shown in Figure
   [82]S2. The R^2Y ≥ 0.97 and Q2 ≥ 0.80 among the different comparison
   groups indicated that the model had good predictive power and a low
   risk of overfitting.

Fig. 1.

   [83]Fig. 1
   [84]Open in a new tab

   Metabolomic analysis of the hypothalamus during postnatal sexual
   development in goats. (a). PCA analysis of metabolomics data in D1, M2,
   M4 and M6 periods. (b). Histogram of up- and down-regulated
   differentially expressed metabolites (DAMs). (c). Venn diagram and
   cluster heatmap of DEGs identified in hypothalamic tissue. (d).
   Classification pie chart of DAMs

Identification and analysis of differential metabolites in the hypothalamus

   A total of 418 DAMs were identified according to the screening criteria
   (VIP > 1, FDR < 0.05, FC ≥ 1.5, or FC ≤ 0.667) (Table [85]S1-6).
   Specifically, pairwise comparisons identified 174 DAMs (80 upregulated,
   94 downregulated) between M2 and D1, 70 DAMs (24 upregulated, 46
   downregulated) between M4 and M2, 71 DAMs (28 upregulated, 43
   downregulated) between M6 and M4, 201 DAMs (70 upregulated, 131
   downregulated) between M4 and D1, 247 DAMs (76 upregulated, 171
   downregulated) between M6 and D1, and 128 DAMs (32 upregulated, 96
   downregulated) between M6 and M2 (Fig. [86]1b). The highest number of
   DAMs was identified in M6 vs. D1, suggesting significant differences in
   the physiological status of goats after birth and sexual maturity. In
   addition, the Venn diagram (Fig. [87]1c) illustrates common and unique
   DAMs across comparative groups. A sequential temporal analysis (M2 vs.
   D1, M4 vs. M2, M6 vs. M4) identified one common DAM, 17-AAG, which
   exhibited significantly higher abundance during the M2 stage.
   Comparative analysis against the D1 stage revealed 87 common DAMs in
   M2, M4, and M6, displaying two distinct expression patterns: higher
   abundance at postnatal D1 and increased abundance during sexual
   development (Fig. [88]1c).

   Further, the identified DAMs were classified (Fig. [89]1d), and after
   removing unannotated metabolites, lipids and lipid-like molecules
   accounted for the highest proportion (33.49%), followed by organic
   acids and derivatives (18.42%). The results of cluster analysis showed
   that lipids and lipid-like molecules were mainly composed of
   glycerophospholipids, fatty acids, steroids and their derivatives
   (Figure [90]S3). Among the glycerol phospholipid metabolites, the
   abundance of Lyso-phosphatidylcholine (LPC),
   Lyso-phosphatidylethanolamine (LPE), and Phosphatidylcholine (PC) was
   significantly higher at D1 stage than at other stages. The abundance of
   most fatty Acyls metabolites such as FAHFA was highly expressed at the
   M2 stage. At the same time, the abundance of several organic acid
   metabolites, such as L-Ornithine, betaine, L-serine, threonine, and
   valine, was also higher at the D1 stage (Figure [91]S4). In addition,
   the expression abundance of metabolites such as anserine, L-Histidine,
   carnosine, and gamma-glutamylmethionine was low in D1 stage, but
   remained high in the subsequent stages of sexual development. These
   results suggest that different metabolites may play different
   biological functions during sexual maturation in goats.

Dynamic changes and functional analysis of differential metabolites in
hypothalamic tissue during sexual development

   To investigate the temporal dynamics of hypothalamic metabolite
   profiles, we conducted K-means clustering and KEGG pathway enrichment
   analysis on the identified DAMs. The results showed that the 418 DAMs
   could be divided into 6 clusters according to their expression patterns
   (Fig. [92]2, Table [93]S7).

Fig. 2.

   [94]Fig. 2
   [95]Open in a new tab

   Heatmap and hierarchical clustering KEGG annotation diagram of DAMs at
   four stages of postnatal sexual development in Jining grey goats. The
   color indicates the expression value of the normalized adjusted DAMs.
   Red and blue indicate upward and downward adjustments, respectively

   Cluster 1 comprised 80 DAMs exhibiting a progressive decline in
   abundance with advancing sexual development. These metabolites were
   enriched in pathways associated with HIF-1 signaling pathway,
   neuroactive ligand-receptor interactions, circadian entrainment,
   tryptophan metabolism, and glutathione metabolism. In contrast, Cluster
   4 contained 58 DAMs that demonstrated an initial increase followed by a
   decline, peaking at the M2 stage. This cluster showed functional
   enrichment in histidine metabolism, beta-Alanine metabolism, and AMPK
   signaling pathway. Cluster 6 (76 DAMs) exhibited an inverse expression
   pattern relative to Cluster 1, with metabolite levels progressively
   increasing during sexual maturation. These DAMs were predominantly
   associated with GABAergic synapse, estrogen signaling pathway,
   glutathione metabolism, steroid hormone biosynthesis, cAMP signaling
   pathway, and histidine metabolism. These results suggest that these
   pathways may be involved in the sexual maturation of goats by
   regulating amino acid synthesis and metabolism, hormone secretion, and
   signaling.

WGCNA analysis identified key metabolites during postnatal sexual development
in goats

   To reduce data dimensionality and minimize interference from background
   noise, thereby facilitating the identification of key DAMs associated
   with phenotypic traits, we integrated previously acquired goat serum
   hormone data with the identified DAMs and conducted a WGCNA (Table
   [96]S8) [[97]41]. Using a soft threshold of β = 7, we identified a
   total of five modules (Figs. [98]3a and [99]S5a). Among these, the
   brown module exhibited significant positive correlations with GnRH
   (R = 0.56, P < 0.01), FSH (R = 0.33, P = 0.04), progesterone (PROG)
   (R = 0.44, P < 0.01), and estradiol (E2) (R = 0.37, P = 0.02)
   (Fig. [100]3b). In contrast, the turquoise module revealed significant
   negative correlations with GnRH (R = -0.45, P < 0.01), FSH (R = -0.58,
   P < 0.01), luteinizing hormone (LH) (R = -0.63, P < 0.01), PROG (R =
   -0.37, P = 0.02), and E2 (R = -0.57, P < 0.01).

Fig. 3.

   [101]Fig. 3
   [102]Open in a new tab

   Weighted gene co-expression network analysis (WGCNA) of DEGs identified
   at different developmental stages of goat hypothalamic tissue. (a).
   Hierarchical cluster tree of co-expression modules identified. (b).
   Heatmap of the correlation between module characteristic metabolites
   and serum hormone content (GnRH, FSH, LH, PROG2, E2) and developmental
   stage (D1, M2, M4, and M6). Red indicates a positive correlation
   between DAMs and serum indices, and blue indicates a negative
   correlation. (c). GS-MM scatterplot of the brown module with GnRH. (d).
   GS-MM scatterplot of the turquoise module with GnRH. (e). KEGG
   enrichment analysis scatterplot of DAMs in the brown module. (f). KEGG
   enrichment analysis scatterplot of DAMs in the urquoise module. (g).
   Expression heatmap of the top 20 DAMs of the brown module. (h).
   Expression heatmap of the top 20 DAMs of the turquoise module

   Based on the gene significance and module membership (GS-MM) results
   (Fig. [103]3c and f), we selected the turquoise and brown modules for
   further analysis. KEGG pathway enrichment analysis indicated that the
   DAMs in the brown module were primarily associated with arginine
   biosynthesis, phenylalanine metabolism, and glutathione metabolism
   pathways (Fig. [104]3d), whereas the DAMs in the turquoise module were
   mostly enriched in tryptophan metabolism, β-alanine metabolism, and
   neuroactive ligand-receptor interaction pathways (Fig. [105]3g).

   Utilizing |KME| for filtering, we identified the top 20 hub DAMs to
   elucidate the key metabolites regulating hormone secretion during goat
   sexual development. Notably, DAMs in the brown module demonstrated the
   highest abundance during the M4 and M6 stages (Fig. [106]3e), whereas
   the DAMs in the turquoise module peaked in expression at the D1 stage
   (Fig. [107]3h). Interestingly, certain DAMs in the brown module, such
   as adrenic acid, DLK, TQH, and FAHFA (22:5/18:2), displayed elevated
   expression levels during both M2 and M4.

Hypothalamic transcriptome WGCNA analysis

   WGCNA analysis was performed to further identify key genes during
   postnatal sexual development in goats using 506 previously identified
   DEGs (Table [108]S9) [[109]42]. Based on a soft threshold of β = 6
   (Figure [110]S5b), a co-expression network modeling was constructed.
   Combining the phenotypic data of serum hormones after birth in Jining
   grey goats (Table [111]S8) [[112]41], 506 genes were further divided
   into 6 modules (Fig. [113]4a). Among them, the blue and green modules
   were significantly positively correlated with serum GnRH and FSH levels
   (p < 0.05) (Fig. [114]4b). The GS-MM correlation results showed that
   there was a significant positive correlation between DEGs and GnRH
   content in the blue and green module (Cor = 0.44, P < 0.01; Cor = 0.41,
   P < 0.01) (Fig. [115]4c and f). The KEGG enrichment analysis results
   indicated that the DEGs within the blue module were significantly
   enriched in several pathways, including the Wnt signaling pathway, TGFβ
   signaling pathway, and IL-17 signaling pathway (p < 0.05)
   (Fig. [116]4d). Interestingly, DEGs in the green module were also
   significantly enriched in the TGFβ signaling pathway (p < 0.05)
   (Fig. [117]4g). The top 20 DEGs in each module were defined as hub
   genes, and the clustered heat map showed that most of the hub genes in
   the blue module were highly expressed during the M2 period
   (Fig. [118]4e), while most hub genes in the green module were highly
   expressed during the M2 and M4 periods (Fig. [119]4h).

Fig. 4.

   [120]Fig. 4
   [121]Open in a new tab

   WGCNA of DEGs identified at different stages of sexual development in
   goat hypothalamic tissue. (a). Hierarchical clustering tree of
   hypothalamic DEGs. (b). Correlation heat map between gene co-expression
   network modules and phenotypic indicators such as serum hormones. The
   red and blue modules indicate positive and negative correlations,
   respectively. (c). GS-MM scatterplot between the blue module and GnRH.
   (d). KEGG enrichment analysis scatterplot of DAMs in the blue module.
   (e). The heatmap of DEGs expression in the top 20 of the blue module.
   (f). GS-MM scatterplot between the green module and GnRH. (g). KEGG
   enrichment analysis scatterplot of differential DAMs in the green
   module. (h). he heatmap of DEGs expression in the top 20 of the green
   module. (i). Correlation heat map between Hub DEGs and Hub DAMs. Red
   signifies a positive correlation, while blue represents a negative
   correlation. Asterisks (*) denote significance levels of p < 0.05

   Additionally, this study computed the Spearman correlation coefficients
   between the hub DAMs identified through metabolomic WGCNA analysis and
   the hub DEGs in the green and blue modules. The results indicated that
   prune homolog 2 with BCH domain (PRUNE2), mitogen-activated protein
   kinase kinase kinase 9 (MAP3K9), potassium voltage-gated channel
   interacting protein 4(KCNIP4), protein kinase cGMP-dependent 2 (PRKG2),
   and PTTG1 regulator of sister chromatid separation, securing (PTTG1)
   were significantly positively correlated with LPC22:5, 2-Arachidonyl
   Glycerol ether, LPE22:5, and Lysops22:5 (P < 0.05). In contrast, WD
   repeat domain 18 (WDR18) exhibited a significant negative correlation
   with these DAMs (P < 0.05). Furthermore, PTTG1 demonstrated a
   significant positive correlation with various metabolites, including
   threonine (P < 0.05). Furthermore, our previous research has shown that
   the transcriptomic sequencing data is consistent with the qRT-PCR data
   (Figure [122]S7), indicating the reliability of our sequencing results
   [[123]42].

Association analysis of metabolome and transcriptome data

   Integrating transcriptomic and metabolomic analyses provides a more
   comprehensive perspective on the regulatory mechanisms underlying
   complex biological processes. Initially, K-means clustering analysis
   and KEGG enrichment analysis were utilized to further explore the key
   genes and pathways involved in goat sexual development (Figure
   [124]S6). Combining these results with metabolomic data, we found that
   the DEGs in cluster 3 of the transcriptomic analysis exhibited a
   similar expression pattern to the DAMs in cluster 1 of the metabolomic
   analysis, while demonstrating an opposing expression pattern compared
   to the DAMs in cluster 6. Both clusters were significantly enriched in
   the glutathione metabolism pathway (p < 0.05) (Figs. [125]2 and
   [126]S6). By integrating these results, we further elucidated the
   characteristics of glutathione metabolism in hypothalamic tissues
   during goat puberty (Fig. [127]5), reflecting the precise regulatory
   patterns of gene expression in the hypothalamus.

Fig. 5.

   [128]Fig. 5
   [129]Open in a new tab

   Schematic diagram of the involvement of DEGs and DAMs in hypothalamic
   tissue in the synthesis of glutathione during postnatal sexual
   development of goats

Discussion

   The hypothalamus is the most upstream tissue of the gonadal axis and
   plays an important role in the sexual development of goats after birth.
   Therefore, this study employed non-targeted metabolomics and
   transcriptomics techniques to analyze the important regulatory
   mechanisms in female goat hypothalamic tissue during sexual maturation.

   In this study, the abundance of lipids and lipid-like molecules, as
   well as organic acids and their derivatives, in hypothalamic tissue
   changed significantly during goat sexual maturation. Among them, most
   lysophosphatidylcholine (LPCs, including LPC14:0, LPC22:3, LPC18:2,
   LPC20:3, LPC16:2, LPC18:3, LPC16:1), lysophosphatidylethanolamine (LPE,
   including LPE16:0, LPE22:3, LPE18:0, LPE20:3, LPE20:2, LPE18:1,
   LPE18:2, LPE14:0), and other metabolites were highly expressed in goats
   after birth. Their expression abundance decreased significantly before
   the first puberty period and then remained stable. Lipids in the
   hypothalamus play a key role in maintaining biological processes such
   as cellular homeostasis, signal transduction, and energy storage
   [[130]43, [131]44, [132]45]. LPCs exert biological functions by
   activating ion channels, promoting the release of inflammatory
   cytokines, and inducing apoptosis and oxidative stress [[133]46]. In
   addition, LPC activates the MAPK signaling pathway, enhances nerve
   growth factor (NGF) production [[134]47], and is also involved in the
   synthesis of acetylcholine in the brain, thereby promoting nerve
   conduction [[135]48, [136]49, [137]50]. However, LPE is structurally
   similar to LPC and is involved in cell signaling and enzyme activation
   [[138]51]. It can also induce neuronal differentiation in PC12 cells
   [[139]52, [140]53]. The highest levels of LPE and LPC were found in the
   D1 stage of Jining goats, possibly related to the regulation of lipid
   metabolism in the hypothalamus after birth, supporting thermoregulation
   and nervous system development [[141]54]. FAHFAs are a new class of
   endogenous lipids that have been shown to promote glucose metabolism
   and have anti-inflammatory and antioxidant effects [[142]55]. Unlike
   most FAHFAs in the serum during goat sexual development, which are
   highly expressed at the M4 stage [[143]56], hypothalamic FAHFA-like
   metabolites exhibit higher expression levels at the D1 and M2 stages.
   This may be related to the role of FAHFAs in regulating glucose
   homeostasis, helping neurons in the hypothalamus perceive energy status
   in vivo, which in turn affects their growth and development [[144]57,
   [145]58, [146]59]. In addition, ACar is produced when L-carnitine
   combines with fatty acids, promoting thermogenesis in brown adipose
   tissue and participating in multiple energy metabolism pathways
   [[147]60, [148]61]. It was found that most ACar metabolites were lowest
   after birth and increased with the gradual maturation of sexual
   development. This suggests that ACar may be involved in regulating
   energy homeostasis during the development of goats after birth. As the
   body develops and body weight reaches a specific threshold, signals are
   transmitted to the hypothalamus, promoting the secretion of GnRH and
   triggering sexual maturation in animals [[149]62, [150]63, [151]64].
   Overall, these studies indicate that lipids play a key role in goat
   sexual maturation and nervous system development.

   Metabolites with different expression patterns may play a similar role
   during sexual development in goats. In this study, six different
   expression patterns of metabolites were identified, which were mainly
   enriched in histidine metabolism, β-alanine metabolism, neuroactive
   ligand-receptor interaction, glutathione metabolism and other pathways.
   All of these, including Anserine, L-histidine, and Carnosine, were
   enriched in the histidine metabolism and β-alanine metabolism pathways
   [[152]65]. In mammals, histaminergic neurons are predominantly found in
   the hypothalamus, and L-Histidine is converted to histamine to regulate
   feeding behavior and energy metabolism, and is involved in the
   synthesis and release of growth hormone and thyroid-stimulating hormone
   (TRH) [[153]66]. Carnosine is a dipeptide composed of alanine and
   L-histidine, and its analogue Anserine is present in skeletal muscle
   and nervous system, and has been identified in the rat hypothalamus
   [[154]67]. Carnosine and Anserine have antioxidant activity and also
   modulate nerve signaling [[155]68, [156]69]. Moreover, neuroactive
   ligand-receptor interactions have been confirmed to play an important
   role in hypothalamic reproductive regulation [[157]70, [158]71].
   Metabolites identified in this study, such as taurine and
   4-aminobutyric acid (GABA), are enriched in this pathway. Taurine is
   one of the most abundant free amino acids in the hypothalamus and
   pituitary tissue of animals [[159]72], with the highest levels found in
   the fetal brain, which decrease with age [[160]73]. This finding aligns
   with the age-dependent reduction of taurine observed in the sexual
   maturation process of goats in this study. Research has shown that
   taurine can inhibit the release of pituitary LH via the hypothalamus
   [[161]74]. Additionally, taurine-mediated GnRH release may activate the
   GABAergic system, significantly suppressing HPG axis activity and
   leading to delayed sexual development [[162]75]. GABA stimulates GnRH
   release in immature male rats while inhibiting GnRH secretion in adult
   rats [[163]76]. In this study, GABA levels increased from D1 to M2,
   potentially related to the rise in GnRH secretion from the
   hypothalamus. Our study reveals dynamic changes in hypothalamic
   metabolites during sexual maturation in goats. During this process,
   changes in neuron development and hormone secretion are accompanied by
   changes in the expression levels of signaling-related metabolites
   (Anserine, L-histidine, Carnosine, taurine, and 4-aminobutyric).

   Through WGCNA analysis by transcriptomics and metabolomics, we
   identified a number of DEGs and DAMs associated with hypothalamic GnRH
   secretion and performed spearman correlation analysis on them. This
   opens up the possibility of further research into the association
   between these DEGs and DEMs. MAP3K9 regulates signal transduction
   through mitogen-activated protein kinase (MAPK) and c-Jun N-terminal
   kinase (JNK) pathways [[164]77]. The MAPK signaling pathway is critical
   to neurodevelopment and GnRH hormone secretion [[165]78, [166]79,
   [167]80]. Studies have shown that PRUNE2 may be involved in
   neurodevelopmental processes in sheep [[168]81]. KCNIP4 is a small
   calcium-binding biomolecule that regulates neuronal excitability in
   response to changes in intracellular calcium levels [[169]82, [170]83],
   and has also been identified as a candidate gene for reproductive
   traits in buffalo [[171]84], poultry [[172]85, [173]86], and sheep
   [[174]87]. 2-Arachidonoylglycerol (2-AG) is a key regulator of
   neurotransmitter release [[175]88] and plays an important role in food
   intake and metabolism [[176]89]. In our study, the expression levels of
   MAP3K9, PRUNE2 and KCNIP4 gradually increased with sexual development
   after birth, reaching a peak at M4 stage, and all of them were
   significantly positively correlated with 2-Arachidonyl Glycerol Ether.
   These results suggest that MAP3K9, PRUNE2, KCNIP4, and 2-Arachidonyl
   Glycerol Ether may be involved in the sexual development of goats after
   birth by regulating hypothalamic GnRH neuronal development and energy
   metabolism.

   The results of joint analysis showed that the glutathione metabolism
   pathway was enriched in both transcriptome and metabolome. Glutathione
   (GSH) is a tripeptide composed of glutamate, cysteine, and glycine that
   is involved in biological processes such as redox, apoptosis, and
   neuromodulation [[177]90, [178]91, [179]92]. Glutamate and its
   receptors are present in all key hypothalamic regions closely related
   to reproduction [[180]93] and have been shown to regulate mammalian
   GnRH neuronal development and neuroendocrine function [[181]94,
   [182]95]. In this study, the elevated L-glutamate content in the
   hypothalamus at the D1 and M4 stages may relate to its regulation of
   GnRH neuronal development and secretion associated with sexual
   maturity. GSH is mainly converted from Glutathione disulfide (GSSG).
   IHD1, on the other hand, is mainly involved in the conversion of GSSG
   to GSH via NADPH from the pentose phosphate pathway [[183]96]. Goats
   have higher levels of hypothalamic GSSH and IDH1 after birth, and the
   ability of GSSG to convert to GSH increases, ultimately leading to
   higher postnatal hypothalamic GSH content, consequently, this may be
   related to postnatal hypothalamic neuronal development [[184]97].
   Therefore, homeostasis of glutathione metabolism may play an important
   role in the sexual maturation of goats.

   However, this study has some limitations. Firstly, molecular biology
   experiments have not yet been conducted on the identified key genes and
   metabolites. Additionally, since this research focuses on Jining grey
   goats, further validation is needed to determine whether these findings
   are applicable to other early-maturing goat breeds.

Conclusion

   In this study, metabolites such as Anserine, L-histidine, Carnosine,
   taurine, and 4-aminobutyric may affect neuron development and hormone
   secretion processes during sexual maturation in goats via histidine and
   phenylalanine metabolic pathways. Furthermore, MAP3K9, PRUNE2, and
   KCNIP4 may jointly regulate the development and energy metabolism of
   hypothalamic GnRH neurons in conjunction with DAMs such as
   2-Arachidonyl Glycerol Ether. Moreover, the homeostasis of glutathione
   metabolism may be involved in the regulation of GnRH hormone secretion
   during sexual maturation. This study offers a new perspective for
   understanding the regulatory mechanisms of the hypothalamus in goat
   sexual development, which may yield important insights for enhancing
   reproductive traits in future animal breeding efforts.

Electronic Supplementary Material

   Below is the link to the electronic supplementary material.
   [185]Supplementary Material 1^ (529.9KB, xlsx)
   [186]Supplementary Material 2^ (1.7MB, docx)

Acknowledgements