Comprehensive Profiling of ceRNA (circRNA-miRNA-mRNA) Networks in Hypothalamic-Pituitary-Mammary Gland Axis of Dairy Cows under Heat Stress

Heat stress (HS) is directly correlated with mammary gland dysfunction and the hypothalamic-pituitary-mammary gland (HPM) axis is involved in regulating stress responses and lactation in dairy cows. Circular RNAs (circRNAs) play major roles in regulating transcription and post-transcription but their expression in the HPM axis of dairy cows under HS is still unclear. In the present study, we performed RNA sequencing to identify diferentially expressed (DE) circRNAs, DE microRNAs(miRNAs) and DEmRNAs, and performed bioinformatics analysis on those in HPM axis-related tissues of heat-stressed and normal cows. A total of 1680, 1112 and 521 DEcircRNAs, 120, 493 and 108 DEmiRNAs, 274, 6475 and 3134 DEmRNAs were identified in the hypothalamic, pituitary, and mammary gland tissues, respectively. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses indicated that the MAPK signaling pathway is potentially a key pathway. Competitive endogenous RNA (ceRNA) networks related to HS response and lactation regulation were established in three tissues. In conclusion, our results indicate that HS induces differential circRNA expression profiles in HPM axis-related tissues, and the predicted ceRNA network provides a molecular basis for regulating the stress response and lactation regulation in heat-stressed dairy cows.


Introduction
Heat stress syndrome is a systemic indication that includes non-specific defense responses and specific disorders of the body systems when cows are stimulated by excessive temperatures that exceed their thermoregulatory ability. Heat stress (HS) affects various physiological and biochemical indexes and reduces feed intake, immunity, and reproductive performance of dairy cows [1]. Importantly, HS inhibits the proliferation of bovine mammary epithelial cells (BMECs), aggravates cellular oxidative stress, and promotes BMECs apoptosis [2] and mammary gland autophagy [3], thereby reducing lactation performance [4]. Regulation of the hypothalamic-pituitary endocrine axis is crucial in response to stimulation of the organism by the external environment [5].
The hypothalamus controls the pituitary endocrine system, and they are thought of as a structural and functional unit, regulating important physiological processes such as thermoregulation and feeding [6,7]. Prolactin (PRL) secreted by the pituitary is essential for maintaining lactation [8]. Under stress, glucocorticoid levels in the blood can rapidly increase under the influence of the hypothalamic-pituitary-adrenal cortex (HPA) axis to enhance the body's resistance to adverse environment [9]. Excitation of hypothalamicpituitary-thyroid (HPT) axis can regulate the oxygen consumption rate of most tissues through the thyroid hormone, altering the thermogenic effect of the body [10]. However, information about HS effects on hormonal changes related to the hypothalamic-

Differences in Endocrine Hormones, Antioxidant Enzymes and Heat Shock Proteins (HSPs) between NHS and HS Groups
HS significantly increased the levels of HPA axis-related hormones CRH, ACTH and COR ( Figure 1A, p < 0.01), in the serum of dairy cows. HPT axis-related hormones TRH and TSH significantly increased (p < 0.01), while those of T3 and T4 significantly decreased ( Figure 1B, p < 0.01) under HS. Compared to the NHS group, PRF, PIF and PRL significantly increased ( Figure 1C, p < 0.01) in the HS group. Meanwhile, lactation related hormones GH and IGF-1 significantly increased (p < 0.01), while INS significantly decreased (p < 0.01) in the serum of heat-stressed dairy cows ( Figure 1D). Moreover, SOD, GSH-Px and LDH activity in the serum of dairy cows significantly decreased (p < 0.01), while MDA content significantly increased (p < 0.01) under HS condition ( Figure 1E). Lastly, the levels of HS marker proteins, HSP70 and HSP90, significantly increased ( Figure 1F, p < 0.01) in the serum of heat-stressed dairy cows.

