Genetic diversity and population patterns of Tibetan sheep from whole-genome analysis, with insights for conservation and breeding

Blood samples were collected from sheep as approved by Qinghai University Animal Ethics Committee (2020/PJ202401-75). The experiment was conducted in accordance with the Chinese Code of Practice for the Care and Use of Animals for Scientific Purposes (1988). All 270 sheep (90 per breed) were healthy adult ewes aged 2–3 years. This age group was selected because it represents the stable adult stage (avoiding juvenile genetic instability or senescent phenotypic changes), and using only ewes eliminated sex-related genetic variation (e.g., Y-chromosome loci) that could confound population structure analysis. Within each breed, the ewes had consistent body weight ranges: Tibetan sheep (38–45 kg), Hu sheep (42–48 kg), and Small Tail Han sheep (48–55 kg). While absolute weights differed (reflecting inherent breed differences), the narrow weight range within each breed ensured no extreme phenotypic outliers (e.g., undernourished or obese individuals) that could affect genomic data quality.

Tibetan sheep (Gonghe County, Qinghai) were raised under traditional extensive pasture-based management—grazing freely on natural alpine meadows year-round, with no supplementary feeding. Hu sheep and Small Tail Han sheep (Minghe County) were raised under semi-intensive management—housed in open-sided sheds at night and grazing on the pastures during the day, with daily supplementary feeding of 0.5 kg/head of concentrated feed (corn–soybean meal mix).

Tibetan sheep (Gonghe County, Qinghai): Sampled from areas with an average altitude of 3200–3500 m—a representative mid-high altitude zone of the Qinghai–Tibet Plateau, characterized by low oxygen (oxygen partial pressure ~65% of sea level), annual average temperature ~2 °C, and strong ultraviolet radiation (UV-B intensity ~1.5 times that of low-altitude areas). Hu sheep and Small Tail Han sheep (Minghe County, Qinghai): Sampled from areas with an average altitude of 2200–2500 m—the eastern low-altitude fringe of the Qinghai–Tibet Plateau, with milder climate (annual average temperature ~6 °C), higher oxygen partial pressure (~85% of sea level), and weaker UV radiation. The altitude difference between the two sampling regions (≈1000 m) aligns with the core contrast between the Tibetan sheep’s native high-altitude habitat and the introduced low-altitude breeds’ adapted environment, ensuring our genetic comparison could effectively capture altitude-related selection signatures.

Blood samples were collected uniformly in mid-July 2024. About 5 mL blood sample per sheep was collected into EDTA-anticoagulated tube and stored at −20 °C until DNA extraction. Genomic DNA was extracted from the blood samples using the magnetic bead method with CWE9600 Magbead Blood DNA Kit (Kangwei Century Biotechnology Co., Ltd., Taizhou, China) to perform DNA extraction. The quality and concentration of the extracted DNA were determined by agarose gel electrophoresis and a NanoDrop spectrophotometer. DNA samples with an OD260/280 ratio between 1.8 and 2.0 and a concentration above 30 ng/μL were used for library construction. The samples were then sent to Compass Agritechnology Co., Ltd. (Tianjin, China) for whole-genome resequencing.

During the library construction process, genomic DNA with a total amount greater than 1 μg was randomly fragmented into approximately 300–350 bp fragments using the universal library construction kit (for MGI)-NadPrep^® DNA and the Covaris™ fragmentation instrument. After end repair, A-tailing, and adapter ligation, DNA fragments of approximately 300–350 bp were selected using NadPrep^® SP Beads. PCR amplification was then performed, followed by purification of the PCR products using NadPrep^® SP Beads. The sequencing library was ultimately obtained. After the library was constructed, an initial quantitative analysis was conducted using Qubit 2.0. Subsequently, the insert fragments of the library were detected using Bioanalyzer^® (Agilent, Santa Clara, CA, USA), and once the results met expectations, the library was then subjected to sequencing on the DNBSEQ-T7 platform; the PE150 sequencing protocol was selected for sequencing.

Quality Control: To ensure the quality of information analysis, fastp was used to conduct a series of quality control (QC) on the raw reads. Remove reads with adapters; remove paired reads when the proportion of N in the sequencing read exceeds 1% of the total base number of the read; and remove paired reads when the number of low-quality (Q ≤ 5) bases in the sequencing read exceeds 50% of the total base number of the read.

