A New Panel of SNP Markers for the Individual Identification of North American Pumas

Pumas Puma concolor are the most widespread terrestrial carnivore in the western hemisphere, ranging from Alaska and northern Canada to the southernmost extent of Argentina and Chile (Sunquist and Sunquist 2002). Regarded as one of the most iconic and well-studied species in the Americas, pumas often attract large amounts of attention as a game species (Logan and Sweanor 2001), for predating upon (at times endangered) wildlife (Turner et al. 1992; Wehausen 1996; Hayes et al. 2000; Schaefer et al. 2000), and for conflicts with humans (Beier 1991; Torres et al. 1996). They are frequently the center of conservation programs targeting entire landscapes (Weber and Rabinowitz 1996) because the removal of pumas, like other large carnivores, can have substantial impacts on the dynamics of both plant and animal communities (Logan and Irwin 1985; Turner et al. 1992; Schmitz et al. 2000; Ripple et al. 2014). In response to the significant ecological impacts of pumas, many wildlife agencies are actively researching and monitoring puma populations to make informed decisions regarding their conservation and management (Culver and Schwartz 2011).

As part of effective wildlife conservation and management programs, genetic studies are often initiated to investigate many aspects of a species' ecology (Frankham 2003). In pumas, genetic studies have concentrated either on the examination of gene flow on continental (Culver et al. 2000) and local scales (e.g., Loxterman [2011]; Andreasen et al. [2012]; Balkenhol et al. [2014]), or on the identification of individuals (Haag et al. 2009; Miotto et al. 2011; Naidu et al. 2011). Microsatellites have been the genetic marker of choice, and numerous loci derived from either domestic cats Felis catus (Menotti-Raymond and O'Brien 1995; Menotti-Raymond et al. 1997, 1999) or directly from pumas (Kurushima et al. 2006; Rodzen et al. 2007) have been developed and used in many studies (Figure 1). However, microsatellite allele sizes can vary among studies as a result of differences in the Taq polymerase, fluorescent dye, genotyping instrument, and analytical software used (Hahn et al. 2001; Vignal et al. 2002). This size variation renders comparing microsatellite genotypes across laboratories challenging and depositing genotypes into databases difficult without calibration of positive control samples. Furthermore, the loci analyzed to date vary across studies, often with little overlap. For example, across 17 genetic studies of western North American pumas published since 2000, 56 total different microsatellite loci were used (Figure 1). Many loci were used once (n = 25), whereas 9 loci were used in at least half the studies and no locus was shared across all studies. As a result of the technical difficulties in comparing microsatellite genotypes and little overlap of loci, genetic data from pumas have often been obtained from smaller geographic regions and rarely integrated into studies across larger landscapes (see references in Figure 1). Because pumas are capable of traversing long distances (>1,000 km; e.g., Thompson and Jenks 2005; Stoner et al. 2008), a consistent panel of genetic markers that can be analyzed in a high-throughput fashion and compared across laboratories without the need for calibrating genotypes will facilitate the long-term and range-wide genetic monitoring of North American puma populations.

Recently, another genetic marker, called a single nucleotide polymorphism (SNP), has become increasingly popular in population genetic studies. A SNP is a single base in the DNA sequence that contains multiple variants (alleles) segregating in the population. There are several features of SNPs that make them effective markers for population genetic studies (see reviews by Vignal et al. 2002; Morin et al. 2004; Garvin et al. 2010). Most importantly, SNPs are numerous, broadly distributed across various genomic regions (coding, noncoding, organellar), and amenable to high-throughput technologies (Garvin et al. 2010). Alleles are encoded simply as bases (i.e., A, C, T, G) rather than amplicon lengths, and are thus independent of the genotyping technology used. Because SNP loci are generally biallelic and have a well-defined mutation model (Garvin et al. 2010), the mathematical descriptions of many parameters using them are simpler than multiallelic loci (e.g., microsatellites). However, the biallelic nature of SNPs often reduces the information content per locus. As a result, more SNPs (2–6×) are theoretically required to match the resolution of microsatellites (Morin et al. 2004), although several empirical studies have demonstrated that <1–3× as many SNPs can perform equally as well (Seddon et al. 2005; Ryynänen et al. 2007; Coates et al. 2009; Fernández et al. 2013). This number can vary, however, depending on a variety of factors, including the parameter in question, per-locus information content, levels of genetic variation, absolute number of markers, and sample size. Furthermore, because SNPs can be genotyped using short fragments of DNA (often <100 base-pairs [bp]), it has been suggested that they may outperform microsatellites in amplification success from noninvasively collected samples (Morin et al. 2004; Seddon et al. 2005; Morin and McCarthy 2007; Ogden 2011).

