Introduction

The Qinghai-Tibet Plateau experiences severe cold, low oxygen levels, high ultraviolet radiation, and strong winds1. Animals have adapted to these harsh environmental factors through long-term natural selection, forming unique groups suited for high-altitude environments2.

Tibetan donkeys in high-altitude regions showed a stronger tolerance to hypoxia and a higher metabolism of cellulose and hemicellulose. These changes contribute to their survival in harsh habitats3. Guo et al. found higher flora diversity and richness in Tibetan donkeys than in Dezhou donkeys at low altitudes4. A study on the high-altitude adaptation of Tibetan donkeys identified the role of the EGLN1 gene in their adaptation5. Zhou et al. conducted second-generation sequencing on six donkey breeds to explore genetic variations related to high-altitude adaptation6. Luo et al. identified 16 novel amino acid variants in proteins encoded by Tibetan donkeys’ mitochondrial DNA. These variants may affect mitochondrial energy production efficiency and enhance adaptation to cold and hypoxic conditions in the Qinghai-Tibet Plateau7. In summary, these studies have revealed adaptation mechanisms in donkeys in extreme environments, focusing mainly on the intestinal flora and genes of Tibetan donkeys. However, compared with other livestock, research on the high-altitude adaptability of Tibetan donkeys is scarce8,9,10.

Metabolites are a key means of communication between the host and the producer microbes11. Metabolomics, through liquid chromatography-tandem mass spectrometry (LC-MS/MS), can achieve high-throughput, high-sensitivity, wide-coverage, and accurate metabolite detection and analysis12,13. It can reveal metabolic changes in organisms during disease development, drug action, environmental exposure, and other processes14,15,16,17. Such studies provide new insights and methods for disease diagnosis, treatment, and prevention.

This study explored the specific mechanism of high-altitude adaptation in Tibetan donkeys using serum non-targeted metabolomics analyses comparing high-altitude Tibetan donkeys and low-altitude Dezhou donkeys (Fig. 1). Our results highlight the adaptive metabolism of donkeys at high altitudes and elucidate its potential mechanism, and these screened biomarkers may provide new targets for donkey breeding.

Fig. 1
figure 1

Workflow overview of the comprehensive analysis of metabolome in Tibetan and Dezhou donkeys.

Full size image

Results

Partial least squares discrimination analysis (PLS-DA)

A total of 678 metabolites were detected, including 342 and 336 metabolites in the positive (ESI+) and negative ion modes (ESI−), respectively. PLS-DA results showed that Component 1 was 31.6%, and Component 2 was 15.2%. The RKZ, CD, and DZ groups were effectively divided into three categories in the combined ion mode, indicating significant differences in their metabolite phenotypes (Fig. 2A). The validation of the PLS-DA model showed that in the combined ion mode, the R2 and Q2 intercepts of the model reached 0.993 and 0.957, respectively. This indicated that the model fitted well, and was highly predictable and suitable for subsequent data analysis (Fig. 2B). The Model Overview and Model parameters are shown in Supplementary Document 1.

Fig. 2
figure 2

PLS-DA model score (A) and model validation (B) for the three groups of serum metabolome samples under the combined ion mode. The abscissa represents the permutation reservation of the permutation test. The ordinate represents the values of the substitution tests for R2 (the red dot) and Q2 (the blue triangle). The two dotted lines represent the R2 and Q2 regression lines.

Full size image

Identification of DMs

According to VIP ≥ 1, p < 0.05, and Fold-Change ≥ 1.218 or ≤ 0.5, 93 DMs were identified in the RKZ_VS_DZ group, including 72 up-regulated and 21 down-regulated (Fig. 3A). In the CD_VS_DZ group, 188 DMs were identified (83 in ESI + and 105 in ESI−), with 163 up-regulated and 25 down-regulated (Fig. 3B). In the RKZ_VS_CD group, 162 DMs were identified (65 in ESI + and 97 in ESI−), including 34 up-regulated and 128 down-regulated (Fig. 3C).

Fig. 3
figure 3

