Whole-genome resequencing to investigate the genetic diversity and mechanisms of plateau adaptation in Tibetan sheep

Abstract

Introduction: Tibetan sheep, economically important animals on the Qinghai-Tibet Plateau, have diversified into numerous local breeds with unique characteristics through prolonged environmental adaptation and selective breeding. However, most current research focuses on one or two breeds, and lacks a comprehensive representation of the genetic diversity across multiple Tibetan sheep breeds. This study aims to fill this gap by investigating the genetic structure, diversity and high-altitude adaptation of 6 Tibetan sheep breeds using whole-genome resequencing data.

Results: Six Tibetan sheep breeds were investigated in this study, and whole-genome resequencing data were used to investigate their genetic structure and population diversity. The results showed that the 6 Tibetan sheep breeds exhibited distinct separation in the phylogenetic tree; however, the levels of differentiation among the breeds were minimal, with extensive gene flow observed. Population structure analysis broadly categorized the 6 breeds into 3 distinct ecological types: plateau-type, valley-type and Euler-type. Analysis of unique single-nucleotide polymorphisms (SNPs) and selective sweeps between Argali and Tibetan sheep revealed that Tibetan sheep domestication was associated primarily with sensory and signal transduction, nutrient absorption and metabolism, and growth and reproductive characteristics. Finally, comprehensive analysis of selective sweep and transcriptome data suggested that Tibetan sheep breeds inhabiting different altitudes on the Qinghai-Tibet Plateau adapt by enhancing cardiopulmonary function, regulating body fluid balance through renal reabsorption, and modifying nutrient digestion and absorption pathways.

Conclusion: In this study, we investigated the genetic diversity and population structure of 6 Tibetan sheep breeds in Qinghai Province, China. Additionally, we analyzed the domestication traits and investigated the unique adaptation mechanisms residing varying altitudes in the plateau region of Tibetan sheep. This study provides valuable insights into the evolutionary processes of Tibetan sheep in extreme environments. These findings will also contribute to the preservation of genetic diversity and offer a foundation for Tibetan sheep diversity preservation and plateau animal environmental adaptation mechanisms.

Keywords: Altitude adaptation; Domestication; Genetic selection; Population genetic structure; Tibetan sheep; Whole-genome sequencing.

Fig. 1

Phylogenetic relationships and population structure of Tibetan sheep. A Geographical distribution map of sample collection. B PCA of Pan sheep and Tibetan sheep. C PCA of the 6 Tibetan sheep breeds. D Utilizing the maximum likelihood algorithm in fast tree software, a phylogenetic tree of the 6 Tibetan sheep breeds was constructed using Argali sheep as the outgroup. E Based on whole-genome SNP data, the population genetic structure of the 6 breeds of Argali and Tibetan sheep were inferred using admixture software (set K = 2–10). PA: Argali sheep; TS: Tibetan sheep; GB: Guinan Black Fur sheep; PT: Plateau Tibetan sheep; ZS: Zhashijia sheep; VT: Valley Tibetan sheep; ZK: Zeku sheep; EL: Euler Tibetan sheep. Each column in the image represents an individual, and the length of each different colored segment indicates the proportion of an ancestor in the individual’s genome. The labels on the image denote the groups

Fig. 2

Gene flow and LD analysis of the 6 breeds of Tibetan and Argali sheep. A Gene flow diagram between Argali sheep and the 6 Tibetan sheep breeds. The arrow direction in the figure corresponds to the direction of gene drift. The x-axis represents the 10 times average standard deviation of the elements in the sample covariance matrix. B Gene flow diagram of the 6 Tibetan sheep breeds. The arrow direction in the figure corresponds to the direction of gene drift. The x-axis represents the 10 times average standard deviation of the elements in the sample covariance matrix. C Distribution map of LD in Argali and Tibetan sheep. The horizontal axis represents the distance between SNPs, and the vertical axis represents the r² value. D Distribution map of LD among the 6 Tibetan sheep breeds. PA: Argali sheep; TS: Tibetan sheep; GB: Guinan Black Fur sheep; PT: Plateau Tibetan sheep; ZS: Zhashijia sheep; VT: Valley Tibetan sheep; ZK: Zeku sheep; EL: Euler Tibetan sheep. The horizontal axis represents the distance between SNPs, and the vertical axis represents the r² value

