The study was conducted according to the guidelines of the Institutional Animal Care and Use Committee of Institute of Animal Science and Veterinary Medicine, Chinese Academy of Sciences (IACUC2021311).

Seven Tibetan sheep breeds exhibiting significant geographical and phenotypic variations were selected from Qinghai Province (plateau-type Tibetan sheep, valley-type Tibetan sheep, Black Tibetan sheep, Qumaari Speckled sheep, Zeku sheep, Oula sheep, and Zashgar sheep). A total of 140 blood samples were collected from unrelated adult ewes (n = 20 per breed). Each sample was meticulously documented with the species name, codes, geographical coordinates of the sampling site (longitude, latitude, and altitude), and phenotypic characteristics (see Figure 1a and Supplementary Table 1 for further details). All tissue samples were preserved in 95% ethanol and stored at −80°C for subsequent genomic analysis. The phenotypes of five traits were measured in 20 different breeds of Tibetan sheep, including three body size traits (body weight, body length, and body height) and two wool quality traits (wool fiber length and fiber fineness).

Total genomic DNA was extracted from the samples, and at least 3 μg of genomic DNA was used to construct paired-end libraries with a read length of 2 × 150 bp for paired-end sequencing. Simultaneously, we obtained the farmer’s consent to collect skin samples. We randomly selected three unrelated plateau Tibetan sheep and three unrelated Oula sheep, and collected 2 square centimeters of skin tissue from the shoulder region using sterile surgical blades. The samples were promptly immersed in liquid nitrogen (−80°C) for subsequent transcriptome sequencing and real-time quantitative polymerase chain reaction (RT-qPCR) analysis. These libraries were sequenced utilizing the Illumina NovaSeq 6000 platform located in Shanghai, China.

ASTP Toolkit v0.18.0 was utilized for stringent quality control of the raw reads based on the following criteria: (1) removal of reads containing ≥ 10% unidentified nucleotides (N); (2) elimination of reads with > 50% bases having Phred quality scores ≤ 20; and (3) discarding reads that aligned to the barcode adapter. The Burrows–Wheeler Aligner (BWA) was employed to map the cleaned reads from each sample to the CAU O.aries_1.0 reference genome¹ using the parameters “mem 4 -k 32 -M”. Here, -k specifies the minimum seed length, and -M marks shorter split alignment hits as secondary alignments (11). Variant calling was conducted across all samples using GATK’s Unified Genotyper (12). Subsequently, SNPs and InDels were filtered through GATK’s Variant Filtration tool, applying specific thresholds (-Window 4, -filter “QD <2.0||FS >60.0||MQ <40.0,” -G_filter “GQ <20”), and variants showing segregation distortion or indicative of sequencing errors were excluded (13). For the determination of the physical positions of each SNP, the ANNOVAR software (14) was applied for alignment and annotation of SNPs or InDels. Structural variations (SVs), encompassing translocations, inversions, and insertions, were identified using the breakdancer software (Max1.1.2.) (15). Copy number variants (CNVs) were classified using CNVnator (0.3.2) (16).

The SNP-only dataset was analyzed utilizing a maximum likelihood algorithm. A phylogenetic tree was constructed using PHYML 3.0 (17) and FastTree, based on the selected optimal nucleotide substitution model, generalized time reversible (GTR). Node support was assessed through 1,000 bootstrap replicates. Principal component analysis (PCA) was conducted using SNP data [excluding SNPs with a minor allele frequency (MAF) less than 0.05] via the GCTA software. Individuals were subsequently clustered into distinct subgroups based on the principal components derived from the analysis. The ADMIXTURE software² was utilized to analyze the genetic structure of the population using SNP data. A mixed model was selected, with K ranging from 2 to 10 (assuming 2 to 10 ancestral populations). All other parameters were set to the software’s default values. The optimal K value, which is closest to the true number of ancestral populations, was determined based on the cross-validation (CV) error values.

Select regions with extremely low or high θ, where the π ratio is in the 5% left and right tails, and those with significantly high Fst values (i.e., the top 5% of Fst), as these are identified as regions that have undergone strong selective sweeps.

