Integrative Genomic, Transcriptomic and Epigenomic Analysis Reveals cis-regulatory Contributions to High-altitude Adaptation in Tibetan Pigs

Abstract

The Qinghai-Tibet Plateau, characterized by its extreme environmental conditions, presents significant challenges to life, making it an ideal region for studying adaptation and evolution. Tibetan pigs, known for their high genetic diversity and exceptional adaptability to high altitudes, serve as excellent models for investigating high-altitude adaptation. While previous studies have extensively identified genetic determinants associated with high-altitude adaptation, the molecular mechanisms, particularly cis-regulatory patterns, remain poorly understood. Here, we conducted a selective sweep analysis using 484 genomes from Chinese and Western pig breeds across various altitudes, revealing 38.56 Mb of genomic regions under selection in Tibetan pigs. Enrichment analysis identified the lung as the primary functional tissue involved in high-altitude adaptation, supported by tissue-specific transcriptional and regulatory patterns observed between Tibetan and Meishan pigs (low altitude). By integrating genomic, RNA-seq, ATAC-seq, and H3K27ac HiChIP data, we constructed comprehensive enhancer-promoter regulatory maps of candidate genes and pinpointed promising genetic determinants associated with high-altitude adaptation, including SNPs in EPAS1, KLF13, SPRED1, and CFD. These loci were predicted to influence chromatin accessibility and the interactions of regulatory elements, with altered binding strength of relevant transcription factors. Further in vitro experiments confirmed that these loci function as allele-specific enhancers, modulating the expression of target genes. Our findings elucidate the regulatory basis of high-altitude adaptation in Tibetan pigs and provide valuable insights for exploring hypoxia-related diseases in livestock and humans.

Keywords: H3K27ac HiChIP; Tibetan pigs; high-altitude adaptation; lung tissue; multi-omics integration.

Fig. 1.

Experimental design and genome-wide analysis of selective signals associated with high-altitude adaptation in Tibetan pigs. a) Schematic representation of the experimental design for multi-omics data integration and analysis. b) Genome-wide distribution of pairwise F_ST values comparing Tibetan pigs (TT) with lowland pigs (lowland). c) Genome-wide selective signals identified through θπ ratios (θπ, lowland/θπ, TT), presented as log₁₀-transformed values. The top 5% significant threshold is indicated by a black dash line, with annotated genes in selected regions. d) Enrichment of 15 chromatin states across 14 tissue types (Adipose, Cecum, Cerebellum, Colon, Cortex, Duodenum, Hypothalamus, Ileum, Jejunum, Liver, Lung, Muscle, Spleen, and Stomach) in candidate genomic regions associated with high-altitude adaptation. Chromatin states included promoters (TssA, TssAHet, and TssBiv), TSS-proximal transcribed regions (TxFlnk, TxFlnkWk, and TxFlnkHet), enhancers (EnhA, EnhAMe, EnhAWk, EnhAHet, and EnhPois), repressed regions (Repr and ReprWk), quiescent regions (Qui), and accessible but did not coincide with any other measured epigenetic marks (ATAC islands) (Pan et al. 2021). e) Tissue-specific promoters (TssA) enrichment across 14 tissue types in candidate genomic regions related to high-altitude adaptation. The enrichment intensity of the All_common promoter serves as a reference for that of each tissue-specific promoter. Statistical significance was determined using two-tailed Fisher's exact test.

Fig. 2.

Tissue-specific transcriptional landscape associated with high-altitude adaptation. a) Principal component analysis (PCA) of RNA-seq data from 12 samples across heart and lung tissues of Tibetan pigs (TT) and Meishan pigs (MS). b) Hierarchical clustering analysis of global gene expression patterns. c) Volcano plot illustrating differentially expressed genes (DEGs) between MS and TT in lung tissue, with significant DEGs highlighted. No sig, no significance; up, upregulated. d) Expression profiles of three high-altitude adaptation-related genes across each group. Statistical significance was determined by unpaired Student's t-test (*P-value < 0.05, **P-value < 0.01, ***P-value < 0.001, ns, non-significant). Pearson's correlations of log₂ (TPM + 1) between the cultured human umbilical vein endothelial cells (HUVECs) of Tibetan highlanders and the lung (e) and heart tissues (f) of Tibetan pigs.

Fig. 3.