Histomorphological Observations of HPM-Axis Tissues and Neurotransmitters in Hypothalamus of NHS and HS Groups
As shown in Figure 2A, the hypothalamic histology was not significantly different between the NHS and HS groups. However, there are histological differences between the NHS group and the HS group in the pituitary and mammary glands. Among them, the percentage of basophils in the pituitary tissue significantly increased under the HS condition ( Figure 2B, p < 0.01). Mammary gland acini in the NHS group were round, but those in the HS group were significantly shrunken with a thicker acinar wall and their round shape was lost ( Figure 2C). Subsequently, we detected neurotransmitters in the hypothalamus and found that DA significantly decreased while β-EP levels significantly increased under HS ( Figure 2F,G p < 0.01). However, GABA, NE and Ach had no significant effect on the hypothalamic tissues of heat-stressed dairy cows ( Figure 2C  tion ( Figure 2B, p < 0.01). Mammary gland acini in the NHS group were round, but those in the HS group were significantly shrunken with a thicker acinar wall and their round shape was lost ( Figure 2C). Subsequently, we detected neurotransmitters in the hypothalamus and found that DA significantly decreased while β-EP levels significantly increased under HS ( Figure 2F,G p < 0.01). However, GABA, NE and Ach had no significant effect on the hypothalamic tissues of heat-stressed dairy cows ( Figure 2C-E, p > 0.05).

Identification and Characterization of circRNAs in HPM Axis Tissues of Dairy Cow under Heat Stress
Expression of circRNAs in the NHS and HS groups were detected and high-quality clean reads obtained to identify circRNAs (Supplementary Table S1). In the 18 sequencing libraries, 81.08%-85.63% of the reads were mapped to the reference genome (Supplementary  Table S2). Approximately 68.97%, 74.58% and 76.59% of the reads in the hypothalamus, pituitary and mammary gland, respectively, were compared to exons ( Figure 3A). Based on the filtering criteria, we conducted length distribution statistics and annotated circRNAs. The percentage and frequence of circRNAs with length between 300-400 nt (11.91%) and greater than 3000 nt (12.96%) are relatively high ( Figure 3B), and the main annotation type was "annot_exons" ( Figure 3C). A violin diagram is generally used for the visualization of gene abundance and expression and can show data density at any location. As shown in Figure 3D, the circRNA expression abundance and data density of each experimental group are appropriate. Therefore, we conducted subsequent analysis.
(11.91%) and greater than 3000 nt (12.96%) are relatively high ( Figure 3B), and the main annotation type was "annot_exons" ( Figure 3C). A violin diagram is generally used for the visualization of gene abundance and expression and can show data density at any location. As shown in Figure 3D, the circRNA expression abundance and data density of each experimental group are appropriate. Therefore, we conducted subsequent analysis.  Furthermore, a total of 1680, 1112 and 521 DE circRNAs were identified in hypothalamus, pituitary and mammary gland, respectively. We found that 1111, 2461 and 250 circRNAs were up-regulated, and 569, 1112 and 271 circRNAs were downregulated in the three tissues ( Figure 4A-C). Moreover, 33 circRNAs were shared in these tissues ( Figure 4D). To verify the RNA-seq results, we randomly selected three circRNAs from HPM axis-related tissues and detected their expression using RT-qPCR. The results were highly consistent with those of RNA-seq ( Figure 4E). circRNAs were up-regulated, and 569, 1112 and 271 circRNAs were downregulated in the three tissues ( Figure 4A-C). Moreover, 33 circRNAs were shared in these tissues ( Figure  4D). To verify the RNA-seq results, we randomly selected three circRNAs from HPM axisrelated tissues and detected their expression using RT-qPCR. The results were highly consistent with those of RNA-seq ( Figure 4E).

Functional Analysis of the Source Genes of DEcircRNAs
To determine the biological processes associated with the source genes of the DE circRNAs in HPM axis tissues, we conducted GO and KEGG pathways enrichment analysis. In the hypothalamus, most source genes were significantly (p < 0.01) enriched in protein binding (GO:0005515), ATP binding (GO:0005524), regulation of GTPase activity (GO:0043087), and cellular response to DNA damage stimulus (GO:0006974) ( Figure 5A). The KEGG pathways were in glutamatergic synapse, MAPK signaling pathway, cAMP signaling pathway, GABAergic synapse, thyroid hormone synthesis, and insulin secretion in the hypothalamus ( Figure 5B, p < 0.05). In the pituitary, cellular component organization or biogenesis (GO:0071840), catalytic activity, acting on a protein (GO:0140096), and cellular protein modification process (GO:0006464) were significantly enriched ( Figure  5C, p < 0.01). The KEGG pathways were mainly enriched in the MAPK signaling pathway, hormone synthesis, secretion, and actin, autophagy ( Figure 5D, p < 0.05). The GO terms

Functional Analysis of the Source Genes of DEcircRNAs
To determine the biological processes associated with the source genes of the DE circRNAs in HPM axis tissues, we conducted GO and KEGG pathways enrichment analysis. In the hypothalamus, most source genes were significantly (p < 0.01) enriched in protein binding (GO:0005515), ATP binding (GO:0005524), regulation of GTPase activity (GO:0043087), and cellular response to DNA damage stimulus (GO:0006974) ( Figure 5A). The KEGG pathways were in glutamatergic synapse, MAPK signaling pathway, cAMP signaling pathway, GABAergic synapse, thyroid hormone synthesis, and insulin secretion in the hypothalamus ( Figure 5B, p < 0.05). In the pituitary, cellular component organization or biogenesis (GO:0071840), catalytic activity, acting on a protein (GO:0140096), and cellular protein modification process (GO:0006464) were significantly enriched ( Figure 5C, p < 0.01). The KEGG pathways were mainly enriched in the MAPK signaling pathway, hormone synthesis, secretion, and actin, autophagy ( Figure 5D, p < 0.05). The GO terms were enriched in cellular protein modification process (GO:0006464), transferase activity (GO:0016740), and regulation of GTPase activity (GO:0043087) in the mammary gland ( Figure 5E, p < 0.01). Additionally, more source genes were enriched in tight junction, MAPK signaling pathway, and hormone synthesis and secretion pathway ( Figure 5F, p < 0.05).

Identification of DE miRNAs in HPM Axis Tissues of Dairy Cows under Heat Stress
To investigate the miRNA expression profiles in the HPM axis tissues of dairy cows under HS, 18 small RNA libraries were constructed and sequenced. The results identified 120, 493 and 108 DEmiRNAs in the hypothalamus, pituitary gland, and mammary gland, respectively, of which 51, 227, and 37 miRNAs were upregulated and 69, 266, and 71 miRNAs were downregulated ( Figure 6A). Furthermore, we verified the expression of miRNAs in the hypothalamus (bta-miR-211, bta-miR-218, and bta-miR-144), pituitary gland (bta-miR-152, bta-miR-148a, and bta-miR-379), and mammary gland (bta-miR-365-3p, bta-miR-196a, and bta-miR-1388-3p) using RT-qPCR. These results were highly consistent with those of RNA-seq ( Figure 6B). Subsequently, we analyzed the function of target genes of DEmiRNAs. As shown in Supplementary Figure S1, 2499, 2649, and 2454 GO terms were enriched in the hypothalamus, pituitary, and mammary gland, respectively. Most target genes of DEmiRNAs were significantly (p < 0.05) enriched in the MAPK signaling pathway and hormone synthesis and secretion in the HPM-axis ( Figure 6C-E).

Identification of DE miRNAs in HPM Axis Tissues of Dairy Cows under Heat Stress
To investigate the miRNA expression profiles in the HPM axis tissues of dairy cows under HS, 18 small RNA libraries were constructed and sequenced. The results identified 120, 493 and 108 DEmiRNAs in the hypothalamus, pituitary gland, and mammary gland, respectively, of which 51, 227, and 37 miRNAs were upregulated and 69, 266, and 71 miR-NAs were downregulated ( Figure 6A). Furthermore, we verified the expression of miR-NAs in the hypothalamus (bta-miR-211, bta-miR-218, and bta-miR-144), pituitary gland (bta-miR-152, bta-miR-148a, and bta-miR-379), and mammary gland (bta-miR-365-3p, bta-miR-196a, and bta-miR-1388-3p) using RT-qPCR. These results were highly consistent with those of RNA-seq ( Figure 6B). Subsequently, we analyzed the function of target genes of DEmiRNAs. As shown in Supplementary Figure S1, 2499, 2649, and 2454 GO terms were enriched in the hypothalamus, pituitary, and mammary gland, respectively. Most target genes of DEmiRNAs were significantly (p < 0.05) enriched in the MAPK signaling pathway and hormone synthesis and secretion in the HPM-axis ( Figure 6C-E).

Discussion
The decline in lactation performance of dairy cows caused by HS should not be underestimated. The hypothalamic-pituitary-endocrine axis is crucial for the stress response and lactation, and the expression of related pathways and key genes to relieve HS is also important information. Consequently, in the present study, we determined the endocrine hormones of the HPM axis responsible, and established a hypothalamic, pituitary, and mammary gland ceRNA (circRNA-miRNA-mRNA) network using RNA-seq data from dairy cows under HS.
Under high temperatures, cows attempt to promote heat dissipation or reduce metabolic heat production to maintain physiological homeostasis. Excessive accumulation of reactive oxygen species (ROS) in cells caused by HS induces antioxidant dysfunction [19]. The scavenging capacity of ROS is mainly reflected in the activity of antioxidant enzymes [20]. In the present study, MDA levels in the serum of dairy cows increased significantly, while that of SOD, GSH-Px and LDH decreased significantly. Therefore, our results confirm that dairy cows experience oxidative stress under HS.
Hormones are highly effective bioactive substances secreted by the endocrine glands or scattered endocrine cells and play a regulatory role through tissue fluid or blood transmission. HS stimulates the HPA axis and promotes the secretion of ACTH, which enhances the activity of the adrenal gland and promotes the secretion of adrenaline and COR to resist the stress response [21]. The secretion of TSH from the pituitary is regulated by TRH from the hypothalamus [22]. T3 and T4 increase oxygen consumption and heat production in the body by promoting the decomposition and oxidation of sugar, fat, and protein [23,24]. The present study found that HS reduces the secretion of T3 and T4 but TSH and TRH secretion, and the expression of TRHR and TSHR significantly increased in heat-stressed cows. This is related to the negative feedback effect of T3 and T4 concentrations in the blood on the release of TSH and TRH. A variety of endocrine hormones regulate the growth, development, and lactation of the mammary gland [8]. The synthesis and secretion of PRL, a polypeptide hormone secreted by the anterior pituitary, is regulated by hypothalamic neurosecretory cells, providing PRF and PIF to the pituitary [24]. In this study, serum PRF and PIF levels in heat-stressed cows increased, which was related to the excitation of the hypothalamic-pituitary axis caused by the stress response. Interestingly, serum GH and IGF-1 levels increased significantly but the expression of GHR and IGF-1 decreased under thermal stimulation, which may affect the growth and function of BMECs. DA released from the hypothalamus binds to the D2 receptor (D2R) in the anterior pituitary, which can inhibit PRL secretion [25]. Moreover, the release of methionine enkephalin from the hypothalamus can also promote the secretion of PRL [26]. Our results showed that HS decreased the secretion of DA in the hypothalamus, but significantly increased the secretion of β-EP. Therefore, the hypothalamic-pituitary axis maintains the normal secretion of PRL in the HS environment through positive and negative feedback.
HSPs are highly conserved stress-response proteins promoting cell survival under various adverse conditions [27]. In our study, 42 DE HSPs were identified in HPM axisrelated tissues by RNA-seq, and the expression levels of HSP40, HSP70, and HSP90 families were different under HS. As molecular chaperones, HSPs are involved in protein folding, transport, and stability [28]. It is worth noting that HSPs with high molecular weight are ATP-dependent proteins [29]. Moreover, HSPs can be used as signaling molecules to inhibit apoptosis and maintain cell homeostasis after adverse stimulation [29,30]. The results of our GO analysis showed significantly enriched terms related to ATPase activation, ATP binding, and cellular response to stress. Consequently, we speculate that HSPs are important regulators in HPM axis-related tissues of heat-stressed dairy cows.
The mechanisms of ncRNAs in lactation and stress response have been studied; they play an important role in regulating RNA transcription, protein translation, and protein interaction [31][32][33]. In our study, a total of 1680, 1112, and 521 DE circRNAs were identified from the in the hypothalamic, pituitary, and mammary gland tissues, respectively. GO analysis identified biological processes related to the stress response of DEcircRNA source genes in three tissues, including cellular response to stress, regulation of cellular response to stress, stress-activated MAPK cascade, and stress-activated protein kinase signaling cascade. Subsequently, KEGG pathway enrichment analysis showed that the DEcircRNAs have roles mainly in MAPK signaling pathway, tight junction, GnRH signaling pathway, and hormone synthesis, secretion and action, which have been shown to be closely related to the stress response and physiological process of lactation. It is worth noting that the MAPK family is a group of serine/threonine kinases, which can be activated by cytokines, neurotransmitters, hormones, and cell stress. The MAPK pathways regulate a variety of cellular activities, including proliferation, differentiation, survival, and death [34]. Previous studies have shown that the MAPK signaling pathway is closely related to HSP70. MAPK can activate and up regulate the expression of heat shock factor 1 (HSF1), and promote the combination of HSF1 and HSP70 heat shock elements, thereby increasing the synthesis of HSP70 [35]. Accumulating evidence has indicated that HS stimulation induces MAPK activation, thereby affecting cell proliferation and apoptosis [36,37]. Additionally, external temperature stimulation leads to loosening of tight junctions in intestinal epithelial cells, causing gastrointestinal dysfunction [38]. A study has shown that HS affects the expression of tight junction proteins (such as ZO-1, claudin-1, and claudin-4) by regulating the ERK1/2-MAPK signaling pathway [39]. The localization of these proteins is maintained by the activation of HSPs by MAPK [38]. Therefore, the MAPK signaling pathway, as a key pathway regulating the HPM axis in heat-stressed dairy cows, should be investigated in subsequent studies.
It is widely known that circRNAs have MREs that can act as ceRNAs to regulate miR-NAs and mRNAs to exert their biological functions [40]. For example, a previous study has shown that circ09863 alleviates the inhibition of FASN expression by combining miR-27a-3p, and regulates TAG synthesis and fatty acid composition in BMECs [41]. Circ08409 relieves the inhibitory effect of miR-133a on TGFB2 expression by binding it and subsequently modulating cell proliferation and apoptosis in BMECs [42]. In the current study, some ncRNAs and key genes in the HPM axis were identified under HS in the ceRNA regulatory network, such as novel_circ_000126-bta-miR-211-HSPH1 and novel_circ_002325-bta-miR-218/bta-miR-33a-GHRHR in the hypothalamus. In the pituitary, the mainly ceRNA network was found in novel_circ_013439-miR-9-z-PRL, novel_circ_006323-bta-miR-124a-GH1, and novel_circ_020252-miR-7864-z -IGF1. Moreover, the subnetworks of novel_circ_011229bta-miR-335-PRLR and novel_circ_016821-bta-miR-129-IGF1 were established in mammary gland. The expression of miR-23a, miR-27a, miR-27b, miR-103 and miR-200a can be upregulated by PRL [43]. One study has shown that miR-15a decreases BMEC viability and lactation, and regulates GHR expression [44]. Additionally, miR-29a-3p inhibits cell proliferation, PRL and GH expression of prolactinoma cells, and promotes apoptosis by inhibiting IGF-1 [45]. Therefore, we hypothesized that dairy cows may respond to the thermal environment by regulating the expression of various RNAs and key pathways in the HPM axis, affecting hormone secretion and lactation performance. The specific mechanisms should be investigated in future studies.
In conclusion, this study indicates that HS may affect the stress response and lactation physiology of dairy cows by regulating the secretion of endocrine-related hormones and release of neurotransmitters, and the expression of circRNAs, miRNAs, and mRNAs in HPM axis-related tissues. The differential gene expression is closely related to the MAPK signaling pathway, which may be a potential regulatory pathway. Moreover, the established ceRNA (circRNA-miRNA-mRNA) regulatory networks provide a theoretical basis for the molecular mechanism of the stress response and lactation regulation of ncRNAs in HPM axis-related tissues under HS.

Animals and Sample Collection
The temperature and humidity index (THI) of the cowshed, and the rectal temperature, and respiratory rate of Holstein cows were continuously monitored under the same management level in summer (August, heat stress group (HS), n = 20) and winter (December, Non heat stress group (NHS), n = 20) in Huai'an Jielong Animal Husbandry Co., Ltd. (Huai'an, China). Information about parity, lactation days and milk yield of the cows in NHS group and HS group are shown in Table 1. There was no significant difference in parity and lactation days between NHS and HS groups (p > 0.05), but the milk yield decreased significantly in HS groups (p < 0.05). The data of THI, rectal temperature, and respiratory rate are shown in Table 2. The temperature (T, • C) and relative humidity (RH, %) of the cowshed were measured daily at 10:00 and 20:00 for a month. Using electronic thermometers at six fixed points in the cowshed. THI criteria were derived from the NRC (1971), and the THI was calculated according to the following equation: THI = (1.8 × T + 32) − (0.55 − 0.55 × RH) × (1.8 × T − 26), where T is the temperature ( • C) and RH is the relative humidity (%) [46]. Rectal temperature and respiratory rate were measured using an animal thermometer and stopwatch, respectively, at 10:00 and 20:00 daily during the experiment. Blood samples (10 mL) from the NHS and HS groups of Holstein cows were obtained from the caudal vein before morning feedings during the trial period and collected into tubes until coagulation. The collected samples were centrifuged at 3000 rpm for 10 min at 4 • C and stored at −20 • C for subsequent serum endocrine hormone and antioxidant index measurements. Blood sample collection was approved and implemented by Huai'an Jielong Animal Husbandry Co., Ltd. (Huai'an, China).

