Abstract

Background

   Heat stress (HS) is a significant challenge in poultry, negatively
   impacting feed efficiency and survival. These adaptive responses could
   lead to disrupted lipid metabolism, impaired immunity, and neural
   damage. We hypothesized that the neuroendocrine system plays a central
   role in HS responses by activating the hypothalamic-pituitary-adrenal
   and sympathetic-adrenal-medullary axes. This study explores the neural
   transcriptomic changes and metabolomic profile subjected to HS, aiming
   to unravel the complexity of nervous adaptive responses, including the
   disrupted lipid metabolism, impaired immunity, and neural damage
   processes.

Results

   Based on the designed chicken HS model, we identified 1,871 serum
   metabolites, with 149 differentially accumulated metabolites (DAMs),
   including ursodeoxycholic acid, niacinamide, and prostaglandin E1.
   These metabolites indicated disrupted lipid homeostasis and potential
   adaptive mechanisms to mitigate stress. Key metabolite biomarkers were
   identified using random forest, including 2-(Methylthio)phenol, and
   Deca-2,4,6-triynedioic acid, enriched in immune and stress-related
   pathways such as glutathione metabolism and nicotinamide metabolism.
   Then, the comparative transcriptome analysis of the hypothalamus and
   pituitary revealed 358 differentially expressed genes (DEGs),
   highlighting tissue-specific adaptive pattern. The hypothalamus DEGs
   enriched in neurotransmitter biosynthesis and Notch signaling pathways,
   while the pituitary DEGs exhibited changes in calcium and insulin
   signaling pathways. Finally, by integrating omics data with machine
   learning, we uncovered key stress-relief mechanisms involving adenosine
   accumulation and adrenaline signaling. Adenosine signaling emerged as a
   critical pathway for thermoregulation through vascular relaxation,
   while adrenaline signaling regulated calcium homeostasis and muscle
   contraction under stress.

Conclusions

   This study provides comprehensive insights into the molecular
   adaptations of HS in broiler, emphasizing the roles of lipid
   metabolism, neuroendocrine regulation, and vascular dynamics. The
   findings highlight the potential of integrating transcriptomics and
   metabolomics with machine learning to identify biomarkers and pathways
   involved in HS adaptation. These results could provide an effective
   theoretical basis of the preventive and treatment measures for HS.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12864-025-11881-7.

   Keywords: Biomarkers, Heat stress, Metabolomics and transcriptomics,
   Neuroendocrine regulation

Background

   The poultry industry is globally essential to food security and
   economic stability. However, with unpredictable climate patterns,
   managing heat stress (HS) has become a challenge, as it can reduce
   growth rates, feed efficiency, and survival rate [[36]1, [37]2]. As a
   consequence, there is a deepening research focus on the biological
   underpinnings of how HS affects broilers [[38]3], for example, heat
   shock proteins (HSPs) [[39]4, [40]5], xanthine dehydrogenase (XDH), and
   periostin (POSTN) [[41]6] are abundant biomarker proteins activated in
   response to HS. While recent studies have increasingly focused on the
   complex responses of blood metabolites and neural tissue as indicators
   of adaptive or maladaptive physiological changes [[42]7, [43]8].

   The neuroendocrine system (NES) played critical role in HS responses.
   HS activates the sympathetic-adrenal-medullary (SAM) axis and the
   hypothalamic-pituitary-adrenal (HPA) axis. The SAM axis releases
   adrenaline, leading to tachycardia and increased respiratory rates in
   chickens, which rapidly accelerates glucose synthesis [[44]9].
   Conversely, the HPA axis releases corticosterone more slowly than
   adrenaline but with longer-lasting effects. Consequently,
   corticosterone is widely recognized as a biomarker of stress. Elevated
   corticosteroid levels during HS can cause hypercholesterolemia,
   cardiovascular issues, gastrointestinal lesions, and impaired immune
   function [[45]10]. Vasodilation, excessive sweating, and rapid panting
   increase blood pH, resulting in respiratory alkalosis. This disrupts
   bicarbonate and ionized calcium balance, impairing egg production
   [[46]11]. Specifically, whole egg weight, eggshell thickness, weight,
   and strength decrease by 3.24%, 1.2%, 9.93%, and 0.66%, respectively
   [[47]12]. The prolonged HS decreases feed intake by 16.4%, reduces body
   weight by 32.6%, and increases feed conversion ratio by 25.6% in
   broilers at 42 days age [[48]13]. Chronic HS disrupts lipid metabolism
   and lowers muscle protein content [[49]14]. High temperatures affect
   broiler production via cecal microbiota and their metabolites, mediated
   by the glucose-sensing receptor–gut peptide–glucolipid metabolism axis
   [[50]15, [51]16]. In parallel, HS responses are primarily managed by
   the hypothalamus and pituitary, which integrate stress signals to
   regulate physiological and behavioral adaptations to mitigate HS
   [[52]17, [53]18].

   The thermoreceptors were primarily found in hypothalamus, eliciting
   responses similar to those of cutaneous thermoreceptor signals [[54]18,
   [55]19]. These hypothalamic signaling provided insights into the
   effects of HS on the nervous system. In addition, HS can lead to
   substantial damage to organs, such as acute degeneration of myocardial
   cells [[56]20], inflammatory cell infiltration in the liver [[57]21],
   and liver injury due to increased autophagy via the Nrf2-Keap1
   signaling pathway [[58]22, [59]23]. In light of these intricate
   biological interactions, traditional analytical methods are often
   insufficient to capture the full scope of interaction network.

   Recently, machine learning (ML) methods are increasingly adapted in
   biomarker discovery because of their substantial promise in handling
   multi-omics datasets and identifying complex interactions [[60]24,
   [61]25]. A recent study integrating non-targeted metabolomics with
   random forest (RF) machine learning analysis identified d-xylose,
   xanthine, and pyruvaldehyde as potential key marker compounds linked to
   freshness deterioration in chilled pork [[62]26]. In this study, we
   employed non-targeted metabolome approach, identifying 1,871
   metabolites across HS and control cohort of chronic heat stress in
   white-feather broiler. We apply RF algorithm to identify hub-elements
   comprising nervous DEGs and key serum metabolites, and reveal
   interaction networks that adenosine accumulation and adrenaline
   signaling. These findings offer profound insights into the metabolomic
   and neural transcriptomic adaptations underlying HS responses in
   chickens.