Use top GO/KEGG to conduct enrichment analysis. In this process, the gene list and gene count for each pathways are derived from the candidate region genes annotated. The p-value is then calculated using the hypergeometric distribution method, with a threshold of p-value <0.05 indicating significant enrichment. This approach identifies the pathways that are significantly enriched in the candidate region genes relative to the whole-genome background, thereby elucidating the primary biological functions of these genes.

Based on inter-population SNPs and linkage disequilibrium, the EMMAX software was utilized to perform association analysis between molecular markers and trait phenotypes,³ identifying markers or candidate genes closely associated with the target traits. These findings were visualized using Manhattan plots and QQ plots.

Total RNA was extracted using the Trizol Reagent (Invitrogen Life Technologies). Subsequently, the concentration, quality, and integrity of the RNA were assessed using a NanoDrop spectrophotometer (Thermo Scientific). To isolate cDNA fragments of the desired 400–500 base pairs in length, the library fragments were purified utilizing the AMPure XP system (Beckman Coulter, Beverly, CA, United States). Products were purified using the AMPure XP system and quantified with the Agilent High Sensitivity DNA Assay on a Bioanalyzer 2100 (Agilent Technologies). Subsequently, the sequencing library was subjected to sequencing on the NovaSeq X Plus platform (Illumina).

Fastp (18) (version 0.18.0) was used to filter adapters or low quality bases from raw reads. The reference genome and gene annotation files were retrieved from the genome database. The filtered reads were aligned to the reference genome using HISAT2 version 2.0.5.

We utilized HTSeq (version 0.9.1) to compare the read count values for each gene as a measure of its original expression, followed by normalization using FPKM. Differential gene expression was analyzed using DESeq (version 1.30.0) with the following screening criteria: an absolute log2 fold change greater than 1 and a significant p-value less than 0.05. Additionally, we employed the R language Pheatmap (version 1.0.8) package to conduct bi-directional clustering analysis of all differentially expressed genes across samples. The heatmap was generated based on the expression levels of the same gene in different samples and the expression patterns of different genes within the same sample, using the Euclidean method to calculate distances and the complete linkage method for clustering.

The SNPs identified through GWAS were mapped to the corresponding genes in the expression dataset to evaluate their impact on gene expression. The significant genes obtained from GWAS were intersected with differentially expressed genes (DEGs), and KEGG/GO enrichment analysis was performed on the overlapping genes. Cluster analysis of the commonly differentially expressed genes was carried out using the Pheatmap package (version 1.0.8) in the R programming language.

Seven DEGs were selected and confirmed by RT-qPCR. The housekeeping gene actin served as an internal control for normalizing mRNA expression levels. Primers were designed using Oligo 6.0 software (Supplementary Table 30). Quantitative real-time PCR was performed in a 20 μL reaction volume containing 10 μL of 2 × Top Green EX-Taq Mix, 2 μL of cDNA, 7 μL of ddH₂O, and 0.5 μL each of forward and reverse primers. The thermocycling conditions were as follows: initial denaturation at 94°C for 30 s; followed by 42 cycles of denaturation at 94°C for 5 s, annealing at 61°C for 35 s; a final extension at 97°C for 10 s, and a dissociation curve analysis stage consisting of 65°C for 60 s and 97°C for 1 s. Relative mRNA expression levels were determined using the 2^−ΔΔCT method. All experiments were conducted with six biological replicates.

Whole-genome resequencing of 140 samples generated 36.4 billion paired-end raw reads with an insert size of 400 bp, and 35.98 billion high-quality reads, achieving an average depth of 14.44 × per sample and an average genome coverage of 99.33%. The GC content varied between 42.77 and 44.30%, the Q20 value was no less than 97%, the Q30 value was at least 93.69%, and the mapping rate against the reference genome surpassed 99.84% (Supplementary Table 2). Furthermore, a total of 46,681,283 SNP sites were identified for purine-pyrimidine transitions and transversions, with a Ts/Tv ratio of 1.86 (Supplementary Table 3). A total of 46.69 million SNPs were identified and utilized for subsequent analyses. The majority of these high-quality SNPs (64.35%) were found in intergenic regions, characterized by T/C and A/G transitions, while only 0.67% were located within exonic regions. The remaining SNPs were distributed as follows: 0.58% upstream, 0.56% downstream, and 33.82% within introns. In total, 314,431 SNPs were detected in exons, of which 46.63% were non-synonymous and 51.09% were synonymous, yielding a non-synonymous to synonymous ratio of 0.913 (Supplementary Table 4). In addition, a total of 1,880,322 insertions and 2,680,046 deletions were identified across seven indigenous sheep breeds (Supplementary Table 5).

