Thirty healthy Diannan small-ear pigs (DNSP) from Xishuangbanna, Yunnan Province, were selected for this study. Blood samples were collected from the ear vein into 5 mL collection tubes containing heparin sodium as an anticoagulant. The samples were stored and transported at −20°C and will be used for genomic DNA extraction. Additionally, the genomic resequencing data of 30 Tibetan pigs were downloaded from the NCBI database (Accession Number: PRJNA186497, Supplementary Table S1).

Moreover, three healthy Tibetan pigs (TP) from Diqing Tibetan Autonomous Prefecture (Altitude about 3,000 m) and three Diannan small-ear pigs (DNSP) from Xishuangbanna (Altitude about 500 m), Yunnan Province, were slaughtered. Heart, liver, spleen, lung, and bone marrow tissues were collected and placed into cryogenic vials pre-filled with RNA preservation solution. These samples were immediately stored in liquid nitrogen for subsequent RNA extraction and transcriptome sequencing. Animal protocol undergoing review and approval by the Experimental Animal Welfare and Animal Experiment Ethics Review Committee of Foshan University (approval number: FOSU202401-28). This process was conducted in accordance with the guidelines established by the Animal Use Committee of the Ministry of Agriculture of China, Beijing, aimed at minimizing animal suffering throughout the study.

An intersection analysis was conducted between the candidate genes identified as being under positive selection in the Tibetan pig population, which are likely involved in adaptation to hypoxic environments, and the differentially expressed genes identified across the five tissues. The intersecting genes, which are both under positive selection in the Tibetan pig population and differentially expressed between breeds, represent the results of the integrated analysis. By calculating Fst values and nucleotide diversity, combined with Tajima’s D, we can visualize the selection of genes in the Tibetan pig population. Tajima’s D is used to detect selection pressure; when one population has a Tajima’s D value below zero while the other is close to zero or positive, it may indicate positive selection in the former’s gene region. Combining Fst values with nucleotide diversity allows for more effective detection of genes under positive selection in the Tibetan pig population compared to Diannan small-ear pigs. The visualization of differentially expressed genes will use gene expression levels normalized by Transcripts Per Million (TPM) and will be presented through box plots.

To validate the results of the sequencing analysis, we performed RT-qPCR analysis on the selected genes. Total RNA was extracted from the corresponding high-expression tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, United States), and then reverse transcribed into cDNA using the TIANGEN Kit KR210831. The 2 × SuperReal PreMix Plus (SYBR Green) was used for the qPCR reaction, with 2 μL of cDNA template and forward and reverse primers added. A reaction volume of 20 μL was used in the qPCR reactions according to the manufacturer’s protocol. qPCR experiments were performed in triplicate, and the average Ct was used for further analysis. The 2^−ΔΔCt method was used to determine relative mRNA abundance. The reaction conditions were as follows: 95°C for 15 min (1 cycle), 95°C for 10 s, and 60°C for 20 s (40 cycles). The sequences of the internal primers used for RT-qPCR are listed in Supplementary Table S4.

The genome resequencing data of 30 Diannan small-eared pigs generated a total of 1.52 Tb of raw data. The whole-genome resequencing of 30 Diannan small-ear pigs generated a total of 1.52 Tb of raw data. After filtering and quality control, the average base quality scores for Q20 and Q30 were 96.22 and 90.33%, respectively. Q20 represents a base quality score where the probability of an incorrect base call is 1 in 100 (99% accuracy), while Q30 indicates a base quality score with a probability of 1 incorrect base call in 1,000 (99.9% accuracy). These metrics reflect the overall sequencing accuracy and the reliability of the data for downstream analysis. The average alignment rate was 99.43%, and the average sequencing depth was 19.33× (Supplementary Table S2). Similarly, the resequencing data of 30 Tibetan pigs downloaded from the NCBI database underwent the same processing, resulting in average base quality scores of 96.64% for Q20 and 90.80% for Q30, an average alignment rate of 99.61%, and an average sequencing depth of 5.81× (Supplementary Table S2). These results indicate that the quality of the whole-genome resequencing data is high and suitable for subsequent analyses.

