Genomic analysis of morphometric traits in bighorn sheep using the Ovine Infinium^® HD SNP BeadChip

Introduction

A goal of molecular ecology is to identify the genetic architecture underlying traits of ecological relevance (Ellegren & Sheldon, 2008; Slate et al., 2009). There is a particular interest in finding regions associated with variation in fitness-related traits, as such traits and the loci underlying them are expected to be subject to strong selection. Under strong directional selection, the genetic variability underlying fitness-related traits should rapidly go to fixation, and yet much phenotypic variation in such traits is observed in the wild (Kruuk, Slate & Wilson, 2008; Chenoweth & McGuigan, 2010). Elucidating the genetic basis of fitness-related traits might help clarify how phenotypic variation is maintained, for example by overdominance or epistatic interactions (e.g., Lappalainen et al., 2011; Johnston et al., 2013).

Sexual selection poses particularly interesting scenarios for fitness-related traits, leading to the so called lek paradox (Borgia, 1979). Under sexual selection, mate choice is either non-random, based on an ornamental trait that may confer benefits to the offspring, or members of one sex compete for access to mates using dimorphic secondary sex characteristics. Over time, under both scenarios evolution should deplete the genetic variation underlying the selected trait, and thereby diminishing offspring benefits. Paradoxically, however, in many systems choice for such traits continues. A classic example of a secondary sexual characteristic often subject to sexual selection is horn size in bovids. Horns are characterized by a keratin sheath around a bony projection from the skull that grows continuously throughout life (Davis, Brakora & Lee, 2011). Across a number of species, horn size in males determines social status and mating access to females (Bro-Jørgensen, 2007). Thus, there is likely selection for variants that confer the ability of individuals to grow large horns, and yet variation remains.

In domestic sheep (Ovis aries) some breeds have horns while others are polled (lacking horns entirely), and from an agronomic production standpoint there is interest in removing horns from certain breeds (Kijas et al., 2012). Soay sheep, a primitive breed now living feral on the islands of St. Kilda, Scotland, have an additional horn “morph”. In females, there are three morphs: normal, scurred (deformed horns composed only of keratin sheaths), and polled. In males there are only two morphs: normal horns and scurs (Johnston et al., 2010; Johnston et al., 2011).

The genetic basis of these differences in horn morphology and development in domestic sheep has been investigated. A single genomic region on chromosome 10 is associated with the presence and absence of horns in domestic breeds (Kijas et al., 2012) including Soay sheep (Johnston et al., 2010; Johnston et al., 2011), and is linked to quantitative differences in horn length of normal-horned male Soay sheep (Johnston et al., 2011). This region contains a single gene, relaxin/insulin-like family peptide receptor 2 (RXFP2). RXFP2 affects osteoporosis and testicular descent in mice and humans (Ferlin et al., 2008; Feng et al., 2009; Yuan et al., 2010); thus, its association with both bone development and secondary sexual characteristics make it an interesting candidate for influencing horn morphology. Furthermore, different genotypes at this locus in male Soay sheep are associated with trade-offs between reproductive success and longevity, which may maintain the different horn morphs through heterozygote overdominance (Johnston et al., 2013). Finally, though not the major QTL underlying horn phenotype, RXFP2 has been implicated in horn development in several association studies in cattle (Bos taurus) (Gautier & Naves, 2011; Allais-Bonnet et al., 2013; Wiedemar & Drögemüller, 2015), indicating that it may have similar function across species.

Bighorn sheep (Ovis canadensis) are iconic North American wild sheep named for their large horns. All individuals have normal horns, though males have much larger horns than females. Previous research has shown that horn size and body mass are important for intrasexual competition among males for reproductive access to females (Coltman et al., 2002; Martin et al., 2013). However, for female bighorn sheep horn length was found to be unrelated to social rank or other life history characteristics, which were more determined by body mass and age (Favre, Martin & Festa-Bianchet, 2008). In addition, horn size determines the trophy status of a male and, in many jurisdictions, whether it can be legally harvested. Such regulations directly influence longevity, as males with fast-growing horns are removed at a younger age (Festa-Bianchet et al., 2008; Festa-Bianchet et al., 2014; Bonenfant et al., 2009; Hengeveld & Festa-Bianchet, 2011; Pigeon et al., 2016).

Previous studies of bighorn sheep have shown that both horn size and body mass are heritable (Coltman, 2005; Coltman et al., 2005; Poissant et al., 2008; Poissant et al., 2012) and quantitative trait locus (QTL) mapping with microsatellite loci highlighted several suggestive regions for different aspects of horn morphology (e.g., volume and base circumference) as well as body mass (Poissant et al., 2012). These regions appear on several chromosomes, but notably included suggestive QTLs for horn volume and base circumference in males on chromosome 10 that spans the region predicted to contain RXFP2. This same region was highlighted by Kardos et al. (2015) who used pooled genome sequencing to search for signatures of selective sweeps in bighorn sheep. Analyses from pooled samples representing 58 bighorn sheep from three populations across Montana (USA) found 21 regions showing signs of selective sweeps. The strongest signal from these regions contained RXFP2, and a comparison to haplotypes in Soay sheep (Johnston et al., 2013) suggested that the selected alleles are associated with large horn size (Kardos et al., 2015). However, there has yet to be an individual-based study of the genomics of horn size, as Kardos et al. (2015) was based on pooled sequencing, and no assessments have been done of other fitness related traits.

