Genomic DNA was extracted using the EasyPure Blood Genomic DNA Kit (TransGen Biotech) from blood samples collected of Northeast Merino rams (NMS, n = 20) from Jilin Qianyang Agriculture and Animal Husbandry Co., Ltd.(Songyuan City, Jilin Province, China). For each individual, 2 × 150 bp paired-end read data were sequenced using DNBSEQ-T7 at Novogene Bioinformatics Institute company (Beijing, China) (Supplementary Table S1). In addition, to better study the population structure and selection signals of NMS, WGS data for 177 published sheep of 11 breeds were obtained from the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/). Including South African Meat Merino (SAM, n = 10), Australian Merino (AMS, n = 25), Rambouillet (RAM, n = 10), Chinese Merino (CMS, n = 20), Dorset (DOR, n = 24), Suffolk (SUF, n = 12), and several Chinese local breeds: Hu (HUS, n = 10), Small-tailed Han (STH, n = 21), Altay (ALT, n = 10), Tibetan (TIB, n = 12), and Mongolian (MON, n = 23) (Supplementary Table S2).

Burrows-Wheeler Aligner (BWA) software (v0.7.13) was used for mapping clean reads from all 197 sheep to the Ovis aries reference genome Oar_rambouillet_v1.0 (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_002742125.1/) using ‘bwa mem’ parameters (Li and Durbin, 2009). And Picard MarkDuplicates tool (v1.115) (https://github.com/broadinstitute/picard) was used to remove duplicate reads from each alignment. GATK (v4.1.4) was used for SNP calling and then the results were filtered with GATK’s “VariantFiltration” module (McKenna et al., 2010). The filtering parameters were ‘QD < 2.0 ||FS > 60.0 || MQ < 40.0 || SOR > 3.0 || MQRankSum < −12.5 ||ReadPosRankSum < −8.0'. In addition, bcftools (v1.8) was used to extract the biallelic loci located on the autosomes from the hard filtered results, and PLINK software (v1.9) for quality control with parameters’--geno 0.05 --mind 0.1 --maf 0.03’ (Purcell et al., 2007; Danecek et al., 2021).

Based on the annotation file of the Oar_rambouillet_v1.0 reference genome, the types of each SNP were annotated by SnpEff software (v5.1d) (Cingolani et al., 2012). Using previously reported methods, the genes with more than five NMS-specific non-synonymous variations were further extracted for analysis (Kawahara-Miki et al., 2011). To better understand the function of these genes, DAVID was used for online Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis (Huang da et al., 2009; Sherman et al., 2022). Pathways with a p-value less than 0.05 were considered significantly enriched.

Nucleotide diversity (pi), expected heterozygosity (H _E), observed heterozygosity (H _O), linkage disequilibrium (LD) decay, and the runs of homozygosity (ROH) for all breeds were calculated and analysed to explore the genomic genetic diversity of NMS. The values of pi were calculated by VCFtools software (v0.1.16) with the parameters ’--window-pi 10,000 --window-pi-step 5,000' (Danecek et al., 2011). PLINK software was used to calculate H _O and H _E with the ‘--hardy’ parameter. The number and length of ROH for each individual were calculated with the ‘--homozyg-density 10 --homozyg-gap 100 --homozyg-kb 100 --homozyg-snp 10 --homozyg-window-het 1 --homozyg-window-missing 5 --homozyg-window-snp 50 --homozyg-window-threshold 0.05’ parameter of PLINK software (v1.9). Based on the analysis of ROH, the genomic inbreeding coefficient was calculated using the following formula: F _ROH = ∑L _ROH/L _AUTO, where L _ROH is the total length of ROH fragments per individual and L _AUTO is the total length of autosomes covered by SNPs sequenced across the genome. LD decay was calculated using the PopLDdecay software (v3.42) with the ‘-MaxDist 1,000’ parameter (Zhang et al., 2019).

After quality control, PLINK software with the ‘--indep-pairwise 50 5 0.2’ parameter was used to remove high LD sites in the dataset. The obtained sites were used for population structure analysis. GCTA software (v1.92.3) was used to perform principal component analysis (PCA) with the parameter ‘grm’ (Yang et al., 2013). Using the pairwise genetic distances matrix calculated by PLINK, a phylogenetic tree was constructed based on the neighbor-joining (NJ) model by MEGAX and visualised with iTOL (https://itol.embl.de/, accessed on 21 June 2023) (Kumar et al., 2018; Letunic and Bork, 2021). ADMIXTURE software (v1.3.0) was used to infer ancestral populations with K = 2–6. For each K, the software was run 10 times with random number seeds and chose the result with the lowest average cross-validation (CV) error (Alexander et al., 2009).

Two methods, integrated haplotype score (iHS) and composite likelihood ratio (CLR) were used to detect regions of the genome subject to selection in the NMS population. For SNPs detected in the NMS population, BEAGLE (v5.4) was used for imputing and phasing genotypes (Browning and Browning, 2009), and selscan software (v1.2) was used to calculate iHS (Szpiech and Hernandez, 2014). The results were normalised by the norm module of selscan with a window size of 50 kb, the final score for each window is calculated based on the number of SNPs in the window with an iHS score greater than 2 in absolute values. In conclusion, the top 5% of the windows with the highest final scores were retained as the candidate area subject to selection. CLR was calculated by SweeD software (v4.0.0) within a non-overlapping 50 kb window, the top 5% of windows with the highest CLR values are regarded as candidate selected areas (Pavlidis et al., 2013). Only candidate regions that were determined by both methods were considered to be under positive selection. To study the unique selection signals of NMS in recent years and their genetic divergence from other breeds, VCFtools software (v0.1.16) was used to calculate the Fixation Index (F _ST) between NMS and two other breeds with a non-overlapping 50 kb window. The MON is an established parental breed of NMS, while the SAM was introduced to improve the meat performance of NMS (Yang et al., 2017). The 5% of the window with the highest value in the F _ST calculation was considered candidate areas, where regions detected by all three methods were considered to be regions with selection differences between breeds. The results of the selected signals are analyzed using SnpEff software (v5.1d) and DAVID in the same way as above.

The sheep quantitative trait locus database (Sheep QTLdb) contains previously reported regions of QTLs and association data associated with important production traits in sheep (Hu et al., 2022). To calculate the main function of the selected area, the results of the selected signals were compared with the sheep QTLdb (published 25 April 2023).

Using the DNBSEQ-T7 platform, 682.4 Gb of raw data was obtained from 20 NMS individuals and the details of the sequencing data are given in Supplementary Table S1. After filtering, 673.01 Gb clean data were retained, and individual genomes of NMS were generated with an average depth of ∼11.90×. After quality control, 27,770,572 high-quality autosomal biallelic SNPs were obtained. In brief, there were 46,860,074 SNPs before quality control, of which a total of 2,423,764 SNPs were removed using the ‘geno 0.05’ parameter, no individuals were removed due to ‘mind 0.1’ and 16,665,738 SNPs removed by ‘maf 0.03’. Among the remaining SNPs, there are a total of 19,833,469 transitions (Ts) and 7,937,103 transversions (Tv) in all SNPs, with a Ts/Tv ratio of 2.50.

In addition, a total of 23,975,257 high-quality SNPs in 20 NMS were detected. Most of the variants were located in intergenic (54.33%) and intronic regions (35.61%), and only 0.62% (including 47,910 non-synonymous variants and 99,926 synonymous variants) were located in exons (Supplementary Table S3). Of the NMS-specific 1,975,256 SNPs, 53.29% and 34.82% of the variants were located in intergenic and intron regions, respectively. The numbers of non-synonymous and synonymous variants were 8,673 and 11,383, respectively, accounting for 0.44% and 0.58% of the total (Supplementary Table S4).

Non-synonymous SNPs specific to NMS were annotated using SnpEff software, resulting in 2,864 genes. Of these, 350 genes containing more than five non-synonymous variants were selected for enrichment analysis. A total of 11 GO terms were significantly enriched (p < 0.05), of which the most significant (p = 0.002653) GO term was “calcium ion binding, GO:0005509”, containing 13 genes. Several biological process terms were related to immunity, such as “antigen processing and presentation of peptide or polysaccharide antigens via MHC class II, GO:0002504”, “antigen processing and presentation, GO:0019882” (Supplementary Table S5). In addition, “heat shock protein binding, GO:0031072”, a molecular function term related to heat shock protein binding activity was enriched (Yue et al., 2020). For KEGG, 13 pathways were significantly enriched (p < 0.05), of which the most significant (p = 2.80E-09) was “Graft-versus-host disease, oas05332”, associated with immune response. Moreover, several pathways related to immunity and disease, such as “Antigen processing and presentation, oas04612”, and “Phagosome, oas04145” were also enriched. Notably, the pathway “Hippo signalling pathway - multiple species, oas04392”, which is associated with a wide range of important production traits was enriched (Supplementary Table S5) (Deng et al., 2016; Deng et al., 2019; Yatsenko et al., 2020; Dos Santos et al., 2022).

To investigate the population relationship between NMS with other breeds, the 2,330,031 SNPs after Linkage pruning, was used for admixture analysis, phylogenetic analysis and PCA.

The results of ADMIXTURE showed that when K = 2, the ancestors of the China local breeds had a single genomic composition, while Merino breeds showed a mixed ancestral component except for CMS. And when K = 5, NMS displayed clear evidence of shared genome ancestry with the Merino breeds (average 55.00%) and China local breeds (average 30.11%) (Figure 1A). For PCA, the genetic data variation was explained by 4.76% and 4.07% of the first and second principal components, respectively. The results showed that Chinese local breeds and the three Merino breeds (RAM, SAM, and AMS) were clustered separately, with CMS forming a distinct cluster. And NMS was located between Chinese local breeds and the three Merino breeds (Figure 1B). The NJ tree also revealed the same pattern of NMS located between Chinese local breeds and Merino breeds (Figure 1C).

To compare the distribution of ROH fragments in different breeds, ROH fragments were classified into five categories according to their length (0–0.5 Mb, 0.5–1 Mb, 1 Mb–2 Mb, 2 Mb–3 Mb, >3 Mb). Most ROH lengths were in the range of 0–0.5 Mb in all breeds, with only TIB, SUF, RAM, and SMA detecting ROH fragments greater than 3 Mb in length (Supplementary Table S6).

The total ROH length of NMS is medium, lower than the four Merino breeds and the two commercial breeds which have been subjected to stronger selection pressure. F _ROH results showed that NMS had a low inbreeding coefficient (0.095732) and ranked 10th in 12 breeds, higher than ALT (0.088212) and STH (0.094908) only, while RAM had the highest inbreeding coefficient (0.213313). (Figure 2A, Supplementary Table S7). In terms of nucleotide diversity, NMS was ranked 4th (0.002971) behind HUS (0.002979), TIB (0.002988), and AMS (0.003066), while SAM had the lowest pi value (0.002590) (Figure 2B and Supplementary Table S7). Similarly, NMS also exhibited a high level of heterozygosity with the H _O (0.262932) and H _E (0.267285) values, ranking fifth and second, respectively (Supplementary Table S7, Figure 2C). Regarding LD decay, the results were similar to those of F _ROH. The r² values of all breeds decreased rapidly with increasing genomic distance, with the fastest decrease in the first 50 kb. For the distance between markers that was greater than 50 kb, the results showed that NMS had a low genome-wide LD, ranking ninth out of 12 breeds, with MON showing the lowest LD and SAM showing the highest (Figure 2D).

Selected regions on the NMS genome were examined using both integrated haplotype score (iHS) and composite likelihood ratio (CLR) methods, and the top 5% of each method was extracted for annotation. A total of 2,039 (Figure 3A, Supplementary Table S8) and 2,221 (Figure 3B, Supplementary Table S9) genes were annotated by iHS and CLR, respectively, and 14.09 Mb of chromosomal regions containing 287 genes were detected by both methods (Supplementary Table S10). Overall, a range of candidate genes associated with economical traits was subject to positive selection, such as wool growth and type regulation (IRF2BP2, GLI2, and PKIG) (Harmon and Nevins, 1997; Mill et al., 2003; Demars et al., 2017; Zhang et al., 2020; Lv et al., 2022), hair follicles development (RBM28, MAP3K7, and WNT3) (Millar et al., 1999; Sayama et al., 2006; Sayama et al., 2010; Warshauer et al., 2015), growth and development-related (NDUFA9, SETBP1, and ZBTB38) (Liu et al., 2010; Khansefid et al., 2018; Zlobin et al., 2021), meat quality and fat deposition (CTCF, PARP4, and USP25) (Rouleau et al., 2010; Hamill et al., 2012; Yuan et al., 2022), cold resistance (DNAJC13, LPGAT1, and PRDM16) (Shi and Manley, 2007; Seale et al., 2008; Lynes et al., 2018; Liu et al., 2022a), and immune-related (PRDM2, GALNT8, and HCAR2) (Al Kalaldeh et al., 2019; Wang et al., 2019; Ghoreishifar et al., 2020), etc. Subsequently, GO and KEGG enrichment analyses were performed on the 287 genes using DAVID. The results showed that KEGG was enriched to only two pathways, but neither was significantly enriched. One of them was “NOD-like receptor signaling pathway, oas04621” (p = 0.052668), containing 6 genes (NLRP12, DNM1L, MAP3K7, LOC101103623, LOC101103376, and LOC101105481) and is associated with immunity (Van Gorp et al., 2014). The other was “Yersinia infection, oas05135”(p = 0.078788), containing 5 genes (MAP3K7, PXN, DOCK1, LOC101103376, and LOC101103623) and is also associated with immunity (Supplementary Table S10). Regarding GO enrichment analysis results, the most significant term (p = 0.005476) was “acute-phase response, GO:0006953” (LOC101120204, LOC101120613, and LOC105601867), which is related to immunity (Pannen and Robotham, 1995). In addition, another two highly significant enrichment (p < 0.01) terms may be related to meat quality. “actomyosin structure organization, GO:0031032” contains 3 genes (EPB41L4B, CDC42BPB, and CDC42BPA), which is related to composition and disassembly of structures made up of actin and myosin or paracrine. Among them, CDC42BPB may affect the synthesis of 3-hydroxybutyric acid and thus the quality of meat (Li et al., 2022). And “high-density lipoprotein particle, GO:0034364” (LOC101120204, LOC101120613, and LOC105601867) (Supplementary Table S11), which is a cellular component term involved in the transport of lipids (Zhao et al., 2022b).

Annotation of candidate regions detected by both methods using the Sheep QTLdb to determine the function of selected genomic regions in the NMS population. The results showed that 361 QTLs were detected in 420 non-overlapping candidate regions (Supplementary Table S12). The largest proportion of these QTLs was related to meat and carcass, with 113 QTLs (31.30%) distributed across 290 candidate regions. Health-related QTLs followed, with a total of 73 QTLs (20.22%) detected in 217 regions. In addition, 19 wool-related QTLs (5.26%) were also detected, distributed across 135 regions. This suggested that NMS was strongly selected for meat and wool production traits during the breeding process.

Based on F _ST, the selection signals between NMS and other breeds were explored. For MON, 2,389 genes (Supplementary Table S13) were annotated and 3.75 Mb of chromosomal regions contained 87 genes (Supplementary Table S14) that were also detected by iHS and CLR. In contrast to within-population selection signals, the genes detected were mainly related to hair production traits and meat-production traits, such as IRF2BP2 (Lv et al., 2022), SETBP1 (Khansefid et al., 2018), and LNX2 (Santana et al., 2015). In addition, PRDM16 (Liu et al., 2022a), which is related to cold resistance, was also detected. Enrichment analyses were performed of these 87 genes using DAVID and 14 terms and 12 pathways were significantly enriched (Supplementary Table S15). The two most significant terms were “D-threo-aldose 1-dehydrogenase activity, GO:0047834” (p = 0.000007) and “synaptic transmission, glutamatergic, GO:0035249” (p = 0.001591), related to lipid accumulation (Poorinmohammad et al., 2022) and wool colour, respectively (Zhang et al., 2023). The most significant pathway was the “Citrate cycle (TCA cycle), oas00020” (p = 0.002265), which contains five genes related to metabolism. In the F _ST results between NMS and SAM, 2,103 genes (Supplementary Table S16) were annotated in the candidate regions. Among them, 1.58 Mb of chromosomal regions contained 23 genes (Supplementary Table S17) that were detected by all three methods. The annotated genes were mainly divided into two categories: immunity and disease resistance, such as GALNT8(Al Kalaldeh et al., 2019), and RASAL2 (Mesure et al., 2010); and body size, such as NDUFA9 (Khansefid et al., 2018), ADGRD1 (Fischer et al., 2016). Enrichment analysis was performed on these 23 genes and one term and three pathways that were significantly enriched (Supplementary Table S18). The most significant pathway was “Chemical carcinogenesis - reactive oxygen species, oas05208” (p = 0.020840), which might be related to environmental adaptation (Li et al., 2017). The only term was “D-threo-aldose 1-dehydrogenase activity, GO:0047834” (p = 0.024574), which was also significantly enriched in the selection signal results with MON.