The principal component analysis (PCA) results showed a clear separation between the two populations, with individuals from the same population clustering together. The PCA results were consistent with the expected classification (Figure 1). The genetic differentiation index (Fst) between the Tibetan pig and Diannan small-ear pig populations, calculated using the sliding window method, was 0.11. In population genetics, higher Fst values indicate greater genetic differentiation between two populations. An Fst value between 0.05 and 0.15 generally indicates moderate genetic differentiation. Therefore, both the PCA results and the Fst value suggest the presence of genetic differences between the two populations, supporting further exploration of these differences.

To identify functional genes at the genomic level that may contribute to the adaptation of Tibetan pigs to hypoxic environments, we selected the Diannan small-ear pig, which lives in low-altitude areas, as the control group. A selection signal analysis was conducted using a combination of the genetic differentiation index Fst and nucleotide diversity θπ. The top 1% of regions were selected as candidate regions under selection, with thresholds of log2(θπ_DNSP/θπ_TP) ≥ −0.27 and Z(Fst) ≥ 1.88 (Figure 2). Gene annotation within the candidate regions revealed that 461 genes were under positive selection in the Tibetan pig population.

The results of the GO function and KEGG pathway enrichment analyses for the 461 genes. The GO enrichment analysis revealed significant enrichment of candidate genes in the following categories (Figure 3A): response to hypoxia, coronary vasculature development, heart morphogenesis, hematopoietic progenitor cell differentiation, and cellular response to UV. These genes may contribute to the adaptation of Tibetan pigs to the hypoxic and high UV environments of high-altitude regions by influencing vascular development, hematopoietic function, and cellular responses to ultraviolet radiation (23–25). In the KEGG pathway enrichment analysis, the candidate genes were predominantly enriched in the following pathways (Figure 3B): Rap1 signaling pathway, Ras signaling pathway, PI3K-Akt signaling pathway, MAPK signaling pathway, and platelet activation pathway. These pathways are primarily involved in vascular formation and repair, suggesting that angiogenesis is a key process in the adaptation of Tibetan pigs to hypoxic conditions (13, 26, 27). By examining the functions of genes associated with hypoxia adaptation within the enriched GO terms and KEGG pathways, and integrating analyses of Fst, nucleotide diversity, and Tajima’s D, we identified 10 candidate genes: HES4, ANGPT1, HIF3A, SPHK2, PCK2, RCN3, HIGD2A, DNM2, IRF9, and SRF. Selection signal analysis indicates that these genes have undergone positive selection in the Tibetan pig population, as evidenced by significantly high Fst values within the selected regions, reduced nucleotide diversity in the selected population, and notable differences in Tajima’s D values between the two populations (Supplementary Figure S1).

RNA from five tissues (heart, liver, spleen, lung, and bone marrow) of both Tibetan pigs and Diannan small-ear pigs was collected for transcriptome sequencing, with three biological replicates per group, resulting in a total of 30 transcriptome sequencing datasets. After filtering and aligning low-quality bases, the average base quality scores of Q20 and Q30 were 97.23 and 94.85%, respectively, and the average alignment rate was 99.23% (Supplementary Table S3). These results indicate that the transcriptome sequencing data are reliable.

To ensure the reliability of the differential expression analysis, we first performed principal component analysis (PCA) on the TPM-normalized gene expression matrix using the PCAtools package in R to detect any outlier samples. The results (Figure 4) showed that tissue samples from the same breed exhibited good clustering across different tissues, with no obvious outliers, indicating high sequencing quality. Subsequently, we used the DESeq2 package to analyze differentially expressed genes (DEGs), with the criteria for selection being |log2Fold Change| > 1 and Padj <0.05. The results of the differential expression analysis for the various tissues are as follows (Figure 5). In heart tissue, we identified 1,762 DEGs (918 upregulated and 844 downregulated). In liver tissue, 2,437 DEGs were identified (1,328 upregulated and 1,109 downregulated). In spleen tissue, 119 DEGs were identified (86 upregulated and 33 downregulated). In lung tissue, 2,504 DEGs were identified (1,070 upregulated and 1,434 downregulated). Finally, in bone marrow tissue, 212 DEGs were identified (159 upregulated and 53 downregulated).