Volcano maps of DMs in the combined ion model. Red, blue, and gray represent upregulation, downregulation, and no change, respectively.

Full size image

Furthermore, the top 30 DMs in each of the three groups were selected for cluster analysis. Each of the three difference groups was clustered into two categories, indicating their similar functional annotations or belonging to the same metabolic pathway (Fig. 4). The expression differences between different samples were obvious.

Fig. 4
figure 4

Heatmaps of the top 30 DMs ranked by VIP values in serum with a cutoff value of VIP > 1 andp < 0.05( t-tests). Each row represents a metabolite and the column corresponds to a sample. The color intensity indicates the relative abundance of respective metabolites, with warmer colors reflecting higher levels and cooler colors indicating lower levels. The hierarchical clustering of both metabolites and samples highlights patterns and relationships among the DMs.

Full size image

KEGG enrichment analysis of DMs

KEGG pathway results showed that DMs mainly clustered in bile secretion, parathyroid hormone synthesis, secretion and action, vitamin digestion and absorption, mTOR signaling pathway, and PI3K-Akt signaling pathway in the RKZ_VS_DZ group (Fig. 5A, Table 1). Metabolites in the CD_VS_DZ group mainly concentrated in vitamin digestion and absorption, tryptophan metabolism, fluid shear stress and atherosclerosis, bile secretion, and HIF-1 signaling pathway metabolism (Fig. 5B, Table 1). Metabolites in the RKZ_VS_CD group mainly concentrated in tryptophan metabolism, vitamin digestion and absorption, parathyroid hormone synthesis, secretion and action, bile secretion, and steroid hormone biosynthesis (Fig. 5C, Table 1).

Table 1 KEGG enrichment statistics.
Full size table
Fig. 5
figure 5

Classification and KEGG pathway analysis of DMs. The A), B), and C) bubble plots depict DMs mapped to their respective metabolic pathways. The dotted lines indicate the significance cutoff with a p-value of 0.05. The size of the bubble corresponds to the number of DMs in the respective KEGG pathway19.

Full size image

Discussion

High-altitude adaptation involves not only genetic regulation but also complex physiological and biochemical responses. We used an untargeted LC-MS metabolomics method to analyze serum metabolism to understand high-altitude adaptation in Tibetan donkeys.

In the three test groups (RKZ, CD, and DZ), we found that the metabolite DCA (deoxycholic acid) was significantly enriched in bile secretion metabolic pathways. Cholic acid (BAs), produced through bile secretion, plays a crucial role as a biological regulator regulating glucose, fat, and energy metabolism20. As a principal constituent of bile, DCA has immunomodulatory and anti-inflammatory properties21,22. It facilitates fat digestion and absorption in the small intestine and also contributes to cholesterol (CHOL) metabolism and induces gallbladder contraction to regulate gallbladder function23. CHOL, a primary component of cell membranes, is crucial for membrane formation, maintenance, and signal transduction24. We found significant enrichment of DCA in the bile secretion pathway across various groups, implying that Tibetan donkeys may have enhanced lipid digestion and absorption due to modulated bile secretion. This potentially improves their energy metabolism and adaptability to high-altitude environments.

Similarly, the DMs of vitamin A, vitamin C, and FMN were primarily enriched in the vitamin digestion and absorption pathways in the three distinct elevation groups. Vitamin A has antioxidant properties that enable direct removal of reactive oxygen species, enhance antioxidant oxidase activity, and promote antioxidant defense mechanisms25. FMN, the active form of vitamin B2, is vital for energy metabolism and participates in various biochemical reactions26,27,28. Similarly, vitamin C inhibits hypoxia-induced mitochondrial damage and apoptosis in human endothelial cells29. Vitamin C potentially contributes to the maintenance of the central vascular system in the Tibetan donkey’s plateau environment. Additionally, studies have shown that vitamin C can promote vasodilation and regulate peripheral vascular responses to local cold stimulation in high-altitude populations30. Thus, vitamins A and C may support the adaptability of Tibetan donkeys in the plateau environment through their antioxidant and vascular regulation effects. Vitamin metabolism might enhance their immunity and production performance.