In this study, we build on these results using individual phenotypic data from a long-term study of marked sheep followed throughout their lives at Ram Mountain, Alberta Canada. Though the costs of genome sequencing continue to decline, large-scale resequencing projects are still out of reach for many study systems. Therefore, we capitalize on the close evolutionary relationship between domestic and bighorn sheep (Bunch et al., 2006) to efficiently genotype many SNPs in many individual bighorn sheep using a genomic technology developed for domestic sheep. By investigating a suite of morphological characters, we aim to examine the genetic architecture underlying complex quantitative traits in wild sheep, and with respect to horn size, assess if the architecture is similar to its domestic relative.

Methods and Materials

Results

Average horn length, horn base circumference, and body mass all showed positive phenotypic correlations, with the magnitude much stronger in males than females (Table S1). All three morphological traits also exhibited significant additive genetic variation, with values on par with other studies of this population (Table 1). In total, 95 individuals were genotyped on the SNP chip and used to filter loci based on GenTrain and GenCall scores. One individual was subsequently removed from further analyses after significant (>5%) pedigree inconsistencies were found. Of the original 606,006 loci on the chip, 474,277 returned genotypes in bighorn sheep. Subsequent filtering removed 8,528 loci based on their levels of missing data, 469,822 based on our minor allele threshold, and 127 loci based on HWE equilibrium. The final dataset contained 3,777 loci, with at least 60 markers on each autosome (average ± SD = 145.3 ± 88.6; Table S2). Such reductions in the number of polymorphic loci are expected in cross-species application of SNP chips (Miller et al., 2012). Of the 94 originally genotyped individuals, 76 had morphological measures and were used in subsequent analyses.

Manhattan plots for each trait are shown in Fig. 1 with corresponding QQ-plots. In all cases genomic inflation (λ) was ≤1, indicating that there was no underlying population structure or other factors which could lead to false positive associations (Freedman et al., 2004; François et al., 2016). No loci were associated at the genome-wide significance level to any of the morphological traits examined. One locus, OAR9_91647990, showed suggestive association with body mass (Fig. 1, indicated with a green arrow).

As expected, there was a general decrease in LD with increasing inter-marker distance, and half-length was estimated to be 412,834 bp (Fig. 2). Based on this half-length estimate we extracted gene names from the Ovis aries gene set (Oar v3.1, genebuild Mar 2015) within a 413,000 bp window on either side of the candidate marker using BioMart (Kinsella et al., 2011) and Ensembl version 89 (Flicek et al., 2014). This returned two genes: U6 spliceosomal RNA, and ENSOARG00000026555, a long intergenic non-coding RNA. No gene ontology (GO) terms were available for either of these genes, and we do not see an immediate connection with body mass.

Discussion

We examined the genetic bases of three fitness-related characteristics in bighorn sheep. To do so, we utilized a new genomic technology originally designed for domestic sheep to rapidly genotype markers in a wild species, then combined these data with phenotypic measures from a long-term individual-based study. We found one locus with suggestive associations to body mass (Fig. 1). Previous QTL mapping with microsatellite loci for these same traits in the RM population highlighted several candidate regions (Poissant et al., 2012); however, our suggestive locus is not near any of the QTLs described in Poissant et al. (2012). In addition, we found no overlap in location between the locus found here and morphological traits in the domestic sheep QTL database (Hu, Fritz & Reecy, 2007; Hu et al., 2013). While it is possible that the sample sizes used in the Poissant et al. (2012) led to an overestimation of effect sizes due to the Beavis Effect (Slate, 2013) we note that the methods underlying QTL mapping and GWAS analyses are different (Slate et al., 2010). Specifically, QTL mapping relies on informative meioses within a pedigree of related individuals, while GWAS uses linkage disequilibrium between loci. In addition, the sample sizes differed, with fewer individuals included in the work presented here. These differences could influence the associations detected.

It is somewhat surprising that we did not see even a suggestive association between horn morphology and the region surrounding RXFP2 on chromosome 10 given the very strong links seen in both domestic sheep and cattle (Gautier & Naves, 2011; Johnston et al., 2011; Johnston et al., 2013; Kijas et al., 2012; Wiedemar et al., 2014) as well as the suggestive QTL for horn volume in bighorn sheep in this same region (Poissant et al., 2012). However, based on the estimate of half-length (412,834 bp) it appears as if we did not have sufficient marker coverage to adequately test for associations in the horns region. Within our set of loci the closest marker to RXFP2 was 698,861 bp away.