To date, the use of SNPs in wildlife conservation and management programs remains relatively uncommon and has often been limited to taxa closely related to species with available genome sequences (e.g., Pertoldi et al. 2010; Cronin et al. 2014; Kraus et al. 2015). There currently exists no set of SNPs specific for use in pumas, but advances in next-generation sequencing technologies have made the process of SNP discovery much more tractable in nonmodel species (Garvin et al. 2010). Here we use Roche 454 pyrosequencing of cDNA prepared from whole blood to assemble the first transcriptome of pumas. We identify SNPs from the transcriptome assembly and use them to develop a high-throughput, 25-SNP genotyping assay for North American pumas Puma concolor couguar. We demonstrate the utility of our assay in 748 samples collected throughout their range and its ability to accurately identify individuals. Furthermore, we compare genotyping success of our SNP assay with that of 12 microsatellites, 8 of which are commonly used in studies of North American pumas, in 46 fecal DNA samples. Finally, we discuss the relevance of our assay as an effective tool for monitoring and managing individual pumas.

We collected whole blood samples from 12 pumas of wild origin via saphenous venipuncture. All pumas were from different locations in Arizona, with the exception of one puma from Bosque del Apache National Wildlife Refuge, New Mexico. Local agency personnel captured wild individual pumas under their live-capture permits. For each individual we filled two ribonucleic acid (RNA) Protect Animal Blood Tubes (Qiagen, Germantown, MD) with approximately 500 μL of whole blood each and immediately placed the tubes in liquid nitrogen, or on ice for <24 h and then in liquid nitrogen, until RNA extraction was performed. Remaining whole blood was stored at −20°C for DNA extraction. We extracted total RNA using the RNA Protect Animal Blood System (Qiagen) according to the manufacturer's recommendations, and combined the two extracted RNA samples per individual for the final sample.

We quantified total RNA in each sample using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and combined the 12 RNA samples in equimolar ratios to obtain a pooled RNA sample. We synthesized cDNA using the SMARTer Pico polymerase chain reaction (PCR) cDNA Synthesis Kit (Clontech Laboratories, Mountain View, CA) following the manufacturer's recommendations. This kit targets polyadenylated RNA and minimizes contamination from ribosomal and other small RNAs. To avoid biases in cDNA library construction and amplification, we prepared six independent cDNA libraries as described above, quantified each on a 2100 Bioanalyzer, and pooled them in equimolar concentrations for the final cDNA library. We sequenced the cDNA library on half of a picotiter plate on a GS-FLX Titanium platform (454 Life Sciences, Branford, CT) at the University of Arizona Genetics Core (http://uagc.arl.arizona.edu/).

We assembled all the sequences into contigs using GSASSEMBLER (454 Life Sciences), employing parameters to screen reads against NCBI's vector-contaminant database, to trim the reads of the SMARTer oligonucleotide primer used in cDNA library preparation, and to remove contigs <100 bp in length. The GSASSEMBLER software also quality-trimmed reads according to default criteria (3′ ends trimmed to approx. Q ≥ 20). We functionally annotated the contigs using FASTANNOTATOR (Chen et al. 2012), which integrates BLAST2GO (Conesa et al. 2005), PRIAM (Claudel-Renard et al. 2003), and RPSBLAST (http://blast.ncbi.nlm.nih.gov) to assign Gene Ontology (GO) terms, protein domains, and Enzyme Commission numbers to contigs.

