Genotyping-in-thousands by sequencing (GT-seq) of noninvasive faecal and degraded samples: A new panel to enable ongoing monitoring of Canadian polar bear populations.

Genetic monitoring using noninvasive samples provides a complement or alternative to traditional population monitoring methods. However, next-generation sequencing approaches to monitoring typically require high quality DNA and the use of noninvasive samples (e.g., scat) is often challenged by poor DNA quality and contamination by nontarget species. One promising solution is a highly multiplexed sequencing approach called genotyping-in-thousands by sequencing (GT-seq), which can enable cost-efficient genomics-based monitoring for populations based on noninvasively collected samples. Here, we develop and validate a GT-seq panel of 324 single nucleotide polymorphisms (SNPs) optimized for genotyping of polar bears based on DNA from noninvasively collected faecal samples. We demonstrate (1) successful GT-seq genotyping of DNA from a range of sample sources, including successful genotyping (>50% loci) of 62.9% of noninvasively collected faecal samples determined to contain polar bear DNA; and (2) that we can reliably differentiate individuals, ascertain sex, assess relatedness, and resolve population structure of Canadian polar bear subpopulations based on a GT-seq panel of 324 SNPs. Our GT-seq data reveal spatial-genetic patterns similar to previous polar bear studies but at lesser cost per sample and through use of noninvasively collected samples, indicating the potential of this approach for population monitoring. This GT-seq panel provides the foundation for a noninvasive toolkit for polar bear monitoring and can contribute to community-based programmes - a framework which may serve as a model for wildlife conservation and management for species worldwide.

INTRODUCTION

Informed wildlife management requires accurate demographic data that include recurring, reliable population estimates and an understanding of the factors that shape population dynamics (Durner et al., 2018; Hamilton & Derocher, 2019; Laidre et al., 2015). Population monitoring is especially urgent for species that are impacted by rapid climate change, are harvested or poached, or reside in habitats that have been heavily altered by human activity (Durner et al., 2018; Laidre et al., 2015; Robinson et al., 2014). Such monitoring can be challenging for species that have large territories, occupy inaccessible habitat, are cryptic or elusive (e.g., nocturnal, fossorial), or are of heightened conservation concern. For such species, capture, direct handling, and invasive sampling may be impractical, inappropriate, or culturally undesirable. Moreover, traditional methods for monitoring of animal populations (e.g., aerial censusing, mark‐recapture, and radiotelemetry) can be expensive, time‐consuming, and stressful for the focal animals (Van Coeverden de Groot et al., 2013; Solberg et al., 2006; Stapleton et al., 2014). When populations have low densities, mark‐recapture and aerial surveys may also be hindered by low probabilities of capture and detection, respectively (Garshelis & Noyce, 2006; Hayward et al., 2002).

Genetic monitoring using noninvasive samples, such as scat, hair, feathers, or skin, affords an alternative that can mitigate some of the challenges of traditional monitoring. Noninvasive genetic monitoring can be deployed on a large scale and with greater frequency, potentially enabling larger sample sizes and improved temporal monitoring. Other benefits include ease of collection, reduced disturbance of study species, potentially decreased spatial or temporal biases, and diminished physical risk to collectors (Carroll et al., 2018; Morin et al., 2018; Steyer et al., 2016; Waits & Paetkau, 2005). Noninvasive genetic monitoring can provide robust and repeatable data for individual identification (including sex), movement, and population trends (e.g. Aziz et al., 2017; Quinn et al.,2019; Schmidt et al., 2020; Schultz et al.,2018). Scat samples in particular can be useful for assessing the health of individuals and populations, as they offer a range of other information, including data on parasite and pathogen presence (Bergner et al., 2019; Cristescu et al., 2019; Weese et al., 2019), diet composition (e.g., via quantitative PCR or DNA metabarcoding – Iversen et al., 2013; Nelms et al., 2019; Ogurtsov, 2018), hormone profiles (Morden et al., 2011; Vynne et al., 2012), and contaminant loads (Lundin et al., 2015, 2016).

Despite the purported advantages of noninvasive genetic monitoring, it poses challenges that have limited its widespread implementation. DNA from noninvasive samples may be degraded due to environmental exposure (Bourgeois et al., 2019; Poinar et al., 1996; Schultz et al., 2018) and contaminated by nontarget species (Carroll et al., 2018; Taberlet et al., 1999). Often less than 5% of the total DNA that scat contains is host DNA, with most DNA coming from pathogens, parasites, commensal bacteria, prey, and off‐target species (Han et al., 2019; Perry et al., 2010; Snyder‐Mackler et al., 2016). Due to low quality and quantity of host DNA, accurate quantification and genotyping of noninvasive samples using next‐generation sequencing (NGS) methods have often proved difficult, as they require high DNA concentration and quality (Graham et al., 2015; Maroso et al., 2018; Taberlet et al., 1999). Such reduced genotyping accuracy due to contaminated and degraded DNA may increase processing efforts and costs, and complicate inferences from collected data.

