Contrasting Patterns of Genomic Diversity Reveal Accelerated Genetic Drift but Reduced Directional Selection on X-Chromosome in Wild and Domestic Sheep Species.

Analyses of genomic diversity along the X chromosome and of its correlation with autosomal diversity can facilitate understanding of evolutionary forces in shaping sex-linked genomic architecture. Strong selective sweeps and accelerated genetic drift on the X-chromosome have been inferred in primates and other model species, but no such insight has yet been gained in domestic animals compared with their wild relatives. Here, we analyzed X-chromosome variability in a large ovine data set, including a BeadChip array for 943 ewes from the world's sheep populations and 110 whole genomes of wild and domestic sheep. Analyzing whole-genome sequences, we observed a substantially reduced X-to-autosome diversity ratio (∼0.6) compared with the value expected under a neutral model (0.75). In particular, one large X-linked segment (43.05-79.25 Mb) was found to show extremely low diversity, most likely due to a high density of coding genes, featuring highly conserved regions. In general, we observed higher nucleotide diversity on the autosomes, but a flat diversity gradient in X-linked segments, as a function of increasing distance from the nearest genes, leading to a decreased X: autosome (X/A) diversity ratio and contrasting to the positive correlation detected in primates and other model animals. Our evidence suggests that accelerated genetic drift but reduced directional selection on X chromosome, as well as sex-biased demographic events, explain low X-chromosome diversity in sheep species. The distinct patterns of X-linked and X/A diversity we observed between Middle Eastern and non-Middle Eastern sheep populations can be explained by multiple migrations, selection, and admixture during the domestic sheep's recent postdomestication demographic expansion, coupled with natural selection for adaptation to new environments. In addition, we identify important novel genes involved in abnormal behavioral phenotypes, metabolism, and immunity, under selection on the sheep X-chromosome.

Introduction

The X chromosome and autosomes differ in several aspects of their molecular evolution and population genetics including their effective population sizes (Ne), mutation rates, susceptibility to genetic drift, patterns of selection and demography (Schaffner 2004; Vicoso and Charlesworth 2006; Hedrick 2007; Keinan and Reich 2010), all of which differentially affect the genetic diversity of the X-chromosome compared with the autosomes (Begun and Whitley 2000; Casto et al. 2010; Heyer and Segurel 2010) and the X-to-autosome (X/A) diversity ratio (Ellegren 2009). As a result, deviations from neutral expectation of 0.75 (X)/1.00 (A) (Charlesworth et al. 1987) due to variable diversity ratios across the chromosomes and in specific genomic regions are expected due to the combined effect of different evolutionary forces (Gottipati et al. 2011; Arbiza et al. 2014). Therefore, comparing patterns of X-chromosome and autosomal diversity can provide insight into the relative contributions of these forces in shaping genome evolution.

Several recent studies have compared genomic diversity on the X-chromosome with the autosomes, including comparative analysis of nucleotide diversity as a function of distance from the near genes in humans (Keinan et al. 2009; Hammer et al. 2010; Gottipati et al. 2011; Arbiza et al. 2014) and great apes (Nam et al. 2015, 2017; Narang and Wilson Sayres 2016). More generally, a variety of evolutionary forces have also been invoked to explain patterns of genetic variability on the X chromosome versus autosomes in model species through empirical data analyses (Andolfatto 2001; Lu and Wu 2005; Singh and Petrov 2007; Cox et al. 2010; Corl and Ellegren 2012; VanBuren et al. 2016) and statistical modeling (Betancourt et al. 2004; Ellegren 2009; Heyer and Segurel 2010; Keinan and Reich 2010; Charlesworth 2012; Evans et al. 2014; Veeramah et al. 2014; Lasne et al. 2017).

Unlike most species studied to date, domestic animals have experienced strong and long-term directional selection during the processes of domestication and breeding and they therefore represent an ideal model to explore the evolution of the X chromosome in comparison both with their wild relatives and among divergently selected domestic populations. Sheep (Ovis aries) were domesticated in the Fertile Crescent ∼ 8000–11000 years ago (Ryder 1983), with diverse native breeds having been developed subsequently across the world with a concomitant demographic expansion, and through adaptation to a diverse range of local environments and varied production systems (Lv et al. 2014). Importantly, their wild relatives (O. orientalis, O. ammon, O. musimon, O. vignei, O. dalli, O. canadensis, and O. nivicola) remain extant (Rezaei et al. 2010), making them an ideal model to disentangle the effects of multiple underlying forces on X chromosome evolution.

Here, we analyzed X-chromosomal genetic variability across species in the genus Ovis, comprising a large data set for domestic sheep and their wild relatives. The data consist of genome-wide BeadChip SNP genotypes of 943 females from 68 domestic populations (fig. 1 and supplementary table S1, Supplementary Material online) and 110 whole genomes of wild and domestic sheep from around the world (fig. 1 and supplementary table S2, Supplementary Material online), representing the largest sample size and population sampling of X-chromosomes in any nonmodel mammal. We compared intra- and interpopulation genetic diversity on the X-chromosome and autosomes. We estimated the ratios of X-to-autosome (X/A) nucleotide diversities (e.g., absolute and relative X/A diversity; Arbiza et al. 2014), and evaluated their patterns as a function of distance from genes. We also compared these patterns in domestic populations of different geographic origin and between Ovis species. Further, we investigated potential explanations for the observed patterns by examining the effects of selection, sex-biased migration, admixture, and genetic drift. Our goals here were 3-fold. First, we aimed to detect general patterns of X-chromosome versus autosomes diversity among species as well as among domestic populations during domestication, demographic expansions, and breeding processes. Second, we sought to disentangle these different forces in shaping the patterns of X-chromosome evolution. Third, we examined selective footprints on the X-chromosome and related these signals to early domestication as well as recent artificial selection in breeding practices and natural selection in different environments.