Result

Identification of metabolites accumulated in the HS

   A total of 1,871 metabolites were detected in serum samples from 86
   chicks (Control: n = 45; HS: n = 41). The largest category was lipids
   and lipid-like molecules, comprising 597 metabolites (31.9%), followed
   by organic acids and derivatives (393 metabolites, 21.0%), and
   Organoheterocyclic compounds (285 metabolites, 15.2%). Organic oxygen
   compounds constituted 161 metabolites (8.6%), while benzenoids also
   contributed 152 metabolites (8.1%). Other notable groups included
   phenylpropanoids and polyketides with 89 metabolites (4.8%), and
   nucleosides, nucleotides, and analogues, comprising 51 metabolites
   (2.7%). A small fraction of metabolites (56 compounds, 3.0%) remained
   unclassified, potentially due to limitations in current annotation
   databases (Additional file 1: Table S1, Fig. 1a). In addition,
   categories with less than 50 metabolites are uniformly classified into
   the Other category (Fig. [63]1b).

Fig. 1.

   [64]Fig. 1
   [65]Open in a new tab

   Categories, numbers and proportions of annotated metabolites. Super
   class with fewer than 50 metabolites are grouped into the Other
   category. A, Taxonomic categories of annotated metabolites. B,
   Categories, numbers and proportions of other class metabolites

   The OPLS-DA analysis for metabolites indicated a significant separation
   between HS and control group (Fig. [66]2a), where the R2Y (96%) and Q2Y
   (80%) of this model were both nearly 1, reflecting strong model
   performance. The permutation test results demonstrated that the R2Y and
   Q2Y values of randomized models (scatter points) were consistently
   lower than those of the original model (solid line), suggesting that
   the original model exhibits robust predictive ability (Additional file
   2: Figure S1). Among these metabolites, 149 DAMs were identified
   between the HS and control groups (Fig. [67]2c; Additional file 3:
   Table S2). Of these, 117 metabolites (78.5%) were significantly reduced
   in the HS group compared to controls (e.g., ursocholic acid and
   niacinamide), while 32 metabolites (21.5%) were elevated in the HS
   group (e.g., prostaglandin E1 (PGE1), polidocanol, and
   3-decenoylcarnitine). The details of the top 50 DAMs is provided in
   Additional file 4: Figure S2. Ursocholic acid, an intermediate in bile
   acid breakdown, generated secondary bile acids, which regulated lipid
   and glucose metabolism, modulate liver inflammation, and support
   productivity, gut health, heat stress relief, and liver protection
   [[68]27]. Prostaglandin E1 is known to regulates blood pressure and
   vascular dilation and can alleviate AngII-induced cardiac hypertrophy
   by activating EP3 receptors and upregulating Ntn1 [[69]28]. Lipids and
   lipid-like molecules were the most abundant differential metabolites
   (69 compounds, 46%), followed by organic acids and derivatives (13%)
   and organic oxygen compounds (11%) (Fig. [70]2d; Additional file 5:
   Table S3). We observed that the result of OPLS-DA analysis was
   consistent with DAMs, where 77.8% DAMs with VIP > 1 (Fig. [71]2b).

Fig. 2.

   [72]Fig. 2
   [73]Open in a new tab

   Identification and characterization of differentially accumulated
   metabolites (DAMs) between heat stress (HS) and control groups. A
   OPLS-DA of untargeted metabolomics data comparing HS group (colored in
   purple) and control group (colored in green). B Distribution of
   variable importance in projection (VIP) values for DAMs. C Volcano plot
   of the detected metabolites between the HS group and control group.
   DAMs are colored in purple (upregulated) and green (downregulated). D
   Histogram of number distributions of DAMs in super class: upregulated
   (purple), downregulated (green)

Machine learning-based selection of signature metabolites

   Here, we trained a RF model to screen metabolite markers based on DAMs
   with VIP > 1. Based on the average rank and cross-validation error, the
   top classification RF (100 trees) recognized 17 metabolites as feature
   indicators, indicating excellent classification performance
   (AUC = 0.97) (Fig. [74]3b). Among these, XBGUTAISSLTGIB-UHFFFAOYSA-N,
   2-(Methylthio)phenol, and Deca-2,4,6-triynedioic acid were considered
   most significant (Fig. [75]3a).