Multiple methods that have been developed to improve NGS of low‐quality samples (e.g., scat, hair, or archaeological samples) use selectively targeted, species‐specific arrays of single nucleotide polymorphism markers or SNPs (Carroll et al., 2018). Traditional NGS methods (e.g., double‐digest restriction‐associated DNA sequencing: ddRADseq) can be used to identify large SNP panels across a focal species genome, from which smaller, informative panels can be selected (Andrews et al., 2018; Blåhed et al., 2018; Hess et al., 2015). DNA capture, SNP genotyping assays (e.g., TaqMan), or amplicon sequencing methods can then be used to genotype the reduced panel with high coverage (reviewed in Meek & Larson, 2019). Indeed, SNP genotyping has been applied to noninvasive samples from a range of wildlife species and has yielded high genotyping success and low genotyping error, reducing the need for systematic replicates that increase cost and effort (Fitak et al., 2016; Kleinman‐Ruiz et al., 2017; Kraus et al., 2015; Schultz et al., 2018; von Thaden et al., 2017).

Genotyping‐in‐thousands by sequencing (GT‐seq) is among the most promising approaches for genotyping noninvasively collected DNA. It uses highly multiplexed PCR to amplify short amplicons, followed by individual barcoding that allows rapid, high‐quality genotyping of targeted SNP panels across thousands of individuals (Campbell et al., 2015). GT‐seq library preparation can be performed with standard laboratory equipment and shows decreased genotyping error and genotyping costs relative to other NGS‐based genotyping methods, including TaqMan assays (Campbell et al., 2015). Because GT‐seq uses suites of multiplexed, species‐specific primers, it may mitigate some of the challenges presented by exogenous DNA and degraded host DNA in noninvasive faecal samples. Combined with new Illumina platforms like NovaSeq, GT‐seq costs can be decreased further by running up to thousands of individuals on a single lane. GT‐seq has been successfully applied to minimally invasive cloacal swab DNA samples collected from western rattlesnakes (Crotalus oreganus) with low rates of genotyping error and discordance relative to RADseq (Schmidt et al., 2020). Natesh et al. (2019) also found high genotyping success using GT‐seq of noninvasive scat samples for Indian Bengal tigers (Panthera tigris tigris). Thus, GT‐seq has the potential to enable efficient and economical genetic monitoring of populations based on noninvasively collected samples, but for implementation, requires a clear guide for development, testing, and validation.

There is a long‐expressed desire by Northern communities in Canada for monitoring practices based on noninvasively collected samples, particularly for polar bears (Ursus maritimus). Current polar bear monitoring is based primarily on surveys of 19 subpopulations (also called management units, MUs), designated largely using mark‐recapture, radio collaring, and aerial surveys (e.g. Stapleton et al., 2014; Taylor et al., 2009). Northern communities, including the Inuit – for whom polar bears are of key cultural and economic importance – have voiced concern about the invasiveness of some of these methods (e.g., mark‐recapture), potential negative impacts on polar bear health and behaviour, lack of inclusion of Traditional Ecological Knowledge in monitoring and management, and lack of collaboration in management activities (Wong et al., 2017; York et al., 2016). Two‐thirds (13) of the 19 polar bear subpopulations are fully or partially found in Canada (Figure 1), highlighting Canada's need to lead on polar bear management. However, surveys for many subpopulations are infrequent due to logistical and economic constraints, and 11 of 19 subpopulations are data deficient (Government of Canada, 2018; Hamilton & Derocher, 2019). As of 2019, only six subpopulations had population estimates <10 years old (Hamilton & Derocher, 2019). These data deficiencies preclude robust estimates of population parameters and have limited implementation of effective management strategies. Thus, as polar bears continue to be impacted by climate change and face limitations in range and prey availability as rapid sea ice decline continues (Fontúrbel et al., 2018; Hamilton & Derocher, 2019; Hunter et al., 2010; Rode et al., 2014), there is an urgent need for new monitoring approaches.

FIGURE 1

Map of the Canadian Arctic showing currently recognized subpopulations that are fully or partially in Canada (solid line polygons). Subpopulation abbreviations are the same as in Table 2. Points correspond to sampled individuals, coloured according to genetic cluster assignment based on structure analysis (pink = Polar Basin, green = M’Clintock Channel, orange = Arctic Archipelago, Blue = Hudson Complex). Black points represent individuals with membership <0.7 to a genetic cluster

TABLE 2

Diversity metrics for 11 Canadian polar bear subpopulations (see Figure S1) based on a combined GT‐seq+ddRADseq data set consisting of 642 individuals genotyped at 322 autosomal loci

Subpopulation	n	H o	H E	GIS	p (GIS>0)	Self‐assignment	Main genetic cluster
Baffin Bay (BB)	50	0.45	0.45	−0.006	.23	0	Arctic Archipelago
Davis Strait (DS)	40	0.45	0.45	0.004	.311	0.01	Arctic Archipelago
Foxe Basin (FB)	82	0.45	0.44	−0.025	.001	0.52	Hudson Complex
Gulf of Boothia (GB)	110	0.45	0.45	0.003	.282	0.73	Arctic Archipelago
Lancaster Sound (LS)	83	0.45	0.45	−0.009	.083	0.03	Arctic Archipelago
M’Clintock Channel (MC)	80	0.45	0.44	−0.023	.001	0.56	M’Clintock Channel
Northern Beaufort Sea (NB)	55	0.44	0.43	–0.019	.01	0.04	Polar Basin
Norwegian Bay (NW)	1	–	–	–	–	–	Arctic Archipelago
Southern Beaufort Sea (SB)	12	0.45	0.42	−0.070	.001	0	Polar Basin
Southern Hudson Bay (SH)	80	0.44	0.43	−0.017	.003	0.65	Hudson Complex
Viscount Melville Sound (VM)	11	0.45	0.44	−0.018	.169	0	M’Clintock Channel/Polar Basin
Western Hudson Bay (WH)	38	0.44	0.43	−0.017	.037	0	Hudson Complex

