Development of a high-density sub-species-specific targeted SNP assay for Rocky Mountain bighorn sheep (Ovis canadensis canadensis)

Assay development

Variant site data were sourced from five studies; an implementation of a domestic sheep SNP-chip (Miller et al., 2012), RAD-seq performed on individuals from Ram Mountain, Alberta and the National Bison Range, Montana (Miller, Hogg & Coltman, 2013) (Dryad: doi: 10.5061/dryad.4qk81), whole genome sequencing performed on one individual from Ram Mountain, Alberta (Miller et al., 2015) (NCBI SRA: SRP052039), pooled whole genome sequencing of 58 individuals from Montana and Wyoming (Kardos et al., 2015) (Dryad: doi: 10.5061/dryad.3f2t2), and the application of a high density domestic sheep SNP-chip (Miller, Festa-Bianchet & Coltman, 2018) (Dryad: doi: 10.5061/dryad.c0p090f). We took the variant call data files (VCFs) from four of these studies; Kardos et al. (2015), Miller et al. (2012), (2015), and Miller, Festa-Bianchet & Coltman (2018), as these studies were already mapped to the Ovis aries 3.1 (OAR3.1) reference genome (Jiang et al., 2014) (NCBI refseq ID: GCF_000298735.1). The Miller, Hogg & Coltman (2013) study was mapped to an older domestic sheep reference genome, therefore we aligned its raw reads to the OAR3.1 reference genome, and then called variants from this alignment to produce a variant call file in a format consistent with that of the other studies. To expedite and reduce computational power required for alignment of these reads, prior to alignment we created a reduced version the OAR3.1 genome. This reduced genome only included sequence data for a 400 bp region, 200 bp upstream to 200 bp downstream, of each of the 50,000 target loci. We produced indexes from this reduced genome using the Bowtie2 v2.3.5.1 (Langmead & Salzberg, 2012) build function and the SAMtools v1.10 (Li et al., 2009) command “faidx”. Then both forward and reverse reads were aligned as pairs to the indexed reference genome using the “—sensitive” pre-set options in Bowtie2. The SAMtools command “view” was then run, with the reference index specified, to convert the paired read alignments from un-indexed .sam files to indexed .bam files. Finally, the SAMtools command “sort” was used to sort the .bam files.

We subsequently filtered the variant site data generated by aligning the Miller, Hogg & Coltman (2013) reads and variant data from the other four studies through a series of steps to reduce the number of variants and obtain optimal loci for our assay. To filter the variant site data we used VCFtools v0.1.15 (Danecek et al., 2011) and BCFtools v1.10.2 (Li et al., 2009). The filtering process is illustrated in Fig. 1. We used the same filtering processes for variants sourced from all five studies with one exception. Since the Kardos et al. (2015) sampling locations were distant from western Canadian populations, to maximize targeting loci polymorphic across populations we filtered out loci with a minor allele frequency <0.3 and only retained loci polymorphic in all three (Kardos et al. 2015) sampling locations.

Our target was to submit 100,000 candidate SNPs to Tecan for probe design, as the probe design process had ~50% success rate in designing paired probes for target loci. To reach this target we retained all the Miller et al. (2012), Miller, Hogg & Coltman (2013), and Miller, Festa-Bianchet & Coltman (2018) SNPs that passed filtering, and subsampled SNPs from the Kardos et al. (2015) and Miller et al. (2015) data which passed filtering. Subsampling of the Kardos et al. (2015) and Miller et al. (2015) SNPs was done by using the VCFtools function “—thin” at increasingly larger increments until we reduced the candidate SNP list to ~100,000 (Fig. 1). We submitted the target SNP list to Tecan as a regions file (bed format), targeting a 5 base pair (bp) region around our target SNPs (2 bp upstream, 2 bp downstream), from which they designed a genotyping assay for 50,000 loci across the Rocky Mountain bighorn sheep genome based upon the OAR3.1 domestic sheep genome. We visualised SNP distribution in a rainfall plot generated using the R (R Core Team, 2013) package karyoploteR (Gel & Serra, 2017). We also plotted chromosomal coverage and inter-SNP distances using ggplot2 v3.3.5 (Wickham, 2011).

Prior to confirming the 50,000 loci assay (the “50k assay”) design, we developed a 10,000 loci trial assay (the “10k assay”) to test the ATG technology. We genotyped a common set of individuals with both assays to assess performance. This comprised DNA samples from 56 individuals from five species and subspecies of mountain sheep. All Rocky Mountain bighorn sheep samples were typed in duplicate to test the accuracy of the assays, yielding 96 samples genotyped (Table 1). To test the performance of the assay on differing biological sample types, Rocky Mountain bighorn sheep DNA was extracted from blood, faeces, and other tissues (i.e., skin and muscle) (Table 2, Fig. 2).