Fig. 3.

   [76]Fig. 3
   [77]Open in a new tab

   Selection of characteristic metabolites and KEGG pathway enrichment
   analysis. A Histogram of Gini importance distribution for
   characteristic metabolites. B Receiver operating characteristic (ROC)
   of the Random Forest (RF) model. C KEGG pathway enrichment analysis of
   metabolites. The bar chart represents the gene ratio for the enrichment
   analysis, while the bubble indicates the significance of the metabolite
   enrichment, with darker colors denoting higher significance. D
   Proteomaps analysis of KEGG pathway enrichment significance

   These metabolites were mapped to the 11 immune and stress related KEGG
   metabolic pathways (Fig. [78]3c; Additional file 6: Table S4),
   including glutathione metabolism, ascorbate, aldarate metabolism,
   nicotinate, and nicotinamide metabolism (Fig. [79]3d). Glutathione, a
   key antioxidant, is known to protected sulfhydryl-containing proteins
   and enzymes in cell membranes from oxidation [[80]29]. Nicotinate and
   nicotinamide, both forms of vitamin B3, can be interconvertible in
   metabolic pathways. Nicotinamide minimizes reactive oxygen species and
   free radicals, shielding cells from oxidative damage and maintaining
   integrity [[81]30].

Functional enrichment analysis of the neural transcriptomic adaptations
underlying HS responses in chickens

   To investigate the molecular basis of neural tissue of the metabolic
   differences detected in HS and control group, transcriptome sequencing
   was performed on pituitary (Control: n = 18; HS: n = 19) and
   hypothalamus (Control: n = 19; HS: n = 19) tissues. Comparative
   transcriptomic analysis between HS and control groups revealed
   significant differential gene expression patterns in both hypothalamic
   and pituitary tissues (Fig. [82]4a). In the hypothalamus, we identified
   160 differentially expressed genes (DEGs), consisting of 74 genes
   significantly upregulated and 86 genes downregulated in HS group
   relative to controls (Additional file 7: Figure S3a; Additional file 8:
   Table S5). The pituitary exhibited an even more pronounced response,
   with 198 DEGs detected, of which 183 genes showed increased expression
   while only 15 genes were suppressed under heat stress conditions
   (Additional file 7: Figure S3b; Additional file 9: Table S6). Notably,
   five genes (ENSGALG00000020895, glutamate metabotropic receptor 2
   (GRM2), nucleoredoxin-like 2 (NXNL2), G protein-coupled receptor 183
   (GPR183), and DExH-box helicase 58 (DHX58)) were upregulated in the
   hypothalamus but downregulated in the pituitary. Additionally, the
   Homocysteine inducible ER protein with ubiquitin like domain 1
   (HERPUD1) gene exhibited upregulated expression in both the
   hypothalamus and pituitary(Fig. [83]4b).

Fig. 4.

   [84]Fig. 4
   [85]Open in a new tab

   Transcriptome analysis between the HS group and control group. A the
   distribution of fold change (FC) values and p-adjust for differentially
   expressed genes (DEGs) in the hypothalamus and pituitary. B The
   comparison of the number of DEGs in the hypothalamus and pituitary. The
   pie chart represents the proportion of upregulated and downregulated
   genes in the hypothalamus and pituitary. The venn diagram illustrates
   the overlap of upregulated genes and downregulated genes between the
   hypothalamus and pituitary. C The interaction network of the Top 10 hub
   genes. Triangles refer to upregulated genes, and circles refer to
   downregulated genes. D The 15 most significantly enriched pathways for
   DEGs in the hypothalamus. The bar chart represents the gene ratio from
   the enrichment analysis, while the bubble chart illustrates the
   significance level of the enrichment analysis. E The 15 most
   significantly enriched pathways for DEGs in the pituitary. The bar
   chart represents the gene ratio from the enrichment analysis, while the
   bubble chart illustrates the significance level of the enrichment
   analysis

   The protein-protein interactions among DEGs were predicted using the
   STRING database and analyzed for network topology using the CytoHubba
   plugin in Cytoscape [[86]31]. The interaction network comprised 168
   genes with a total of 910 interactions, and the topological parameters
   for each node (Additional file 10: Figure S4; Additional file 11: Table
   S7). The interactions network of top 10 genes with the highest Maximum
   Clique Centrality (MCC) identified that the fibrinogen alpha chain
   (FGA) gene as the core gene, which enhanced immune responses through
   innate immunity and T cell-mediated pathways(Fig. [87]4c) [[88]32].

   Functional enrichment analysis of DEGs in neural tissues highlighted
   distinct biological roles for the hypothalamus and pituitary tissues
   (Fig. [89]4d-e, Additional file 12: Figure S5, Additional file 13–14:
   Table S8−9). DEGs in the hypothalamus were enriched in the
   Phenylalanine, Tyrosine, and Tryptophan metabolism pathways, suggesting
   an association with the generating of neurotransmitters like dopamine,
   norepinephrine, and serotonin, which were known to serve as pivotal
   modulators of basic metabolic processes [[90]33, [91]34]. The Notch
   signaling pathway was highly enriched, indicating a potential role for
   the hypothalamus’s involvement in neural adaptability [[92]35, [93]36].
   In contrast, DEGs in the pituitary were enriched in calcium signaling
   and insulin signaling pathways, which were essential for hormone
   release, metabolic regulation, and response to external signals
   [[94]37]. The insulin signaling pathway further emphasizes the
   pituitary’s role in maintaining metabolic balance by regulating cell
   growth, hormone secretion, and metabolic processes [[95]38].
   Additionally, the enrichment of glutathione metabolism in the pituitary
   suggests a protective role against HS, ensuring proper function under
   high metabolic demands [[96]39].