. 1.

—Geographic origins of the samples and distribution of the six wild Ovis species (O. orientalis, O. ammon, O. musimon, O. vignei, O. dalli, and O. canadensis). Geographic distributions of the six wild Ovis species were adapted from figure 1 in Rezaei et al. (2010).

Materials and Methods

The Ovine SNP50K BeadChip and Whole-Genome Sequence Data

We used publically available Illumina Ovine 50 K (54,241 SNPs) BeadChip SNP and whole-genome sequence data of domestic sheep (O. aries) and their wild relatives (e.g., O. orientalis, O. ammon, O. musimon, O. vignei, O. dalli, and O. canadensis). The ovine SNP BeadChip data consisted of genotypes of 943 females representing 68 worldwide domestic populations/breeds in previous investigations (Kijas et al. 2012; Zhao et al. 2017), in addition to 6 males and 2 females of O. orientalis obtained from the NextGen Project (http://projects.ensembl.org/nextgen/) (fig. 1 and supplementary table S1, Supplementary Material online). Details on quality control for the Ovine SNP50K BeadChip genotypes were provided in Supplementary Material.

In addition, we obtained a total of 110 whole genome sequences of domestic and wild sheep with a sequencing depth of 5–14× (fig. 1 and supplementary table S2, Supplementary Material online). Among the genomes, sequences of 63 domestic and 26 wild sheep (17 O. orientalis, 4 O. vignei, 2 O. dalli, and 3 O. canadensis) were from the NextGen Project, and 19 additional domestic and wild (1 O. ammon and 1 O. musimon) sheep sequences were from our previous work (Yang et al. 2016) (fig. 1 and supplementary table S2, Supplementary Material online). Quality control, filtering, and alignment for sequencing reads were detailed in Supplementary Material.

SNP Calling and Validation

We conducted population-scale SNP calling using a Bayesian approach as implemented in the SAMtools “mpileup” command (Li et al. 2009) with the parameters “-pup -q 1 -C 50 -S -D -m 2 -F 0.002” “-popf –Q 20 –d 1 –D 10000000.” Only the bases with quality scores >20 (Q ≥ 20) were considered, and then we kept the positions in which genotypes are reported from at least one individual for the nucleotide diversity calculation. After SNP calling, we filtered raw variants at the individual level using the following criteria, and SNPs meeting either of the criteria were removed: 1) the SNPs are not in biallelic sites; 2) low-quality SNPs with the coverage depth < 5 or > 100 and RMS (root mean square) mapping quality < 20. To minimize the number of false positive calls, we also applied more strict filtering criteria in domestic and wild sheep, separately: 1) MAF < 0.05; 2) the proportion of missing genotypes > 50%; and 3) P value of the chi-squared test for Hardy–Weinberg equilibrium (HWE) test > 0.00001 (P value > 0.00001). SNPs meet any of the above criteria were filtered. These high-quality SNPs were referred as whole-genome sequence (WGS) SNPs.

Further, we examine the accuracy of our high-quality SNPs using the method as described in Yang et al. (2016). In summary, we compared the SNPs with those from Build 143 of the sheep dbSNP database in the National Center for Biotechnology Information (NCBI; http://www.ncbi.nlm.nih.gov/SNP). We annotated all the high-quality SNPs using the ANNOVAR software (Wang et al. 2010) and categorized them into three groups according to their genomic locations such as intronic, exonic, and intergenic.

Estimates of Nucleotide Diversity and Its Pattern along X Chromosome

We calculated within-population genomic diversity metrics and estimated the linkage disequilibrium (LD) pattern between pairwise SNPs (Supplementary Material online). In addition, genetic differentiation and population structure were examined based on the Ovine BeadChip and WGS data sets (Supplementary Material online).

The nucleotide diversity (π) on per-site basis was calculated using the formula by Nei and Li (1979) in VCFtools (Danecek et al. 2011). For SNPs in intronic, exonic, and intergenic locations, we calculated π separately. To correct for the effect of potential difference in mutation rates between X-chromosome and autosomes, we calculated mean interspecies divergence (D) (Jukes and Cantor 1969) for X-chromosome and autosomes between sheep and cattle (Bos taurus) using the software MEGAv.7 (Kumar et al. 2016) under the Jukes–Cantor model (Jukes and Cantor 1969). We then normalized the nucleotide diversity (π) using the D estimates following the method in Gottipati et al. (2011). In summary, we aligned sequences of domestic sheep and cattle from UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/oviAri1/vsBosTau4/), and computed the mean sheep–cattle divergence as the fraction of fixed differences using a Jukes–Cantor correction for recurrent mutations (Jukes and Cantor 1969).

In order to examine the impact of selection on linked sites via affecting estimates of nucleotide diversity and absolute X/A diversity (X-to-autosome nucleotide diversity ratio), we evaluated the above two diversity indices as a function of physical distances from the nearest genes in 1 Mb windows, and sorted these windows into 20 bins at an increasing distance. Physical distance from each SNP site to the nearest gene was estimated following the method of Nam et al. (2017).The 95% confidence intervals for the estimates of π and absolute X/A diversity were obtained using a nonparametric bootstrapping procedure of resampling 1000 random data set from each bin.