In the selection signal analysis, we identified 10 candidate genes—HES4, ANGPT1, HIF3A, SPHK2, PCK2, RCN3, HIGD2A, DNM2, IRF9, and SRF—through functional annotation of genes under selection in the Tibetan pig population and a review of relevant literature (12, 28–36). By examining the expression of these 10 genes across five tissues in the two populations, we found that HIF3A showed significant expression in heart tissue (Figures 6A,B), HIGD2A and PCK2exhibited significant differential expression in liver tissue (Figures 6C,D,G,H), IRF9 showed significant differential expression in bone marrow (Figures 6E,F), and RCN3 exhibited significant expression in lung tissue (Figures 6I,J). The results from the selection signal analysis show that the regions in the genome under selection exhibit high Fst values. In population genetics studies, it is generally considered that when Fst values range between 0.05 and 0.15, there is moderate genetic differentiation between populations. In this joint analysis, the Fst values of the five selected genes were all greater than 0.05, indicating moderate genetic differentiation in the regions where these genes are located. When populations undergo positive selection, the selected population typically displays lower nucleotide diversity, and the Tajima’s D values between populations differ, with the Tajima’s D of the selected population usually being less than 0. In the regions where these genes are located, the nucleotide diversity of the Tibetan pig population is lower than that of the Diannan small-ear pig, and its Tajima’s D value is not only lower than that of the Diannan small-ear pig population but also less than 0. This indicates that these five genes are under positive selection in the Tibetan pig population. Combined with the results of transcriptome analysis, we found that there were also significant differences in the expression of these five genes in the tissues. Therefore, HIF3A, RCN3, HIGD2A, PCK2, and IRF9 have been identified as key functional genes associated with hypoxia adaptation through the integrated analysis of selection signals and transcriptome data.

We validated the RNA-seq data for the five selected genes, HIGD2A, HIF3A, RCN3, PCK2, and IRF9, and found significantly higher expression levels in Tibetan pigs living in high-altitude areas compared to Diannan small-ear pigs from lower altitudes (Figure 7). These results suggest that identifying DEGs through transcriptome data is effective and reliable.

Archaeological and genetic data provide significant evidence that human migration and cultural expansion during the late Pleistocene (22,000–18,000 years ago) and the early Neolithic period (10,000–7,000 years ago) influenced the ecological and agricultural development of the plateau regions. This process introduced agriculture to the Himalayas and facilitated the establishment of animal husbandry on the plateau, impacting the domestication and gene flow of high-altitude animals (37–39). The origin of Tibetan pigs is closely linked to the migration of the Tibetan people, as evidenced by 5,000-year-old pig bones from the Karuo archaeological site, which suggest that the Tibetans domesticated Tibetan pigs (11). Due to their long-term adaptation to the plateau environment, Tibetan pigs have evolved unique physiological traits, such as a small body size, long hair, and thick subcutaneous fat, which help them withstand the cold climate. Additionally, Tibetan pigs have a higher red blood cell count, increased hemoglobin levels, and more efficient metabolism, allowing them to utilize oxygen more effectively and survive at high altitudes (40–42). In contrast, Diannan small-ear pigs primarily inhabit the subtropical regions of southern Yunnan Province and exhibit genetic adaptations to the hot and humid climate of tropical and subtropical environments. These adaptations include smaller ears, a slightly longer snout, and a short, thick neck, enabling them to thrive in low-altitude, humid conditions (43–45). Overall, the significant differences in environmental adaptability and biological characteristics between Tibetan pigs and Diannan small-ear pigs make them ideal models for studying the adaptability and traits of pig breeds in different environments. In this study, we first used principal component analysis (PCA) and the genetic differentiation index (Fst) to investigate the two populations. The PCA results revealed clear genetic differences between the two populations, while the Fst value indicated a moderate level of genetic differentiation between them, consistent with our expectations.