Characterization of key genes involved in adenosine accumulation induced by
HS

   Based on metabolite markers between HS and control group, comparative
   transcriptome analysis of two tissues, and function enrichment
   analysis, we constructed stress-relieving networks of adenosine
   accumulation induced by HS (Fig. [97]5a). During heat stress,
   vasodilation occurred alongside increased blood flow, facilitating heat
   dissipation to maintain normal body temperature. The vascular tone
   regulation involves multiple factors, including adenosine-mediated
   signaling. Adenosine binds to adenosine A2 receptors (ADORA2), which
   activated G protein-coupled signaling pathways. The stimulatory G
   protein (Gs) enhances adenylate cyclase (AC) activity, encoded by the
   ADCY1 gene (p-adjust < 0.05, upregulated 1.12-fold) (Fig. [98]5b),
   thereby increasing cyclic adenosine monophosphate (cAMP) production.
   Elevated cAMP activates protein kinase A (PKA), which phosphorylates
   and inactivates myosin light chain kinase (MLCK) [[99]40].

Fig. 5.

   [100]Fig. 5
   [101]Open in a new tab

   Analysis of key genes associated with adenosine accumulation induced by
   heat stress (HS). A The regulatory pathways of vascular contraction and
   relaxation (created with biorender). B Boxplot of relevant genes
   (ENSGALG00000019233, ENSGALG00000015203, MYH11, ACAT2, MYL3, ADCY1) and
   metabolites (Adenosine, Sphingosine 1-phosphate (SIP)). In the figure,
   the green box plot represents the Control group (Control, C), and the
   purple box plot represents the heat stress group (HS, T)

   This inhibition reduced the phosphorylation of myosin light chains
   (MLC), encoded by the MYL3 (Myosin, light chain 3) gene
   (p-adjust < 0.05, downregulated by 0.57-fold)(Fig. [102]5b), and
   attenuated the interaction between myosin heavy chains (MHC, encoded by
   the MYH11 gene) and actin filaments. Consequently, vascular smooth
   muscle relaxed, promoting vasodilation and contributing to
   thermoregulation during HS [[103]41].

HS-induced metabolite action pathways of adrenaline signaling

   HS activated the HPA and SAM axes [[104]17], leading to a marked
   increase in the secretion of glucocorticoids, adrenaline, and
   isoproterenol (Fig. [105]6a). Adrenaline play a central role in
   modulating muscle contraction and relaxation through various signaling
   pathways [[106]42], aiding poultry in coping with acute stress
   conditions.

Fig. 6.

   [107]Fig. 6
   [108]Open in a new tab

   Metabolic pathways of adrenaline signaling induced by heat stress (HS).
   A Mechanism of Action of Adrenaline(created with biorender). B, Boxplot
   of relevant genes (ADCY1, ATP1B1, PLN, MYH15, MYL3, CACNG3, CAMK2B) and
   metabolites (Isoproterenol). In the figure, the green box plot
   represents the Control group (Control, C), and the purple box plot
   represents the heat stress group (HS, T)

   Specifically, adrenaline binds to β₁ and β₂ adrenergic receptors,
   thereby activating adenylate cyclase (AC). This enzyme converts ATP to
   cyclic adenosine monophosphate (cAMP), subsequently elevating
   intracellular calcium concentrations via downstream signaling
   [[109]43]. The expression of myosin-related genes was significantly
   altered under stress conditions. MYH15 (Myosin, heavy chain 15)
   (p-adjust < 0.05; 0.76-fold down regulation; Fig. [110]6b) and MYL3
   (p-adjust < 0.05; 0.57-fold down regulation; Fig. [111]6b) encode the
   heavy and light chains of myosin, respectively, both of which are
   essential structural components of the contractile apparatus.

   Concurrently, the calcium regulatory network was perturbed: the
   phospholamban (PLN) gene (p-adjust < 0.05, downregulated by 0.29-fold)
   (Fig. [112]6b) encoded phosphorylated phospholamban (PLB), a
   sarcoplasmic reticulum protein that controls calcium reuptake and
   release by modulating SERCA (sarco/endoplasmic reticulum Ca²⁺-ATPase)
   activity [[113]44]. In contrast, CACNG3 (p-adjust < 0.05, upregulated
   by 1.10-fold) (Fig. [114]6b) encoded the γ subunit of dihydropyridine
   receptor (DHPR), a voltage-gated calcium channel critical for
   excitation-contraction coupling in muscle cells. Furthermore, ATP1B1
   (p-adjust < 0.05, upregulated by 1.02-fold) (Fig. [115]6b) encoded the
   β-subunit of sodium-potassium ATPase (INaK), which maintains
   cardiomyocyte ion homeostasis by actively transporting sodium and
   potassium ions across the membrane [[116]45].