Read Mapping: After data filtering, an index was built based on the provided reference genome. Then, Clean Reads were aligned to the reference genome using the BWA0.7.17 (men) software. The Clean Reads were sorted and indexed using samtools 1.7. The GATK 4.1.8.0 software’s built-in modules were used to remove duplicates from the Bam files. Based on the Bam files, the sequencing depth and genomic coverage of each sample were then statistically analyzed to prepare for subsequent variant detection.

Variant Calling: Single nucleotide polymorphism (SNP) mainly refers to the DNA sequence polymorphism caused by the variation in a single nucleotide at the genomic level, including single base transitions and transversions. HaplotypeCaller module of the software GATK (4.18) was used to perform variant detection on multiple samples with the processed alignment files. The detected variants were filtered using VariantFiltration, and the SNPs were annotated using ANNOVAR (v2016Feb01), which performs gene-based annotation, region-based annotations, filter-based annotation, and other functionalities annotations. To reduce the error rate of SNP detection, SNPs were filtered as follows: QD < 2.0, FS > 60.0, MQ < 40.0, SOR > 3.0, MQRankSum < −12.5, ReadPosRankSum < −8.0; then, SNPs were strictly filtered: SNP cluster filtering (no more than 2 SNPs within 5 bp); SNP filtering near Indel (SNPs within 5 bp of Indel are filtered out); and for sites with GQ (Genotype Quality) less than 20.0, the genotyping quality tag of the sample was marked as low GQ.

Population Structure Analysis: Principal component analysis (PCA) was carried out using PLINK software (2.00a3.3) to visualize the genetic relationships among samples. The population structure was further analyzed using STRUCTURE 2.3 software, with the number of assumed populations (K). The multiple K values were determined based on the given population information (typically, K equals the number of populations), and population structure analysis was conducted using the Admixture 1.3.0 software. The smaller the cv value was, the better the corresponding K value was, and thus the optimal number of clusters (K) was determined to be the number of populations. The genetic distance matrix of IBS was constructed using Plink (version 1.9), and based on this, a phylogenetic tree was constructed using the Neighbor-Joining (NJ) method. From the analysis results, we can roughly infer which samples have more similar genetic characteristics.

Linkage Disequilibrium Analysis: Linkage Disequilibrium (LD) refers to the non-random combination relationship between different loci’ alleles within a population. That is, when the probability of two alleles (A, B) located on the same chromosome existing simultaneously is greater than the probability of their simultaneous occurrence due to random distribution within the population, these two loci are said to be in LD state. Usually, it is represented by the values of D’ and r². Generally speaking, in linkage disequilibrium analysis, the LD value of wild species is relatively low, while that of domesticated species will be larger due to the positive selection effect. This is because in a DNA sequence, there is a linkage relationship between loci. Different loci’ linkages form “haplotypes”, and with the accumulation of recombination, specific haplotypes will be weakened and gradually disappear. Due to the relationship between recombination rate and linkage distance, the intensity of linkage disequilibrium between loci and their surrounding loci continuously decreases as the distance increases. For a newly generated haplotype, since recombination has not had time to disrupt the linkage between loci, the distance of linkage disequilibrium between them is often relatively far. Under neutral conditions, if a haplotype is relatively new, its frequency is often low, while a haplotype with a higher frequency will need a long period of time to reach a relatively high frequency due to the influence of random drift. If a population undergoes positive selection, the surrounding loci linked to favorable loci will rapidly increase in frequency due to the hitchhiking effect. Therefore, the haplotype containing the favorable locus has a higher frequency and, due to the short time it has experienced, also has a longer LD influence range. The PopLDdecay software 3.41 was used to calculate the r² values of linkage disequilibrium for multiple populations. The r² values of SNPs within different distances were calculated, and the attenuation of population LD r² values with the increase in the distance between loci was observed. Generally, the r² values of wild species decay faster than those of local varieties and cultivated varieties.

Genetic Diversity Analysis: Genetic diversity analysis includes the analysis of the proportion of polymorphic markers in populations (P_N), the expected heterozygosity (H_E), the observed heterozygosity (H_O), and nucleotide diversity (P_i). H_O and H_E were analyzed using the methods proposed by Sun et al. [10] and Nei [11]. The Plink software was used to calculate the minimum allele frequency (MAF) for each locus; then, a self-written R script was used to calculate P_N. P_i was calculated using VCFtools 0.1.16 to evaluate the genetic diversity within each population. Polymorphism information content (PIC) is a quantitative measure used to assess the degree of polymorphism in genetic markers, primarily employed to evaluate the effectiveness of genetic markers in linkage analysis. PIC is calculated using the following formula:PIC = 1 − ∑(p_i²) − ∑∑(2p_i²p_j²) where p_i and p_j represent the frequencies of the i-th and j-th alleles, respectively.