In this study, 461 candidate genes were identified through genomic selection signal analysis. After conducting KEGG pathway and GO functional enrichment analyses, and reviewing relevant literature, 10 potential hypoxia adaptation genes were identified: HES4, ANGPT1, HIF3A, SPHK2, PCK2, RCN3, HIGD2A, DNM2, IRF9, and SRF. Chen et al. (28) performed single-cell transcriptome analysis on yak heart tissue and, through cell type annotation and differential expression analysis, found that altitude changes primarily affected the expression patterns of smooth muscle cells and vascular endothelial cells. Four key transcription factors (MEF2B, FOXP4, ARID5A, and HES4) were identified in smooth muscle cells. Among them, HES4 can help cells cope with environmental stress, such as hypoxia, by regulating the expression of heat shock proteins (HSPs). Therefore, HES4 may be a potential gene for high-altitude hypoxia adaptation in Tibetan pigs. The ANGPT1 gene, a crucial angiogenic factor, can promote tissue angiogenesis and enhance oxygen supply (46). It has been identified as a positively selected gene in the Tibetan population (47). Previous research has also identified ANGPT1 as a candidate gene for high-altitude adaptation by comparing the whole-genome sequences of wolves from different altitudes (29). HIF3A has been reported to promote high-altitude adaptation in both Tibetan pigs (12) and yaks (48), although some studies suggest it is also associated with obesity (49). Hypoxia activates hypoxia-inducible factors (HIFs), leading to changes in the expression of downstream genes. HIF-3α can directly bind to erythropoietin (EPO) and induce its transcription (50). EPO is a key growth factor that regulates red blood cell count and hematocrit, and is controlled by hypoxia-inducible factor-2α (HIF-2α). It is produced in the liver during early development and later in the kidneys (51), playing a crucial role in hypoxia adaptation (52–54). In mouse kidney fibroblast cell lines, knockdown of the sphingosine kinase 2 (SPHK2) gene was found to increase the synthesis of HIF-2α and EPO (55). Other studies have shown that SPHK2 is associated with hypoxia-inducible factor-1α (HIF-1α) within protein complexes, and is enriched in the promoter regions of HIF target genes such as vascular endothelial growth factor (VEGF), enhancing local histone H3 acetylation and transcription (30). Under various internal and external stress factors such as hypoxia, glucose deprivation, ATP depletion, calcium overload, and reduced protein degradation, the endoplasmic reticulum can lose its homeostasis (56). Phosphoenolpyruvate carboxykinase 2 (PCK2) is involved in cellular metabolism, including glucose metabolism and the tricarboxylic acid cycle (TCA), and can promote endoplasmic reticulum stress, playing a vital role in maintaining cellular homeostasis and activity (57). RCN3, a calcium-binding protein, plays a role in calcium homeostasis and protein folding within the cell. Calcium ions have multiple functions in cardiac and pulmonary tissues, including promoting myocardial contraction and relaxation, maintaining heart rhythm, regulating alveolar ventilation and gas exchange, and participating in cell signaling and metabolic regulation. Calcium ion homeostasis is crucial for the normal function of cardiac and pulmonary tissues (58, 59). HIGD2A, which encodes a hypoxia-induced domain protein, is involved in the assembly of the mitochondrial respiratory chain and is essential for cell survival under hypoxic conditions (33). DNM2 is a key regulator of intracellular membrane dynamics, including endocytosis and vesicle transport. Research by Joshi et al. has found a regulatory relationship between DNM2 and hypoxia-inducible factor-1α (HIF-1α) (34). Interferon regulatory factor 9 (IRF9) is an important transcription factor in the interferon signaling pathway. Studies have shown that interferon signaling can regulate the expression or activity of HIF-1α, thereby influencing cellular adaptation to hypoxia. Conversely, HIF-1α can regulate some interferon-related genes, affecting immune responses (60, 61). Serum response factor (SRF) plays a critical role in regulating cell growth, differentiation, and survival. Research has shown that SRF is important in hypoxia adaptation, as it can regulate the expression of various genes related to the hypoxic response (36, 62, 63). In summary, through selection signal analysis, we have identified 10 key candidate genes potentially involved in the hypoxia adaptation of Tibetan pigs.