Abbreviations: G IS, inbreeding coefficient; H E, expected heterozygosity; H o, observed heterozygosity; n, sample size.

Noninvasive scat surveys may enable more direct community participation, and provide a cost‐effective complement to traditional polar bear monitoring methods. Noninvasive scat samples could be obtained through community‐level monitoring programmes, in which Inuit hunters are remunerated for field sampling efforts and surveys are repeated regularly to better track the trajectory of polar bear subpopulations. Noninvasive scat surveys have already been used in brown bears (Ursus arctos) as an alternative source of DNA to high quality samples, such as muscle (e.g. Giangregorio et al., 2019). We have also established that sufficient DNA for GT‐seq protocols can be obtained and quantified from field‐collected polar bear scat (Hayward et al., 2020). Thus, polar bears present an opportunity to demonstrate GT‐seq panel development and validation, application, and usefulness for population monitoring. Importantly, this application of GT‐seq will be in collaboration with Northern Canadian communities and are predicted to have real socioeconomic benefits for the communities involved.

In this study, we test the practicality of using GT‐seq for SNP genotyping of noninvasively collected faecal samples from polar bears and apply GT‐seq to population monitoring, using degraded samples that heretofore could not be genotyped with NGS‐based methods, to expand our understanding of Canadian polar bear genetic structure. We show that our optimized GT‐seq panel can be used to distinguish among individuals (including noninvasively collected samples), and to characterize genetic structure of the Canadian polar bear population.

MATERIALS AND METHODS

GT‐seq panel development

To identify potential SNPs for our GT‐seq panel, we screened 411,094 SNP loci identified from ddRADseq data generated from 327 polar bears in Jensen et al. (2020). The distribution of these samples is depicted in Figure 1 and Table 1 of Jensen et al. (2020). We filtered this data set to retain loci that were successfully genotyped in ≥85% of individuals and had a minor allele frequency of at least 0.25. Retained loci were thinned to one SNP per 50,000 bp using the polar bear reference genome (GenBank GCA_000687225.1; Liu et al., 2014) to reduce possibility of linkage (Table S1 presents the full filtering workflow). This resulted in 442 retained SNP loci. To allow for sex identification, we added two SNPs known to be sex‐linked in polar bears (Pagès et al., 2009). We validated our GT‐seq panel of 442 SNPs and designed primers as described in Supporting Information: GT‐seq Panel Development and Validation. To ensure our potential panel could distinguish individuals with high reliability, validation included calculating probability of identity (PID) and probability of identity (siblings) for an increasing number of loci from our ddRADseq data. We could reliably identify individuals with as few as 34 SNPs, but included a panel of >300 SNPs to achieve both individual identification and sufficient resolution for population genetic analyses.

TABLE 1

GT‐seq genotyping success for five different sources of polar bear DNA, including 365 total samples collected across 10 subpopulations. Percent individuals successfully genotyped before and after samples with no detectable polar bear DNA were removed (based on qPCR results) are shown, and average percent missing data by locus for each sample type (individuals with >50% missing data removed)

Sample type	n	% individuals >50% data	% individuals >50% data after removal of samples qPCR = 0 ng/µl	Average % missing data by locus (range)
Set muscle (MS)	101	96.2	–	2.2 (0.31–45.3)
Colon faeces (CF)	69	88.5	90.1	2.1 (0.31–23.3)
Biopsy (BP)	134	95.7	–	3.1 (0.31–39.1)
Harvest muscle (HV)	38	97.4	–	1.1 (0.31–2.5)
Field faeces (FF)	23	30.6	62.9	14.9 (0.31–50.0)

Set muscle, tissue from corresponding muscle and colon sets; Colon faeces, faeces removed from the colon of corresponding muscle and colon sets; Biopsy, biopsy tissue sample; Harvest muscle, tissue from annual polar bear harvest; Field faeces, noninvasively collected scat from the field.

Sample collection

Our work draws on multiple sample sources (Figure S1), including archived harvest tissue samples (HV; n = 38) from bears taken by Inuit hunters, and biopsy samples (BP; n = 138) that are housed in collections of the Nunavut and Northwest Territories governments in accordance with Government of the Northwest Territories and Government of Nunavut research permits. We also received “sample sets” from hunted bears, which contained fat, tissue, liver, and lower intestine with feces, and were collected in accordance with wildlife research permits ARI #WL 500540 to MB and WL‐2019–061 to SCL. These were used to estimate faecal genotyping error by comparing genotypes from set tissue (MS; n = 108) and faeces from the colon (CF; n = 78), with expectation that these estimates will be conservative. Noninvasive, field‐collected faecal samples (FF; n = 72) were also located and collected by Inuit hunters under wildlife research permit WL‐2018–006 to SCL. All sample types were stored in –20°C or –80°C at Queen's University until subsampling. Many of these samples were too degraded for use in previous NGS‐based studies and thus, represent opportunities to expand sampling coverage and draw new inferences regarding population structure. As the ultimate goal is to use our GT‐seq method on field‐collected scat samples, we took duplicate subsamples from 21 FF samples to assess within‐sample genotyping error and variation in genotype quality. Subsamples of MS, HV, and BP were stored in 100% ethanol at –20°C or –80°C and CF and FF subsamples were stored without ethanol at –20°C until DNA extraction.