Selection Signature Analysis: The fixation index (FST) and genetic diversity (θπ) have been proven to be an extremely effective method for detecting selection elimination regions. Specifically, when exploring functional areas closely related to the living environment, it often leads to strong population evolution. The two methods jointly screen for strong population evolution, facilitating the selection of target genes. Here, the results of the two analysis methods were combined and selected the positions that appeared simultaneously in the top (0.01%) of the FST sliding window analysis and the top (0.01%) of the θπ Ratio analysis. These positions were taken as candidate target sites. The reference genome information of the corresponding species was downloaded from the ENSEMBL planform. The ANNOVAR software was used to annotate the candidate target SNPs onto the corresponding genes. Then, the gene function enrichment analysis of the annotated genes was conducted based on the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.

After initial variant detection and filtering, we implemented quality filtering (to ensure SNP reliability) and functional prioritization (to identify biologically relevant SNPs) to identify high-confidence, functionally relevant SNPs for downstream analysis (e.g., population structure, selection signatures). A total of 51,378 SNPs were screened out (Table S1), including their genomic location, associated gene, variant type, and functional impact. These SNPs are linked to core traits: to high-altitude adaptation (e.g., HDAC5, BMP2) and economic traits (e.g., SLC27A2 for meat quality, KRTAP13-1 for wool).

Genetic diversity analysis includes the analysis of the proportion of polymorphic markers (P_n) in the population, the expected heterozygosity (H_e) of the population, the observed heterozygosity (H_O) analysis, and nucleotide diversity (P_i) and polymorphism information content (PIC). The results of genetic diversity analysis for each population are presented in Table 1. The results indicated that Tibetan sheep (ZY) had an average Pn of 0.6399, which was not significantly different from Small Tail Han sheep (XWHY, 0.5975) or Hu sheep (HY, 0.638) (p > 0.05), indicating similar levels of polymorphic loci across breeds. Tibetan sheep exhibited the lowest He (0.2489) and Pi (0.2533) among the three breeds. Specifically, He of ZY was 8.29% lower than XWHY (0.2714) and 4.56% lower than HY (0.2608), while Pi of ZY was 8.22% lower than XWHY (0.2760) and 4.60% lower than HY (0.2655) (p < 0.05 for all comparisons). The polymorphic information content (PIC) of SNP loci in Small Tail Han sheep, Hu sheep, and Tibetan sheep were 0.1661, 0.1725 and 0.1731, respectively, indicating low polymorphism (PIC < 0.25) in all three populations.

Principal component analysis (PCA) revealed that the first three principal components (PCA1, PCA2, and PCA3) could explain for 20.69%, 12.26%, and 14.18% of the genetic variation, respectively (Figure 1A–C). Hu sheep and Tibetan sheep were relatively independent and clustered together, respectively. In contrast, Small Tail Han sheep exhibited more dispersed, indicating greater heterogeneity within this breed. This pattern suggests that after being introduced to the region, some Small Tail Han sheep individuals may have undergone hybridization, resulting in non-purebred descendants. A population phylogenetic tree was constructed by the Neighbor-Joining method using MEGA 12 software, and the result was shown in Figure 1D. This result was consistent with the PCA results. The three groups were clustered into groups, and the groups were relatively independent. The ZY clade was the most compact, reflecting low within-breed genetic diversity; the XWHY clade was the most divergent (consistent with PCA results); and the HY clade was intermediate in compactness. No inter-clade mixing was observed. Population structure analysis revealed that when K = 2, the population structure was best explained, with clear genetic admixture patterns among different populations. ZY and HY populations showed unique genetic components, while XWHY and the others had a degree of gene exchange (Figure 1E). The LD decay analysis Figure 1F) revealed the genetic diversity of different populations and the rate of decay for each population. Populations with rapid decay generally have higher genetic diversity, while those with slow decay generally have lower genetic diversity. In this study, the decay of ZY was relatively fast, indicating relatively higher genetic diversity, while the decay of HY was relatively slow, indicating relatively lower genetic diversity.