Discussion

   HS undeniably impacts broiler growth and physiological function in ways
   that are as multifaceted as they are concerning, particularly when
   viewed in the larger context of global warming, which causes economic
   losses and welfare concerns in poultry production. The stress response
   mechanisms in broilers, especially the activation of the SAM and HPA
   axes, with adrenaline rapidly inducing effects like tachycardia and
   increased respiratory rates [[117]17]. Although corticosterone can lead
   to more insidious outcomes such as hypercholesterolemia, cardiovascular
   stress, gastrointestinal distress, and immune suppression [[118]42]. In
   these processes, the hypothalamus and pituitary gland, are central
   players orchestrating these responses, and our comparative analysis
   between stressed and non-stressed groups uncovered a series of dramatic
   shifts in gene expression profiles of the nervous system and the serum
   metabolite profiles. Specifically, we identified 149 DAMs and 358 DEGs
   across tissues. Integration with ML model revealed their intricate
   connections. Our analysis specifically highlighted metabolites
   associated with important biological pathways, particularly those
   involved in immune response and stress adaptation, such as glutathione
   metabolism and nicotinate/nicotinamide metabolism.

   For metabolites accumulation, we observed both key indicators of lipid
   homeostasis, free fatty acids and cholesterol, showed a marked
   increase, further highlighting how HS disrupts lipid metabolism—an
   observation that aligns well with earlier studies demonstrating similar
   patterns of lipid metabolism disorder in heat-stressed broilers
   [[119]46]. Interestingly, OPLS-DA analysis validated these DAMs, such
   as ursodeoxycholic acid, niacinamide, and prostaglandin E1, which not
   only play critical roles in lipid regulation [[120]47], but also are
   implicated in blood pressure control and the alleviation of
   stress-related damage [[121]48]. Prostaglandin E1, in particular, was
   previously noted by Hii et al. to ameliorate intestinal and systemic
   inflammation as well as multi-organ dysfunction in animals under HS,
   further underscoring its relevance here [[122]49]. Alongside these
   insights, ML models have been increasingly utilized for metabolite
   prioritization, which effectively narrow down the vast pool of
   potential metabolites, enhancing the efficiency and accuracy of marker
   metabolites detection [[123]50, [124]51]. Metabolites like
   2-mercaptophenol, nicotinamide, and nigellidine emerged as key players,
   all closely tied to the intricate mechanisms underlying HS adaptation
   [[125]29, [126]52].

   The 3-hour post-heat stress (HS) sampling timepoint was selected based
   on its alignment with the widely adopted 2–4 h acute HS exposure
   paradigm established in prior research [[127]53, [128]54]. Substantial
   evidence indicates that stress-responsive genes, including heat shock
   transcription factor 1 (HSF1) [[129]55–[130]57], nuclear receptor
   subfamily 3 group C member 1 (NR3C1) [[131]58, [132]59], and
   corticotropin releasing hormone (CRH) [[133]60, [134]61] gene, exhibit
   peak expression dynamics synchronously with cortisol surges and
   chromatin remodeling at glucocorticoid receptor-binding sites. Although
   this experimental design effectively captures early-phase
   transcriptional adaptations such as Notch signaling activation and
   neurotransmitter biosynthesis, it presents inherent limitations in
   resolving delayed metabolic feedback mechanisms like nicotinamide
   salvage pathways [[135]62]. Furthermore, the single timepoint approach
   impedes differentiation between transient compensatory responses (e.g.,
   acute cholesterol elevation) and pathological dysregulation. Future
   investigations should implement serial sampling protocols (0.5 h, 3 h,
   6 h) to elucidate kinetic relationships within adenosine signaling
   pathways and vascular adaptation processes, while maintaining the
   3-hour benchmark to facilitate cross-study comparisons.

   The observed differential expression patterns of key genes in the
   hypothalamus and pituitary under heat stress provide valuable insights
   into the potential molecular mechanisms underlying the neuroendocrine
   response to thermal challenge. Notably, the contrasting regulation of
   GRM2, NXNL2, GPR183, and DHX58 between the hypothalamus (upregulated)
   and pituitary (downregulated) suggests tissue-specific roles in
   mediating heat stress adaptation. The upregulation of GRM2 (glutamate
   metabotropic receptor 2) gene in the hypothalamus aligns with its role
   as a modulator of excitatory neurotransmission. Glutamate signaling is
   critical for stress-responsive pathways [[136]63], and increased GRM2
   expression may reflect heightened neuronal activity to coordinate
   systemic responses to heat stress. Conversely, its down regulation in
   the pituitary could indicate a feedback mechanism to attenuate
   excessive stress hormone release [[137]64], such as ACTH or prolactin,
   to prevent overactivation of the hypothalamic-pituitary-adrenal (HPA)
   axis [[138]65]. The consistent up regulation of HERPUD1 in both tissues
   underscores its critical role in endoplasmic reticulum (ER) stress
   response [[139]66, [140]67]. Heat stress likely initiates protein
   misfolding cascades, with HERPUD1 activation serving as a regulatory
   node that coordinates the unfolded protein response (UPR) to
   reestablish proteostatic equilibrium [[141]68]. Mechanistically, this
   proteostatic mechanism is directly governed by the transcription factor
   Nrf1 through its binding to antioxidant response elements in the
   HERPUD1 promoter—a regulatory relationship evidenced by abolished
   HERPUD1 induction in Nrf1 knockout cells [[142]69]. This ubiquitous
   upregulation suggests ER stress as a central feature of heat stress
   pathophysiology across neuroendocrine tissues.

   Hypothalamic DEGs were highly enriched in pathways related to
   neurotransmitter synthesis (e.g., glutamate and GABA pathways) and
   neural plasticity, underscoring the hypothalamus’ central role in
   orchestrating neuroendocrine responses to physiological and
   environmental stressors [[143]18, [144]19]. Meanwhile, pituitary DEGs
   highlighted its involvement in calcium and insulin signaling pathways,
   suggesting its pivotal role in maintaining metabolic homeostasis
   through hormone secretion [[145]70], a link further evidenced by the up
   regulation of NF-κB signaling genes, which hints at the activation of
   inflammatory responses in the context of HS [[146]71]. Furthermore,
   genes tied to lipid metabolism, such as PPARα, showed altered
   expression patterns that echoed the metabolomic data, adding layer of
   corroboration to the story of disrupted lipid regulation [[147]46].
   However, while these tissue-specific findings elucidate HPA axis
   engagement, the cross-tissue signaling networks connecting peripheral
   organs (e.g., gut, liver) to central neuroendocrine regulation remain
   unresolved. Emerging evidence suggests that peripheral metabolites may
   coordinate systemic adaptation to heat stress through interorgan
   signaling, though the precise molecular conduits remain
   uncharacterized. For instance, gut-derived ursodeoxycholic acid (UDCA),
   which is significantly elevated under heat stress, potentially crosses
   the blood-brain barrier via organic anion transporters (e.g., OATP1A4
   [[148]72]) and binds to the farnesoid X receptor (FXR) in the
   hypothalamus [[149]73, [150]74], thereby modulating the expression of
   glucocorticoid sensitivity-related genes (e.g., NR3C1) [[151]75]. This
   mechanism aligns with hepatic observations, where UDCA suppresses
   CYP7A1 expression to mitigate cholesterol accumulation, synchronously
   alleviating lipid metabolic dysregulation in both central and
   peripheral tissues [[152]76]. Another key metabolite, PGE1, may exhibit
   dual regulatory roles: (1) activating hypothalamic EP3 receptors to
   stimulate corticotropin-releasing hormone (CRH) release [[153]28,
   [154]77], and (2) enhancing hepatic gluconeogenesis via the cAMP-PKA
   pathway to fuel thermoregulatory energy demands [[155]78, [156]79].
   Furthermore, nicotinamide trafficking across tissues may integrate
   metabolic and epigenetic regulation—the hypothalamus salvages
   nicotinamide via nicotinamide phosphoribosyltransferase (NAMPT) to
   sustain NAD + synthesis for hormonal stability [[157]80], while the
   heart methylates excess nicotinamide via nicotinamide
   N-methyltransferase (NNMT) [[158]81], providing substrates for
   hypothalamic CRH gene demethylation [[159]82]. Although these
   mechanisms require further validation, they suggest that metabolites
   serve as molecular bridges, translating environmental stress signals
   into epigenetic and metabolic responses. Future studies could employ
   spatially resolved transcriptomics to map metabolite gradients within
   vascular networks, coupled with conditional knockout models targeting
   gut-liver-brain signaling nodes to dissect multi-organ crosstalk
   mechanisms.

   The vascular adaptations during HS also deserve attention, as
   vasodilation and enhanced blood flow are vital for dissipating heat and
   stabilizing body temperature, mechanisms strongly influenced by
   adenosine signaling [[160]83]. This pathway mediates vasodilation
   through activation of adenosine A2 receptors and associated G
   protein-coupled signaling [[161]84], as highlighted by both our
   findings and prior studies like Ritchie et al., which explored how
   caffeine, an adenosine receptor antagonist, altered thermoregulation
   under heat stress [[162]85]. Building on this, our exploration of
   adrenaline signaling pathways revealed nuanced regulation of genes such
   as MYH15, PLN, CACNG3, and ATP1B1, all of which play roles in muscle
   contraction, calcium handling, and ion homeostasis.

   These results align well with existing research, showing how
   β-adrenergic signaling governs calcium dynamics, including the
   phosphorylation of skeletal muscle CaV1.1 channels to boost calcium
   influx during stress, a mechanism critical for muscle contractility
   [[163]86]. The role of phospholamban, too, emerged as central; its
   regulation of the SERCA pump was shown to impact calcium reuptake,
   influencing muscle strength and recovery [[164]87]. Altogether, these
   results shed light on the complex, interconnected processes that
   underpin the adaptation of broilers to HS, encompassing neuroendocrine
   control, metabolite accumulation, and transcriptomic shifts.