All samples were collected across 11 of the 13 Canadian subpopulations between 1998 and 2019, with mean year of sample collection withing subpopulations ranging from 1999 (M’Clintock Channel) to 2018 (Southern Hudson Bay; Figure 1). For nine subpopulations, we had at least 10 sampled individuals (range: 10–95). However, there were only four samples from Davis Strait and 3 samples from Western Hudson Bay, and no samples from either Norwegian Bay or Kane Basin. While this represents a sampling limitation, we supplement our GT‐seq data with ddRADseq data (Jensen et al., 2020) for our final assessment of Canadian polar bear structure and diversity, and note that census population size estimates for the two subpopulations for which we have no samples are small (KB = 357 individuals, NW = 203 individuals; Hamilton & Derocher, 2019).

DNA extraction

Whole genomic DNA was extracted from faecal samples using the QIAamp Fast DNA Stool Mini Kit (Qiagen) according to the manufacturer's protocols. For BP, HV, and MS samples, genomic DNA was extracted using a modified salt extraction protocol (Aljanabi & Martinez, 1997), with an RNaseA (Thermofisher Scientific) step included. Once extracted, DNA extracts from all tissue samples were run on 1.5% agarose gel stained with RedSafe Nucleic Acid Staining Solution (iNtRON Biotechnology) to assess quality and quantified using a Nanodrop ND_1000 spectrophotometer (Nanodrop Technologies Inc.). We used a polar bear‐specific qPCR assay targeting the F2 gene to quantify the amount of polar bear DNA in both CF and FF samples (Hayward et al., 2020). To gauge the value of running faecal samples through this qPCR assay as a screening tool, we devised a small double‐blind experiment in which we randomly divided FF samples into two subsets: 1. 8 samples for which we assayed DNA quantity before GT‐seq library construction and sequencing; and 2. 89 samples (+21 duplicates) for which we assessed DNA quantity only after sequencing had already been performed (see Supporting Information: qPCR Experiment for qPCR experiment methods and results).

GT‐seq genotyping and genotype calling methods

A full description of GT‐seq panel optimization can be found in Supporting Information: Panel Optimization. Using optimized GT‐seq conditions and primers for our final GT‐seq panel, we prepared two libraries for faecal samples (CF, FF) and another three for tissue samples (HV, BP, MS). The faecal libraries contained all 150 faecal samples, as well as the 21 FF duplicates to assess subsample variation in genotyping error. We included 284 tissue samples in the tissue libraries and two technical replicates (BP) to assess genotyping error for GT‐seq. Genotyping error was measured conservatively as percent discordance between technical replicates. Within our total of 457 samples, we included 65 paired sets of MS and CF samples collected from the same individuals allowing us to compare GT‐seq genotyping error between sample types and obtain a “best case scenario” estimate of genotyping error for FF samples. Here, we assume that the muscle genotype is “correct” and calculate the genotyping error (percent discordance) for CF as a proxy for FF samples. Prior to calculating genotyping error within and among sample types, samples with >50% missing data were removed from the GT‐seq data set, as we considered samples with >50% data (i.e. genotypes for at least half of panel loci) to have been successfully genotyped. We considered a sample of 161 SNPs (50%) genotyped to be a “success” because it enabled sufficient resolution for both individual identification and analysis of population structure, although far fewer SNPs can be used to reliably distinguish among individuals (Supporting Information Results: GT‐seq panel validation). To assess the validity of the sex identities provided by our two GT‐seq SNPs, we compared the GT‐seq determined polar bear sex to hunter‐provided sex for 293 samples that were successfully genotyped at >50% loci and for which we had field data.

Library preparation followed the original protocols of Campbell et al. (2015), modified based on our pilot tests (see Supporting Information: Panel Optimization), and using only primers for a final optimized panel of 327 SNPs. Libraries were sequenced using an Illumina MiSeq at Queen's University. We used the GT‐seq pipeline available on github (https://github.com/GTseq/GTseq‐Pipeline) for filtering to a minimum depth of 10 and genotype calling, as suggested by Campbell et al. (2015). However, we were also interested in examining genotype discordance between different SNP calling models, as some discrepancy has been found between the GT‐seq pipeline and other workflows used to call RADseq data (Schmidt et al., 2020). As bcftools (Li, 2011) is a common genotyping tool Qnd was used to call genotypes for our original ddRADseq data set (Jensen et al., 2020), we used this same workflow to process our raw GT‐seq sequencing reads. Briefly, reads from GT‐seq were aligned to the polar bear reference genome (assembly version ursmar_1.0, PMID: 24813606) using the bwa‐mem v0.7.17 aligner (Li & Durbin, 2009). Alignments were sorted, indexed, and read pairs were fixed using tools from the samtools v1.9 suite (Li et al., 2009). A target file of the GT‐seq assay SNP positions was used with bcftools mpileup to call and produce a VCF file of the targeted sites. The VCF was filtered in vcftools (Danecek et al., 2011) for a minimum depth of either six (data set BCF‐6) or 10 (data set BCF‐10). Output files for each workflow (GT‐seq, BCF‐6, BCF‐10) were converted to genpop format to compare percent missing data by locus. Average mismatch rates for all 457 samples were also calculated between the GT‐seq pipeline and the bcftools workflows. The data set with the least missing data (BCF‐6) was used to calculate genotyping error and percent missing data, and for all subsequent population genetic analyses.