The majority of insertions and deletions were situated in intergenic regions (63.71%). Additionally, exon InDels predominantly comprised in-frame deletions or insertions (56.62%), leading to alterations in the reading frame of protein-coding genes, which often exhibited a multiple of the triplet codon length (Supplementary Table 6). After quality control, the following copy number variants (CNVs) were detected: 24,803 in GY, 238,655 in HZ, 24,984 in OL, 24,960 in QM, 23,435 in SG, 22,029 in ZK, and 24,868 in ZS (Supplementary Table 7). A total of 337,491 reliable structural variants (SVs) were identified. Specifically, this included 3,469 insertions, 186,018 deletions, 15,204 inversions, 46,226 intrachromosomal translocations, and 86,574 interchromosomal translocations (Supplementary Table 8).

The observed heterozygosity (Ho) of seven Tibetan sheep populations ranged from 0.358 to 0.372, while the expected heterozygosity (He) ranged from 0.349 to 0.358 (Supplementary Table 9). The observed heterozygosity was consistently higher than the expected heterozygosity, Fis values for ZS, SG, and GY were all less than 0. To gain a more comprehensive understanding of the distribution of chromosomal variations in Tibetan sheep, we established a 1 Mb sliding window and generated density maps for genes, SNPs, CNVs, and InDels within each window (Figure 1b). Notably, the genomic variation exhibits a high degree of homogeneity. The constructed weighted phylogenetic tree elucidates the relationships among seven Tibetan sheep breeds, and the resulting neighbor-joining (NJ) tree provides evidence of distinct separations between breeds, with each breed forming its own branch, although some overlap is observed in Tibetan sheep populations (Figure 1d). On the other hand, PCA cluster analysis offers further empirical support for the delineation of these groups, thereby demonstrating the robust consistency of the selected samples (Figure 1c). The mixed model was employed to assess the genetic structure of the population. At K = 2, the cross-validation error reached its minimum, resulting in the division of 140 germplasm resources into two distinct groups. When K = 3, SG tends to diverge in a different direction from the main group, while the remaining breeds are distributed across the other two clusters. When K = 7, the majority of the genetic information of GY, HZ, OL, QM, ZK, and ZS appears to stem from a shared ancestral population (red). When K = 10, the Tibetan sheep population exhibits extremely rich genetic diversity, which may be attributed to its highly complex population history, extensive geographical distribution, or diverse adaptive traits (Figures 1e,f). Linkage disequilibrium (LD, r²) decreases to half of its maximum value within less than 10 kb, with ZK populations exhibiting the fastest decline and SG populations showing the slowest decline. These findings indicate that SG species exhibit strong LD and are more prone to inbreeding, which may be a result of long-term artificial selection (Figure 1g). The maximum likelihood tree, based on the number of optimal migration events, provided evidence of gene flow between distinct Tibetan sheep populations (Figure 1h).

Through the analysis of genomic regions, we identified loci with high levels of Fst and nucleotide diversity ratio (θπ), which are associated with hypoxic adaptability. Specifically, by comparing varieties below 3,000 meters and those above 4,000 meters, we identified 2,900 and 2,618 candidate positively selected regions in the SG (Fst >0.082, θπ >0.496) (Supplementary Table 9) and HZ (Fst >0.044, θπ >0.365) (Supplementary Table 10) varieties, respectively. These two varieties exhibited 607 and 586 candidate genes, respectively, during variety-specific selection events, with a total of 218 co-selected genes (Figure 2a). These shared genes were selected through critical adaptive responses to hypoxic conditions in high-altitude environments, exhibiting higher levels of differentiation. For instance, several positively selected genes associated with hypoxic adaptation were identified in one or more of the studied populations, including DOCK8, IGF1R, JAK1, SLC4A7, TMTC2, and VPS13A. KEGG enrichment analysis revealed that these genes were predominantly enriched in signal transduction, development and regeneration, and endocrine system (Figure 2b and Supplementary Table 11).