Conclusions

   This study systematically investigates the physiological and molecular
   adaptations of broilers to HS, integrating metabolomic and
   transcriptomic profiling. HS disrupts lipid metabolism, as shown by key
   metabolites like ursodeoxycholic acid and niacinamide, which aid stress
   mitigation and vascular stability. These metabolites are functionally
   implicated in stress adaptation and vascular function regulation
   through their antioxidant and anti-inflammatory activities.
   Transcriptomic analysis highlights the roles of the hypothalamus in
   neurotransmitter production and the pituitary in hormone regulation and
   metabolic balance. Mechanistic insights reveal adenosine signaling
   supports thermoregulation via vasodilation, while adrenaline signaling
   maintains calcium homeostasis and muscle function. Machine learning
   identified biomarkers such as 2-(Methylthio) phenol and
   Deca-2,4,6-triynedioic acid, linked to glutathione and nicotinamide
   metabolism pathways, offering strategies to mitigate HS. These findings
   reveal the specific mechanisms by which heat stress affects the
   neuroendocrine system in poultry, highlighting the interconnections
   between neuroendocrine axes and organ damage. The molecular signatures
   and pathway interactions identified in this study propose testable
   targets for dietary intervention trials, particularly through
   modulation of nicotinamide salvage pathways and adenosine receptor
   signaling, pending validation in controlled feeding studies.

Materials and methods

Chick materials and sampling

   All chickens were obtained from the fast-growing white-feathered pure
   line, produced by Xinguang Agricultural and Animal Industrials Co.,
   Ltd. (Mile, China). The 86 healthy broilers at 42 days, raised under
   standard environmental and nutritional conditions (NY/T 33-2004), were
   randomly divided into a HS group (41 chicks) and a control group (45
   chicks). While, the HS was induced by raising the lab temperature from
   24℃ to 38 ℃ within 0.5 h and continued to be treated for 3 h [[165]88].
   The temperature monitoring showed at Additional file 15: Tables S10.
   After the treatment, blood samples from the wing veins of all
   individuals were collected (serum was isolated for metabolomics
   sequencing). In this study, all experimental subjects were anesthetized
   via intramuscular injection of ketamine hydrochloride (50 mg/mL) at 60
   mg/kg body weight. Surgical anesthesia depth was confirmed by loss of
   palpebral reflex and absence of corneal reflexes. Following achievement
   of unconsciousness, broilers were euthanized by cervical dislocation.
   Following the standard of the Poultry Performance Terminology and
   Poultry Determination (China, NY/T 823–2020), pituitary and
   hypothalamus tissues were obtained for subsequent RNA extraction.