In the CD_VS_DZ and RKZ_VS_CD groups, DMs such as NFK, n-acetylserotonin, and melatonin were enriched in tryptophan metabolism. Tryptophan metabolism is closely linked to carbohydrate, protein, and fat metabolism, as well as oxidative stress and inflammation regulation31,32. Moreover, tryptophan regulates immunity, neuronal function, and intestinal homeostasis through kynurenine metabolism33. Previous research has shown that Tibetan donkeys have significantly higher heme levels than plain donkeys under hypoxic conditions, with heme dioxygenase content increasing with heme levels at high altitudes. Heme dioxygenase catalyzes L-tryptophan conversion to N-formyl kynurenine (NFK), initiating tryptophan catabolism34. At high altitudes, light, temperature, and oxygen levels influence animal circadian rhythms. Particularly, melatonin plays a crucial role in the regulation and promotion of sleep. Melatonin is synthesized in the pineal gland from tryptophan conversion to serotonin35. It plays essential roles in circadian rhythm regulation, sleep cycles, and immune function36. Typically, less than 5% of tryptophan is utilized for serotonin and melatonin synthesis, and the remaining 95% is metabolized via the kynurenine pathway in the liver37,38. High-altitude living can elevate pro-inflammatory cytokines and glucocorticoids, enhancing tryptophan metabolism through the kynurenine pathway39,40. Consequently, animals at high altitudes may have increased energy production due to upregulated tryptophan biosynthesis in response to hypoxic stress35. This suggests that robust tryptophan metabolism in high-altitude donkeys potentially aids basal metabolism in such environments, enhancing hypoxia endurance.

In CD_VS_DZ and RKZ_VS_CD groups, calcidiol and AMP DMs were both enriched in parathyroid hormone (PTH) synthesis and secretion pathways. PTH synthesis and secretion are primarily regulated by serum calcium levels, with additional influence from factors such as vitamin D and neurotransmitters41. PTH regulates calcium levels42, and calcidiol is a hydroxylated form of vitamin D343,44,45,46. During hypoxia, diminished cellular adenosine triphosphate (ATP) levels activate adenosine monophosphate-activated protein kinase (AMPK), which indirectly fosters autophagy via nodular sclerosis complex 2 phosphorylation47. AMPK activity responds to intracellular energy status, particularly changes in the ATP/AMP ratio. Reduced ATP and elevated AMP levels, which indicate increased cellular energy expenditure, augment AMPK activity, modulate cellular energy metabolism, and maintain energy balance48,49. AMPK activation is influenced by energy status and pathological conditions such as glucose deficiency, ischemia, hypoxia, and oxidative stress50. Overall, AMP plays an important role in Tibetan donkeys’ adaptation to high altitudes.

Materials and methods

Sample collection

A total of 30 blood samples, 10 from each region, were collected from donkeys in three regions: Shigatse (RKZ, 5,000 m above sea level) and Changdu (CD, 3,500 m above sea level) in Tibet, as well as Dezhou (DZ, 38 m above sea level) in Shandong. All subjects were adult healthy donkeys aged 5–10 years, in medium condition, and fed with local feed. Blood samples were collected using disposable vacuum coagulant tubes and left at room temperature for 20 min. They were then centrifuged at 12,000 rpm for 10 min to separate the serum. The serum samples were promptly frozen in liquid nitrogen for subsequent non-targeted metabolomics detection. This study was approved by the Institutional Animal Care and Use Committee (IACUC) of Qingdao Agricultural University.

Sample preparation