Chromatin accessibility landscape underlying high-altitude adaptation. a) Average ATAC-seq peak signals relative to transcription start site (TSS) in lung tissues across replicates. b) Principal component analysis of 12 ATAC-seq samples. c) Genomic distribution of ATAC-seq peaks. d) Correlations analysis between chromatin accessibility and expression levels of the shared genes identified as both differentially expressed genes (DEGs) and differentially accessible genes (DAGs) in lung tissues. Gene-peak associations were determined based on nearest TSS (±2 kb). Statistical significance was assessed using linear regression (two-sided P-values). e) Representative example of regions with differential chromatin accessibility at the EPAS1 locus between MS and TT lung tissues. RefSeq, reference sequence. Accessibility differences were evaluated using unpaired Student's t-test (*P-value <0.05, **P-value <0.01). f) Heatmaps visualization of transcription factors (TF) binding preferences in differentially open regions (DORs) from TT-MS comparisons in each group. Color intensity represents normalized TF expression levels. P-values were calculated using hypergeometric tests and adjusted for multiple comparisons with the Benjamini–Hochberg method. g) Enrichment analyses of candidate selected regions within general ATAC open regions, with control regions as reference. Significance was determined by two-tailed Fisher's exact test. h) Comparative analysis of normalized accessibility (ATAC signals) surrounding SNPs (±200 bp) within versus outside selected regions. Data distributions are shown as violin plots with embedded box plots (median line, interquartile range box, and full data range whiskers). n = 31,515 (SNPs within selected regions), and 30,000 (SNPs outside selected regions). P-values were obtained using unpaired two-tailed Student's t-test.

Fig. 4.

Functional analysis of H3K27ac-mediated chromatin interactions related to high-altitude adaptation. a) Venn diagram showing the overlap between differential loop genes (DLGs) and differentially expressed genes (DEGs) in lung and heart tissues. P-values were calculated using the hypergeometric test. b) IGV (Integrative Genomics Viewer) tracks and normalized H3K27ac HiChIP signals demonstrating that EGLN3 is regulated by stronger enhancer-promoter interactions in MS compared to TT in lung tissues. c) Significantly enriched GO terms (biological processes) for genes associated with differential chromatin loops between TT and MS in lung and heart tissues. d) Heatmaps displaying transcription factors (TF) binding preferences at differential loop anchor regions from TT and MS comparison in each group. P-values were calculated using hypergeometric tests and adjusted for multiple comparisons with the Benjamini–Hochberg method. The color scale represents the normalized expression levels of the corresponding TFs.

Fig. 5.

Integrated analysis for identifying candidate variants and target genes. a) Flowchart illustrating the model for identifying candidate variants and their corresponding target genes associated with high-altitude adaptation. Firstly, we kept SNPs in the top 5% selected regions (F_ST ≥ 0.207 and log₁₀ θπ ratio (θπ_lowland/θπ_TT) ≥ 0.157). Secondly, we filtered for SNPs showing significant allele frequency differences (|ΔAF| ≥ 0.5) between the Tibetan pigs and other lowland breeds. Thirdly, the resultant SNPs located in ATAC peaks regions were further selected. Fourthly, we screened for SNPs residing in H3K27ac HiChIP loop anchors (enhancer-promoter contact regions). Fifthly, SNPs corresponding to expressed target genes were finally retained. b) Functional annotations of candidate variants. c) Cumulative fractions of the number of candidate SNPs associated with each target gene (mean = 2.66). d) Cumulative fractions of the number of target genes linked to each candidate SNP (mean = 3.47). e) Cumulative fractions of genomic distances for identified variant-target gene (TG) connections (mean = 108.62 kb). f) Bar chart displaying the number of genes skipped by SNPs before interacting with their target genes. The accompanying pie chart shows the proportion of SNP that skip at least one gene before their target gene.

Fig. 6.

Role of the candidate SNPs in regulating EPAS1 in lung tissue. a) IGV plot showing the distribution of candidate SNPs potentially regulating EPAS1 expression in lung tissue. b) LD block illustrating the linkage of variants in the region surrounding the candidate SNPs. c) Normalized H3K27ac HiChIP signals indicating the contact frequency between the region surrounding the candidate SNPs and the EPAS1 promoter in lung tissues. Comparison of chromatin accessibility in the region around (±200 bp) SNP rs342893663 (d), and rs323786058 (e) between TT and MS in lung tissue. Luciferase assays in HEK-293T cells comparing enhancer activity between the two alleles of SNP rs342893663 (f), and rs323786058 (g). h) SNP rs323786058 is predicted to alter the binding strength of the transcriptional factor EGR1. i) Luciferase-based promoter assay evaluating the transcriptional effect of TF EGR1 on the EPAS1 promoter region in PIEC cells. P-values were calculated using an unpaired Student's t-test, *P-value <0.05, **P-value <0.01.

Fig. 7.

Candidate SNPs for the novel adaptive target gene KLF13 in lung tissue. a) IGV plot and normalized HiChIP signals demonstrating the regulatory role of candidate SNPs potentially involved in KLF13 transcription in lung tissue. b) LD block showing the linkage of candidate variants within a 19.697 kb region. c) Comparison of chromatin accessibility in the region around (±200 bp) rs337794133 between TT and MS in lung tissue. d) Luciferase-based enhancer assay results comparing the enhancer activity of the two alleles of SNP rs337794133 in HEK-293T cells. P-values were calculated using an unpaired Student's t-test, *P-value <0.05, **P-value <0.01, ***P-value <0.001. e) Competitive EMSA of rs337794133. Protein-binding affinity to DNA oligos was evaluated by adding unlabeled rs337794133-C or rs337794133-T to the reaction. f) SNP rs337794133 is predicted to alter the binding strength of the transcriptional factor ETS1.