Transcriptomic Profiling of Hypoxia-Adaptive Responses in Tibetan Goat Fibroblasts

Abstract

The Tibetan goat (Capra hircus) exhibits remarkable adaptations to high-altitude hypoxia, yet the molecular mechanisms remain unclear. This study integrates RNA-seq, WGCNA, and machine learning to explore gene-environment interactions (G × E) in hypoxia adaptation. Fibroblasts from the Tibetan goat and Yunling goat were cultured under hypoxic (1% O₂) and normoxic (21% O₂) conditions, respectively. This identified 68 breed-specific (G), 100 oxygen-responsive (E), and 620 interaction-driven (I) Differentially Expressed Genes (DEGs). The notably higher number of interaction-driven DEGs compared to other effects highlights transcriptional plasticity. We defined two gene sets: Environmental Stress Genes (n = 632, E ∪ I) and Genetic Adaptation Genes (n = 659, G ∪ I). The former were significantly enriched in pathways related to oxidative stress defense and metabolic adaptation, while the latter showed prominent enrichment in pathways associated with vascular remodeling and transcriptional regulation. CTNNB1 emerged as a key regulatory factor in both gene sets, interacting with CASP3 and MMP2 to form the core of the protein-protein interaction (PPI) network. Machine learning identified MAP3K5, TGFBR2, RSPO1 and ITGB5 as critical genes. WGCNA identified key modules in hypoxia adaptation, where FOXO3, HEXIM1, and PPARD promote the stabilization of HIF-1α and metabolic adaptation through the HIF-1 signaling pathway and glycolysis. These findings underscore the pivotal role of gene-environment interactions in hypoxic adaptation, offering novel perspectives for both livestock breeding programs and biomedical research initiatives.

Keywords: DEGs; HIF-1 signaling pathway; WGCNA; gene–environment interactions; machine learning.

Figure 1

Immunofluorescence of goat fibroblast cells. (A) Immunofluorescence of D group goat fibroblast cells showing red fluorescence and blue cell nuclei (DAPI); (B) Immunofluorescence of Y group goat fibroblast cells showing red fluorescence and blue cell nuclei (DAPI).

Figure 1

Immunofluorescence of goat fibroblast cells. (A) Immunofluorescence of D group goat fibroblast cells showing red fluorescence and blue cell nuclei (DAPI); (B) Immunofluorescence of Y group goat fibroblast cells showing red fluorescence and blue cell nuclei (DAPI).

Figure 2

Sample structure and analysis. (A) Heatmap of sample correlations; (B) Principal component analysis (PCA) of samples.

Figure 3

Analysis of differentially expressed genes (DEGs). (A) Volcano plot of DEGs with different effects; (B) Venn diagram of DEGs.

Figure 4

Functional enrichment and protein–protein interaction (PPI) networks. (A,B) GO biological process (BP) enrichment for environmental stress genes and genetic adaptation genes, highlighting the top 10 terms with the lowest p-values. (C,D) KEGG pathway enrichment, displaying the top 20 enriched pathways for each gene set. (E,F) PPI network diagrams, identifying key hub genes involved in hypoxia adaptation.

Figure 5

Identification of hypoxia-related genes using Lasso regression and Random Forest. (A) Optimization of Lasso regression parameters for hypoxia-related gene selection; (B) Random Forest-based ranking of key genes contributing to hypoxia adaptation, as determined by mean decrease in Gini scores. (C) The top 20 candidate genes identified by Random Forest in the environmental stress gene set. (D) The top 20 candidate genes identified by Random Forest in the genetic adaptation gene set. (E) The overlapping genes identified by both machine learning methods among the environmental stress genes were MAP3K5, TGFBR2 and ITGB5. (F) The overlapping genes identified by both machine learning methods among the genetic adaptation genes were RSPO1, TGFBR2, and ITGB5.

Figure 6

Module–trait relationships and gene significance analysis. (A) Heatmap showing the correlation between gene modules and traits under different conditions; (B) Scatter plot of MEblack module membership vs. gene significance; (C) Scatter plot of MElightcyan module membership vs. gene significance; (D) Choosing the best soft-threshold power.

Figure 6

Module–trait relationships and gene significance analysis. (A) Heatmap showing the correlation between gene modules and traits under different conditions; (B) Scatter plot of MEblack module membership vs. gene significance; (C) Scatter plot of MElightcyan module membership vs. gene significance; (D) Choosing the best soft-threshold power.

Figure 7

Enrichment and PPI network analysis of black and light cyan modules. (A,B) GO_BP enrichment in the black (A) and light cyan (B) modules; (C,D) KEGG pathway enrichment in the black (C) and light cyan (D) modules; (E,F) PPI networks of the black (E) and light cyan (F) modules.

Figure 8

Different letters (a, b, c) indicate significant differences between groups (p < 0.05, Duncan’s multiple range test), while the same letter indicates no significant difference.