Gene expression regulation plays a crucial role in species adaptation to environmental changes. This complex biological process allows organisms to respond to and adapt to various environmental shifts by modulating the activity of specific genes (64). Therefore, in this study, we conducted transcriptomic analysis on five tissues—heart, liver, spleen, lung, and bone marrow—from Tibetan pigs and Diannan small-ear pigs. These tissues encompass key functions such as blood circulation, oxygen transport, and hematopoiesis, enabling us to identify important functional genes involved in hypoxia adaptation in Tibetan pigs through differential gene expression analysis. After identifying differentially expressed genes across various tissues, we conducted a combined analysis of the 10 candidate genes for hypoxia adaptation, previously identified through selection signal screening, and the differentially expressed genes. The results revealed that five genes—HIF3A, RCN3, HIGD2A, PCK2, and IRF9—not only underwent positive selection in the Tibetan pig population but also showed significant differential expression in various tissues. Animals living in high-altitude regions often exhibit enhanced cardiopulmonary function to cope with hypoxic environments by increasing blood flow rate and oxygen content (5, 6, 65). We found that the HIF3A gene was differentially expressed in heart tissue. As a member of the HIF family, HIF3A, along with other members like HIF1A, can promote angiogenesis under hypoxic conditions (66). The high expression of HIF3A in heart tissue may enhance the formation of blood vessels within the heart, thereby improving its blood supply capacity. The RCN3 gene was differentially expressed in lung tissue. Hypoxic environments can trigger endoplasmic reticulum stress in cells, potentially leading to protein misfolding (67). As an endoplasmic reticulum calcium-binding protein, the high expression of RCN3 in the lungs of Tibetan pigs may help enhance the cell’s ability to cope with hypoxia-induced endoplasmic reticulum stress, thus protecting lung cells from hypoxic damage. PCK2 and HIGD2A were both highly expressed in the liver tissue of Tibetan pigs. PCK2 encodes an enzyme involved in gluconeogenesis, and its high expression may enhance this pathway, helping to maintain blood glucose levels and ensuring sufficient energy supply to glucose-dependent tissues such as the brain (31). HIGD2A is involved in regulating mitochondrial function, and its high expression may help maintain mitochondrial respiratory efficiency in the liver, ensuring energy balance under hypoxic conditions (33, 68). The IRF9 gene was highly expressed in the bone marrow tissue of Tibetan pigs. IRF9 is a key regulator in the interferon signaling pathway, playing a role in the immune response. The IRF gene family is also important in the differentiation of hematopoietic cells (69). While IRF9 is primarily considered an immune regulator, its potential role in hematopoiesis should not be overlooked. In summary, through combined analysis, we identified five key genes that contribute to the hypoxia adaptation of Tibetan pigs.

In this study, although several important candidate genes related to the high-altitude adaptation of Tibetan pigs were identified through genomic selection signal analysis and transcriptomic analysis, there are still some limitations that need to be considered. This research primarily focused on two specific pig populations—Tibetan pigs and Diannan small-ear pigs. While these pig breeds exhibit significant differences in adapting to different environments, the differences in their genomes and transcriptomes may not represent the adaptive evolution of all high-altitude and low-altitude pig breeds. Therefore, the generalizability of the research findings requires further validation. This study experimental validation was limited to qPCR, which somewhat restricts the in-depth understanding of the functions of these genes. While qPCR can verify the expression differences of candidate genes in various tissues, this method does not directly reveal the specific functional mechanisms of these genes in low-oxygen adaptation. Therefore, future research should adopt additional experimental methods for functional validation to enhance the depth and reliability of the study, thereby providing a comprehensive understanding of the molecular mechanisms underlying Tibetan pigs’ adaptation to the low-oxygen environment of high altitudes.