Metabolite extraction

   The serum samples were pre-processed before metabolomics sequencing
   (Additional file 16: Table S11). Here, we adapted non-targeted
   metabolomics sequencing using an ultra-performance liquid
   chromatography-tandem mass spectrometry system. Separation was achieved
   with a C18 column (1.8 μm, 2.1 mm × 100 mm). The mobile phases were
   0.1% formic acid in water (A) and 0.1% formic acid in acetonitrile (B),
   with a flow rate of 0.35 mL/min, column temperature at 45 °C, and
   injection volume of 2 µL. Samples were kept at 4 °C during analysis.
   After separation, the samples were analyzed using high-resolution ion
   mobility mass spectrometry (parameters listed in Additional file 17:
   Table S12). Data were collected in full-scan mode (m/z 50–1000) with
   MSE, alternating between low-energy (CE 4 eV) and high-energy (CE 20–45
   eV) scans for ion fragmentation. Argon (99.99%) served as the collision
   gas. The quality control samples were prepared by mixing equal volumes
   of extraction solutions from all individual samples. The correlation of
   the quality control samples was 0.99 (Additional file 18: Figure S6),
   indicating stable instrument operation and reliable results. To ensure
   data reliability, four internal standards were employed: [²H₄]-succinic
   acid (succinic acid-D4), L-2-chlorophenylalanine, D-luciferin free
   acid, and [²H₄]-cholic acid (cholic acid-D4). These standards were
   selected to cover diverse chemical classes (organic acids, amino acid
   analogues, and bile acids), thereby monitoring extraction efficiency
   and instrumental performance across metabolite categories. The relative
   standard deviations (RSDs) of all internal standards were below 30%,
   indicating stable instrument response. Additionally, the retention time
   shifts of these standards across all samples were confined to ± 18 s,
   further demonstrating the robustness of the chromatographic method.
   These metrics collectively confirm the stability of the analytical
   platform and the reliability of the detection process. While the raw
   LC-MS data were processed using Progenesis QI V2.3 for baseline
   filtering, peak detection, alignment, and normalization (key
   parameters: 5 ppm precursor tolerance, 10 ppm product tolerance, 5% ion
   threshold). Compounds were identified using databases (HMDB, LipidMaps,
   Metlin, etc.) based on charge-mass ratio, fragments, and isotopes.
   Peaks with > 50% missing values were removed, missing data were imputed
   with half the minimum value, and positive/negative ion data were merged
   into a single matrix.

Identification of signature metabolites accumulated in HS and control cohort

   Using A two-tailed Student’s t-test to determine the significance of
   metabolites’ differences between groups. The Orthogonal Partial Least
   Squares Discriminant Analysis (OPLS-DA) was performed to rank each
   metabolites’ contribution (variable importance in projection, VIP) to
   phenotypic discrimination. Differentially accumulated metabolites
   (DAMs) were identified by comparing metabolite abundance between HS
   (n = 41) and control (n = 45) groups, with FC calculated as HS/Control.
   The screening criteria such as p-value < 0.05,|log[2](FC)| ≥ 0.26 were
   applied.

   Random Forest is an efficient feature selection technique and has shown
   promising applications in metabolome data analysis [[166]89]. We used
   the RF function of caret package in R (v4.3.2). The Random Forest model
   was implemented with 100 decision trees. The feature selection process
   employed all available metabolomic data to ensure a comprehensive
   assessment of potential candidates. Following feature selection, the
   dataset was partitioned into training (65%) and validation (35%)
   subsets using stratified random sampling to maintain balanced class
   distributions. The predictive model was then constructed using the
   selected features, and to enhance the model’s robustness, we employed
   10-fold cross-validation during the training phase. This process
   involved dividing the training data into ten subsets, training the
   model on nine of these subsets, and validating it on the remaining
   subset. This procedure was repeated ten times, ensuring that each data
   point was used for both training and validation. Feature importance was
   quantified using permutation importance (z-scores) derived from 1000
   Monte Carlo permutations, with significance assessed at p-value < 0.05
   based on null distributions generated through permutation. Model
   performance was rigorously evaluated on the validation set, with
   classification performance quantitatively assessed by calculating the
   area under the receiver operating characteristic curve (AUC). The
   selection was guided by the significance of these metabolites in
   classification, accounting for variations across metabolome samples.
   Next, we deduced the promising metabolites that overlapped among DAMs,
   metabolites with VIP ≥ 1 and the featured metabolites.