Fig. 3

Analysis of SNPs unique to Argali and Tibetan sheep. A Venn diagram of unique and shared SNPs in Argali and Tibetan sheep. PA: Argali sheep; TS: Tibetan sheep. B KEGG enrichment analysis of Argali-specific SNPs. C GO enrichment analysis of Argali-specific SNPs. D Enrichment analysis of KEGG pathways for Tibetan sheep-specific SNPs. E Enrichment analysis of GO terms for Tibetan sheep-specific SNPs

Fig. 4

Selective signatures of Tibetan sheep during the introduction of improvements. A Manhattan plot of genome-wide selective sweep signals of Argali and Tibetan sheep. The horizontal axis represents the chromosomes, and the vertical axis represents the xpclr_norm values. Taking the first 1% of the xpclr_norm values (pclr_norm ≥ 3.726865) as the selection criterion, the value is represented by a dashed line, which indicates the selective elimination region during the domestication process of Tibetan sheep. The text in the figure indicates genes related to domestication with significant selection intensity. B KEGG functional enrichment of selectively eliminated genes. C GO functional enrichment of selectively downregulated genes

Fig. 5

Selective elimination analysis of adaptation to Tibetan sheep. A Manhattan plot of genome-wide selective sweep signals in Zhashijia sheep (ZS) and Valley Tibetan sheep (VT). ZS distributed at an altitude of 4,300 m were used as the experimental group, and VT distributed at an altitude of 1,800 m were used as the control group. The horizontal axis represents the chromosomes, and the vertical axis represents the xpclr_norm values. Taking the first 1% of the xpclr_norm values (xpclr_norm ≥ 3.842473) as the selection criterion, the value is represented by a dashed line, which indicates the selective elimination region during the domestication process of Tibetan sheep. The text in the figure indicates genes associated with plateau adaptation with significant selection intensity. B KEGG functional enrichment of selectively eliminated genes. C GO functional enrichment of selectively downregulated genes

Fig. 6

Heart transcriptome results in the high- and low-altitude groups Heart transcriptome results for the high-altitude (Zhashijia sheep, ZS) and low-altitude (Valley Tibetan sheep, VT) groups. A Volcano map of DEGs. HZ: Heart transcriptome of Zhashijia sheep; HV: Heart transcriptome of Valley Tibetan sheep. The X-axis represents the fold change in the difference after conversion to log₂, and the Y-axis represents the significance value after conversion to log₁₀. Red represents upregulated DEGs, blue represents downregulated DEGs, and gray represents non-DEGs. B Bubble map of enriched KEGG pathways. C Bar graph graphs of enriched GO pathways. D PPI network of DEGs (|log₂fold-change| ≥ 2). The color of the circles represents the number of node connections, with redder colors indicating higher numbers of node connections. The color of the connection lines represents the interaction score, with redder colors indicating higher scores

Fig. 7

Lung transcriptome results in the high- and low-altitude groups. A Volcano plot of DEGs. LZ: Lung transcriptome of Zhashijia sheep; HV: Lung transcriptome of Valley Tibetan sheep.The X-axis represents the fold change in the difference after conversion to log₂, and the Y-axis represents the significance value after conversion to log₁₀. Red represents upregulated DEGs, blue represents downregulated DEGs, and gray represents non-DEGs. B Bubble map of enriched KEGG pathways. C Bar graph of enriched GO pathways. D PPI network of DEGs (|log₂fold-change| ≥ 2). The color of the circles represents the number of node connections, with redder colors indicating higher numbers of node connections. The color of the connection lines represents the interaction score, with redder colors indicating higher scores

Fig. 8

Functional enrichment analysis of common differentially expressed genes in selecte sweep regions, the heart transcriptome, and the lung transcriptome. A Bubble map of enriched KEGG pathways. B Bar graphs of enriched GO terms