Population genetic analyses

To evaluate the usefulness and power of GT‐seq for population genetic analysis and further extend our understanding of polar bear population structure, we estimated several metrics of genetic diversity and structure using a combined data set called and filtered with our BCF‐6 workflow. This data set comprised individuals genotyped at 322 autosomal SNPs using GT‐seq, as well as individuals previously genotyped for these same SNPs using ddRADseq (Jensen et al., 2020). After combining GT‐seq and ddRADseq data, we removed all replicates (similarity >0.8, known duplicate subsamples, tissue from tissue‐colon sets), and all individuals with >50% missing data. For individuals that had been included in both data sets, only GT‐seq replicates were retained. Our final GT‐seq+ddRADseq data set contained 642 individuals genotyped at 322 autosomal loci and 2 sex‐linked loci.

For the combined GT‐seq + ddRADseq data set, we estimated observed and expected heterozygosity (H o/H E) and inbreeding coefficients (G is) for each subpopulation, as implemented in genodive v3.04 (Meirmans, 2020). Using the related package in r (Pew et al., 2015), we used our empirical allele frequencies to simulate pairs of individuals with known relatedness. One hundred pairs were simulated for each of the following categories: unrelated individuals, half siblings, full siblings, and parent‐offspring pairs. We used this simulated data set to test the ability of multiple relatedness estimators to distinguish among relatedness categories. We also calculated and plotted pairwise relatedness for all 642 individuals in our combined data set (Figure S4). To assess population substructure, we used Bayesian clustering analysis as implemented in structure 2.3.4 (Pritchard et al., 2000). We evaluated the number of clusters in the data (K) from 1 to 10, with 20 iterations of each, and with a run length of 300,000 MCMC following a burnin period of 100,000 MCMC. The most likely value of K was identified using the DeltaK method of Evanno et al. (2005) in structure harvester (Earl & vonHoldt, 2012). To complement the structure analysis, we also used discriminant analysis of principal components (DAPC; Jombart et al., 2010), implemented in adegenet in r (Jombart & Ahmed, 2011) to evaluate the number of genetic clusters in the data. The best‐fit value of K was tested using the find.clusters function, retaining all PCs, and the Bayesian information criterion (BIC). The chosen value of K was selected based on the lowest BIC value, and a DAPC plot was generated based on this clustering, retaining sufficient PCs to capture 80% of the variance. Principal component analysis (PCA) was also performed using the r package ade4 v. 1.7–16 (Dray & Dufour, 2007) with a priori groups assigned based on sample type (MS, HV, BP, CF, FF) and data type (GT‐seq, ddRADseq) to confirm structuring was independent of sample and data types (Figure S5).

To estimate assignment accuracy to each subpopulation and structure‐derived genetic cluster, we used principal components analyses and Monte Carlo cross‐validation procedures implemented in the assignpop package in r (Chen et al., 2018). For the subpopulation version of analysis, we removed subpopulations with <10 samples (Norwegian Bay), whereas for the genetic cluster version, we only retained individuals with >0.7 membership to a single genetic cluster. From these data we built a predictive model using a support vector machine (model svm) classification based on training sets composed of the most informative 75% loci and a random sample of 75% of individuals in the data set. The remaining 25% of individuals were used to test the rate of assignment, which was then averaged across 30 iterations.

RESULTS

Calling workflow comparisons

Our final GT‐seq data set, called using the published GT‐seq pipeline (Campbell et al., 2015), included 325 autosomal SNPs and two sex‐linked markers. An additional three loci were removed from the bcftools workflow data sets after filtering for minimum depth (depth 6, 10), leaving each with 322 autosomal loci and two sex‐linked markers. We removed the same three loci from our GT‐seq pipeline data set to enable direct comparison of genotypes and missing data by locus across the calling methods. Based on all 457 samples, there was an average of 25.4% missing data for the GT‐seq calling pipeline, whereas missing data were 23.9% and 21.3% for BCF‐10 and BCF‐6 calling workflows, respectively. Regardless of bcftools calling workflow, genotype mismatch with the GT‐seq workflow was 1.1% (SD = ±5.3%) on average. Based on these results and the potential for easy comparison with existing ddRADseq data, we chose to use the data set generated from the BCF‐6 calling workflow to assess genotyping error and analyse population structure.

Genotyping success