Histological Evaluation
To prepare slides for histological evaluation, fresh hypothalamus, pituitary, and mammary gland samples were fixed in 4% paraformaldehyde and tissue sections were subsequently dehydrated through steps of graded alcohol, cleared in xylene, and embedded in paraffin blocks. Sections of 6 µm thickness were cut from each block. The deparaffinized sections were stained with hematoxylin and eosin. The tissue sections were then observed using a microscope (Olympus Corporation, Tokyo, Japan).

Detection of Neurotransmitters
High performance liquid chromatography-tandem mass spectrometry (HPLC-MS/MS) was used to detect γ-Aminobutyric acid (GABA), acetylcholine (ACh), dopamine (DA), and norepinephrine (NE) content. β-Endorphin (β-EP) ELISA Assay Kit (Nanjing Ruiyuan Biotechnology Co., Ltd., Nanjing, China) was used to determine it according to the manufacturer's protocol. The absorbance was read at 450 nm using a microplate reader and β-EP content was calculated based on absorbance.

RNA Extraction and Sequencing
RNA-seq was performed using hypothalamus, pituitary, and mammary gland tissues sample. First, total RNA was extracted from tissue samples. Ribosomal RNAs (rRNAs) were depleted to construct a total RNA-seq library. RNAs were fragmented into short fragments using fragmentation buffer and reverse-transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP (dUTP instead of dTTP) and buffer. The cDNA fragments were purified using the QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end-repaired, poly(A) added, and ligated to Illumina sequencing adapters. Uracil-N-Glycosylase (UNG) was used to digest the second-strand cDNA. The digested products were size-selected using agarose gel electrophoresis, amplified by PCR, and sequenced using Illumina Nova-Seq 6000 by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China).
After total RNA was extracted, Small RNA molecules in a size range of 18-30 nt were enriched using polyacrylamide gel electrophoresis (PAGE). Then, the 3 -adapters were added and the 36-44 nt RNAs were enriched. The 5 -adapters were then ligated to the RNAs and the ligation products were reverse-transcribed using PCR amplification. The 140-160 bp size PCR products were enriched to generate a small RNA library and sequenced using Illumina Nova-Seq 6000.

