From promise to practice: pairing non-invasive sampling with genomics in conservation

Sample collection

This study was conducted in the North Cascades National Park, Washington, USA (Fig. 1). Sites within this national park present the opportunity to sample American pika along steep elevational transects where climates change rapidly over short linear distances, while controlling for other environmental and historical factors. Additionally, while pika are currently abundant in the park, this area has been disproportionally affected by climate change (Karl et al., 2009). Pika populations were sampled along two elevational transects (Pyramid Peak (PP) and Thornton Lakes (TL)) between July and August 2013. Sites within transects ranged from 450 m to 1,700 m, representing an approximate 6 °C gradient in mean annual temperature (Briggs et al., 1997) over less than 6.5 km distance (Fig. 1).

Non-invasive snares were used to obtain hair samples from 12 individuals at four sites along each of the two transects (n = 96) following Henry & Russello (2011). To minimize resampling the same animal, snares were set a minimum of 15 m apart and only one sample from each snare was used. Subsequent genetic data were used to test the assumption that each sample possessed a unique genotype (see below). All samples were collected under United States Department of Interior National Park Service permit #NOCA-2014-SCI-0022 and in accordance with animal care protocol (A11-0371) as approved by the University of British Columbia’s Animal Care & Biosafety Committee.

DNA isolation, genomic data collection and SNP discovery

Total genomic DNA was extracted using the DNA IQ Tissue and Hair Extraction Kit (Promega, Madison, WI, USA) following the manufacturer’s protocol. Each sample contained 60 hair follicles with the majority of the hair shaft removed under a dissecting microscope to reduce protein and other contamination. All DNA extractions were conducted in a separate laboratory free of concentrated PCR products. Negative controls were included in each extraction to monitor contamination. DNA quantifications were conducted using real-time PCR fluorescence measurements of double stranded DNA (Blotta et al., 2005) and the Quant-it kit (Life Technologies, Foster City, CA).

Genomic DNA was converted into nextRAD genotyping-by-sequencing libraries (SNPsaurus, LLC). The nextRAD method uses selective PCR primers to amplify genomic loci consistently between samples. Genomic DNA (10 ng or less depending upon extraction yield) was first fragmented with Nextera reagent (Illumina, Inc), which also ligates short adapter sequences to the ends of the fragments. Fragmented DNA was then amplified, with one of the primers matching the adapter and extending 9 arbitrary nucleotides into the genomic DNA with the selective sequence. Thus, only fragments starting with a sequence that can be hybridized by the selective sequence of the primer will be efficiently amplified. The resulting fragments are fixed at the selective end, and have random lengths depending on the initial Nextera fragmentation. Because of this, amplified DNA from a particular locus is present at many different sizes and careful size selection of the library is not needed. For this project, an arbitrary 9-mer was chosen from those previously validated in smaller genomes, which did not appear to target repeat-masked regions in the publically available American pika genomic scaffolds (Ensembl, release 74, Ochotona_princeps.74.dna_sm.toplevel.fa) and that would approximate the results of standard RAD sequencing projects using SbfI (Baird et al., 2008).

Since these samples were collected non-invasively, it was important to assess the proportion of sequence reads in each sample that originated from the target organism relative to other environmental sources prior to conducting genotyping analysis. This was done using a custom script (SNPsaurus, LLC) that randomly sampled 1,000 high-quality reads from each sample and aligned those to the publically available American pika genomic scaffolds as well as subjected them to a blastn (Altschul et al., 1997) search of all sequences in the NCBI non-redundant database. Only samples that had greater than 50% sequencing reads that mapped to Ochotona princeps were retained for genotyping analysis.

The genotyping analysis used custom scripts (SNPsaurus, LLC) that created a reference from abundant reads present between 500 and 2,000 times across the combined set of samples, mapping all of the reads to the reference allowing two mismatches. The identified variants were then filtered by removing loci that had more than the expected maximum of two alleles and those that were present in less than 10% of all samples.

Following assembly, mapping and variant detection, the data were further filtered to maximize data quality. We retained only those loci that were genotyped in ≥50% of individuals from each transect, had a minor allele frequency ≥0.05, and a minimum coverage of 6X for homozygotes (affording 95% confidence in the genotype) while heterozygotes were required to have a minimum of 2X coverage per allele for each individual. These values were chosen to minimize null alleles and sequencing errors from biasing homozygote and heterozygote genotype calls, respectively. We then removed loci that displayed significant deviation from Hardy-Weinberg equilibrium (HWE) in more than two sites per transect as assessed using the method of Guo & Thompson (1992) as implemented in Genepop 4.3 (Raymond & Rousset, 1995; Rousset, 2008).

To ensure that only non-redundant samples were included in subsequent analyses, we conducted genotype matching across a random subset of 100 loci. We conducted the match analysis and calculated the multi-locus probability of identity (Waits, Luikart & Taberlet, 2001) for the 100 randomly chosen loci using GenAlEx (Peakall & Smouse, 2006). Only samples with unique genotypes were retained.