We extracted DNA from differing biological material using different protocols for each material type. Biological materials were collected from wild animals under protocols approved by the University of Alberta (certificate No. 610901), University of Calgary (protocol No. BI11R-14), and Université de Sherbrooke Animal Care Committee (Protocol MFB01). Tissue samples were digested and extracted using the DNeasy Blood and Tissue kit (Qiagen N.V., Venlo, the Netherlands), following the manufacturer’s recommended procedure. Blood was treated with ammonium–chloride–potassium (ACK) (Brown, Hu & Athanasiou, 2016) prior to undergoing the same procedure as other tissue samples. For faecal samples, three pellets from each sample were soaked in 1X phosphate-buffered saline (PBS) for 20 min then swabbed with a cotton-tip applicator. The applicator was washed with Aquastool (MoBiTec GmbH, Goettingen, Germany) then DNA extracted following the Aquastool protocol. Each sample was quantified and normalized to a concentration of 8 ng/uL, resulting in a total of 80 ng of DNA processed per sample during library preparation. For the library preparation we used two Allegro Targeted Genotyping V1 kits from Tecan Genomics. The first kit targeted the 10k assay, which was a subset of SNP loci from the full 50k assay, and the second kit targeted the 50k assay. We followed the manufacturer’s protocol provided by Tecan with the fragmentation digest time increased from 15 to 22.5 min. Following library preparation, we sequenced a total of 96 samples: 80 Rocky Mountain bighorn sheep (40 individuals in duplicate), four desert bighorn sheep, four Dall sheep, four Sierra Nevada bighorn sheep, and four Stone sheep using both kits. The 10k assay was sequenced on an Illumina (San Diego, CA, USA) NextSeq using a NextSeq V2.5 mid-output 300 cycle kit which provided 120,000,000 two-directional 150 bp reads. Thus, each sample had an average coverage of 125X (120,000,000 reads ÷ 10k loci ÷ 96 samples = 125X coverage). The 50k assay samples were sequenced as part of a 384-individual library on a NextSeq V2.5 high-output 300 cycle kit providing 400,000,000 two-directional 150 bp reads. Thus, each sample had an average coverage of 20.83X (400,000,000 reads ÷ 50k loci ÷ 384 samples = 20.83X coverage).

The raw reads generated by both sequencing runs were analysed using FastQC v0.11.9 (Andrews, 2010). Raw reads were trimmed using TrimGalore v0.6.5 (Martin, 2011) with the commands “—paired” and “—phred33”. Additionally, the command “—adapter2” followed by adaptor sequence removed adaptors from the R2 read. The commands “—three_prime_clip_R1/R2 5” and “—quality 30” were used to trim five bp from the three-prime end of both the R1 and R2 reads, and trim anything with a quality less than 30. We divided our reads by sequencing run (for the 10 or 50k assay) and by subspecies prior to calling variant sites.

To align the reads to the OAR3.1 genome and call variant sites, we used Bowtie2, SAMtools, and BCFtools following the same process and parameters as previously described for realigning the Miller, Hogg & Coltman (2013) reads. Subsequently, we filtered the 10 VCF files individually. We removed indels using VCFtools, then removed loci with a quality score lower than 30 and a depth of less than six using BCFtools. We excluded loci with more than 30% missing data and loci with a minor allele frequency of less than 0.05 in Rocky Mountain bighorn sheep and 0.15 in all other species using VCFtools. Samples with more than 25% missing data were removed, and the resulting dataset was re-filtered based on the earlier minor allele frequency requirements to ensure all retained SNPs still had suitable minor allele frequencies using VCFtools. Finally, we filtered to retain loci with no more than 10% missing data using VCFtools.

Assay validation

To evaluate the efficiency of the two assays, we examined reads from the 56 genotyped individuals, 96 samples in total including duplicates, sequenced on the 10 and 50k assays. First, we examined the overall efficiency and on-target efficiency of each assay. We defined the overall efficiency as the total number of SNPs recovered relative to the number of loci targeted, which, due to the nature of the SPET technology, could be higher than 100% as novel SNPs may be discovered in the flanking regions of the targeted loci. We defined on-target efficiency as the number of targeted SNPs recovered relative to the number of loci targeted.

