Genomic analyses of Asiatic Mouflon in Iran provide insights into the domestication and evolution of sheep

Abstract

Background: Asiatic mouflon (Ovis gmelini) consists of several subspecies mainly distributed in Armenia, southern Azerbaijan, Cyprus, northern, southern, and western regions of Iran, and eastern and central regions of Turkey nowadays. Genome analyses of Asiatic mouflon in Iran revealed that they could have diverged from the direct ancestor of domestic sheep, and showed genetic introgression into domestic sheep after domestication. However, the impact of the Asiatic mouflon subspecies in Iran on sheep domestication remains unclear.

Results: Here, we conducted a comprehensive population genomics analysis of Asiatic mouflon in Iran with 788 whole-genome sequences (including 40 from Asiatic mouflon), 1104 whole mitogenomes (105 from Asiatic mouflon), and 239 Y chromosomes (21 from Asiatic mouflon). Whole-genome sequence analyses revealed two subpopulations of Asiatic mouflon in Iran: O. gmelini_2 limited on Kaboodan Island in Urmia Lake National Park and O. gmelini_1 over a wide geographic area. Phylogenetic analyses of Asiatic mouflon in Iran based on uniparental variants revealed a monophyletic lineage with the mitochondrial haplogroups C/E, and clustered into a monophyletic with Y-chromosomal lineage HY2 of sheep. Additionally, introgression tests detected significant signals of genetic introgression from O. gmelini_2 to four sheep populations (e.g., Garut, Bangladeshi, Nellore, and Sumatra) in South and Southeast Asia. In the four sheep populations, selective tests and introgression signals revealed that the wild introgression could have contributed to their body size, fat metabolism and local adaptation to the hot and humid environments in the Indian Peninsula.

Conclusions: Our results clarified subpopulation structure of Asiatic mouflon in Iran, identifying two distinct groups: O. gmelini_1 and O. gmelini_2. Additionally, we suggest a potential genetic contribution to domestic sheep by introgression, with maternal haplogroup C and paternal lineage HY2 likely originating from the Asiatic mouflon populations in Iran. Our findings offer new insights into domestication of sheep and subsequent introgressions events from wild relatives to domestic populations.

Fig. 1

Population genetics structure of Asiatic mouflon in Iran (O. gmelini). a Geographic distribution of Asiatic mouflon used in this study. Purple, light blue, green and dark blue dots indicate the location of O. gmelini_1, O. gmelini_2, Cyprian mouflon and Anatolian mouflon; b Principal components 1 and 2 for Asiatic mouflon; c Neighbor-joining (NJ) trees of Asiatic mouflon based on whole-genome SNPs, using bighorn sheep as outgroup; d Population genetic structure of the Asiatic mouflon and domestic sheep inferred from the program the sNMF v1.2 (K = 7) using whole-genome SNPs

Fig. 2

Phylogeny of Asiatic mouflon, urial, and domestic sheep inferred from mitogenomes and Y-chromosomal genetic variants. a Phylogeny of Asiatic mouflon, urial, and domestic sheep inferred from 767 mitogenomes; b Phylogeny of Asiatic mouflon, urial, and domestic sheep inferred from 239 Y-chromosomal SNPs; c Network of mitogenome haplotypes; d Network of Y-chromosome haplotypes. Bootstrapping values above 50% were labeled on the phylogenetic tree. Bighorn sheep were used as the outgroup in phylogenetic analyses based on mitogenome and Y-chromosome genetic variants

Fig. 3

Demographic history of Iranian Asiatic mouflon and domestic sheep. a Effective population sizes (N_e) across six domestic sheep populations and the two subpopulations of Asiatic mouflon in Iran inferred from SMC++; b Inference of split times between O. gmelini_1 and the other wild and domestic sheep populations. The vertical lines represent the split time between each pair of populations.; c Effective population size inference of O. gmelini_1 and O. gmelini_2 using GONe; d Maximum Likelihood Network Orientation reconstruction for the relationship between O. gmelini_1 and O. gmelini_2 using OrientAGraph

Fig. 4

The phylogenetic relationship among Asiatic mouflon subspecies and domestic sheep. a Geographical distribution of Asiatic mouflon and samples information of 32 Asiatic mouflon in Iran and 301 domestic sheep samples; b Phylogenetic tree inferred based on complete Cyt-b sequences using the maximum likelihood algorithms. Bootstrap values above 50% were labeled on the phylogenetic tree. c Network of Asiatic mouflon in Iran and Turkey and domestic sheep lineage of C and E built from complete Cyt-b sequences

Fig. 5

Genetic introgression of Asiatic mouflon in SSA domestic populations. a Admixture proportions of Asiatic mouflon in Iran and SSA domestic populations inferred from the program the sNMF v1.2 (K = 3) using whole-genome SNPs. b Heatmap of D-statistics and f₄-ratio calculated by Dsuite v0.5 applied to data including O. gmelini_1, O. gmelini_2, and SSA domestic populations; c Venn diagram for the shared introgression regions among four SSA domestic populations; d Manhattan plot showing signals of selection in the four introgressed SSA populations (BGE/NES/GUR/SUM) detected by the PBS with the dashed line marking 99.9% threshold. The sheep QTLs overlapped with the top 1% PBS regions are shown, e PBS and f_dm values around the significant introgressed genomic region on Chr.10. The selected introgression region (SIR), QTL, and adjacent gene PCDH9 were labeled with coral, purple, and green bars, respectively

Fig. 6

Candidate gene PCDH9 associated with the significant introgression signal on Chr. 10. a Linkage disequilibrium (LD) blocks in the significant introgression region and gene PCDH9 (chr10:41,104,271-43,073,231). The purple and pink bars represent the location of the significant introgression region and gene PCDH9; b Pattern of LD decay in the introgressed populations. The red line is calculated from the whole 10 chromosomes, the green line is calculated from the beginning of the significant selected introgression region (SIR) (chr10:41,104,271) to the next 800 kb; c D_xy and F_ST around the significant introgression region on Chr. 10 (chr10:41,104,271–41,313,804), calculated by pixy v1.2.7.beta1; d Haplotype pattern of the significant selected introgression region (SIR) on Chr. 10 (chr10:41,204,271–41,213,804) (left) and exons of PCDH9 (right)