After removing individuals with >50% missing data, a total of 365 samples, including three duplicate samples, remained in our GT‐seq data set. Average missing data dropped to 3.2% from 21.3% when these individuals were removed. Percent individuals successfully genotyped ranged from 30.6% to 97.4% depending on sample type (Table 1), recalling that most of the FF samples were not screened before sequencing to allow for testing of the qPCR screening tool (Hayward et al., 2020). When considering only the FF samples that would have passed screening (>0 ng/µl polar bear DNA), percent individuals for which we had sufficient SNPs to estimate genetic identity (>34 SNPs) was 80.0% and percent individuals “successfully genotyped’” as per our cutoff (>50% loci genotyped) was 62.9%. See Supporting Information: qPCR Experiment for detailed results.

For retained individuals with >50% data, average missing data per locus was similar across sample types, with FF samples missing the most at 14.9% of loci (Table 1). Both BP replicates were successfully genotyped at >50% loci giving a conservative mean genotyping error rate of 0.17%. Of the 21 FF samples for which duplicate subsamples were included in the GT‐seq runs without qPCR prescreening, one or both subsamples failed to amplify in 20 of 21 cases. In most cases (15/21), both subsamples failed to amplify, although there were some cases with large differences in genotyping success between FF duplicates. For example, one FF subsample was genotyped at only eight loci, whereas its corresponding subsample was genotyped at all 322 loci. For the one FF sample for which both subsamples successfully genotyped, the genotyping error was 10.7%. Fifty‐six of the 65 muscle‐colon faeces sets had both sample types from the same individual successfully genotyped, giving a mean genotyping error estimate of 6.8% for faecal samples. Of the 322 retained autosomal loci, 2 had >50% missing data across individuals. Removal of these loci in future GT‐seq runs may increase genotyping success, particularly for field‐collected faeces.

We had a total of 293 bear samples for which hunters provided a sex identity. Of these, 78.6% were genotyped successfully at one or both sex‐linked loci. There were 18 instances where a sample had only a genotype at one of the sex‐linked loci, and only one case of allelic dropout. For samples with genotypes at both sex loci, genotype concordance was 99.7%. Hunter identification of bear sex matched the sex genotype provided by the first sex locus for 95.9% of samples, matched the second sex locus for 91.1% of samples, and matched both sex loci for 90.1% of samples.

Population structure

Based on the final combined data set of 642 individuals (GT‐seq+ddRADseq data at 322 GT‐seq loci), all subpopulations display similar levels of genetic diversity (Table 2). Expected heterozygosity ranged from 0.42 to 0.45 and GIS values were mostly negative, ranging from −0.006 to 0.004.

We discovered 13 true recaptures in our combined data set (relatedness >0.80), which were removed from all subsequent analysis. Recaptures spanned sample types, with pairs existing between or within sample types: FF (n = 8), CF (n = 4), MS (n = 5), BP (n = 4), and HV (n = 4). For example, one bear identified from a BP sample (nonlethal) was later detected using a CF sample (lethal). Based on comparison of multiple relatedness estimators, we determined that the quellergt correlation coefficient was best for our combined data set (correlation between observed and expected values = 0.97; Queller & Goodnight, 1989). Results from simulations with the quellergt coefficient yielded pairwise relatedness density plots with marked separation for full siblings, half siblings and parent‐offspring pairs, and unrelated individuals (Figure 2). Although relatedness among individuals in our data set was mostly unknown prior to analysis, we can see >50 empirical relatedness values of approximately 0.5 (Figure S4), as would be expected for parent‐ offspring or sibling pairs.

FIGURE 2

Density plot of pairwise quellergt relatedness values generated from simulations (100 per relatedness category) based on allele frequencies from our combined GT‐seq+ddRADseq data set of 642 individuals genotyped at 322 autosomal loci. Each colour represents a different relatedness category (Full sibs, full siblings; Half sibs, half siblings; P‐O, parent‐offspring; Unrelated, unrelated pair)

Our structure and dapc analyses revealed similar patterns in population genetic structure. structure analysis suggested an optimal value of K = 4, with ln Pr(X|K) plateauing around K = 4 (Figure 3d) and deltaK greatest at K = 4 (Figure 3e). For the find.clusters dapc analysis, the lowest BIC scores occurred at K = 3 and K = 4 (Figure 3c). For both these analyses, the genetic clusters at K = 4 correspond to three geographic regions typically referred to as the Hudson Complex, the Arctic Archipelago, and the Polar Basin, with an additional cluster corresponding to the subpopulation of M’Clintock Channel (Figure 3a,b, Figure 1). Barplots showing K = 3 and K = 5 are provided in Figure S6.

FIGURE 3

Results of two clustering analyses performed with combined GT‐seq+ddRADseq data set of 642 bears genotyped at 322 autosomal loci. (a) find.clusters discriminant analysis of principal components (DAPC) displaying inferred clustering at K = 4. Each point represents an individual bear and inertial ellipses surround each genetic cluster. (b) structure barplot showing inferred clustering at K = 4. Each colour corresponds to a distinct genetic cluster and each vertical bar represents an individual and their proportional membership in each cluster. (c) Plot of Bayesian information criterion (BIC) from DAPC analysis for each number of clusters evaluated. (d) Plot of ln P(K) for each number of clusters evaluated in the structure analysis. (e) Plot of deltaK for each number of clusters evaluated in the structure analysis