Further, we estimated relative X/A diversity, the ratio of the estimates of absolute X/A diversity (Gottipati et al. 2011), to unravel differential forces of demographic history and selection among populations from different continental origins (e.g., Middle East, Central and East Asia, South Asia, Africa, America, and Europe). Different from the absolute X/A diversity, relative X/A diversity overcomes the effects stemming from differences in mutation rates between X-chromosome and the autosomes and other potential issues (e.g., difference in SNP ascertainment biases between X-chromosome and autosomes) (Arbiza et al. 2014). We assumed the initial sheep domestication in the Middle East (Zeder 2008; Lv et al. 2015), and calculated the relative X/A diversity between non-Middle Eastern populations to the Middle Eastern populations following the method of Gottipati et al. (2011). We also examine the relationship between relative X/A diversity and physical distances from the closest gene in the comparisons as detailed above.

Analysis of Sex-Biased Migration

To explore the potential effect of sex-biased demographic events on the observed pattern (Keinan and Reich 2010), we calculated Q statistics, a test for sex-biased migration (Goldberg et al. 2017). Under a simple demographic model, that is, two isolated populations originated from one common ancestor, and have been panmictic and constant in population size since their split, autosome-to-X genetic drift ratio (Q) is equal to 0.75 (Keinan et al. 2009). We measured sex-biased genetic differentiation on autosomes and X chromosome, that is, autosomal () and X-chromosomal () FST values, and then estimated the Q statistic in the formula, following the methods of Keinan and Reich (2010) and Waldman et al. (2016). Based on the assumption of initial sheep domestication in the Middle East (Zeder 2008; Lv et al. 2015), we calculated the Q statistic between different continental groupings using the X-linked and autosomal SNPs. From the theoretical prediction, the Q statistic is equal to 0.75, while deviation from the expected value suggests possible occurrence of sex-biased demographic events (Ramachandran et al. 2008; Emery et al. 2010).

Simulations of Genetic Drift

We simulate the genetic drift on X-chromosome and autosomes in the WGS data using an updated version of the POWSIM software (Ryman and Palm 2006) as detailed in Lamichhaney et al. (2012) and Lamichhaney et al. (2017). Simulations of the expected distribution of FST among domestic populations were conducted under a selective neutrality model. We then compare the genetic drift on X and autosomes using the distributions of simulated and observed FST values. To obtain the selectively neutral SNP loci, the following criteria were applied, and SNPs met any of the criteria were removed: 1) SNPs are within genic regions; 2) SNPs had less than five reads from each sample; 3) SNPs had >20% missing rates in a population; 4) SNPs had negative observed FST values. After the filtering, we obtained 27,167 autosomal and 1,795 X-linked neutral SNPs. For X-linked SNPs, we simulated an infinitely large base population with a major frequency of 0.863, and the base population was split into six subpopulations. The simulation was run for 1,795 repetitions with the effective size of Ne = 8,350, generation time of t = 840 (3 years/generation), and expected FST = 0.049. For autosomal SNPs, simulation with 27,167 repetitions was conducted with a major frequency of 0.899 for the base population and six subpopulations. The simulation was run with the effective sizes of Ne = 10,000, generation time of t = 840 (3 years/generation), and expected FST = 0.0412.

Tests for Selective Sweeps on the X-Chromosome and Autosomes

To characterize the selective forces shaping patterns of genetic variability on the X chromosome and autosomes, we applied three statistical tests to detect signatures of selective sweeps. We calculated FST values (Weir and Cockerham 1984), integrated haplotype scores (iHS) (Voight et al. 2006), and cross-population extended haplotype homozogysity scores (XP-EHH) (Sabeti et al. 2007) between domestic sheep and their putative wild ancestors Asian mouflon as well as between different continental groups of sheep populations.