Next, we cleaned the raw sequence reads by removing the cDNA library primer, trimming the 3′ ends of each read to a minimum phred-scaled base quality score of 20, and filtering reads shorter than 100 bp using CUTADAPT v1.0 (Martin 2011). We mapped the cleaned reads to the assembly using BWA v0.6.2 (Li and Durbin 2009). We identified and scored polymorphisms using FREEBAYES v0.9.4 (Garrison and Marth 2012) and filtered the resulting polymorphisms to include only biallelic SNPs with a phred-scaled quality score ≥20 and a minimum allele count ≥2. We performed a reciprocal BLAST (http://blast.ncbi.nlm.nih.gov) search to exclude potential false-positive SNPs as a result of the alignment of sequence paralogs. Briefly, we selected a fragment spanning approximately 100 bp on either side of each SNP and performed a nucleotide BLAST search against the RefSeq database. We omitted a SNP from further analyses if its fragment had multiple matches within the same genome with an e-value <1 × 10⁻³⁰ and similarity >90%. Finally, we removed all SNPs whose sequence had a best match to mitochondrial or ribosomal RNA genes.

In addition to the 12 samples described above, we collected 520 tissue samples from hunter-harvested pumas as part of a monitoring program by the Arizona Game and Fish Department (AZ; n = 532; Figure 2). We also used tissue samples previously collected by Culver et al. (2000), which represent a majority of the range of P. c. couguar (NAM; n = 216; Figure 2). We extracted DNA using a BioSprint 96 automated system (Qiagen) with Dynabeads magnetic bead technology (Invitrogen, Grand Island, NY) and quantified samples using the Quant-iT PicoGreen reagent (Invitrogen).

Using the candidate SNPs identified above, we designed four primer multiplexes for simultaneous analyses of 28, 29, 36, and 36 (n = 129 total) SNPs using the Assay Design v3.1 (Sequenom, San Diego, CA). This software designs PCR primers and a single, internal “extend primer” for each locus. Next, the software combines the loci into a multiplex that minimizes interactions between oligonucleotides and maximizes their amplification and genotyping potential following the instrument's requirements. We screened each multiplex for polymorphic SNPs in a subset of DNA samples (n = 191) on the MassARRAY system (Sequenom) at the University of Arizona Genetics Core. The 191 samples included the 12 pumas used to construct the cDNA library and 179 of the AZ pumas. We developed a final multiplex, PumaPlex 1.0 (hereafter referred to as PumaPlex), from the polymorphic SNPs using the Assay Design v3.1 as described above. Using PumaPlex, we genotyped the remaining samples (n = 748 total). All PCR reagents and cycling conditions for the PumaPlex runs can be found in Table S1 (Supplemental Material). We replicated 24 samples twice within and between the different Sequenom analyses and included negative controls for each assay. The replicate with the least number of called genotypes was removed from subsequent analyses. We annotated the SNPs on PumaPlex using a nucleotide BLAST search and assigned GO terms using the genomic resources available for the domestic cat (www.ensembl.org).

To further validate genotypes, we selected 8 SNP loci at random for verification using traditional PCR and Sanger sequencing in the 12 original puma samples. For each SNP we designed primers using PRIMER3 (Rozen and Skaletsky 2000) to target a fragment between 150 and 250 bp in length and to include ≥50 bp of flanking sequence on either side of the SNP. We amplified each SNP using approximately 20 ng of template DNA, 0.5 μM each of forward and reverse primer described above, 1× PCR buffer (USB, Cleveland, OH), 1.5 mM MgCl₂, 0.25 mM each dNTP, 0.05% BSA, and 0.5 U Taq DNA polymerase (USB) in a final volume of 20 μL. Cycling conditions consisted of an initial melting step of 5 min at 95°C, then 40 cycles of 95°C for 45 s, 52–62°C for 45 s, and 72°C for 1 min, with a final extension step at 72°C for 7 min. We purified all PCR products using the ExoSAP-IT PCR Clean-up Kit (USB) according to manufacturer's specifications. The PCR products were sequenced on an ABI3730 DNA Sequencer (Applied Biosystems, Foster City, CA) at the University of Arizona Genetics Core in both forward and reverse directions. We visually inspected the sequences for errors and assembled them using Sequencher v5.1 (Gene Codes, Ann Arbor, MI). We manually identified genotypes at each SNP locus after visual inspection of the chromatograms.

Using the results from PumaPlex, we removed individuals missing >20% of genotypes (n = 21; see Results). For each locus we calculated summary statistics in GENALEX v6.5 (Peakall and Smouse 2012), including allele frequencies, observed heterozygosity, and expected heterozygosity. Using only AZ pumas, we tested each locus for significant deviations from Hardy–Weinberg equilibrium with a chi-square test implemented in GENALEX and screened for the presence of null alleles following the procedure described by Girard (2011) using the median F_IS (calculated in GENALEX) and 100,000 simulations. This method is effective at detecting locus-specific deviations from neutrality caused by null alleles in both panmictic and inbred populations using any codominant marker. We calculated the observed, theoretical, and sibling probability of identity (P_ID) according to Waits et al. (2001) also in GENALEX.

We genotyped 46 DNA samples extracted from feces previously collected in southwestern AZ and identified as puma in origin (Naidu et al. 2011, 2014a), using both PumaPlex and 12 microsatellite loci: FCA026, FCA035, FCA043, FCA052, FCA057, FCA077, FCA090, FCA096, FCA144, FCA176, FCA221, FCA229 (Menotti-Raymond et al. 1999). The loci, eight of which are commonly used in studies of North American pumas, were selected specifically for their high polymorphism and previous use in pumas from the same region (McRae et al. 2005; Naidu et al. 2011; Figure 1). We followed PCR conditions for microsatellite loci as reported by Menotti-Raymond et al. (1999) and fragment analysis as described in Naidu et al. (2011). We replicated genotypes three times for both marker sets. We calculated the overall proportion of called genotypes as the number of called genotypes divided by the product of the total number of reactions performed (3 × 25 × 46 = 3,450 for SNPs and 3 × 12 × 46 = 1,656 for microsatellites). We then constructed consensus genotypes (25 × 46 = 1,150 for SNPs and 12 × 46 = 552 for microsatellites) and estimated the frequencies of false alleles (genotyping errors between replicates) and allelic dropout in GIMLET v1.3.3 (Valière 2002) using two degrees of stringency: 1) at least two replicates must produce a genotype at a given locus, and 2) all three replicates must produce a genotype at a given locus. We compared results between the SNP and microsatellite data sets using a Wilcoxon signed-rank test in R v3.0.0 (www.r-project.org).

