Population genetic structure of Texas horned lizards: implications for reintroduction and captive breeding

Sampling and DNA extraction

A number of volunteers collected 542 Texas horned lizard tissue samples across Texas, New Mexico, and Colorado between 2009 and 2017 (Fig. 1). While most of these samples were collected using the cloacal swab method described in Williams et al. (2012), some were from tissues collected from road kill and toe clips collected as part of other population studies. Field activities were approved by Texas Parks and Wildlife Department (SPR-1006-763). Texas Christian University Institutional Animal Care and Use Committee (IACUC) provided approval for this research (protocol 01/08). Samples were mapped onto EPA level II and III ecoregions using ArcGIS Pro and each sample was classified as belonging to a specific ecoregion (Omernick, 1987; Omernick, 1995; Omernik & Griffith, 2014) (https://www.epa.gov/eco-research/ecoregions).

We extracted DNA by incubating swabs and tail or toe clips overnight in 300 µl lysis buffer and 15 µl of Proteinase K (20 mg/ml) at 55 °C. The following day, a half volume of 7.5 M ammonium acetate was added to precipitate proteins. Then 0.7 volume of isopropanol was added to the supernatant and the samples were placed at −20 °C overnight to precipitate the DNA. Finally, DNA was pelleted and washed in 70% ethanol, air dried, and resuspended in 100 µl 10 mM Tris pH 8.5.

Genotyping and sequencing

We amplified 10 microsatellite loci in three multiplexes using 10 µl polymerase chain reactions (PCR; Williams et al., 2012). Four loci (PcGATA49, PcGATA61, PcGATA60, PcGATA31) reported in Williams et al. (2012) gave evidence of high levels of null alleles in multiple populations, or problems with large allele dropout (PcGATA 49) and so are not included in this study. PCR reactions contained 10–50 ng DNA, 0.2 µM of each primer, 1× Qiagen Multiplex PCR Master Mix with HotStarTaq, Multiplex PCR buffer with 3 mM MgCl₂ pH 8.7, and dNTPs. Reactions were cycled in an ABI 2720 thermal cycler (ThermoFisher Scientific, Waltham, MA, USA). The cycling parameters were one cycle at 95 °C for 15 min; followed by 35 cycles of 30 s at 94 °C, 90 s at 60 °C, 90 s at 72 °C; then a final extension at 60 °C for 30 min. Following amplification, all reactions were diluted with 200 µl dH₂O. For each sample, 0.5 µl of diluted product was loaded in 10 µl HIDI formamide with 0.1 µl LIZ-600 size standard (ThermoFisher Scientific, Waltham, Massachusetts, USA) and electrophoresed on an ABI 3130XL Genetic Analyzer (ThermoFisher Scientific, Waltham, MA, USA). Genotypes were scored and binned using Genemapper 5.0 (ThermoFisher Scientific, Waltham, Massachusetts, USA). We reamplified a subset of samples (∼10%, n = 55 individuals) to estimate the genotyping error rate.

We amplified a 353 bp fragment near the 5′ end of the mitochondrial control region using primers PcCR F—CTTATGATGGCGGGTTGCT and PcCR R—GGCTGTTAAATTTAT CCTCTGGTG for all 542 individuals. We also amplified the mitochondrial NADH dehydrogenase subunit 4 (ND4) and the tRNAs Histidine, Serine, and Leucine region using the primers ND4—ACCTATGACTACCAAAAGCTCATGTAGAAGC and Leu—CATTACTTTTACTTGGATTTGCACCA from Arèvalo, Davis & Sites (1994) in 49 individuals from across the sampled range. PCR reactions contained 10–50 ng DNA, 0.25 µM of each primer, 1× Qiagen Multiplex PCR Master Mix with HotStarTaq, Multiplex PCR buffer with 3 mM MgCl₂ pH 8.7, and dNTPs. Reactions were cycled in an ABI 2720 thermal cycler. The cycling parameters were one cycle at 95 °C for 15 min; followed by 35 cycles of 30 s at 94 °C, 15s at 55 °C, 30 s at 72 °C; then a final extension at 72 °C for 5 min. Reactions were cleaned enzymatically with ExoI and rSAP using the manufacturer’s protocols (New England Biolabs Ipswich, Massachusetts, USA). Products were sequenced in both directions using PCR primers and internal primer ND4#2-TACGACAAACAGACCTAAAATC from Arèvalo, Davis & Sites (1994) for ND4 and electrophoresed on an ABI 3130XL Genetic Analyzer. Sequences were trimmed, edited, and put into contigs using Sequencher 4.8 (Gene Codes Corporation, Ann Arbor, MI USA). Sequences were aligned in MEGA 10 (Kumar et al., 2018) using Clustal W (Thompson, Higgins & Gibson, 1994). All unique sequences have been deposited in GenBank (accession numbers: MK100594 –MK100710).

Microsatellite genotype analyses

The cryptic nature and low density of Texas horned lizards in most places made it difficult to find large numbers of individuals from single localities. To calculate traditional genetic diversity metrics, however, we grouped 449 of the 542 horned lizards into 16 sampling sites with ≥10 individuals in each (Table 1, Fig. 2). We grouped lizards into these sites because we were interested in estimating the genetic diversity present in protected areas such as Wildlife Management Areas and State Parks where individuals might be captured for future captive breeding purposes. In addition to these nine protected areas, we also chose seven other areas such as counties with ≥10 sampled lizards, to increase our number of sampling sites to test for isolation-by-distance. The coordinates for these sites were the centroids of the sampled lizard locations. For these 16 sites, we used GenAlEx 6.5 (Peakall & Smouse, 2006; Peakall & Smouse, 2012) to calculate observed (H_O) and expected heterozygosity (H_E) and F_IS. We used HP-RARE (Kalinowski, 2005) to calculate allelic richness (A_R) at each site to standardize comparisons across sample sizes. We tested for heterozygote excess and deficiencies and genotypic linkage disequilibrium using GENEPOP 4.0 (Rousset, 2008). We used sequential Bonferroni correction to determine significance for these tests. MICRO-CHECKER 2.2.3 (Van Oosterhout et al., 2004) was used to determine the presence of null alleles, large allele dropout, or issues with scoring due to stuttering and to calculate null allele frequencies (Brookfield, 1996 eq1). As there was some evidence of null alleles, we used the ENA correction method from Chapuis & Estoup (2007) implemented in the software FreeNA to calculate global and pairwise F_ST. We used these corrected values to test for isolation-by-distance using the Mantel test in GenAlEx. Because the magnitude of F_ST is influenced by heterozygosity, we also present the standardized measure F′_ST developed by Meirmans & Hedrick (2011). We also used a principal component analysis (PCA) in GenAlEx to visualize the pattern of pairwise F′_ST between sampling sites and compared these to the STRUCTURE results (see below) which utilized all samples.

We used STRUCTURE 2.3.4 (Pritchard, Stephens & Donnelly, 2000) to cluster multilocus genotypes from all 542 samples, including those not found within the 16 sampling locations mentioned above. We assumed admixture and correlated allele frequencies with no location prior and set the burn-in to 10⁴ iterations and ran the MCMC (Monte Carlo Markov Chain) for 10⁶ iterations. STRUCTURE can give misleading results both for the number of populations and individual ancestry if there is uneven sampling across clusters (K; Puechmaille, 2016; Wang, 2017). We used the recommendations of Wang (2017) and set the prior for admixture to allow α to vary between clusters and we decreased the initial α from 1.0 to 0.2. We ran ten independent runs for K = 1–10. The most likely K was identified using the method of Evanno, Regnaut & Goudet (2005) and by determining the K with the highest LnP(D) before values started to plateau (Pritchard, Stephens & Donnelly, 2000). We then used CLUMPP 1.1.1 (Jakobsson & Rosenberg, 2007) to average across the ten runs for the most likely K. We considered individuals to be admixed between clusters when their ancestry (q) was ≥ 0.10 in each of two or more clusters. We chose this value since similar cut-off values have been used in a number of other studies (Barilani et al., 2005; Vähä & Primmer, 2006; Sanz et al., 2009; Bohling, Adams & Waits, 2013; Johnson et al., 2015). The null alleles we found in this study are not expected to impact the accuracy of these assignments given the very small difference in F_ST values after correction (see ‘Results’) and because simulation studies have found only a slight reduction in the power of assignment tests in the presence of null alleles (Carlsson, 2008; Marsh et al., 2008).

Mitochondrial sequence analyses

We used GenAlEx to calculate the frequency of haplotypes and haplotype diversity (h) of the control region for the 16 sampling sites mentioned above. We used CONTRIB 1.02 (Petit, Mousadik & Pons, 1998) to calculate haplotype richness (H_R) at each site. We compared mtDNA population structure for these 16 sites to the nuclear microsatellite structure by calculating ϕ_PT for both data sets and conducting AMOVAs in GenAlEx.

We constructed single gene trees for both ND4 and the control region due to the difference in sampling for each marker. The model for nucleotide substitution was that with the highest value of the Bayesian information criterion (BIC) in MEGA. The resulting models were HKY + G for ND4 and HKY + G + I for the control region. These models were then used to create maximum likelihood trees in MEGA using the nearest-neighbor-interchange and a strong branch swap filter. Bootstrap values were calculated using 1,000 replicates. Phrynosoma asio (JN809342.1) was used to root the ND4 gene tree and Phrynosoma blainvillii (NC_036492.1) was used to root the control region tree. Bayesian analyses using MrBayes 3.2 (Huelsenbeck & Ronquist, 2001; Ronquist et al., 2012) were also conducted for the ND4 region. The analysis began with two parallel runs with random starting trees and were run for 2 × 10⁷ generations and sampled every 10³ generations. Burn-in was the first 25% of generations and a 50% majority rule consensus tree was calculated.

We evaluated whether well-supported Texas horned lizard clades (bootstrap ≥ 80) identified in the gene trees had experienced past population expansion or were stable, using all samples sequenced for the ND4 region and a subset of samples for the control region. We used DnaSP 6.11.01 (Rozas et al., 2017) to calculate the frequency of haplotypes, haplotype diversity (h), nucleotide diversity (π), and θ (Watterson’s estimator). We also conducted neutrality tests in DnaSP, including Tajima’s D (Tajima, 1989) and Fu’s F_S (Fu, 1997) to compare the number of rare and common mutations to the null hypothesis of a stable population, in contrast to an expectation of an excess of low frequency mutations following rapid population growth. When D and F_S were significantly different from neutral expectations for the ND4 data, we used the mismatch distribution to estimate the time since population growth, tau (τ) as τ = 2 µkt, where t is the time in generations, µis the mutation rate per site per generation, and k is the sequence length (Rogers, 1995). We calculated tau in Arlequin 3.5.1.3 (Excoffier & Lischer, 2010). We assumed a generation time of two years (Ballinger, 1974), and used an online calculator to calculate time since population growth (Schenekar & Weiss, 2011). We used a substitution rate of 0.00805 (substitutions/site/million years) which was estimated for ND4 in geckos (Macey et al., 1999) and used by Blair & Bryson (2017) in their study of species delimitation in Phrynosoma.

Microsatellite data

Most loci across the 16 sampling sites were in Hardy-Weinberg equilibrium and all loci were in genotypic linkage equilibrium. There was one locus/site comparison (out of 160) that had a significant heterozygote deficit after sequential Bonferroni correction. In addition to this one comparison, MICRO-CHECKER also identified five other loci/site comparisons that had null alleles (Table S1). There were four loci that gave evidence of null alleles at one site each and another locus that had evidence for null alleles at two sites (Table S1). There were a total of 19 individual/locus non-amplifications in the entire dataset, spread across six of the loci. We found eight allele errors across 1,140 alleles in the 55 individuals we re-genotyped for an error rate of 0.007. These errors occurred across five loci (per locus error rates for these five loci ranged from 0.009 to 0.018). A total of seven samples (12%) had one allele error each.

Genetic diversity was generally high (mean ± SE: H_E = 0.83 ± 0.006, A_R = 8.47 ± 0.18, n = 16 sampling sites) and similar across sites, with the exception of the population on Matagorda Island which had lower heterozygosity and allelic richness than mainland sites (Table 1). We corrected for null alleles when calculating F_ST; however, this made virtually no difference in the observed patterns (F_ST values generally only decreased by 0.001). Population differentiation was modest but significant with a global F_ST of 0.044 (95% CI [0.033–0.061]) and pairwise F_ST ranging from 0.003 to 0.115 (Tables S2 and S3). Standardized values were much higher with a global F′_ST of 0.302 and pairwise values ranging from 0.011 to 0.665. The first two axes of the PCA of pairwise F′_ST explained 52.4% of the variance and revealed a tight cluster of 11 sites (RPQRR, Midland Co., CMA, Yoakum Dunes WMA, Matador WMA, Seminole Canyon SP, Grey Co., Mitchell Co., Camp Bowie, Colorado, E. New Mexico) and two other more separated sets of sites (Chaparral WMA, Starr Co., Matagorda Island WMA) and (Brewster Co., Hueco Tanks SP; Fig. 3). Most pairwise F_ST comparisons were significantly different from zero except for 23 comparisons which were all between sites within the cluster of 11 sites in the PCA (Table S2). There was a significant pattern of isolation-by-distance for all 16 sampling sites (R² = 0.33, P = 0.003) and for the 11 sampling sites that formed the tight cluster in the PCA (R² = 0.24, P = 0.02).

The STRUCTURE analysis using all samples revealed three major genetic clusters or populations, both using the Evanno, Regnaut & Goudet (2005) method and the LnP(D) method of Pritchard, Stephens & Donnelly (2000; Fig. S1). At K = 3, samples were partitioned into west, north, and south populations (Fig. 4). At increasing levels of K (4–10), new clusters were simply added as admixture in the large north population. Sub-clustering the south population produced a split between Matagorda Island WMA and the mainland (Fig. S2). Sub-clustering the west population revealed two clusters, one composed of west individuals and another population from Brewster Co. around the Elephant Mountain WMA (Site #1 in Fig. 2, Fig. S2). Most of the individuals (83%; 25 of 30) assigned to the west sub-population had western clade haplotypes (see mitochondrial section below), whereas most of the individuals (95%; 21 of 22) assigned to the Elephant Mountain WMA sub-population had eastern clade mitochondrial haplotypes. Sub-clustering the north population did not reveal more populations. The north, south, and west populations found by the STRUCTURE analyses were concordant with the pattern of pairwise population F′_ST seen in the PCA. The 11 sites that clustered close together all belonged to the north population, whereas the Chaparral WMA, Starr Co. and Matagorda Island WMA belonged to the south population. Hueco Tanks SP was assigned to the west population and Brewster Co., which fell between the northern population and Hueco Tanks SP in the PCA, had individuals assigned to north and west populations and individuals that were admixed between the populations.

The west population is found within the Chihuahua Deserts and Madrean Archipelago ecoregions and the south population is found south of the Balcones Escarpment in the Southern Texas Plains, East-Central Texas Plains, and Western Gulf Coastal Plain ecoregions (Table 2, Fig. 1). The north population is found north and north-west of the Balcones Escarpment within a number of level III ecoregions including the High Plains, Central Great Plains, Southwestern Tablelands, Edwards Plateau, and Cross Timbers which collectively belong to the level II South-Central Semi-Arid Prairies ecoregion (Table 2, Fig. 1). Some individuals with genotypes from the north population (n = 28 individuals) also occurred in the Chihuahua Deserts ecoregion in Texas. There were two individuals that had high (q > 0.90) ancestry assignment in the north population but were found in the South Texas Plains, and one individual with high ancestry assignment in the south population that was found in the Central Great Plains ecoregion.

Most individuals (90%; 486 of 542) were strongly assigned to a single population. The remaining individuals (10%; 56 of 542) had evidence of admixture. Admixture between the south and north populations (n = 41 of 488 individuals) was more common than between west and north populations (n = 13 of 421 individuals; Fisher exact test, P = 0.001). The remaining two individuals were admixed between all three populations; one was found in southern Texas whereas the other was found in northern Texas. Admixed individuals with north and west ancestry were found in Brewster Co. where the north and west populations meet (n = 4 individuals) and in widely separated areas including Colorado (n = 4), Seminole Canyon SP (n = 3), Matador WMA (n = 1), Yoakum Dunes WMA (n = 1), CMA (n = 1), and Starr Co. (n = 1). Admixed individuals with north and south ancestry were found south of the Balcones Escarpment (n = 14) in the south population and far north of the Escarpment (n = 27) in the north population.

Mitochondrial data

Haplotype diversity at the control region was high with a total of 86 haplotypes found across all 542 individuals. Haplotype diversity (h) and richness was more variable between the 16 sites than microsatellite diversity (Tables 1 and 3). For instance, Midland, Yoakum Dunes WMA, and RPQRR had relatively low haplotype diversity and richness compared to their microsatellite diversities. Despite the variability in mitochondrial diversity there was a positive correlation between haplotype richness and allelic richness (r_S = 0.57, P = 0.02). Population subdivision was considerably higher for the mtDNA control region than for the microsatellite loci with 32% of the variance between sites for the control region (ϕ_PT = 0.32, AMOVA, P = 0.001) and 8% of the variance between sites for the microsatellite loci (ϕ_PT = 0.08, AMOVA, P = 0.001).

The control region gene tree revealed two clades between haplotypes found in western areas (New Mexico and far western Texas) and more eastern areas (Fig. 5). The western clade was well supported in this analysis (bs = 80). For the control region, there were a total of ten unique haplotypes in the western clade and 76 haplotypes in the eastern clade. There was an average of 3.91% divergence between the clades (range: 2.01–5.31%), 1.43% divergence within the eastern clade (range: 0.28–2.92%), and 0.08% divergence within the western clade (range: 0.28–1.72%). P. cornutum differed from P. blainvillii by 12.40% (range: 11.31–13.74%). The western haplotypes mainly belonged to individuals with high microsatellite ancestry in the west population (mean ± SE q: 0.97 ± 0.003, n = 28), with the exception of two individuals that had high north ancestry (q = 0.99) and were found in Brewster Co. Within Brewster Co., there were also 21 individuals with high west ancestry (q > 0.90) that had eastern haplotypes, indicating that this region is an area of admixture between the western and eastern mitochondrial clades. We found no evidence for separate mitochondrial clades corresponding to the north and south populations detected by STRUCTURE. Individuals with the eastern haplotypes had microsatellite ancestry within the north, south, and west populations. Control region haplotypes were, however, largely unique to the geographic regions encompassed by the three main nuclear microsatellite populations (west, north, and south) as indicated by the colored symbols in Fig. 5. Only seven of the 86 haplotypes were shared between regions. The south region had the highest haplotype diversity (h = 0.90), followed by the west (h = 0.88), and then the north region (h = 0.72). To compare the three regions, we used the standardized estimate of ϕPT using the method of Hedrick to ensure that ϕPT = 1.0 when populations have non-overlapping sets of haplotypes (Meirmans & Hedrick, 2011). ϕPT was 0.97 (P = 0.001) among the three regions, reflecting the low sharing of haplotypes between them.

The ND4 gene tree also revealed two well-supported clades between haplotypes found in western areas (New Mexico and far western Texas) and more eastern areas (Fig. 6). Both the maximum likelihood and Bayesian analyses recovered the same tree topology. There were a total of ten unique haplotypes in the western clade and 21 unique haplotypes in the eastern clade. There was an average of 7.07% divergence between the clades (range: 6.50–7.80%), 0.07% divergence within the eastern clade (range: 0.10–1.70%), and 0.04% divergence within the western clade (range: 0.10–0.60%). P. cornutum differed from P. asio by 16.40% (range: 15.80–16.80%). Similar to the control region, the western clade comprised individuals with high microsatellite ancestry in the west population. Similar to the control region, we did not find separate clades for the ND4 gene that corresponded to the north and south populations detected by STRUCTURE.

Using all samples sequenced for the ND4 region and a subset of samples sequenced for the control region (same samples as the ND4 region plus all western clade samples), haplotype diversity and sequence diversity were higher for the eastern clade than the western clade (Table 4). This result is consistent with the larger geographic range of the eastern clade and the higher effective population size of the eastern clade as indicated by larger theta values (Table 4). We found evidence for past population expansions for both clades; Fu’s F_S was significantly negative (P < 0.05) relative to neutral expectations for both western and eastern clades of ND4 and the control region. Tajima’s D was significantly negative (P < 0.05) only for the ND4 data (Table 4). Using the data from the ND4 region, the time of expansion was estimated to be 151,288 years ago for the western clade and 467,196 years ago for the eastern clade.