Tissue Sample Preparation and Loci Development

Small quantities of blood were collected by venal puncture at the Phoenix Zoo and transferred onto Gensaver 2.0 cards (GenTegra, Pleasanton, California). Blood cards were stored in envelopes at room temperature with silica beads until extraction at the Western Ecological Research Center San Diego Field Station genetics laboratory. All samples were extracted with the Gentra Puregene Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol, with minor modifications, including cell lysis in the presence of Proteinase K and dithiothreitol, an overnight deoxyribonucleic acid (DNA) precipitation step, and 5–30 minutes centrifugation steps at 21,194 relative centrifugal force and 4 degrees Celsius (°C). We quantified DNA concentration on a Qubit fluorometer (Life Technologies, Carlsbad, Calif.).

Microsatellites used for this study were developed by Wood and others (2022) and more information can be found in the U.S. Geological Survey (USGS) data release (Wood and others, 2022). We provide some details in the following sections. For initial marker discovery, a tetrameric enriched genomic library was created and sequenced on an Illumina MiSeq with paired 250 base-pair (bp) reads (Nali and others, 2014). SeqMan NGen (version 11, DNASTAR, Inc., Madison, Wisconsin) was used to generate a de novo assembly from the paired fastq files (raw data). Next, the fasta files generated from the previous step were used to scan for microsatellite loci and primer pairs with the program msatcommander 1.0.8_beta (Rozen and Skaletsky, 2000; Faircloth, 2008). To develop the microsatellite multiplex panels, we selected approximately 500 unique loci, appended the Nextera tags to the 5-prime end of the forward and reverse primers (5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG was appended to each forward primer, 5’-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG to each reverse primer), and generated lists of primer pairs forming primer-dimers using a Multiple Primer Analyzer (MPA; https://www.thermofisher.com/us/en/home/brands/thermo-scientific/molecular-biology/molecular-biology-learning-center/molecular-biology-resource-librar y/thermo-scientific-web-tools/multiple-primer-analyzer.html). We used the following process to ensure that the designed multiplexes would have minimal interference from primer-dimer formation. First, we ran a custom R script to cluster the loci into five groups that would minimize primer-dimer formation (using primer-dimer lists generated by MPA). Next, we selected 32 loci per group that had the lowest number of primer-dimers, and we tested each locus for successful amplification using 2 DNA samples. Finally, we selected 25 successfully amplifying loci from each of the 5 groups (herein referred to as multiplex) for a total of 125 loci genotyped for T. rufipunctatus (app. 1; Wood and others, 2022).

Microsatellite loci were genotyped using a multiplex polymerase chain reaction (PCR) protocol for targeted amplicon sequencing, which is briefly outlined in the next section (see D’Aloia and others [2017] and Wood and others [2022] for more details on the multiplex PCR and Nextera barcoding methods). We amplified the set of 125 loci in 5 multiplex PCR reactions across all samples, 25 loci per mix, with QIAGEN Multiplex PCR kits. Multiplexes were pooled by sample, and the samples were barcoded using Illumina’s S5 and N7 Nextera primers (Illumina, Inc., San Diego, Calif.) for the barcoding PCR. Before sequencing, all barcoded individual libraries were pooled into a single library and cleaned with Ampure XP beads (Beckman Coulter, Brea, Calif.). Finally, each library was diluted to a 2-nanomolar (nM) concentration and sequenced on an Illumina MiSeq with paired 150 bp reads.

Loci Verification

To evaluate the microsatellite loci, we used a subset of genetic samples from a previous study (Wood and others, 2018) that were representative of wild sites occupied by T. rufipunctatus throughout Arizona and New Mexico. We selected 142 individuals from 8 sites throughout the Verde River, Salt River, San Francisco River, and Gila River Basins.

We used a custom Python script called “amplicon” (https://bitbucket.org/cornell_bioinformatics/amplicon) to extract reads from the Illumina runs, assign them to the appropriate locus and individual, and identify the top two haplotypes for individuals at all loci (in other words, their diploid genotypes). We used the default options except for the following parameters: -c 1 (minimum number of samples per haplotype), -a 0.001 (minimum minor allele frequency), -l 75 (minimum haplotype length), and -r 5 (maximum read count ratio between the two alleles in each sample). We also applied an additional missing data filter after running the script to remove loci with greater than 20-percent missing data and individuals with greater than 10-percent missing data. We used the R packages adegenet v.2.1.3 (Jombart and Ahmed, 2011) and poppr v.2.9.1 (Kamvar and others, 2014) to estimate the number of alleles, expected heterozygosity, and to determine if any loci were out of Hardy-Weinberg equilibrium (HWE). We used the method of Brookfield (1996) to estimate the frequency of null alleles for each locus with the R package popgenreport (Adamack and Gruber, 2014), and we removed any locus with an estimated null allele frequency greater than 0.3, following the recommendations of Dakin and Avise (2004). We also tested for linkage disequilibrium among loci using the R package genepop v.1.2.2 (Rousset, 2008) with 5,000 iterations.

Paternity and Sibship Analysis

We assigned paternity and relatedness among eight litters composed of multiple female and male snakes born between 2014 and 2023 at the Phoenix Zoo breeding facility. We also completed a sibship analysis for six wild-caught gartersnakes from Canyon Creek that were brought into the Phoenix Zoo breeding facility in 2017 and 2018. We describe all snakes by their identification numbers that were individually assigned by the Phoenix Zoo for tracking purposes.

We used the program “Colony” (Wang, 2004; Jones and Wang, 2010) to estimate paternity assignments and relatedness among snakes. This program estimates the likelihood that offspring are related to each other, whether full siblings or half siblings, and estimates the likelihood that any genotyped adult male included in the analysis is the father. A total of eight breeding groups were analyzed separately. Each analysis included genotypes of known mother-offspring relationships and genotypes for as many as four adult males that were strategically grouped with a female in the same enclosure to allow free engagement in courtship and breeding behavior. Using the allele frequencies and known mother-offspring relationships, “Colony” generates expected genotypes of the males contributing to each offspring in a litter with a 95-percent confidence level. “Colony” settings were dioecious and diploid species, no inbreeding, polygamous mating system for males and females, full sibship size scaling, weak maternal and paternal likelihood priors that parents are included in the dataset, and five runs each with the full-likelihood analysis method with high precision. We also completed analyses using varying genotyping error rates (0–2 percent) to assess the accuracy of our results.

We estimated pairwise relatedness for all snakes (past and present) in the Phoenix Zoo using the package “related” v1.0 (Pew and others, 2015) accessed through R version 4.2.2 (R Core Team, 2022). For pairwise estimates of relatedness, we used two estimators. The first was the pairwise estimate of relatedness (R_qg) from Queller and Goodnight (1989), which was devised to estimate the average degree of relatedness within groups. We also used the pairwise estimate of relatedness (R_w) from Wang (2002) because this estimator was reported to handle better the uncertainty in allele frequency estimates than the R_qg estimator. We computed 95-percent confidence intervals for each of these point estimates. We used the simdata function in “related” to compare estimates of relatedness for R_qg and R_w using 100 simulated individuals based on our data. Using the simulated individuals, we computed estimates for each relatedness type (fullsib [FS], halfsib [HS], parent-offspring, unrelated) and calculated the correlation coefficient between the observed and expected values for each estimator.

Loci Verification and Relatedness Comparisons

Using the wild dataset for loci evaluation, we identified and removed 37 loci with greater than 20-percent missing data and 1 additional locus that was monomorphic. We also identified 4 loci that indicated significant deviation from Hardy-Weinberg equilibrium, 10 loci with null allele frequencies greater than 0.3, and 3 loci with significant linkage disequilibrium. These 55 loci were removed before further analysis with the captive dataset. Locus evaluation and filtering steps resulted in a total of 70 loci used for sibship, relatedness, and individual heterozygosity estimation (app. 1).

Comparisons between the relatedness estimators (R_qg and R_w) using 100 simulated individuals for each relatedness type (fullsib, halfsib, parent-offspring, unrelated) were similar (fig. 1), and the correlation coefficient between the expected and observed values for each estimator were high (R_qg=0.935 and R_w=0.938). For all comparisons of relatedness reported, we only report R_w (Wang, 2002). We provide pairwise relatedness for all breeders (past and present) in the Phoenix Zoo in appendix 2.

Relatedness of Wild Caught Snakes from Canyon Creek

“Colony” inferred fullsib and halfsib relationships among the snakes collected from Canyon Creek between 2016 and 2018. One pair of snakes has an inferred full-sibling relationship with greater than 0.95 probability and pairwise estimates of relatedness greater than 0.5 (table 1). This full-sibling relationship consisted of a male (13547) and a female (13548) that were collected from Canyon Creek in 2017. One pair of snakes has an inferred half-sibling relationship with greater than 0.95 probability and pairwise estimates of relatedness greater than 0.28 (table 1). This half-sibling relationship consisted of two males (13943 and 13946) collected from Canyon Creek in 2018. The remaining snakes collected from Canyon Creek indicated high relatedness estimates, and “Colony” results indicated half-sibling relationships, but the relationship could not be inferred with confidence (p<0.95; table 1).

Paternity and Relatedness

We completed paternity analyses of known mother-offspring relationships that resulted from cohabitation, with as many as four males in 2014, 2019, 2022, and 2023. In 2014, female 11345 was paired with two males (12303 and 10809). “Colony” inferred male 12303 as the sire of the six offspring genotyped, and pairwise R_w values between the dam and inferred sire were consistent with this inference (fig. 2). Pairwise R_w values between male 10809 and offspring were low, indicating that there was no relationship with the other male housed with dam 11345 and the offspring.

In 2019, the paternity of seven offspring born from a single female (12627) were inferred to be from two sires (13546 and 13547) with high probability (fig. 3). Because of higher-than-average relatedness (R_w=0.29) between the two sires, we ran additional paternity analyses where we removed one of the inferred sires from the dataset and reran “Colony” again, each with only one or the other sire present in the dataset to assess confidence in the number and probability of paternity inference. In both instances, two sires were still inferred: one sampled from the data and one inferred from the data but not sampled, and the probability of multiple paternity was still greater than 0.95. Two other males also were excluded from paternity (13943 and 10809).

In 2022, four females that were housed together with three males, each gave birth to eight or more offspring. For female 12641, “Colony” inferred eight of nine offspring to have arisen from two sires. One individual was not included in the analysis because of high missing data (14841), so the paternal contribution for this snake is unknown. The two inferred sires (13546 and 13943) indicated high relatedness (R_w=0.49) and were assigned as possible half siblings in the “Colony” analysis but with low confidence (table 1). To assess confidence in the multiple paternity inference, we removed one of the sires and reran “Colony,” each with only one or the other sire present in the dataset to assess if the number of sires and probability of paternity would change. In both instances, “Colony” inferred two sires, one sampled and one inferred from the data with high probability (>0.95). Figure 4 provides the inferred pedigree information and pairwise relatedness among the offspring and all adults in the breeding group.

For female 12627, all 12 of the offspring born had an inferred paternity to male 13546 (fig. 5). Although male 13943 was not inferred as a sire, related estimates between this male and the offspring born to dam 12627 were relatively high (mean=0.210), which was likely because of the high relatedness that was estimated between male 13943 and the inferred sire 13546 (R_w=0.49; table 1). For female 14477, all eight of the offspring born in 2022 had an inferred paternity to male 13946 (fig. 6). Similarly, the 11 offspring born to female 14479 in 2022 had an inferred paternity to male 13946 (fig. 7). Based on these analyses, a total of three sires were inferred across the four breeding groups from 2022, with one sire contributing to two breeding groups, and two sires contributing to multiple paternity of a single breeding group.

In 2023, two females (14477 and 14479) were paired with the same three males as in 2022. For female 14477, 9 of the 12 offspring born in 2023 had an inferred paternity to male 13943 (fig. 8). The remaining three offspring failed to sequence, so no inference was made. Similarly, the seven offspring born to female 14479 in 2023 had an inferred paternity to male 13943 (fig. 9).