To evaluate accuracy, we compared genotypes between duplicate samples within and between the two assays. We used PLINK (Purcell et al., 2007) to calculate the identity by state (IBS) for all genotypes, calculated as a proportion of identical loci between two samples excluding missing data. Rocky Mountain bighorn sheep duplicate samples were compared within and between the two assays. Other subspecies were only compared between the two assays which shared common genotypes. To compare the accuracy on DNA extracted from differing biological materials, we separated the Rocky Mountain bighorn sheep samples by type; blood, faecal, and other tissue derived DNA. To further evaluate the accuracy of our assays we examined the concordance between genotypes generated from our 10K and 50K panels and the genotypes produced in Miller, Hogg & Coltman (2013) and Miller, Festa-Bianchet & Coltman (2018). To do this we identified individuals genotyped in both our and Miller, Hogg & Coltman (2013) and/or Miller, Festa-Bianchet & Coltman (2018) and then calculated the identity by state between the pairs of genotypes using PLINK.

We used PLINK to calculate observed (H_obs) and expected heterozygosity (H_exp) within sampling locations of the Rocky Mountain bighorn sheep, as well as the total sample set, using only one of the duplicated samples. We then calculated Wright’s inbreeding coefficient f = (H_exp − H_obs)/H_exp (Wright, 1922). Additionally, we calculated minor allele frequencies for all SNPs using PLINK. Nei’s (1972) standard genetic distances between Rocky Mountain bighorn sheep sampling locations were calculated in adegenet (Jombart, 2008; Jombart & Ahmed, 2011). We then used the R package vegan v2.5.7 (Oksanen et al., 2013) to perform simple Mantel tests on the genetic distances and geographic distances between sampling locations. Finally, to illustrate the ability of each assay to resolve population structure, we plotted the first two principal components extracted from principal component analyses (PCA) of Rocky Mountain bighorn sheep data using the R package adegenet. As samples were genotyped in duplicate, we expect individuals successfully genotyped in duplicate to overlap or be in close proximity on the PCA, depending on the error rate between duplicates.

Assay development

From the five studies we sourced a total of ~39.7 million variant sites. After filtering, we retained 100,000 SNP loci and submitted them to Tecan for assay development. From this list Tecan designed an assay which could target 50,000 loci. For more information on the numbers of loci sourced from each study and included in the 50k assay see Table 3, for full details of the assay designed by Tecan see electronic material 1. The 50,000 loci were relatively evenly distributed within and between chromosomes (Figs. 3 and 4) and had an average inter-loci distance of 51,663 bp ± 60,107 (Fig. 5) (based upon inter-loci distances from the Jiang et al. (2014) genome assembly).

Post-sequencing filtering

After trimming and alignment, we identified 163,176 and 769,835 variant sites from the 10 and 50k assays, respectively, in Rocky Mountain bighorn sheep. After filtering we retained 10,215 and 45,367 SNPs in the 10 and 50k datasets, respectively, of which 8,082 and 38,499 SNPs were “on-target” in the 10 and 50k SNP datasets, respectively (Table 4). For further details of SNPs retained throughout filtering steps in the other mountain sheep species and subspecies, see Table 4.

Assay validation

The 10 and 50k assays had overall efficiencies 102% and 91%, respectively, for Rocky Mountain bighorn sheep samples. On-target efficiency was lower at 81% and 77% for the 10 and 50k loci assays, respectively, based on Rocky Mountain bighorn sheep samples. For other subspecies, overall efficiency ranged from 32% to 43% and 17% to 42%, respectively, and on-target efficiency ranged from 14% to 29% and 12% to 31% for the 10 and 50k SNP assays, respectively.

IBS between Rocky Mountain bighorn sheep duplicate samples ranged from 96% to 97% within and between the 10 and 50k SNP assays (Table 4). The values of IBS between assays for other species and subspecies of mountain sheep ranged from 93% to 96% (Table 5). The IBS of genotypes from DNA extracted from blood and other tissue exceeded those for DNA from faeces (Table 6). The only sample type to lose sufficient numbers of genotypes during filtering to cause samples to fail quality control were faecal material, with the 10k and 50k assays losing 21% and 56% of samples, respectively (Table 6).

Rocky Mountain bighorn sheep samples had observed heterozygosity of 28% and 29% and expected heterozygosity of 33% and 33% for the 10 and 50k SNP assays, respectively. Wright’s inbreeding coefficient for the pooled Rocky Mountain sheep sample set was 0.133 and 0.135 for the 10k SNP and 50k SNP assays, respectively (Table 7). There was a strong pattern of isolation-by-distance among the five sampling locations of Rocky Mountain bighorn sheep in both SNP assays (Fig. 6), with geographic distance accounting for 88% and 84% of genetic distance between populations (10k assay r² = 0.88, p = 0.008; 50k assay r² = 0.84, p = 0.017), respectively. Population structure was also reflected in the distribution of samples in PCA plots (Fig. 7). Notably, the two points positioned between Cadomin Mountain and Ram Mountain are genotypes from the same individual from Ram Mountain that had Cadomin-admixed ancestry according to pedigree data (Coltman et al., 2002; Poirier et al., 2019).Furthermore, as we expected points from replicated genotypes appeared as overlapping points or in extremely close proximity.