A total of 713 genomic regions were identified as potential regions under positive selection. After annotating 713 candidate regions, a total of 207 genes were obtained (Table S2). Functional enrichment analysis of the genes located in these regions showed significant enrichment in several biological processes (Figure 2A–F), revealing significant enrichment in pathways such as ascorbate and aldarate metabolism, Pentose and glucoronate interconversions, retinol metabolism, chemical carcinogenesis—receptor activation, bile secretion, porphyrin metabolism, drug metabolism—cytochrome P450, chemical carcinogenesis—DNA adducts, metabolism of xenobiotics by cytochrome P450, steroid hormone biosynthesis, drug metabolism—other enzymes, oxytocin signaling pathway, cGMP-PKG signaling pathway, calcium signaling pathway, long-term potentiation, osteoclast differentiation, apelin signaling pathway, biosynthesis of cofactors, gap junction, cellular senescence, circadian rhythm, C-type lectin receptor signaling pathway, MAPK signaling pathway, type II diabetes mellitus, growth hormone synthesis, secretion and action, and oocyte meiosis. Functional enrichment analysis indicated that these candidate genes mainly play crucial roles in biological processes and cellular components. The enrichment results showed that the candidate genes are primarily involved in energy and antioxidation, blood vessels and circulation, hormones and growth, detoxification and metabolism, immunity and bone metabolism, and nervous system and reproduction. The genomic analysis revealed strong selection signals in key genes associated with hypoxia adaptation across multiple biological pathways. Within the hypoxia sensing and HIF signaling pathway, HDAC5 was identified as stabilizing HIF-1α through deacetylation, thereby upregulating EPO and other target genes, while BMP2 demonstrated angiogenic properties via VEGF pathway activation, enhancing skeletal muscle capillary density. BMPR1B emerged as a crucial mediator of BMP2 signaling, regulating pulmonary artery smooth muscle proliferation to mitigate high-altitude pulmonary hypertension. In the reprogramming of energy metabolism, SLC25A46, a mitochondrial carrier protein that enhances fatty acid β-oxidation efficiency under hypoxic conditions, regulates cellular metabolic adaptations. NDUFS7, a critical subunit of ETC Complex I, maintains ATP synthesis stability during oxygen deprivation, and COQ5 is the coenzyme Q10 synthase that reduces ROS leakage while improving ETC efficiency. Vascular and hematological adaptation mechanisms were reflected in LTBP2’s regulation of TGF-β-mediated vascular remodeling to increase pulmonary artery elasticity, preventing high-altitude pulmonary hypertension, while GJA3 maintained gap junction communication to coordinate cardiomyocyte electrophysiological synchronization under hypoxia, showing strong association with cardiac hypoxia tolerance. In the dual oxidase system, DUOX2 played a pivotal role in balancing hypoxia-induced ROS production between signaling and oxidative damage. The oxidative stress defense system featured by the DUOX2/DUOXA2 complex regulates H₂O₂ production as a key component of Tibetan sheep’s antioxidant capacity, alongside SORD’s maintenance of NAD⁺/NADH redox balance through its sorbitol dehydrogenase activity. Concurrently, the study identified multiple candidate genes linked to economically important traits. Growth-related genes included BMP2 and BMPR1B regulating skeletal muscle differentiation and bone development, FGF9 promoting muscle precursor cell proliferation, IGF2BP2 stabilizing IGF2 mRNA to enhance growth, and MYCBP2 driving muscle fiber hypertrophy through cell cycle regulation. Meat quality traits were associated with SLC27A2’s modulation of intramuscular fat deposition, LIPH’s influence on flavor compound release through phospholipid metabolism, PDGFD’s promotion of adipocyte differentiation affecting subcutaneous fat thickness, and HMBS’s role in meat color stability via heme synthesis. Wool production characteristics involved KRTAP family genes determining fiber strength and fineness, LTBP2 affecting wool crimp and elasticity through collagen regulation, and GJA3 participating in follicular cell communication during fiber growth cycles. Reproductive performance was linked to BMPR1B’s effect on ovarian granulosa cell function through the FecB mutation, the HOXA family’s (particularly HOXA10) regulation of uterine receptivity, IZUMO1R’s mediation of sperm–egg binding, and CYP19’s control of estrous cycle regularity via estrogen synthesis catalysis. These candidate genes involved in the signaling pathways collectively form the comprehensive adaptation mechanism and production performance of Tibetan sheep in extreme environments such as high-altitude hypoxia, intense ultraviolet radiation, and cold conditions.