In the FST-based test, we calculated the average FST value for 1 Mb nonoverlapping windows based on the Beadchip using GENEPOP v4.3 (Rousset 2008), and for 200-kb nonoverlapping windows based on the WGS SNPs using VCFtools (Danecek et al. 2011). In the iHS and XP-EHH analyses, we first inferred ancestral alleles by comparing the BeadChip and the WGS SNPs with sheep.oarv3_1.ancestral_alleles (BeadChip SNPs: http://www.sheephapmap.org/download/results/AncestralAlleles; WGS SNPs: ftp://ftp.ebi.ac.uk/pub/databases/nextgen/ovis/variants/supplemtary_info/; Jiang et al. 2014). Only SNPs which possess ancestral alleles (BeadChip: 1,132 X-chromosomal and 44,031 autosomal SNPs; WGS: 52,325 X-chromosomal and 6,073,533 autosomal SNPs) were kept for the following analyses. We then inferred haplotypes for each sample using the SHAPEIT v2 program (Delaneau et al. 2013). Finally, we computed the fraction of extreme iHS (|iHS| > 2) or XP-EHH (XP-EHH > 2) and its empirical P values for each pair of comparison using 1 Mb and 200-kb nonoverlapping windows for the BeadChip and the WGS SNPs, respectively. Both the analyses were implemented in Selscan (Szpiech and Hernandez 2014) assuming 1 Mb = 1 cM across the sheep genome. In addition, we quantified and compared the selective intensities between X-chromosome and autosomes by estimating the ratio of selective SNPs (i.e., |iHS| > 2 or XP-EHH > 2) to all the X-chromosomal or autosomal SNPs.

Further, we conducted functional annotation for genes in selective regions using Oar_v4.0 (https://www.ncbi.nlm.nih.gov/genome/?term=sheep). For those candidate genes, we performed the phenotype enrichment analysis in the MGI database (Mouse Genome Informatics; http://www.informatics.jax.org/) using the Enrichr program (Chen et al. 2013; Kuleshov et al. 2016). In addition, we implemented Gene Ontology (GO) enrichment and pathway analyses using the DAVID Bioinformatics Resources 6.8 (available from http://david.abcc.ncifcrf.gov) (Huang et al. 2009a, 2009b). We set all known sheep genes as the background gene set and assessed the statistical significance (P value) by the Binomial distribution test.

Results and Discussion

Genome Assembly and SNPs

We analyzed a total of 3,760 Gb raw pair-ended sequence reads from 110 whole genomes (82 domestic and 28 wild individuals) from seven Ovis species at an average coverage of∼ 8.4× for both X-chromosome and autosomes (supplementary table S3, Supplementary Material online). After strict filtering, we obtained a total of 11,743,108 (autosome: 11,420,115; X: 322,993) high-quality SNPs in sheep and 28,926,374 (autosome: 28,358,469; X: 567,905) high-quality SNPs in their wild relatives (supplementary table S4, Supplementary Material online), with 7,605,294 SNPs shared between domestic and wild species. As expected, most of the high-quality SNPs were located in intergenic regions, followed by intronic and exonic regions (supplementary table S5, Supplementary Material online). Overall, 94.43% of the SNPs identified in sheep and 79.76% in wild ovines were validated using the sheep dbSNP database, indicating a high reliability of the called SNP variants (SNP validation rate in sheep: 89.46–96.21%; Yang et al. 2016; compared with, for examples, cross-species SNP validation for domestic horse Equus caballus in the Przewalski’s horse Equus ferus przewalskii: 50.8%; McCue et al. 2012).

For the WGS data, 280,304 and 442,000 X-linked SNPs in domestic and wild sheep were analyzed, excluding those in the pseudoautosomal region (PAR: 42,689 and 125,905 SNPs in domestic and wild species, respectively; see the Results below). After the quality control and removing SNPs in the PAR (7 SNPs; see below), the BeadChip data set consisted of 45,822 autosomal and 1,187 X-linked SNPs, from a wide collection of native sheep breeds from around the world with varied phenotypic traits such as the types of main product (e.g., wool, meat, and milk), growth rate, body size, tail type, coat color, and disease resistance.

X-Chromosome Pseudoautosomal Region

To locate the PAR, the observed heterozygosity (HO) for each SNP in the BeadChip was calculated within females and an additional set of 924 males from the world’s sheep populations, separately. In the BeadChip, the first 7 SNPs, located between the positions 1.12–5.55 Mb (SNPs: rs400492871–rs401737908), showed an average HO = 0.25 in males (supplementary fig. S1, Supplementary Material online), whereas almost all the remainder 7.14–134.93 Mb showed estimates of Ho = 0 (supplementary fig. S1B, Supplementary Material online). In contrast, an average estimate of HO = 0.31 was observed for those in females (supplementary fig. S1A, Supplementary Material online). Therefore, these observation suggested that the PAR is likely located between 0 and 7.10 Mb on the X-chromosome in sheep. This observation is roughly congruent with the recombining region (i.e., PAR) in the sheep reference genome assembly (Jiang et al. 2014), in which the genomic region 0–7.05 Mb was suggested for the PAR.

In the WGS data, we observed significantly lower (Wilcoxon P < 0.01) mean π estimates in both the autosomes and X-linked region in sheep (πA = 0.0031, πX = 0.0016) than counterpart values in Asian mouflon (πA = 0.0038, πX = 0.0019) (fig. 2 and supplementary table S6, Supplementary Material online), which is consistent with a loss of genetic diversity due to a population bottleneck during domestication (VanBuren et al. 2016). In both sheep and Asian mouflon, we also observed a pronounced shift in π on both sides of the PAR boundary: average estimate of π in the PAR was ∼3.6-fold higher than in X-linked regions (fig. 2 and supplementary table S6, Supplementary Material online). This observation could be explained by a higher recombination rate in the PAR than the X-specific portion during male meiosis (Galtier 2004; White et al. 2012). Similar patterns of variation in π were also observed in all the other wild sheep (supplementary table S6, Supplementary Material online). Taken together previous inferences about the location of the PAR in sheep and other ruminants (Van Laere et al. 2008; Das et al. 2009; Jiang et al. 2014;), patterns of both genetic diversity indices (i.e., HO and π) supported for the first time a clear and consistent PAR boundary (∼7.10 Mb) across the genus Ovis.

. 2.

—Nucleotide diversity and divergence along the X chromosome. (A) Nucleotide diversity (π) values of domestic sheep (Ovis aries) are showed in 50 kb windows and 10 kb steps along chromosome X. (B) Nucleotide diversity (π) values of Asian mouflon (O. orientalis) are showed in 50 kb windows and 10 kb steps along chromosome X. (C) FST dynamic variation among domestic sheep (O. aries) populations are plotted in 100 kb windows and 50 kb steps along chromosome X. (D) FST dynamic variation between sheep (O. aries) and Asian mouflon (O. orientalis) are plotted in 100 kb windows and 50 kb steps along the X chromosome. The pseudoautosomal region (PAR) is represented in red dots. A sequential and conserved domain with lower nucleotide diversity and higher divergence is marked by the red dotted lines.

We observed significantly lower (Wilcoxon P < 0.01) genetic divergence (i.e., mean FST values across loci) in the PAR region than the X-linked regions between sheep populations as well as between sheep and Asian mouflon (fig. 2). A similar pattern of lower mean FST in the PAR than the X-linked regions has also been observed among worldwide pig populations and between domestic pigs and wild boar (Burgos-Paz et al. 2013). Similar patterns of genetic differentiation for PAR versus X-linked between domestic animals and their wild species could be partially attributed to some early genetic introgression (e.g., via interbreeding) from the wild relatives into domestic species following domestication, which is known to have occurred commonly (Larson and Burger 2013; Lv et al. 2015; Barbato et al. 2017; Zhao et al. 2017). However, an elevated divergence in PAR relative the X-linked region has been found between human and other primate species such as chimpanzee and macaque, potentially owing to complex speciation events (Patterson et al. 2006; Cotter et al. 2016).

Within-Population Genetic Diversity and between-Population Genetic Differentiation for X-Chromosome versus Autosome

In general, within-population genetic diversity indices (HO, HE, Pn, and Ne) estimated from X-linked SNPs were lower than those from autosomal SNPs in the BeadChip data (supplementary table S7, Supplementary Material online). Meanwhile, we observed higher level of LD in the X-linked region than chromosome 6, similar in the size and the number of SNPs to those in the X-linked region, within populations (supplementary fig. S2, Supplementary Material online). Additionally, the X-linked region showed more ROHs, a larger average size of ROHs and a higher coverage of ROHs than autosomes in sheep and the wild species (supplementary figs. S3 and S4, Supplementary Material online). The BeadChip data showed a reduction in genetic diversity for the X-linked versus autosomal SNPs at a population and continental level, but this was much smaller than the observed reduction in π of the X-linked region relative to autosomes in the WGS data (supplementary tables S6 and S7, Supplementary Material online) and the 1/4 reduction under assumption of equilibrium (Charlesworth et al. 1987; Ellegren 2009). Similar differences have been observed in X chromosome analyses of worldwide pig populations using SNP and sequence data (Amaral et al. 2011; Esteve-Codina et al. 2011; Burgos-Paz et al. 2013). The most likely explanation for our results is SNP ascertainment bias in the BeadChip as indicated by the longer branch lengths in the phylogenetic trees built based on the BeadChip data than those based on the WGS data (Supplementary Material and supplementary fig. S5, Supplementary Material online), in which SNPs with the highest minor allele frequencies (MAFs) were selected for both the autosomes and the X-chromosome.

In the analyses of population genetic structure using BeadChip SNPs, both the autosomal and X-linked SNPs revealed apparent geographic differentiation (i.e., Middle East, Central and East Asia, South Asia, Africa, America, and Europe) in the STRUCTURE (with the optimal K = 6; see Supplementary Material and supplementary figs. S6 and S7, Supplementary Material online) and AMOVA analyses (supplementary table S8, Supplementary Material online). Nevertheless, the X-linked SNPs showed significantly (Wilcoxon P < 0.01) higher percentage of molecular variance than autosomal SNPs between continental groups (X: 8.42%; autosome: 4.42%) as well as between populations within continental groups (X: 9.99%; autosome: 8.35%; supplementary table S8, Supplementary Material online). Similarly, higher levels of genetic divergence for X-linked SNPs relative to autosomal loci were observed in the WGS data (supplementary table S6, Supplementary Material online), while the SNPs in the PAR failed to separate the sheep population into continental groups in the PCA analyses (supplementary fig. S8, Supplementary Material online).

Patterns of Nucleotide Diversity along the X-Chromosome

In Ovis species, we observed the lowest mean π in exons followed by intronic and intergenic locations in the X-linked regions and autosomes (fig. 3). We found one conserved domain (43.05–79.25 Mb) featuring lower π in the X-linked region across all species (fig. 4). Nam et al. (2015) proposed three potential forces for the extremely low X-chromosomal nucleotide diversity in great apes, such as close inbreeding, selective sweeps, and low mutation rates in large genomic regions.

. 3.

—Nucleotide diversity (π) and the X-linked region/autosome diversity ratio (X/A) in exons, introns, and intergenic regions in the seven Ovis species. The phylogenetic relationships among the seven Ovis species investigated and the divergence times (in Ma) at the top were adapted from figure 1 in Rezaei et al. (2010) The top, middle (within the box), and bottom boundary lines of the boxes represent 25%, 50% (median value), and 75% of π.

. 4.

—Nucleotide diversity along the X-linked regions for the seven species of Ovis. Black bars at the top of each panel indicate the regions where the nucleotide diversity is 20% less than the mean diversity across chromosome X in the species. A sequential and conserved domain with low nucleotide diversity level is marked by the red dotted lines. Distribution of protein coding genes on chromosome X was plotted at the top.

Here, inbreeding is unlikely to explain the pattern because no conserved region of low diversity was observed on the autosomes (supplementary fig. S9, Supplementary Material online). In order to clarify the selective effect in shaping these genomic variations, we identified selected genomic regions using the criteria of |iHS| (integrated haplotype score; Voight et al. 2006) > 2 and XP-EHH (cross-population extended haplotype homozogysity; Sabeti et al. 2007) > 2 in tests between domestic sheep and Asian mouflon, and further compared selection intensities between X-linked regions and autosomes. We found significantly (P < 0.01) larger numbers of selective SNPs per Mb in autosomes (|iHS|: 106.99 SNPs/Mb; XP-EHH: 87.72 SNPs/Mb) than in the X-linked region (|iHS|: 16.74 SNPs/Mb; XP-EHH: 12.65 SNPs/Mb) based on the WGS data (table 1). In addition, similar selective sweeps analyses in the BeadChip data showed lower intensities in the X-linked region than in autosomes (table 1). Different from previous evidence for stronger selective sweeps on X-chromosome than autosomes in the primates (Hammer et al. 2010; Nam et al. 2015), less intensive selection on the X-linked region observed here could be due to the fact that the majority of QTLs (https://www.animalgenome.org/cgi-bin/QTLdb/OA/summary? summ= chro&qtl=1, 658&pub=132&trait=225) (supplementary fig. S10A, Supplementary Material online) and genes (supplementary fig. S10B, Supplementary Material online) (Xu and Li 2017) for production and functional traits under strong and long-term directional selection are located on the autosomes instead of the X-chromosome in domestic animals (supplementary fig. S10, Supplementary Material online). However, we did not detect a large conserved region of highly reduced nucleotide diversity on the autosomes. Thus, these results suggested that selective sweeps, at least alone, cannot explain the observation. In addition, this domain showed higher levels of genetic divergence between sheep and Asian mouflon (fig. 2). Together with the evidence of low π values in the conserved domain in the ancestral species of sheep, our evidence implied that this pattern could not be ascribed to reduced recombination rates (Nam et al. 2015). Interestingly, we observed an extremely high density of coding genes in this domain (fig. 4), and, thus, a large number of conserved coding regions in these genes could partly account for the large diversity reduction observed in the X-linked region.

Table 1

Comparison of Selection Intensity (number of SNPs/Mb) between the X-Linked Region and Autosomes Based on the Selection Tests between Domestic Sheep and Asian Mouflon using the iHS (Voight et al. 2006) and XP-EHH (Sabeti et al. 2007) Approaches

	Illumina Ovine 50K BeadChip		Whole-Genome Sequences
	iHSa	XP-EHHb	iHSa	XP-EHHb
X-linked region	0.34	0.22	16.74	12.65
Autosome	0.83	0.38	106.99	87.72

Number of SNPs selected based on the threshold of |iHS| > 2.

Number of SNPs selected based on the threshold of XP-EHH > 2.

In five of the species (sheep, Asian mouflon, European mouflon, argali, and urial), mean estimates of π in the X-linked regions was ∼60% of those in autosomes for the exonic, intronic, and intergenic locations (fig. 3), which is substantially lower than the expected X/A ratio of 3/4 under neutral evolutionary models (Charlesworth et al. 1987). Nevertheless, the ratios were in the range (X/A ratio = 0.32–1.47) observed in previous studies in primates and other model species (Ellegren 2009; Hammer et al. 2010; Arbiza et al. 2014; Nam et al. 2015). Higher X/A diversity ratios (> 1) were found in O. canadensis and O. dallis, which are more distantly related to domestic sheep (fig. 3). Similar cases have rarely been observed in other species or populations except for African populations of D. melanogaster reported by Andolfatto (2001). Andolfatto’s argument is that autosomal inversions work as recombination suppressors and will therefore lead to reduced autosomal variability when backgrounds selection or the frequency of selective sweeps is large (Andolfatto 2001). Mammals usually do not have the genomes littered with inversion like Drosophila, so it is an unlikely explanation. Here, the most likely explanation for the observation could be that the species are more divergent from domestic sheep and, thus, we can only map reliably in the least variable regions.

Patterns of Diversity Related to Distances from Genic Regions

In all the six continental groups of populations, the SNP estimates of π increased with physical distance from genes (in 1 Mb windows and the windows were sorted into 20 bins at an increasing distances) in autosomes and X-linked regions (fig. 5), consistent with increased levels of selection on genes and linked regions (McVicker et al. 2009; Arbiza et al. 2014). However, a much lower and flatter trend in π values was observed at SNPs within 250 kb from genes in X-linked regions (π ≤ 0.002) compared with autosomes (π ≤ 0.004) (fig. 5). In addition, we observed a mean normalized X/A ratio (0.58) that was significantly lower than the expected values of 0.75 (3/4). The X/A ratio showed a slight decrease with an increase in distance in all the geographic groupings (fig. 5). This differs from a positive correlation between X/A diversity and distance to genes in human (McVicker et al. 2009; Gottipati et al. 2011; Arbiza et al. 2014) and the great apes (Nam et al. 2015; Narang and Wilson Sayres 2016), which has been mainly explained by stronger directional selection at genes and their linked sites on the X chromosome (Arbiza et al. 2014; Nam et al. 2015). Therefore, the lower impact on diversity-reduction in the X-linked region than the autosomes could be due to the more intensive recent directional selection on autosomes than in the X-linked region inferred above. An alternative explanation is the evolutionary inertness of the X-chromosome, which by the hemizygosity in male and the X-inactivation in females are subject to an intensive purifying selection and thus do not follow the normal evolutionary changes as observed for autosomes.

. 5.

—Nucleotide diversity (π) and the X-chromosome/autosome diversity ratio (πX/πA) with the physical distance to the nearest gene in the Ovis species. (A) π and πX/πA in the continental groupings of domestic populations. (B) π and πX/πA in the seven Ovis species. 95% confidence intervals from 1,000 bootstrapping iterations are showed in each bin.

As in domestic sheep, we also observed similar overall trends of π and X/A diversity with distance from genes in wild sheep (fig. 5). In particular, we found that π was nearly constant with the distance from genes on the autosomes (O. canadensis: π = 0.00127–0.00145; O. dallis: π = 0.00115–0.00123), while π decreased with distance from genes in the X-linked regions (fig. 5) in O. canadensis and O. dallis. In addition, the X/A diversity ratio (O. canadensis: 1.29–1.75; O. dallis: 1.47–1.99) were >1 and showed a decrease with distances from genes in O. canadensis and O. dallis (fig. 5). Similar patterns have rarely been revealed in primates and other model species (but see Evans et al. 2014). Forces including sex differences in reproductive success (Charlesworth 2012), varying population size and natural selection on GC content (Evans et al. 2014) have been invoked to explain such observations in different scenarios. However, we cannot differentiate the relative contribution of difference forces to the pattern observed here due to the limitation of data (e.g., the field data and de novo genomic data) in O. canadensis and O. dallis.

Relative X/A Diversity in Continental Population Groupings

In all five pairwise comparisons based on the WGS data, relative autosomal diversity was nearly constant within 250 kb distance from genes (fig. 6). Contrasting with the general assumption of higher diversity at the domestication center, we found significantly higher (∼3%) autosomal diversity in African, American, and European than in the Middle Eastern populations (fig. 6). A possible explanation could be multiple introgression from the Middle Eastern domestication center and admixture in these regions (Chessa et al. 2009). In addition, we noted that the Middle Eastern populations are fat-tailed and probably do not represent wholly the original thin-tailed domesticates that spread to Europe, Africa, and America. For the five pairwise groups, the relative X-linked diversity was approximately equal between the Middle Eastern and the Central and East Asian populations within 250 kb from the genes (fig. 6), with no significant correlation with distance from genes (the Pearson’s r = 0.45, one-tailed test P = 0.132). In contrast, we observed slight lower X-linked diversity in the Middle Eastern populations than populations of the other four continental groupings (fig. 6), and the relative X-linked diversity showed a significant positive correlation (the Pearson’s r = 0.711–0.925, one-tailed P < 0.05) with distance from genes for South Asian and African populations within 250 kb.

. 6.

—Relative nucleotide diversity (π) and their ratio (X/A) with the physical distance to the nearest gene for the six continental groupings of sheep populations. (A) The ratio of normalized diversity in the groupings of populations from the Middle East (ME) domestication center to those in the other five groupings of populations with the physical distance to the nearest gene for X-linked region. (B) The ratio of normalized diversity in the groupings of populations from the Middle East (ME) domestication to those in the other five groupings of populations with the physical distance to the nearest gene for autosome. (C) Relative X/A diversity in each pair groupings of populations with the physical distance to the nearest gene. 95% confidence intervals from 1,000 bootstrapping iterations are showed in each bin. CEA, Central, and East Asian populations; ME, Middle Eastern populations; SA, South Asian populations; AF, African populations; AM, American populations; EU, European populations. Non-ME, Central and East Asian/South Asian/African/American/European populations.

We also estimated relative X/A diversity among different continents groups. Four pairs of the relative X/A diversity (SA/ME, AF/ME, AM/ME, and EU/ME) showed significant positive (the Pearson’s r = 0.61–0.77, one-tailed P < 0.05) correlation with distances from genes (fig. 6), which could be due to the moderate differences in X-linked diversity among populations from different continents groups. Remarkably, within the domestic sheep the X-linked diversity is more variable across different continents than autosomal diversity (fig. 6 supplementary tables S6, S8, and S9, Supplementary Material online). It may be noted that purifying selection of mtDNA diversity has been postulated to act on the long term (Ho et al. 2011), and, thus, similar long-term purifying selection might have acted on the X-linked region as well. The remaining X/A comparison (CEA/ME) showed a lack of correlation with distances from genes (the Pearson’s r = 0.206, one-tailed P = 0.312). The relative X/A diversity between Central and East Asian and the Middle Eastern populations was nearly constant (X/A = 1.043) across all the scales of distances (≤ 250 kb) from genes (fig. 6), which indicated genomic similarity and similar effect of gene proximity between populations from the two regions. This is most likely due to frequent genetic exchange between populations from the two geographic groups in the past (Zhao et al. 2017). In addition, we observed nearly constant pattern of autosomal diversity, but apparently increasing pattern of X-linked diversity and relative X/A diversity with increased genetic distances for the comparison between Middle Eastern and non-Middle Eastern populations (fig. 6). In human, an increase in both autosomal and X-linked nucleotide diversity, but an invariant X/A diversity, with distance from genes was observed in the comparison between African and non-African populations (Arbiza et al. 2014). The difference in these patterns between sheep and human could be partly due to different demographic forces acted on them.

In addition, we noted that most of the genome sequences have the coverage of 10–20× and, consequently low or rare variants were missed. In addition, the number of false calls, despite the statistical treatments and quality control (QC) filtering, might be relatively higher. Thus, the allele callings might influence the final estimates. Nevertheless, we ran a site frequency spectrum (SFS) using the program ANGSD v.0.902 (Analysis of Next Generation Sequencing Data; Korneliussen et al. 2014) and did not find significant (Student’s t-test p = 0.503) difference in allele frequencies between the ones observed and those expected under the mutation-drift equilibrium (supplementary fig. S11, Supplementary Material online). In this sense, the influence of missed low or rare variants and false calls could be very slight, if there were.

Genetic Drift on X and Autosomes

To further investigate whether the patterns of nucleotide diversity and X/A diversity related to distance to genes were primarily driven by the genetic drift, we performed the simulations of genetic drift on X-chromosome and autosomes. For X chromosome, there was a good agreement between the observed and simulated FST data (fig. 7). Around 60% of the loci, both in the simulated and observed data sets, showed low FST values in the range of 0–0.049 (fig. 7). Under the selectively neutral model, average value of simulated FST (0.049) is roughly equal to the average value of observed FST (0.050). The slightly higher average value of observed FST was derived from a bit longer tail of FST values in the observed data set, which could be caused by selection (Lamichhaney et al. 2012). For autosomes, similar results were observed in the distributions of simulated and observed FST values (fig. 7). As a result, comparison of the simulated distributions of FST values indicated an accelerated genetic drift on X chromosome (average FST = 0.049) than on autosome (average FST = 0.040). Thus, our results suggest that genetic drift is one of the primary forces shaping the differential genomic patterns on X-chromosome and autosomes. Nevertheless, only the effect of genetic drift, but not other complex demographic scenarios (e.g., domestication bottlenecks), has been considered in the coalescent simulation, Thus, caution is needed when interpreting the results, because it still cannot fully disentangle the effects of demography and reduced or weak selection.

. 7.

—Distribution of observed and simulated FST values. (A) Density distribution of FST values in simulated and observed data sets on X chromosome. (B) Cumulative distribution of the observed and simulated FST values on X chromosome. (C) Density distribution of FST values in the simulated and observed data sets on autosome. (D) Cumulative distribution of observed and simulated FST values on autosome.

Scan for Selective Signatures

In the selection tests between sheep and Asian mouflon using the iHS (integrated haplotype score), FST, and XP-EHH (cross-population extended haplotype homozogysity) statistics (Weir and Cockerham 1984; Voight et al. 2006; Sabeti et al. 2007) we identified a total of 54 and 45 extreme windows in the X-linked region and autosomes in the WGS data, respectively (supplementary fig. S12, Supplementary Material online). In these selective windows, we annotated 60 genes on the X-chromosome and 93 genes in the autosomes (supplementary table S10, Supplementary Material online). In the phenotype enrichment analysis using the MGI database, we found a set of functional genes associated with behavioral response to stress, nervous system development, lethality, pigmentation, and infertility (fig. 8). In particular, a number of selective genes were detected in the X-linked region related to the abnormal behavioral phenotypes, such as abnormal response to novelty, hypoactivity, reduced long-term potentiation, abnormal hair growth, and abnormal heart and ventricular development (fig. 8 and supplementary tables S11 and S12, Supplementary Material online). These phenotypes are the most likely traits that have been under strong selection during domestication (Kerje et al. 2004; Carneiro et al. 2014; Jensen 2014; Medugorac et al. 2017).

. 8.

—Significantly associated phenotypes for the selective genes in the phenotype enrichment analysis through the MGI database (Mouse Genome Informatics; http://www.informatics.jax.org/) using the Enrichr program. (A) The X-linked candidate genes. (B) The autosomal candidate genes.

Overall, we detected four common regions across the three tests, and seven genes (TRIM24, WDR36, CAMK4, STARD4, HAO1, ADRA1D, and SMOX) were annotated (supplementary fig. S13, Supplementary Material online). Of these, three (WDR36, CAMK4, and ADRA1D) function in the behavior and nervous system. In particular, CAMK4, found in the neuropathic pain signaling pathway, affects nervous system development and function in yak and relates to reduced fear or increased stress activity (Medugorac et al. 2017). In addition, ADRA1D is enriched in the neuroactive ligand-receptor interaction pathway and potentially regulates abnormal behavior in the horse to adapt to new environments (Librado et al. 2015).

In order to detect genomic regions under recent selection, similar tests were also implemented between continental population groupings using the Beadchip data, which consisted of a worldwide collection of native sheep (supplementary table S1, Supplementary Material online). We detected a total of 72 selective windows in the comparisons between the five continental groupings and the Middle East populations (supplementary fig. S14, Supplementary Material online). Of these selective windows, 19 were overlapping in two of the three different tests: only two windows (73–74 Mb and 74–75 Mb) were detected in at least two different continental regions, while the rest 17 were in specific geographic regions (supplementary table S13, Supplementary Material online). In the two windows, we annotated four genes, and two of them (DACH2 and POF1B) are related to the ovarian development, and could be under similar recent artificial selection for high fertility. We also annotated a total of 53 genes (supplementary table S13, Supplementary Material online) in the other seventeen windows identified in specific geographic regions. We found that a majority (64%, 34/53; e.g., C1GALT1C1, IL1RAPL2, and LAMP2) of genes are associated with metabolism, immunity, and diseases. Nevertheless, we did not found any enriched GO terms and pathways for the candidate genes, which could result from natural selection for local adaptation and resistance to various pathogens in different environments (Lv et al. 2014; Yang et al. 2016).

Conclusions

In conclusion, we present a comprehensive analysis of X-chromosome and autosomal diversity across the Ovis genus. We extend previous efforts in quantifying the patterns of X-chromosome diversity and X/A ratio to domestic animals and their wild relatives. We observed greater reduction in X-linked diversity relative to autosomal diversity in these species. In contrast to the typically positive correlation observed in the primates and other model animals, we detected decreased X/A diversity as a function of distances from genes. Our evidence suggested that reduced directional selection but accelerated genetic drift in the X-linked region than autosomes as well as sex-biased migration could have led to the X/A diversity pattern. In addition, our results provide evidence that demographic events including multiple migrations and admixture during the expansion process are likely to have played a role in shaping patterns of the relative X/A diversity between the Middle Eastern and non- Middle Eastern populations. Finally, we identified a set of novel candidate targets of artificial and natural selection on the sheep X-chromosome. Thus, we find a pattern of X/A diversity and provided evidence implicating forces such as selection, genetic drift, migration, and admixture, which need to be invoked for explaining them in domestic species and their wild relatives.

Data Availability

The Illumina Ovine 50K SNP BeadChip used in analyses can be accessed in Kijas et al. (2012), Zhao et al. (2017), and the NextGen Project (http://projects.ensembl.org/nextgen/). The whole-genome sequence data are obtained from the NextGen Project and Yang et al. (2016).

Author Contributions

M.-H.L. designed and supervised the study. Z.-H.C., M.Z., and F.-H.L. analyzed the data with contribution from K.N. M.W.B. contributed to the data. M.-H.L., F.-H.L., M.Z., and Z.-H.C. wrote the paper with contribution from K.N. and M.W.B. All the authors reviewed and approved the final manuscript.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.