Population Genomic History of the Endangered Anatolian and Cyprian Mouflons in Relation to Worldwide Wild, Feral, and Domestic Sheep Lineages

Abstract

Once widespread in their homelands, the Anatolian mouflon (Ovis gmelini anatolica) and the Cyprian mouflon (Ovis gmelini ophion) were driven to near extinction during the 20th century and are currently listed as endangered populations by the International Union for Conservation of Nature. While the exact origins of these lineages remain unclear, they have been suggested to be close relatives of domestic sheep or remnants of proto-domestic sheep. Here, we study whole genome sequences of n = 5 Anatolian mouflons and n = 10 Cyprian mouflons in terms of population history and diversity, comparing them with eight other extant sheep lineages. We find reciprocal genetic affinity between Anatolian and Cyprian mouflons and domestic sheep, higher than all other studied wild sheep genomes, including the Iranian mouflon (O. gmelini). Studying diversity indices, we detect a considerable load of short runs of homozygosity blocks (<2 Mb) in both Anatolian and Cyprian mouflons, reflecting small effective population size (Ne). Meanwhile, Ne and mutation load estimates are lower in Cyprian compared with Anatolian mouflons, suggesting the purging of recessive deleterious variants in Cyprian sheep under a small long-term Ne, possibly attributable to founder effects, island isolation, introgression from domestic lineages, or differences in their bottleneck dynamics. Expanding our analyses to worldwide wild and feral Ovis genomes, we observe varying viability metrics among different lineages and a limited consistency between viability metrics and International Union for Conservation of Nature conservation status. Factors such as recent inbreeding, introgression, and unique population dynamics may have contributed to the observed disparities.

Keywords: conservation; domestication; genomics; sheep.

Fig. 1.

The geographic distribution and phylogeny of the studied sheep lineages. A) Geographic distribution of the sheep species considered wild by IUCN (data from IUCN). The distribution ranges of the CYM and ANM are shown separately as points. The two points with the additional symbols in the middle denote the sampling locations. The primary populations of the feral-considered EUM are also shown separately as points in the smaller panel. B) MDS analysis of the studied sheep lineages, using 1-outgroup-f3 statistics as distance proxy. C) NJ tree of the studied sheep lineages, using (1-outgroup-f3) as distance proxies and using goat as outgroup. Bootstrap support was calculated using 500 replicates and all branches have 100% support.

Fig. 2.

Population size changes among sheep lineages. PSMC analysis of high-coverage individuals from each lineage, A) for mouflons and domestic sheep and B) for N. American, Siberian, and Asian wild sheep. The x axis shows time in a log scale, and the y axis shows the estimated effective population size. We assumed a generation time of 3 years and a mutation rate of 1.5 × 10⁻⁸.

Fig. 3.

Heterozygosity and diversity estimates among sheep lineages. A) Genome-wide heterozygosity values estimated using genotype likelihoods. Only the high-coverage genomes cym008 and oga018 were included to represent CYM and ANM, respectively. B) Within-population autosomal diversity values estimated using pairwise 1-outgroup-f3 statistics per lineage. C) Comparison of autosomal versus chrX diversities, each estimated using pairwise 1-outgroup-f3 statistics. The regression line was generated with the loess algorithm in the R stats package.

Fig. 4.

ROH in sheep genomes. A) Number of ROH segments >500 kb plotted against the total length of the segments found in each individual. The Anatolian and Cyprian populations are only represented by the high-coverage individuals oga018 and cym008, respectively. B) Size distribution of ROH segments divided into four classes (0.5 to 1, 1 to 2, 2 to 3, and 3 to 5 Mb). Inbreeding times corresponding to each size class were estimated assuming a generation time of 3 years and a recombination rate of 1.5 cM/Mb (Thompson 2013; Kardos et al. 2018). The x axis is given in log scale. C) Proportion of ROH segments >500 kb (F_ROH) in each individual’s genome.

Fig. 5.

Mutation load estimates using GERP scores. RMLs were calculated as average GERP scores weighted by the number of derived variants, only including variants found in conserved regions (GERP > 4) (von Seth et al. 2021). We used goat alleles to infer the derived state. Only the high-coverage individuals oga018 and cym008 from the Anatolian and Cyprian populations, respectively, were included.

Fig. 6.

Coevaluation of genetic viability metrics and IUCN status. A) Correlations between viability metrics per individual compared across sheep lineages assessed by the IUCN. Regression lines were calculated using the method “loess” in the R stats package. B) The IUCN status compared with the genetic viability metrics heterozygosity (π), F_ROH, and RML.