Faecal genotyping

Both assays exhibited lower success using DNA extracted from faecal samples than DNA extracted from other sample types. The 10 and 50k assays only retained ~79% and ~43% of faecal samples, respectively, and had accuracies of 93% and 92%, respectively. These values are low relative to the overall genotyping success of 85% and accuracy of 99.9% reported by Gavriliuc et al. (2022) for ATG genotyping faecal samples from domestic horses (Equus caballus). However, there are methodological differences between the studies. First, Gavriliuc et al. (2022) used swabs from fresh faecal samples, whereas our faecal samples were in the environment for an undetermined amount of time prior to collection. Secondly, Gavriliuc et al. (2022) typed 48 samples at 279 SNPs using a sequencing kit which provided ~800,000 reads, and thus had a sequencing coverage almost three times greater than the 50k assay. With greater sequencing coverage we may have observed higher success and accuracy.

Unlike Gavriliuc et al. (2022), our samples were exposed to the environment for an undetermined amount of time. During this time samples would have been exposed to ultraviolet light and, despite most of our samples being collected in the winter, potentially above optimal temperatures. Thus, our samples may have experienced more environmental degradation than those of Gavriliuc et al. (2022). We also previously found a high rate of microsatellite genotyping failures using the same samples (Deakin et al., 2020), necessitating typing each sample in triplicate to generate reliable results. However, despite these issues, sampling fresher material is not easily done, given that bighorn sheep inhabit difficult to access terrain in remote locations.

There are two options for improving genotyping success and accuracy. The first is to increase sequencing depth. Samples in the 10 and 50k assay were sequenced with 125X and 20.8X, respectively, and faecal samples in the 10k assay had elevated accuracy (+1.8%) and genotyping success (+26%) over the 50k assay. This trend is also reflected in the accuracies and genotyping success for the two other sample types (Table 6) and relative to Gavriliuc et al. (2022). This suggests increasing sequencing coverage also increases genotyping rate and accuracy in faecal samples, thus we recommend sequencing samples prepared by the ATG protocol with as much depth as feasible. Secondly, a challenge may be the presence of polymerase chain reaction (PCR) inhibitors in DNA samples extracted from faecal samples. PCR inhibitors may not be removed in DNA purification (Morin et al., 2001) and thus may inhibit efficient amplification. The ATG protocol used in this paper uses PCR, so it is possible PCR inhibitors have contributed to the poor performance of our assays on faecal material. Rectifying this by using different methodological approaches such as DNA extraction using magnetic beads (Flagstad & Stacy, 1999), or a post extraction clean-up such as ethanol precipitation (Green & Sambrook, 2016) may further also improve genotyping success and accuracy.

Population structure

To test the utility of the assays for population genomic studies of Rocky Mountain bighorn sheep, we visualized the genetic similarity among individuals using Principal Component Analysis and by plotting isolation by distance. As expected, Rocky Mountain bighorn sheep samples grouped together by sampling location, and sampling locations separated relative to geographic separation. We also observed a strong pattern of isolation-by-distance among sampling locations using both assays, with over 80% of variance explained by geographic distance alone. This high degree of isolation-by-distance was expected, as previous analyses of bighorn sheep in this region identified a strong pattern of isolation-by-distance (Deakin et al., 2020).

Application to other mountain sheep

Our assays had varying efficiencies with other mountain sheep species and subspecies. The variation in efficiencies appears to be explained by the phylogenetics of North American wild sheep. Efficiency decreases with time since divergence between the test and target taxon. This was expected, Miller, Hogg & Coltman (2013) found that cross species applications of SNP assays yield exponentially less polymorphic SNPs as time to divergence increases. Of the sub-species of bighorn sheep, desert bighorn had the highest efficiency, which is expected given desert and Rocky Mountain bighorn sheep diverged more recently than Rocky Mountain and Sierra Nevada bighorn sheep (Buchalski et al., 2016). We see a further reduction in efficiencies for Dall and stone sheep, with efficiencies of ~33% and ~32% for the 10k SNP assay and ~21% and ~17% for the 50k SNP assay, respectively, which is expected given that bighorn and thinhorn sheep diverged ~1 mya (Rezaei et al., 2010). Despite reduced efficiency, genotype accuracy for all taxa were high, in the range of 93–96%.