Considering the physical properties of coat color, we employed a comprehensive comparative analysis, evaluating the OL, GY, HZ, ZS, and ZK varieties against the black-coated HZ variety. For OL (Fst >0.037, θπ <−0.305), GY (Fst >0.039, θπ >0.273), QM (Fst >0.039, θπ <−0.250), and ZS (Fst >0.044, θπ <−0.257) and ZK (Fst >0.039, θπ <−0.363), a total of 2,136, 2,625, 2,193, 2,376, and 2,034 selection regions were identified (Supplementary Tables 12–16), encompassing 37 co-selection genes (Figure 2c). Among the shared genes, including TCF25, MITF, and MC1R, which are associated with hair color, these genes are involved in transport and catabolism, regulation of the endocrine system, development and regeneration processes, as well as signaling molecules and interaction pathways (Figure 2d and Supplementary Table 17).

GWAS serve as a crucial tool for elucidating the genetic underpinnings of traits, revealing the relationship between genetic variation and these traits. To account for known confounding factors, we incorporated gender and age as fixed effects into the mixed linear model using a general linear modeling approach. Based on the significance criterion of a p-value with a −log10 value exceeding 4, we identified several significant SNP loci associated with body height (BH), body length (BL), and body weight (BW) at the genomic level (Supplementary Tables 18–20). The p-value Manhattan plot revealed significant associations between SNPs and BH, BL, and BW. Through gene annotation, 30 significant loci related to body shape traits were identified. For the BH trait, five significant SNPs were located on OAR1, while one SNP was found on OAR2, OAR5, OAR21, and OAR6 (Figures 3a,d). For BL, five SNPs showed significant associations, with candidate regions identified on OAR1, OAR3, and OAR23. Notably, the five most significant SNPs were located on OAR2 (Figures 3b, e). For BW, 10 significant SNPs were identified, with the most prominent associations occurring in the EVC and EVC2 regions on OAR6 (Figures 3c,f). These SNPs demonstrated a strong positive correlation with body size traits (Table 1), exceeding the predefined threshold for significance.

Similarly, through GWAS analysis, several key genes associated with wool traits and their genetic variations were identified (Supplementary Tables 21, 22). Further annotation analysis revealed 19 significant loci linked to these traits. For wool fiber length (WL), the two most significant SNPs were located on OAR1, while one SNP was identified on each of OAR26, OAR11, OAR15, OAR24, OAR19, and OAR4, and two SNPs were found on OAR13 (Figures 4a,c). Regarding wool fiber fineness (WF), five SNPs were significantly associated with candidate regions on OAR8, OAR11, and OAR2, with the two most significant SNPs located on OAR6 (Figures 4b,d). We identified a set of significantly correlated SNP loci, with detailed information regarding these significant SNPs presented in Table 2.

To comprehensively assess the characteristics of the candidate genes annotated by SNPs, we conducted enrichment analyses to identify the functional categories associated with these genes. In terms of body size traits, the enriched KEGG pathways are primarily associated with the cytoskeleton in muscle cells, cGMP-PKG signaling pathway, cAMP signaling pathway and hedgehog signaling pathway (Figures 5a–c and Supplementary Table 23). For wool fiber traits, the focus is mainly on the melanogenesis, TNF signaling pathway, Wnt signaling pathway, EGFR tyrosine kinase inhibitor resistance, hedgehog signaling pathway and Rap1 signaling pathway (Figures 5d,e and Supplementary Table 24).