mRNA sequencing and differentially expressed genes

   Total RNA extraction, transcriptome libraries construct, and paired end
   sequencing were consistent with previous research [[167]90]. Total RNA
   was extracted using TRIzol reagent according to the instructions, the
   purity and quantification of RNA were determined using a NanoDrop 2000
   spectrophotometer (Thermo Scientific, USA), and the RNA integrity was
   evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies,
   Santa Clara, CA, USA) (Additional file 19: Table S13). Transcriptome
   libraries were constructed using the VAHTS Universal V6 RNA-Seq Library
   Prep Kit according to the instructions. The libraries were sequenced
   using the Illumina NovaSeq 6000 sequencing platform, and 150 bp
   paired-end reads were generated. Approximately 50 million raw reads
   were obtained for each sample (Additional file 19: Table S13). Raw
   reads of fastq format were firstly processed using Fastp [[168]91] and
   the low quality reads were removed to obtain the clean reads. Then
   about 47 million clean reads for each sample were retained for
   subsequent analyses. The clean reads were mapped to the reference
   genome (GRCg6a) using HISAT2 (v2.1.0) [[169]92]. FPKM [[170]93] of each
   gene was calculated and the read counts of each gene were obtained by
   HTSeq-count (v0.11.2). Genes with zero expression across all samples
   were excluded from subsequent analysis. Differentially expressed genes
   (DEGs) between HS and control groups were identified using the DESeq2
   package in R (v4.3.2), with gene expression differences quantified by
   fold change (FC). The analysis compared HS samples directly against
   controls, using a p-adjust threshold of < 0.05 (Benjamini-Hochberg
   method) and an absolute|log2(fold change)| >0.59 to identify
   differentially expressed genes (DEGs). Then, we brined the STRING
   website ([171]https://cn.string-db.org) and Cytoscape software (v
   3.9.1) to construct the DEGs interaction network.

Enrichment analysis

   Metabolite Enrichment Analysis was performed using the MetaboAnalyst
   5.0 web-based platform ([172]https://www.metaboanalyst.ca). Metabolites
   were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG)
   metabolic pathway database. Gene Functional Enrichment Analysis was
   conducted using the R package clusterProfiler (v4.12.1). Ensembl gene
   IDs were converted to NCBI Entrez Gene IDs via the Omicshare platform
   ([173]https://www.omicshare.com) to ensure compatibility with
   annotation resources, prior to enrichment analysis with the Gallus
   database org.Gg.eg.db. Differentially expressed genes were annotated
   against Gene Ontology (GO) biological processes and KEGG pathways. GO
   term and KEGG pathways with p-value < 0.05 were considered
   significantly enriched, adjusted for multiple testing using the
   Benjamini-Hochberg (BH) method.

Supplementary Information

   [174]12864_2025_11881_MOESM1_ESM.xlsx^ (10.6KB, xlsx)

   Additional file 1: Table S1. The metabolites categories
   [175]12864_2025_11881_MOESM2_ESM.pdf^ (950.8KB, pdf)

   Additional file 2: Figure S1. OPLS-DA permutation test plot
   [176]12864_2025_11881_MOESM3_ESM.xlsx^ (26.7KB, xlsx)

   Additional file 3: Table S2. The details of 149 DAMs
   [177]12864_2025_11881_MOESM4_ESM.pdf^ (459.5KB, pdf)

   Additional file 4: Figure S2. Characteristic distribution of the Top 50
   DAMs
   [178]12864_2025_11881_MOESM5_ESM.xlsx^ (9.7KB, xlsx)

   Additional file 5: Table S3. The class proportion of DAMs
   [179]12864_2025_11881_MOESM6_ESM.xlsx^ (11KB, xlsx)

   Additional file 6: Table S4. The enrichment details of DAMs
   [180]12864_2025_11881_MOESM7_ESM.pdf^ (138.5KB, pdf)

   Additional file 7: Figure S3. Volcano plot of the differentially
   expression genes in hypothalamus (A) and pituitary (B)
   [181]12864_2025_11881_MOESM8_ESM.xlsx^ (22KB, xlsx)

   Additional file 8: Table S5. The details of 160 hypothalamic DEGs
   [182]12864_2025_11881_MOESM9_ESM.xlsx^ (24.7KB, xlsx)

   Additional file 9: Table S6. The details of 198 pituitary DEGs
   [183]12864_2025_11881_MOESM10_ESM.pdf^ (1.2MB, pdf)

   Additional file 10: Figure S4. Interaction network of differentially
   expressed genes
   [184]12864_2025_11881_MOESM11_ESM.xlsx^ (24.9KB, xlsx)

   Additional file 11: Table S7. The details of PPI network topology
   [185]12864_2025_11881_MOESM12_ESM.pdf^ (949.9KB, pdf)

   Additional file 12: Figure S5. Bar plot of GO terms (Top 60) enriched
   for hypothalamus (A) and pituitary (B) differentially expressed genes
   [186]12864_2025_11881_MOESM13_ESM.xlsx^ (14KB, xlsx)

   Additional file 13: Table S8. The enrichment details of hypothalamic
   DEGs
   [187]12864_2025_11881_MOESM14_ESM.xlsx^ (14KB, xlsx)

   Additional file 14: Table S9. The enrichment details of pituitary DEGs.
   [188]12864_2025_11881_MOESM15_ESM.xlsx^ (9.7KB, xlsx)

   Additional file 15: Table S10. The temperature monitoring records of HS
   [189]12864_2025_11881_MOESM16_ESM.xlsx^ (9.9KB, xlsx)

   Additional file 16: Table S11. Pre-processing of serum samples before
   metabolome sequence
   [190]12864_2025_11881_MOESM17_ESM.xlsx^ (10.3KB, xlsx)

   Additional file 17: Table S12. The gradient program of
   ultra-performance liquid chromatography -tandem mass spectrometry
   system (comprised an ultra-performance liquid chromatography system and
   a high-resolution mass spectrometer)
   [191]12864_2025_11881_MOESM18_ESM.pdf^ (5.7MB, pdf)

   Additional file 18: Figure S6. Correlation matrix diagram of quality
   control samples
   [192]12864_2025_11881_MOESM19_ESM.xlsx^ (28.4KB, xlsx)

   Additional file 19: Table S13. RNA-seq Data Production and Quality
   Assessment Metrics

Acknowledgements