Despite 90 individuals having values of less than 0.7 membership to a particular genetic cluster (black dots in Figure 1), most of the 642 individuals in our combined ddRADseq and GT‐seq data set were assigned to a single cluster. One exception was Viscount Melville (VM), which comprised individuals with high assignment probabilities to both the Polar Basin and M’Clintock Channel clusters. The self‐assignment test performed by subpopulation (Table 2) suggests that few subpopulations are highly genetically distinguishable. The Gulf of Boothia (GB) and Southern Hudson Bay (SH) subpopulations displayed the highest self‐assignment rates at 0.73 and 0.65, respectively. Our second self‐assignment test shows that the genetic clusters suggested by our structure analysis are typically more highly distinguishable than the subpopulations, with self‐assignment rates >0.80 for three of the clusters (Table 3). Self‐assignment was lowest for the Polar Basin cluster (0.35).

TABLE 3

Results for a self‐assignment test performed in AssignPOP in R for the four genetic clusters suggested by our structure analysis. Individuals with <0.7 membership to one cluster were removed from our combined GT‐seq+ddRADseq data set of 642 bears genotyped at 322 autosomal loci

Genetic cluster	n	Self‐assignment
Arctic Archipelago	230	0.98
Polar Basin	66	0.35
Hudson Complex	188	1.00
M’Clintock Channel	68	0.83

n, sample size.

DISCUSSION

Noninvasive samples (e.g., faeces) are viable sources of DNA for new genetic methods of population monitoring, which may help mitigate limitations associated with traditional approaches. Here, we developed a final GT‐seq assay of 324 SNPs optimized to genotype polar bears based on DNA from noninvasively collected faecal samples. We demonstrated (1) successful GT‐seq genotyping of DNA from field‐collected polar bear faeces, and (2) the practicality of GT‐seq for distinguishing individuals, assessing relatedness, and expanding our understanding of genetic population structure and diversity.

Individual identification

We determined that all existing unrelated bears could be distinguished with as few as 34 SNPs (Figure S2), and that 80% of field‐collected scat with detectable polar bear DNA could be assigned a genetic identity (i.e., >34 SNPs genotyped). Given the current global polar bear population estimate of 23,315 (Hamilton & Derocher, 2019), our panel of 322 autosomal and two sex‐linked SNPs is sufficient to reliably discern bears among most relationship categories (e.g., siblings, unrelated) with high certainty (Figure 2). Given this ability to differentiate individuals based on both faecal samples and tissue, our GT‐seq assay may be particularly useful for genetic mark‐recapture studies using data from traditional sources and noninvasive samples (e.g., scat, urine, hair snags).

Comparing workflows, error rates, and genotyping success

GT‐seq uses a reduced panel of SNPs to genotype samples with degraded DNA, such as scat or archived tissue samples, and as long as there are loci in common, GT‐seq data can be combined with data from other methods (e.g., ddRADseq). However, differences in genotyping calling methods between GT‐seq and RADseq data sets may contribute to discordance among data sets (Schmidt et al., 2020). To address this, we called genotypes from our GT‐seq data using both the original calling pipeline for GT‐seq (Campbell et al., 2015) and a bcftools workflow previously used to call genotypes for our ddRADseq data (Jensen et al., 2020). As genotypic discordance between the two calling workflows was low (1.1%) and levels of missing data were moderate (21.3%–25.4%), we suggest that workflow choice is primarily a matter of preference. Regardless, we recommend using consistent variant calling protocols when working with both GT‐seq and ddRADseq data.

GT‐seq genotyping success was high across all sample types. In particular, field‐collected scat samples with qPCR‐detectable polar bear DNA were genotyped successfully (using our criterion of >50% SNPs genotyped) for 62.9% of individuals and demonstrated an average of 14.9% missing data. Based on MS:CF comparisons and FF sample replicates not screened before sequencing, we estimated genotyping error for field scat to be around 10%. Genotyping success using GT‐seq was expected to be relatively low for field feces compared to tissue samples, as there is an unknown length of environmental exposure, and each sample contains a unique mix of DNA from nontarget species. Thus, to minimize effort and money expended on faecal samples unlikely to yield GT‐seq genotypes, we evaluated a qPCR screening assay (Hayward et al., 2020). The inclusion of a qPCR screening step substantially increased the percentage of field faecal samples successfully genotyped using GT‐seq from 30.6% to 62.9%, with 22 of the 23 (95.6%) field faecal samples that passed this step being successfully genotyped. In one case, duplicate subsamples from the same field scat exhibited a large difference in genotyping success (1 failure, 1 success), suggesting that samples can be heterogenous in the amount of host DNA. Thus, we recommend using repeated subsampling for scat samples for which genotyping is critical, and using a sample quality screening step such as our qPCR protocol for any future research making use of GT‐seq and noninvasive samples. Additional testing of DNA extraction protocols for scat may also help to improve DNA recovery.

Assessment of population structure

By combining our new GT‐seq data with ddRADseq data (Jensen et al., 2020), we were able to reassess Canadian polar bear population structure using a larger data set (642 vs. 358 individuals from Jensen et al., 2020) and make use of samples that were too degraded for genotyping with ddRADseq (e.g., faecal, degraded biopsy). Diversity metrics estimated using these combined data, including new genotypes from noninvasive samples, do not vary markedly among subpopulations and resemble those of previous studies based on SNPs, mtDNA, or microsatellites (e.g. Jensen et al., 2020; Malenfant et al., 2016), although our estimates may be higher due to selection of high diversity SNPs during panel design. Despite many subpopulations having only minor differences in genetic diversity, high assignment rates for a few subpopulations (e.g., GB, SH, MC) may have implications for conservation and management planning.