It is interesting that the extent of LD reported here (∼400,000 bp) is an order of magnitude less than found in a previous assessment of LD in bighorn sheep from RM (∼4,000,000 by Miller et al., 2011) using an order of magnitude fewer markers (308 vs. 3,777 loci). Analogous decreases in LD with the addition of markers have been seen in other species including cattle (McKay et al., 2007; Porto-Neto, Kijas & Reverter, 2014), domestic sheep (García-Gámez et al., 2012; Kijas et al., 2014), and flycatchers (Ficedula albicollis; Backström et al., 2006; Kawakami et al., 2014).

In light of our failure to detect genome-wide significant associations, we more formally quantified the expected power of a marker to detect a hypothetical causal QTL given the average minor allele frequency and genome wide critical p-value for the loci in this study. To do so we used an R script developed by Minikel (2012) which implements the QTL association feature of the Genetic Power Calculator (Sham et al., 2000; Purcell, Cherny & Sham, 2003). Specifically, this script estimates the expected power to detect an association given an estimate of the QTL effect size, the number of samples genotyped, and the average level of linkage disequilibrium among markers. For our analyses we varied effect sizes from 0–1.0, sample sizes between 50–500 individuals, and three levels of linkage disequilibrium (0.75, 0.50, and 0.25). This exploration showed that even at extreme effect sizes for the QTL and levels of LD well above what was seen at the half-length estimate (∼0.23; Fig. 2), the number of samples used in our GWAS analyses was likely not enough to have the power to detect all associations (Fig. 3). Note that these simulations assume that unrelated individuals were used in the GWAS, so the presence of related individuals in our test set will boost power slightly. In general, the simulations indicate that our marker coverage likely increased the chance of Type II errors (missing true associations). Similar results were found with simulations and whole genome sequences of collared flycatchers (Kardos et al., 2016). However, we do not believe this diminishes the association observed, as it has no effect on Type I errors (detecting false associations).

The power of our association analyses was likely also weakened by the cross-species application of a SNP chip originally derived for domestic sheep. While the two species are closely related (Bunch et al., 2006), and have a highly syntenic karyotype (Poissant et al., 2010), loci were selected for inclusion on the chip based on variability in domestic sheep, leading to ascertainment bias when applied to bighorn sheep (Lachance & Tishkoff, 2013). This bias would also increase the chance of Type II errors, as we are unable to assess bighorn sheep specific variants.

Recent research has suggested that a selective sweep occurred around the RXFP2 region in bighorn sheep (Kardos et al., 2015). In this scenario, multiple generations of sexual selection for large horns led to the fixation of genetic variation in the RXFP2 regions. If true, that fixation would preclude detection of associations in the current study. The region described by Kardos et al. (2015) spans ∼350,000 bp and while the Ovine Infinium^® HD SNP BeadChip contains 57 SNPs in this region, none of these loci were polymorphic in our sample of sheep from RM.

Finally, the lack of strong associations could be due to the fact that in this species these complex phenotypes are not single-locus traits. Instead, there may be many loci of small effect that jointly contribute to the phenotype, similar to the “missing heritability” phenomenon seen in many quantitative traits (Manolio et al., 2009; Yang et al., 2010). Other studies of the genetic architecture of complex phenotypes in wild populations have also found that they tend not to be controlled by single loci of large effect, but rather are polygenic (Husby et al., 2015; Bérénos et al., 2015; Kardos et al., 2016; Silva et al., 2017). New methods, such as chromosome partitioning, can now investigate this possibility (Yang et al., 2011; Robinson et al., 2013; Santure et al., 2013). Unfortunately, we cannot utilize chromosome partitioning at this time due to the small number of individuals typed on the 700k SNP chip. Attempts to use this method with our data produced unstable estimates of per-chromosome heritability (results not shown). More broadly, if these traits are truly polygenic it helps to explain how their variation is maintained despite strong directional selection (Rowe & Houle, 1996).

Conclusion

The lack of associations found here highlights the challenges of identifying genes underlying traits in non-model systems. While cross-species application of this SNP chip provided a rapid and affordable way to genotype many loci across a large number of individuals, as high-throughput sequencing costs continue to decline we expect this method to be superseded by those that allow for simultaneous marker discovery and genotyping in the species of interest (e.g., Andrews et al., 2016). Future studies could build on our findings by using high-throughput sequencing to increase the number of loci, individuals, and populations used. Improved genomic resources for bighorn sheep (Coltman, Hogg & Miller, 2013; Kardos et al., 2015; Miller et al., 2015) including whole genome sequence will enable fine mapping of associations, as well as detection of novel associations. Consideration of additional populations will allow for assessing the consistency of associations observed. In addition, haplotype-based analyses (Browning & Browning, 2011) or chromosome partitioning methods (Yang et al., 2011; Robinson et al., 2013; Santure et al., 2013) can detect novel associations and highlight if the traits fit a polygenic framework.

Supplemental Information