Outlier locus detection and annotation

Polymorphic loci were screened for statistical outliers using the Bayesian simulation method of Beaumont & Balding (2004) as implemented in Bayescan 2.1 (Foll & Gaggiotti, 2008). This analysis was run independently for each transect, with all samples coded by site (PP1-PP4, TL1-TL4). We used a prior odds value of 10, with 100,000 iterations and a burn-in of 50,000 iterations. We identified loci that were significant outliers at a q-value of 0.20. A q-value is a false discovery rate (FDR) analogue of the p-value, with the former only defined in the context of multiple testing, whereas the latter is defined on a single test. Consequently, a 20% threshold for q-values is much more stringent than a 20% threshold for p-values in classical statistics. To test for non-random association of genotypes, linkage disequilibrium was assessed between all pairs of outlier loci in each population using the exact test of Guo & Thompson (1992) and 10,000 dememorization steps, 100 batches, and 10,000 iterations per batch as implemented in Genepop 4.3 (Raymond & Rousset, 1995; Rousset, 2008). In addition, each haplotype from all nextRAD-tags that contained outlier loci were subject to a blastn (Altschul et al., 1997) search of all sequences in the NCBI non-redundant database (word size = 11; mismatch scores = 2, −3; maximum e-value = 10–15). To reduce annotations to repetitive sequences in the database, we required either a unique blastn hit or a top hit with an e-value that was at least an order of magnitude lower than the next closest hit.

Population genetic analyses

We segregated loci into two datasets including: (1) all loci identified as an outlier (“outlier dataset”); and (2) all loci not identified as an outlier (“neutral dataset”). The neutral dataset was used to conduct standard population genetic analyses for quantifying the extent and distribution of variation within and among sites. Within sites, proportion of polymorphic loci, observed (H_o) and expected (H_e) heterozygosity, and gene diversity (N_g) were calculated using Arlequin 3.5 (Excoffier & Lischer, 2010). Global tests for heterozygote deficit were conducted using Fisher’s method and 10,000 dememorization steps, 100 batches, and 10,000 iterations per batch as implemented in Genepop 4.3 (Raymond & Rousset, 1995; Rousset, 2008). The inbreeding coefficient, F_is, was calculated for each site as implemented in Genetix (Belkhir et al., 2004), with significance assessed using 1,000 permutations. To evaluate whether site-level genetic diversity was correlated with elevation and sample size, we conducted linear regression analyses implemented in R v. 3.1 (R Development Core Team, 2011).

Levels of genetic differentiation among groups were estimated by pairwise comparisons of θ (Weir & Cockerham, 1984), as calculated in Genetix (Belkhir et al., 2004), and evaluated using 1,000 permutations. The hierarchical organization of genetic variation within and among transects was calculated using an analysis of molecular variance (amova) as implemented in Arlequin 3.5 (Excoffier & Lischer, 2010), with significance assessed using 1,000 permutations. In addition, the model-based clustering method implemented in Structure 2.3.4 (Pritchard, Stephens & Donnelly, 2000) was used to infer the number of discrete genetic units across both transects. Run length was set to 100,000 MCMC replicates after a burn-in period of 100,000 using correlated allele frequencies under a straight admixture model. We varied the number of clusters (K) from 1 to 10, with 10 replicates for each value of K. The most likely number of clusters was determined by plotting the log probability of the data (ln Pr(X|K)) (Pritchard, Stephens & Donnelly, 2000) across the range of K values tested and selecting the K where the value of ln Pr(X|K) plateaued as suggested in the Structure manual. We also employed the ΔK method (Evanno, Regnaut & Goudet, 2005) as calculated in Structure Harvester (Earl, 2011). Results for the identified optimal values of K were summarized using clumpp (Jakobsson & Rosenberg, 2007) and plotted using distruct (Rosenberg, 2004). In order to test for unrecognized substructure in the broader Structure analysis, we repeated the above analysis for each transect separately using neutral and outlier loci.

Data quality

The mean starting DNA concentration recovered from the non-invasively collected hair samples was 0.55 ng/µl with as little as 1 ng total for some samples. The mean number of sequencing reads per sample was 1,863,634. Ten samples yielded less than 100,000 sequencing reads, likely due to the degraded quality and very low quantity of starting DNA. Nineteen additional samples had less than 50% of their sequencing reads mapping to O. princeps. Sixteen of these samples had high proportions of sequence reads matching with two small mammals that likely co-occur in the sampling area (Mus musculus (n = 13) and Spermophilus (n = 3)), with others matching Homo sapiens (n = 2) and Zea mays (n = 1). The above samples (n = 29) with low overall sequence reads or a low proportion mapping to O. princeps (or both) were removed, leaving 67 samples from eight sites across two elevational transects that were subject to all downstream analyses (Fig. 1 and Table 2).