Sample clustering patterns and geographic distributions based on combined GT‐seq and ddRADseq data were also largely consistent with genetic groups called the “Hudson Complex,” “Arctic Archipelago,” and “Polar Basin,” which have been previously described using microsatellite and SNP data (Jensen et al., 2020; Malenfant et al., 2015, 2016; Paetkau et al., 1999). However, a new genetic cluster coincident with the M’Clintock Channel subpopulation emerged in both our structure and dapc analyses (Figures 1 and 3), a pattern not evident from analysis of ddRADseq data alone (Jensen et al., 2020). Rather than M’Clintock Channel comprising a major genetic cluster, it seems more likely that this new grouping reflects subtle genetic differentiation within the Arctic Archipelago (see Malenfant et al., 2016, 2020).

Considering the mean sample collection year for M’Clintock Channel of 1999 and the rapidness of environmental changes in the Canadian Arctic, our data may not represent contemporary genetic structuring of the Canadian polar bear population. Further, our data set contains samples collected between 1998 and 2018 (a 20‐year span), which presents a potential temporal confound that is not uncommon in studies investigating polar bear population structure (e.g. Malenfant et al., 2016; Paetkau et al., 1995, 1999). Regardless, we show here that our GT‐seq panel provides sufficient power to adequately capture known patterns that will be relevant to northern and federal governments for polar bear management. Importantly, we can detect these patterns through use of noninvasive and degraded samples, and GT‐seq genotyping can be done with greater cost‐efficiency than other genotyping methods irrespective of the number of samples available to be assayed at a time (Campbell et al., 2015). Thus, GT‐seq presents an opportunity to iteratively assess the genetic structure of Canada's polar bear populations – a flexible, cost‐effective means to continuously update our understanding as regular sampling is performed and governments shift towards noninvasive methods of population monitoring.

Other applications of GT‐seq and considerations

Further considerations during GT‐seq panel design and optimization may be required for other species and applications. For example, certain characteristics of a species may influence the number of markers required to discern individuals and populations, including mating system and linkage patterns within a genome. In the future, for polar bears we may wish to optimize a panel for assessment of adaptation in real time by including SNPs that are flagged as potentially under selection, or perhaps a panel to distinguish between bear species (e.g., grizzly and polar bears), similar to the application of GT‐seq to invasive brown (Rattus norvegicus) and black rats (R. rattus; Sjodin et al., 2020).

Although we intended to include at least 30 individuals from each Canadian subpopulation in the ddRADseq data set (Jensen et al., 2020) that was used to develop our GT‐seq panel, some subpopulations had few samples available. Ascertainment bias may thus be of concern as our GT‐seq panel may not fully represent all Canadian subpopulations, nor was it designed with the inclusion of data from other Arctic nations. Thus, application to the global polar bear population may require some panel modifications. Establishment of baseline global polar bear population structure through GT‐seq of noninvasive samples will be especially important for examining trends in dispersal patterns, diversity and structure, and behaviour in the context of a rapidly changing climate.

CONCLUSIONS

The need for monitoring using noninvasive genetic methods in the midst of rapid environmental changes, as well as the desire to integrate Indigenous ways of knowing and western science, is not unique to polar bear monitoring. Here, we respond to a long‐expressed desire by northern communities for a means to monitor polar bears noninvasively. We demonstrate that we can use noninvasively collected samples and GT‐seq to distinguish among individuals and to quantify population structure at levels comparable to other methods of genetic monitoring (e.g. Malenfant et al., 2016; Paetkau et al., 1995, 1999), but with greater efficiency even at small batch sizes (Campbell et al., 2015). The GT‐seq assay that we developed here provides a foundation for new community‐based programmes that can use noninvasive methods to improve temporal monitoring of polar bear populations and directly inform conservation efforts and government policy. The envisioned programmes will incorporate the perspectives of Indigenous communities throughout the planning and monitoring processes, and provide both social and economic benefits to them. This GT‐seq protocol is also intended to serve as the basis for a comprehensive toolkit to assess important aspects of polar bear and ecosystem health (e.g., contaminants, parasite load, diet). With GT‐seq at the core of the toolkit, a suite of data can be provided to communities and territorial governments annually. As this framework and our GT‐seq protocol can be adapted to other species and for other research questions, the monitoring methodology we propose here can be adapted and applied as a model for inclusive wildlife monitoring worldwide.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

Kristen M. Hayward participated in study design, conducted laboratory work, performed bioinformatic processing and analysis, and drafted the manuscript. Evelyn L. Jensen participated in study design, genomic data collection, bioinformatic processing, and analysis. Rute B.G. Clemente‐Carvalho participated in genomic data collection and was crucial for panel optimization. Marsha Branigan, Markus Dyck, and Peter V.C. de Groot managed access to samples and provided key insights to polar bear ecology and management implications. Zhengxin Sun performed laboratory work and Christina Tschritter aided with analysis. Stephen C. Lougheed conceived and coordinated study design and helped draft the initial manuscript. All authors contributed to manuscript revisions.

Supplementary Material

Click here for additional data file.