In fulfillment of data archiving guidelines (Wenburg 2011; Whitlock 2011), we deposited raw sequence data in GenBank (Data A1, Archived Material, http://www.ncbi.nlm.nih.gov/sra/SRX633288) and made the transcriptome assembly, the raw SNPs identified, and all individual genotypes accessible in Data A2 (Archived Material, http://doi.pangaea.de/10.1594/PANGAEA.835154). We have also made the first SNP data from PumaPlex available as a spatially intuitive genetic database for pumas called the Puma Genetic Database (Naidu et al. 2014b), hosted on Environmental Systems Research Institute's ArcGIS Online through the University of Arizona (http://www.arcgis.com/home/item.html?id=4d9e04e504bb453691fbff736df49b3b). This online database links each individual sample with SNP genotypes and a variety of metadata including collection date, sample type, Global Positioning System location, etc., and is displayed as an interactive map. Users can filter for criteria of interest and subsequently download the respective information.

Our cDNA library sequencing produced 358,049 sequence reads with a mean length of 509.6 (SD = 135.2) bp (Table 1). We assembled the reads into 2,109 contigs with an N50 of 671 bp (N50 indicates 50% of all bases were found in contigs of this length or larger; Table 1). We assigned 12,231 GO terms to 1,244 contigs (59%) and were unable to assign a functional annotation to 865 (41%) of contigs (Table S2, Supplemental Material). The GO term assigned most often in the category of Biological Process was “oxygen transport” (n = 70), and the other commonly assigned terms included a variety of housekeeping and immune-related processes (Figure S1, Supplemental Material). Genes that localize to the “cytosol” or are “integral to the plasma membrane” were most often represented in the category of Cellular Component (Figure S2, Supplemental Material), and genes involved in “protein binding” were the most frequent Molecular Function (Figure S3, Supplemental Material). We ascribed 73 Enzyme Commission codes to 64 contigs and annotated 707 contigs with at least one protein domain (Table S2, Supplemental Material).

After cleaning and quality control, 268,744 (75.1%) sequence reads and approximately 92 million bp (50.5%) remained (Table 1). We were able to accurately map 34,675 (12.9%) of the cleaned reads to our transcriptome assembly. We identified 434 candidate SNPs that passed our quality criteria. We screened 129 of the candidate SNPs using four different Sequenom multiplexes and verified 30 (23.3%) as polymorphic. We were able to incorporate 25 of the verified SNPs into PumaPlex (the information necessary to construct PumaPlex using Sequenom or other methods can be found in Table S3 and Data S1, Supplemental Material). We annotated each of the 25 SNPs using either the domestic cat Felis catus or tiger Panthera tigris genomes, of which 16 SNPs were found in 3′ untranslated regions, 7 were coding but silent (synonymous) changes, and only 2 were missense (nonsynonymous) changes (Table 2). A functional description of each gene represented on PumaPlex and their associated GO terms is available in Table S4 (Supplemental Material).

We genotyped 772 total samples (748 pumas + 24 replicated samples) and our overall genotyping rate was 95.2%. Across 24 replicated samples, we found 100% agreement between called genotypes. We omitted 14 AZ and 7 NAM individuals missing genotypes at >20% of loci. Our final genotyping rate was 97.4% in AZ and 95.3% in NAM. In 12 individuals tested at 8 loci, we found 100% agreement between Sequenom- and Sanger-called genotypes. In AZ, minimum allele frequency ranged between 0.069 and 0.490, observed heterozygosity between 0.066 and 0.534, and expected heterozygosity between 0.128 and 0.500 (Table 3). Only three loci (PP05, PP08, and PP19) showed significant deviations from Hardy–Weinberg equilibrium after Bonferroni correction (P < 0.05; Table 3). Using the median F_IS of 0.061 we observed no evidence for the presence of null alleles (Figure 3). The theoretical P_ID in AZ and NAM was 7.7 × 10⁻⁹ and 1.6 × 10⁻⁷, respectively, and no genotype matches (observed P_ID) occurred when >17 or >22 loci, respectively, were included (Figure 4). Sibling P_ID was approximately one order of magnitude greater than the observed P_ID in both AZ and NAM at 7.1 × 10⁻⁵ and 3.6 × 10⁻⁴, respectively (Figure 4).

We genotyped 46 fecal DNA samples using both PumaPlex and 12 microsatellite loci. The overall proportion of called genotypes was significantly greater for SNPs (59.8%) than microsatellites (39.9%; Wilcoxon signed-rank test, v = 802.5, P = 7.1 × 10⁻⁵; Figure 5a). Using a cutoff of at least two called genotypes per locus, SNPs produced significantly more consensus genotypes (Wilcoxon signed-rank test, v = 709, P = 3.1 × 10⁻⁴; Figure 5a) and fewer false alleles (Wilcoxon signed-rank test, v = 64, P = 0.025; Figure 5b) than microsatellites. The frequency of allelic dropout was also less but not statistically significant (Wilcoxon signed-rank test, v = 29.5, P = 0.123; Figure 5b). After applying a more stringent cutoff of at least three genotypes called per locus, the proportion of genotypes successfully called was approximately 2.5 times greater for SNPs than microsatellites (Wilcoxon signed-rank test, v = 597.5, P = 3.9 × 10⁻⁶; Figure 5a). The frequency of false alleles did not differ (Wilcoxon signed-rank test, v = 55.5, P = 0.88; Figure 5b) between SNPs and microsatellites, and we observed no allelic dropout in either data set.

In this study we sequenced and assembled the transcriptome of a pool of 12 pumas in order to identify candidate SNPs for use in population genetic studies. Despite the fact that mammalian peripheral blood is primarily composed of nonnucleated erythrocytes and considered to contain little transcript abundance, we were able to produce >2,000 contigs and annotate 59% of them. The most common functional annotations were consistent with those reported from human erythrocytes (Kabanova et al. 2009) and with well-known properties of blood such as oxygen transport. We also identified a variety of genes implicated in immune-related processes, suggesting that future comparison of transcriptome profiles between healthy and diseased individuals could improve our understanding of wildlife health and diseases (Burczynski and Dorner 2006). Furthermore, the use of our sequence reads and transcriptome assembly can greatly improve the identification and annotation of genome sequences from pumas or other closely related species (Wolfsberg and Landsman 1997; Kan et al. 2001; Denton et al. 2014).

The PumaPlex assay is currently set up to genotype 25 SNPs per sample in either a 384- or 96-sample format on the Sequenom MassARRAY system. This technology is available at many laboratories and genetic core service centers across the country, where the workflow, from sample to results, can be completed in as little as 1 d (www.sequenom.com). This type of automated workflow is much faster than traditional microsatellite analysis, saving time and costs for wildlife practitioners. PumaPlex also requires very small amounts of DNA (approx. 10 ng/sample), thus providing an efficient and economical use of samples containing limited DNA quantity such as in noninvasive and museum-collected samples. We demonstrated that the MassARRAY system generates high call rates and reproducible genotypes, consistent with that reported elsewhere for the technology (Pusch et al. 2002; Cronin et al. 2014). Our genotypes were also congruent with those verified using traditional Sanger sequencing. Unlike other high-throughput SNP panels designed for wildlife species with available genome sequences (Cronin et al. 2014; Kraus et al. 2015), there were no SNPs previously known to be polymorphic in pumas. The MassARRAY system provided an effective method to screen candidate SNPs identified from small amounts of next-generation sequencing data.

With PumaPlex we were able to accurately identify individuals. Both theoretical and observed P_ID surpassed the recommended threshold (0.01–0.0001) proposed by Waits et al. (2001) using a minimum of 11 and 14 loci, respectively. Sethi et al. (2014) recommended ≥32 SNPs with a mean minimum allele frequency of 0.2 to reduce errors in mark–recapture studies, albeit less is sufficient when allele dropout is low (<0.6) and the number of SNPs can adequately identify individuals. The authors also demonstrated methods to minimize mark–recapture errors through the use of repeated genotyping or error-tolerant matching protocols (e.g., Creel et al. 2003). With PumaPlex, we observed no allelic dropout and low P_ID, demonstrating that the assay is useful for studies that require individual identification, such as estimating and monitoring population sizes, home range sizes, paternity, kinship, and forensic applications. The P_ID we calculated includes all 25 loci, despite 3 loci (PP05, PP08, and PP19) deviating significantly from Hardy–Weinberg proportions. Even after removal of these three loci, the P_ID (9.5 × 10⁻⁸) still remained well below the recommended threshold. We reported the P_ID across all loci to estimate the minimum possible, and we recommend PumaPlex users exclude markers that significantly deviate from the expected proportions (e.g., using tests of Hardy–Weinberg equilibrium, null alleles, or selection) prior to estimating P_ID. All loci, however, will provide useful information when determining sibship or excluding parents.

It is also important to understand the potential effects of ascertainment bias, or the systematic deviations from theoretical expectations between the population used for SNP discovery and the population under study. In our case, we used a small SNP discovery panel (n = 12) from AZ, which may have biased loci to those with alleles at intermediate frequencies in the population (Morin et al. 2004). This is of relatively little concern when used to analyze AZ, but must be taken into account when examining other, more distantly related populations. This bias can account for the higher mean minimum allele frequency and smaller P_ID in AZ compared with NAM, in addition to population structure and other demographic factors. For this reason we refrained from reporting several statistics in the NAM population (i.e., Hardy–Weinberg equilibrium, test for null alleles) and provided the data solely to demonstrate polymorphism across the subspecies. We suggest that researchers who use PumaPlex are cautious when analyzing populations genetically distinct from AZ for parameters sensitive to the allele frequency spectrum, such as genetic variation, effective population size, and demographic history, or account for ascertainment bias using statistical methods described elsewhere (Wakeley et al. 2001; Nielsen and Signorovitch 2003). However, PumaPlex still provides many useful, individual-based applications less sensitive to this bias, including individual identification, paternity, and assignment tests (Morin et al. 2004).

Although we did not investigate population structure in this study, work by Morin et al. (2009) suggested that approximately 30 SNPs would only be able to detect moderate differentiation (F_ST = 0.01), whereas ≥80 SNPs would be necessary to characterize demographic independence (F_ST < 0.005). This argues that PumaPlex, in its current form, is better suited to monitoring and identification of individuals rather than detecting population structure, although strong structure (F_ST >> 0.01) is likely detectable. In 22/25 (88%) of loci, we observed a heterozygosity deficit, which can often be attributed to demographic effects such as inbreeding or population structure (Wahlund effect). Technical artifacts such as null alleles may also produce this effect, yet our results suggest that null alleles are not present (Figure 3). In Arizona, McRae et al. (2005) observed differentiation between pumas collected north and south of Interstate 40, implicating the highway as a possible barrier to gene flow. Because our Arizona data set (AZ) includes the entire state, population structure is a plausible explanation for the heterozygosity deficit.

Ongoing work in our laboratory is testing additional candidate SNPs to expand the number of loci on future versions of PumaPlex. Although a majority of the SNPs on PumaPlex are noncoding, it is generally considered that variants in coding or other regulatory elements (such as 3′and 5′untranslated regions) are more likely to cause phenotypic effects and be under selection. Only three loci differed from Hardy–Weinberg proportions; therefore, the assumption of neutrality cannot be excluded for the remaining loci. However, there are special situations where strong selection can have no effect on Hardy–Weinberg proportions (see Lewontin and Cockerham 1959), and, in contrast with overdominance and underdominance, directional selection has relatively little impact on genotype frequencies (Waples 2015). This assumption of neutrality does not differ from that of other markers, such as microsatellites, which also may be linked to selected loci or themselves be the targets of selection (reviewed by Haasl and Payseur 2013). Future versions of PumaPlex will target an increasing number of intergenic and synonymous SNPs in addition to those ascertained from all six subspecies of pumas to expand its potential applications while minimizing potential impacts of ascertainment biases and selection.

The ability to identify individuals has become popular with the expansion of population surveys using noninvasive sampling (Waits and Paetkau 2005). In pumas, noninvasive genetic studies have focused on the amplification of microsatellites from fecal DNA, although often encountering poor genotyping success (e.g., Ernest et al. 2002; Miotto et al. 2007; Naidu et al. 2011). Our genotyping rate of SNPs in tissue samples was much higher (95.2%) than from fecal samples (59.8%); this result is consistent with the harsh desert conditions in southwestern Arizona and northern Mexico, including high temperatures and extensive exposure to ultraviolet light that can cause substantial damage to DNA (Friedberg 2003; Hofreiter et al. 2015).

Several studies have proposed that SNP genotyping, which often utilizes short DNA fragments, may be more conducive than microsatellites to the conditions of DNA often encountered in noninvasively collected specimens (Morin et al. 2004; Seddon et al. 2005; Morin and McCarthy 2007; Ogden 2011). Our results support this hypothesis, because PumaPlex produced more overall and consensus genotypes with fewer false alleles than microsatellites. A similar result has also been reported from wolf Canis lupus scats (Fabbri et al. 2012); however, our study is the first to corroborate this using the MassARRAY system. It is possible that dissimilar PCR amplification and cycling conditions explain the difference between SNPs and microsatellites. However, this is unlikely because conditions for both marker sets have been optimized either in silico (SNPs) or through extensive empirical use in puma scats (microsatellites; Ernest et al. 2000, 2002, 2003; Miotto et al. 2007; Naidu et al. 2011). The optimized PCR conditions for the SNPs on PumaPlex are standardized across all MassArray genotyping systems, demonstrating the advantage of such an assay relative to the time and resources necessary to develop study-specific conditions often required for microsatellite analysis. Nevertheless, the reduction in genotype error is important because errors can substantially inflate population size estimates using mark–recapture methods (Waits and Leberg 2000). However, users must also take the information content of the loci into account, because microsatellites are more informative per locus than SNPs on account of the increased number of alleles (Morin et al. 2004).

We created the Puma Genetic Database to have this SNP data set stored in a secure location and make it available for download in multiple standardized formats (e.g., CSV, File Geodatabase, KML, Shapefile) offered through the “Perform Analysis” tools in ArcGIS Online. This service can facilitate the exchange of data among researchers and laboratories with authorized access to the Puma Genetic Database. This database is now live at (http://www.arcgis.com/home/item.html?id=4d9e04e504bb453691fbff736df49b3b) and can be viewed publicly, but interested persons should contact the corresponding author to request authorization to upload and download data from the Puma Genetic Database. Using the options available in ArcGIS Online, this database facilitates easy filtering and selection of individual genotypes based on a variety of criteria (e.g., any user-defined geographic region, sex, date, sample type, country, etc.). The Puma Genetic Database will be continually updated as more genotyping is performed using PumaPlex. We encourage future users of PumaPlex to share their data and to contact the corresponding author to obtain instructions on formatting and submitting data for the database. If researchers intend to design independent genotyping assays using other techniques, the contig sequences and SNP information necessary are provided in Data A2 (Archived Material, http://doi.pangaea.de/10.1594/PANGAEA.835154) and Data S1 (Supplemental Material). Genotypes from these alternative genotyping technologies can also be submitted to the Puma Genetic Database.

In this study we developed PumaPlex, the first set of SNPs specific for use in pumas, and demonstrated its effectiveness in genotyping DNA extracted from both tissue and fecal material. Furthermore, we created a database to host the genetic and other information collected from pumas to promote data-sharing among researchers. Both tools are valuable resources to facilitate range-wide conservation genetic assessments and management of North American pumas.

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any supplemental material. Queries should be directed to the corresponding author for the article.

Data S1. Sequences of the 25 contigs containing the single nucleotide polymorphisms (SNPs) on PumaPlex. Each SNP and its two alleles are indicated within brackets (e.g., “[A/G]”). These SNPs were developed from a transcriptome sequenced from a pooled cDNA library of 12 pumas Puma concolor couguar collected in Arizona and New Mexico, USA, between 2009 and 2010. The file is in FASTA format and can be viewed in any text editor.

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S1 (26 KB TXT).

Table S1. The concentration and volume of polymerase chain reaction (PCR) reagents and cycling conditions for each step of the PumaPlex genotyping using the MassARRAY system (Sequenom). These PCR reaction (rxn) conditions were used to genotype 25 single-nucleotide polymorphisms in all puma Puma concolor couguar tissue and scat samples collected in North America from 1983 to 2013. The first step (“30 Plex, Full Enzyme”) is a traditional PCR, the second step (“SAP Reaction”) is a clean-up protocol using shrimp alkaline phosphatase (SAP), and the third step (“iPLEX Extension 30 Plex, Full Enzyme/Term”) is a single-base extension PCR using the internal, extension primer. All reactions were carried out at the University of Arizona Genetics Core using the manufacturer's specifications and supplied reagents. The file is formatted for use in Microsoft Excel. Other abbreviations: HPLC = high-performance liquid chromatography,

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S2 (12 KB XLSX).

Table S2. The annotations for each contig of the transcriptome assembly produced by FASTANNOTATOR (Chen et al. 2012). The transcriptome was sequenced from a pooled cDNA library of 12 pumas Puma concolor couguar collected in Arizona and New Mexico, USA, between 2009 and 2010. The annotations include a summary of the best BLAST result to GenBank's nonredundant database (“Top BLAST Hit,” “Hit Length,” “E-value,” “Bit Score”), associated gene ontology (GO) terms (see www.geneontology.org), Enzyme Commission (EC) numbers (“Enzyme,” see http://enzyme.expasy.org/), and protein domains identified (“Domain,” see http://pfam.xfam.org/). The dashes indicate contigs unable to be annotated. The file is in a tab-delimited, text format.

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S3 (704 KB TXT).

Table S3. The primers and information necessary to build PumaPlex on the Sequenom MassARRAY system (see Sequenom MassARRAY Assay Design v3.1 user manual; table 3 p. 70). The oligonucleotides shown here were used to amplify the 25 loci on PumaPlex in all puma Puma concolor couguar tissue and scat samples collected in North America from 1983 to 2013. WELL = the well number assigned to the assay; multiplexed assays are assigned the same well number, TERM = termination mix, SNP_ID = name of the input sequence, 2nd–PCRP = secondary amplification primer (includes secondary tag), 1st-PCRP = primary amplification primer (includes primer tag), AMP_LEN = amplicon length (in bases; includes primer tags and maximum SNP sequence length), UP_CONF = uniplex amplification score; this score indicates how well the multiplexed amplicon design meets the design criteria, MP_CONF = multiplex amplification score; this score indicates how well the multiplexed amplicon design meets the design criteria, Tm(NN) = extend primer melting temperature, calculated by Nearest Neighbor method, PcGC = percent GC content of the extend primer, PWARN = assay design warning codes (D/d = primer dimer potential between primers of multiplexed assays; G/g = primer subsequence of contiguous G nucleotides; H/h = primer hairpin or self-dimer potential; S/s = false priming potential), UEP_DIR = direction of MassEXTEND® (F = Forward, R = Reverse), UEP_MASS = extend primer mass, UEP_SEQ = extend primer sequence, EXT1_CALL = name given to the analyte 1 mass peak in the mass spectrum, EXT1_MASS = mass of analyte 1, EXT1_SEQ = sequence of analyte 1, EXT2_CALL = name given to the analyte 2 mass peak in the mass spectrum, EXT2_MASS = mass of analyte 2, EXT2_SEQ = sequence of analyte 2. The file is formatted for use in Microsoft Excel.

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S4 (46 KB XLS).

Table S4. Description and Gene Ontology (GO) terms (see www.geneontology.org) assigned to the genes of the 25 single nucleotide polymorphisms (SNPs) on PumaPlex. These 25 SNPs were developed from a transcriptome sequenced from a pooled cDNA library of 12 pumas Puma concolor couguar collected in Arizona and New Mexico, USA, between 2009 and 2010. All data were retrieved from the genomic resources available for the domestic cat Felis catus (www.ensembl.org). The file is formatted for use in Microsoft Excel.

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S5 (30 KB XLSX).

Figure S1. The frequency of the 50 most common Biological Process Gene Ontology (GO) terms assigned to the transcriptome sequenced from a pooled cDNA library of 12 pumas Puma concolor couguar collected in Arizona and New Mexico, USA, between 2009 and 2010. A complete list of GO terms and detailed descriptions is maintained by the Gene Ontology Consortium (see www.geneontology.org).

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S6 (136 KB PDF).

Figure S2. The frequency of the 50 most common Cellular Component Gene Ontology (GO) terms assigned to the transcriptome sequenced from a pooled cDNA library of 12 pumas Puma concolor couguar collected in Arizona and New Mexico, USA, between 2009 and 2010. A complete list of GO terms and detailed descriptions is maintained by the Gene Ontology Consortium (see www.geneontology.org).

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S7 (125 KB PDF).

Figure S3. The frequency of the 50 most common Molecular Function Gene Ontology (GO) terms assigned to the transcriptome sequenced from a pooled cDNA library of 12 pumas Puma concolor couguar collected in Arizona and New Mexico, USA, between 2009 and 2010. A complete list of GO terms and detailed descriptions is maintained by the Gene Ontology Consortium (see www.geneontology.org).

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S8 (133 KB PDF).

Reference S1. Naidu A, Fitak R, Culver M. 2014a. Landscape genetics of mountain lions (Puma concolor) in southwestern Arizona. Final report to the Arizona Game and Fish Department Habitat Partnership Committee, Project number HPC-09-406, Tucson, Arizona.

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S9 (1723 KB PDF).

Reference S2. Naidu A, Fitak R, Culver M. 2014b. Data sharing for wildlife management: the puma genetic database. Final report to the Arizona Game and Fish Department Habitat Partnership Committee, Project number HPC-10-705, Tucson, Arizona.

Found at DOI: http://dx.doi.org/10.3996/112014-JFWM-080.S10 (1020 KB PDF).

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any archived material. Queries should be addressed to the corresponding author for the article.

To cite this archived material, please cite both the journal article (formatting found in the Abstract section of this article) and the following recommended format for the archived material.

Data A1. We deposited the raw sequence data into the Sequence Read Archive in GenBank under accession number SRX633288, titled “A transcriptome resource for pumas Puma concolor couguar” (available at: http://www.ncbi.nlm.nih.gov/sra/SRX633288). The transcriptome was sequenced from a pooled cDNA library of 12 pumas collected in Arizona and New Mexico, USA, between 2009 and 2010.

Data A2. Fitak, R et al. (2014): PumaPlex 1.0: Data generated during the development of the SNP markers and SNP genotypes of pumas. doi:10.1594/PANGAEA.835154.

We would like to thank those who provided samples and expertise, especially S. Bender, the Navajo Nation Zoo, L. Searles of the Southwest Wildlife Conservation Center, the late Andy Salazar, and the Arizona Game and Fish Department. We thank the NSF-IGERT program in comparative genomics and the Habitat Partnership Committee grants HPC 09-406 and 10-705 for funding. We are deeply indebted to C. Ramsower at the University of Arizona Genetics Core for assistance with Sequenom genotyping. We are grateful to T. Edwards, S. Blount, C. Voirin, E. Mohandesan, the Subject Editor, and three anonymous reviewers for their useful comments on earlier drafts of this manuscript.

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.