We identified 9,825 SNPs that met the minimum parameters for recovering genotypes. To minimize linkage, we retained the highest coverage SNP from each contig, resulting in 3,830 SNPs. Twenty-seven SNPs deviated from HWE in two or more sites per transect and were removed from the dataset. Consequently, all downstream analyses were based on genotypic data at 3,803 SNPs. All 67 of the retained samples possessed unique genotypes at a random subset of 100 loci (average probability of identity within each sampling site = 1.1 × 10⁻²³), suggestive of unique individuals.

Outlier detection and annotation

Outlier detection identified 37 loci along the TL transect and 18 loci along the PP transect, none of which were shared. There was no evidence of significant deviation from linkage equilibrium for any pairwise comparison of outlier loci across populations. Fourteen outlier loci unambiguously matched sequences from the NCBI nr database, five of which annotated to genes of known functions (Table 1). Locus 57863_76 identified from the PP transect mapped to the receptor tyrosine kinase-like orphan receptor 2 (ROR2) gene that is part of a conserved family that function in developmental processes including skeletal and neuronal development, cell movement and cell polarity (Green, Kuntz & Sternberg, 2008). Likewise, locus 108547_114 identified from the TL transect annotated to another gene encoding a cell surface tyrosine kinase receptor (beta-type platelet-derived growth factor receptor), but for members of the platelet-derived growth factor family (Shim et al., 2010). Locus 28594_45 was similar to the laminin alpha 3 gene in humans that codes for a protein that is essential for formation and function of the basement membrane, with additional functions in regulating cell migration and mechanical signal transduction (Hamill, Paller & Jones, 2010). Lastly, locus 23486_75 was annotated to the hephaestin-like 1 (HEPHL1) gene that may function as a ferroxidase and may be involved in copper transport and homeostasis, while locus 33398_46 mapped to thioredoxin-related transmembrane protein 4 (TMX4) that may act as a reductase in the calnexin folding complex (Sugiura et al., 2010).

Population genetic analyses

The proportion of polymorphic loci varied across the sampling sites, with the lower elevation sites (PP1, PP2, TL1) exhibiting substantially lower numbers (P = 0.661–0.777) than found at the mid- and high-elevation sites (P = 0.837–0.947) for both transects (Table 2). Similar trends were seen for gene diversity along both transects where the low and mid- elevation sites recorded the lowest values relative to higher elevation sites (Table 2). Indeed, both measures of site-level genetic variation were significantly correlated with elevation (P: r² = 0.557, p = 0.034; N_g: r² = 0.738 p = 0.006; Fig. 2). Although P significantly correlated with sample size (r² = 0.635, p = 0.018), this was not the case for gene diversity (r² = 0.0813, p = 0.493) or elevation (r² = 0.184, p = 0.289). This general trend of increasing variation with elevation seemed to hold for observed heterozygosity along the PP transect, but were stable across TL sites (H_o: 0.336–0.368; Table 2). Yet, all low (PP1, TL1) and mid-low (PP2, TL2) sites exhibited significant, genome-wide evidence of heterozygote deficit. Interestingly, significant inbreeding was also detected at the low and mid-low elevation sites along PP, with no such evidence at the higher elevation sites (PP3, PP4; Table 2). All sites along the TL transect exhibited evidence of inbreeding (Table 2).

The AMOVA revealed that a significant amount of variation (p < 0.0001) was exhibited both among transect (4.14%, d.f. = 1) and among sites within transect (2.01%, d.f. = 6), with the remaining found within populations (93.05%, d.f. = 126). These patterns were congruent with those from pairwise θ estimates, with the highest values generally displayed by among transect comparisons, but where all comparisons were significant (Table 3).

The Bayesian clustering analyses based on 3,748 neutral loci revealed strong evidence for two clusters within the dataset (ΔK₂ = 249.0), corresponding to the two transects (Fig. 1). When analyzing the PP transect separately, additional substructure (K = 2; Fig. 3) was found using neutral (ΔK₂ = 33.1) and outlier (ΔK₂ = 123.1) loci. In both cases, the low elevation site (PP1) represents a largely distinct cluster relative to all other sites. Similarly, substructure was found along the TL transect when using outlier loci, with strong evidence for two (ΔK₂ = 473.3) and three (ΔK₃ = 314.6) clusters. In the K = 2 plot, TL2 represented a distinct cluster, while in the K = 3 plot, the lower elevation sites (TL1, TL2) were each separate clusters relative to the high elevation sites. Although the Evanno, Regnaut & Goudet (2005) method would favor K = 2, the method described by Pritchard, Stephens & Donnelly (2000) for inferring the optimal number of clusters would suggest K = 3 given that ln Pr(X|K) clearly plateaus at this value (Table S1). No substructure was found along the TL transect based on neutral loci.