Identification of circRNAs, miRNAs and mRNAs
To obtain high quality clean reads, raw reads were filtered using fastp (v0.18.0). Short reads alignment tool Bowtie2 (v.2.2.8) was used for mapping reads to ribosome RNA (rRNA) database. The rRNA removed reads from each sample were then mapped to the cow reference genome by HISAT2 (v2.1.1). After alignment with reference genome, the reads that could be mapped to the genomes were discarded, and the unmapped reads were collected for circRNA identification. From both ends of the unmapped reads, 20 mers were extracted and aligned to the reference genome to identify unique anchor positions within the splice site. Anchor reads that aligned in the reversed orientation (head-to tail) indicated circRNA splicing and were subjected to find_circ to identify circRNAs. The identified circRNAs were subjected to statistical analyses for type, and chromosome and length distribution.
All clean tags were searched against the miRBase database to identify known miRNAs. The unannotated tags were aligned with the cow reference genome according to their genome positions and hairpin structures predicted by the Mireap_v0.2 software to identify novel miRNA candidates.
Transcript reconstruction was carried out with the Stringtie software (v1.3.4), which, together with HISAT2, allow biologists to identify new genes and new splice variants of known genes.

Differentially Expressed of circRNAs, miRNAs and mRNAs
To identify DEcircRNAs, DEmiRNAs and DEmRNAs in the HPM axis tissues between the NHS and HS groups, the edgeR package (http://www.rproject.org/ (accessed on 15 September 2021)) was used. We identified circRNAs and miRNAs with fold change (FC) ≥ 2 and p-value < 0.05 in a comparison between groups as significantly DEcircRNAs and DEmiRNAs, respectively. Significantly DEmRNAs were identified with a FC ≥ 2 and false discovery rate (FDR) < 0.05.

Function Enrichment Analysis
To understand the potential roles of DEcircRNAs, DEmiRNAs and DEmRNAs between NHS and HS groups, Gene ontology (GO) enrichment analysis and the Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses were performed to investigate their biological function using DAVID (http://david.abcc.ncifcrf.gov/ (accessed on 24 September 2021)).

Construction of the Competing Endogenous RNA (DEcircRNA-DEmiRNA-DEmRNA) Regulatory Network
The ceRNA network was constructed based on the ceRNA theory as follows: (1) The expression correlation between mRNA-miRNA or circRNA-miRNA was evaluated using the Spearman Rank correlation coefficient (SCC). Pairs with SCC < −0.7 were selected as negatively coexpressed circRNA-miRNA pairs or mRNA-miRNA pairs, and both mRNA and circRNA were miRNA target genes, and all RNAs were differentially expressed. (2) The expression correlation between circRNA-mRNA was evaluated using the Pearson correlation coefficient (PCC). Pairs with PCC > 0.9 were selected as coexpressed circRNA-mRNA pairs, and both mRNA and circRNA in these pairs were targeted and negatively coexpressed with a common miRNA. (3) A hypergeometric cumulative distribution function test was used to test whether the common miRNA sponges between the two genes were significant. As a result, only the gene pairs with a p-value less than 0.05 were selected.

Quantitative Real-Time Polymerase Chain Reaction (RT-qPCR) Validation
Total RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) and the density of sample RNA was quantified spectrophoto-metrically at 260/280 nm. The RNA was reverse-transcribed using Prime Script TM RT Master Mix (TaKaRa, Otsu, Japan) according to the manufacturer's protocol. The sequences of the primers are listed in (Supplementary Tables S7 and S8). All RT-qPCR runs were performed on a StepOne Plus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) using the AceQ qPCR SYBR Green Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China) in a reaction volume of 20 µL. The relative quantification of miRNAs, circRNAs, and mRNAs were normalized to that of U6 or β-actin using the 2 −∆∆Ct method.

Statistical Analysis
The statistical analyses were performed in SPSS version 20.0 (SPSS Inc., Chicago, IL, USA) and the graphs were generated using GraphPad Prism 8.4.2 (La Jolla, CA, USA). The significance of differences in physiological and biochemical parameters were determined using One-way ANOVA with Duncan's multiple range test, and the data are expressed as mean and standard error of mean (SEM). p-values < 0.05 were considered significant.