The serum samples were analyzed using the LC-MS (Thermo Scientific, Vanquish Horizon UHPLC System-Q Exactive HF-X) platform. After thawing at room temperature, 100 µL of each sample was transferred into 1.5 mL centrifuge tubes, and 300 µL of methanol was added. An internal standard solution of 2-chloro-l-phenylalanine (0.02 mg/mL), prepared by mixing equal volumes of methanol and acetonitrile, was added. The solution was vortex-mixed for 30 s and subjected to low-temperature ultrasound extraction for 30 min at 5 °C and 40 kHz. Next, the ultrasonic extract was chilled at -20 °C for 30 min, centrifuged at 4 °C for 15 min at 13,000 × g, and the supernatant was collected and dried under nitrogen. Afterward, a reconstitution solution was prepared by mixing 100 µL of acetonitrile and water in a 1:1 ratio. The dried extract was reconstituted with the reconstitution solution using an ultrasound probe for 5 min at 5 °C and 40 kHz. It was then centrifuged for 10 min at 4 °C and 13,000 × g. The supernatant was transferred to a lined injection bottle for LC-MS analysis.

UPLC-Orbitrap MS

Chromatographic separation was performed using an ACQUITY UPLC HSS T3 column (100 mm × 2.1 mm i.d., 1.8 μm; Waters, Milford, USA) preheated to 40 °C. A 2 µL serum sample was injected and analyzed at 4 °C. Samples were eluted using a mobile phase for ESI + and ESI– consisting of 95% water and 5% acetonitrile with 0.1% formic acid as solvent A, and 47.5% acetonitrile, 47.5% isopropanol, and 5% water with 0.1% formic acid as solvent B. The flow rate was 0.35 mL/min. The mobile phase (A: B) elution gradient was as follows: 0% for 0–1.5 min, 80%: 20% at 1.5 min, and 0%: 100% at 9.5 min, followed by 3 min of re-equilibration.

Serum metabolite analysis

Raw data were performed by ProgenesisQI (Waters Corporation, Milford, USA) for processing, and a three-dimensional data matrix in CSV format was exported. The data matrix was generated from the pre-processing results, which were composed of the retention time (RT), the mass-to-charge ratio (m/z) values, and peak intensities.

Statistical analysis

The data matrix obtained by searching the database was uploaded to the Majorbio cloud platform (https://cloud.majorbio.com) for data analysis. Fistly, data preprocessing was performed using the MetaboAnalystR package of R Studio (Version 3.3.0, https://github.com/xia-lab/MetaboAnalystR) on the Majorbio Cloud Platform. Metabolic features detected at least 80% in any set of samples were retained. To reduce the errors caused by sample preparation and instrument instability, the response intensities of the sample mass spectrometry peaks were normalized using the sum normalization method, to obtain the normalized data matrix. Meanwhile, the variables of QC samples with relative standard deviation (RSD) > 30% were excluded and log10 logarithmicized, to obtain the final data matrix for subsequent analysis.

A multivariate statistical analysis was performed using ropls package of R Studio (Version 1.6.2, http://bioconductor.org/packages/release/bioc/html/ropls.html) from Bioconductor on the Majorbio Cloud Platform. PLS-DA was used to identify metabolic variations between comparable groups. PLS-DA is a supervised learning model. To verify model reliability, group labels were randomly permuted for each sample, followed by modeling and prediction. All metabolite variables were scaled to Pareto Scaling prior to OPLS-DA analysis. Model validity was evaluated from model parameters R2 and Q2, which provided information about interpretability and predictability, respectively, of the model and ensured that over-fitting is avoided. Each run produced R2 and Q2 values. After 200 random permutations, regression lines for Q2 and R2 were generated. Metabolites were considered DMs if they had a P-value less than 0.05 in the Student’s t-test and a Variable Importance in Projection (VIP) score above 1.0 within the PLS-DA model. Results were graphically represented in volcano plots.

For further analysis, KEGG pathway enrichment, heatmap generation, and cluster analysis were performed using Scipy in Python (version 1.0.0, https://docs.scipy.org/doc/scipy/). Additionally, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/, accessed on 20 October 2023) was used for metabolite pathway analysis.

Conclusion

We analyzed differences in serum metabolomics between high-altitude Tibetan and low-altitude Dezhou donkeys. Four significantly enriched metabolic pathways and nine related DMs were identified. These DMs are associated with nutrient digestion and absorption, anti-inflammatory responses, disease resistance, anti-aging, and immune function in Tibetan donkeys at high altitudes. This study provides new insights into the adaptive regulatory mechanisms of Tibetan donkeys.