To elucidate the molecular mechanisms underlying the distinctive long hair traits of Tibetan sheep, we conducted transcriptomic analyses on skin samples from two Tibetan sheep breeds, ZY and OL, which exhibit significant differences in wool length. On average, 143,236,464 and 126,767,572 original reads were produced in the ZY and OL cDNA libraries, respectively. The GC contents of ZY and OL were 47.09 to 48.30% and 48.05 to 48.25%, respectively (Supplementary Table 25). Based on the significance threshold (|FC| >2 and p-value <0.05), a total of 1,065 DEGs were identified, including 434 up-regulated and 631 down-regulated genes (Figure 6a). The hierarchical clustering showed that there were significant differences in the expression profiles of DEGs between the treat group (GY) and the control group (OL) (Figure 6b). GO analysis revealed significant enrichment in several biological processes, including fibrillar collagen trimer, melanin metabolic process, peptidase regulator activity, and cellular lipid metabolism process (Figure 6c and Supplementary Table 26). A total of 39 significant pathways were identified through KEGG pathway analysis. The top 20 enriched pathways are presented in Figure 6d, indicating that DEGs (73 up-regulated and 155 down-regulated) were predominantly enriched in ECM-receptor interaction, PPAR signaling pathway, focal adhesion, fatty acid elongation and biosynthesis of unsaturated fatty acids these pathways (Supplementary Table 27). Based on the threshold for DEGs analysis, we identified and screened protein pairs with a node score greater than 0.95 in the STRING database to construct the PPI network (Figure 6e). According to the interaction results, there are four up-regulated gene modules (HOXC6, HOXA6, HOXA7, and HOXB6) and four down-regulated gene modules (COL1A1, COL1A2, COL3A1, and SPARC). Additionally, some modules contain both up-regulated genes (such as ATP6, COX2) and down-regulated genes (such as CYTB). Notably, the HOX homeobox genes, which encodes transcription factors (TFs), exhibited a significantly high degree of connectivity within this network. To further investigate the distribution characteristics of differentially expressed TFs, we conducted a statistical analysis of the number of differentially expressed TFs in each family within the control group. The results indicate that the homeobox family exhibits the most significant differential expression, with 10 upregulated TFs identified, including HOXC6, HOXC8, HOXA6, HOXA10, HOXC9, HOXA1, HOXB6, and HOXC10, as well as six downregulated factors (Figure 6f and Supplementary Table 28). Notably, the Hox subfamily demonstrates the most extensive distribution pattern within this family. This finding highlights the critical role of the HOX homeobox family, particularly its Hox subfamily, in gene regulatory networks.

To further evaluate the expression levels of candidate genes identified through GWAS, we examined the intersection between these candidate genes and DEGs associated with two distinct wool traits. SNPs were identified in critical regions of 10 key overlapping DEGs. For instance, a non-synonymous SNP in the coding region of WNT16 was predicted to alter the protein sequence, whereas synonymous mutations were observed in the coding regions of NES and TIAM1. SNPs were also detected in intergenic regions of PRKCA, MAP3K8, MMP7, and OVOL2, as well as within introns of COL7A1, KDM8, and ZNF385D (Table 3). The expression levels of 10 overlapping DEGs exhibited differential expression and distinct clustering patterns, with 8 genes up-regulated and 2 genes down-regulated (Figure 7a). Subsequently, KEGG enrichment analysis was conducted on these overlapping genes to elucidate their biological functions. The results indicate that these genes are involved in multiple metabolic pathways, including the Wnt signaling pathway, Ras signaling pathway, protein digestion and absorption, TNF signaling pathway, melanogenesis and chemokine signaling pathway (Figure 7b and Supplementary Table 29). Heat maps were utilized to visualize the genotypes of overlapping genes, and the results demonstrated that the polymorphism of the GY genotype was significantly higher compared to that of the OL genotype (Figure 7c).

To validate the accuracy of the DEGs assay at the transcriptome level, we selected seven genes and evaluated their expression levels using six replicates for each sample. The results demonstrate a high degree of consistency between the two methods at the expression level, with a Pearson correlation coefficient of 0.75 (Figure 8). This indicates that the RNA-seq data exhibit strong reliability. The selected candidate genes exhibited strong concordance at both the transcriptional and actual expression levels, further substantiating their potential significance in subsequent functional analyses.