We describe the only known tetraploid parthenogenetic species of amniote that has haploid genomes from four distinct ancestral bisexual (gonochoristic) species. These genomes were brought together through three hybridization events that occurred over a time span of hundreds or thousands of years, the last of which occurred recently in captivity without any experimental manipulation. Inheritance of alleles at seven microsatellite deoxyribonucleic acid loci and DNA sequence data for adenosine deaminase intron 9 in the new tetraploid species confirms its parentage and, together with DNA quantification, demonstrates that tetraploidy and high heterozygosity are maintained parthenogenetically generation after generation. Comparisons of univariate and multivariate variation in scalation between the tetraploids and their parental taxa reveal a strong maternal similarity in morphological characters, with little to no significant differences, suggesting that such tetraploid females may exist in old museum samples misidentified as the maternal ancestor. This research on specimens of known parentage adds support to earlier papers that provide important insights for taxonomic treatment of many parthenogenetic clones of teiid lizards that occur in the Western Hemisphere.
INTRODUCTION
We have previously discussed the importance of laboratory colonies for understanding important aspects of the biology of unisexual species of whiptail lizards that occur in nature (e.g., Cole et al., 2014, 2017), including the following: 1) interspecific hybridization of gonochoristic (bisexual) species often results in F1 hybrids; 2) some F1 hybrid females produce all-female lineages by parthenogenetic cloning, dispensing with the need for spermatozoa in just one generation; 3) although diploid parthenogens produce diploid clones, they can also hybridize with a male of another species, which results in triploid offspring; 4) some of these triploids produce triploid all-female clones, although many clones (diploid and triploid) develop genetic modifications over time (e.g., Dessauer and Cole, 1989; Cole et al., 2019); 5) triploid parthenogens can both produce triploid clones and hybridize with a male of another species, resulting in tetraploid offspring; and 6) some of the tetraploids produce tetraploid all-female clones.
The varieties of diploid, polyploid, unisexual, bisexual, parthenogenetic, and clonal species mentioned above are being used for continuing basic research in the molecular and cellular mechanisms involved in their biological systems (e.g., Lutes et al., 2010, 2011; Newton et al., 2016). Promising areas of research include processes of oogenesis, gene regulation, the switch from reproduction based on fertilization to cloning of unfertilized eggs in one generation, and the molecular distinction between successful parthenogenesis in hybrids versus hybrid sterility, among other things. This research involves development of extensive sequence data that are being deposited with online data information services (e.g., GenBank, RefSeq, and the Sequence Read Archive [SRA]) and the need for scientific names to be applied to unique clonal unisexual species to provide clear and unambiguous communication in the literature. The need for such clarity is the foundation for having scientific names recognized in the international literature (ICZN, 1999). Although it is rare to name a taxon that resulted from hybridization in captivity by normal means without any experimental manipulation, these organisms fit all the requirements for recognizing species. Indeed, scientific names are also applied to domesticated animals, having an extensive literature in the veterinary sciences (e.g., donkeys, Equus asinus; see Todd et al., 2022), and to triploid parthenogenetic organisms for which it is not known whether they originated in the field or captivity (e.g., crayfish, Procambarus virginalis; see Vogt et al., 2018; Vogt, 2021).
Here we describe a unique unisexual species of hybrid origin that clones itself parthenogenetically and is reproductively isolated from all other species. Tetraploid all-female clonal lineages (up to the F7 generation to date) were produced after hybridization of females of Aspidoscelis exsanguis (a triploid parthenogen) with males of Aspidoscelis gularis (a diploid gonochoristic species; Fig. 1). We provide molecular evidence that the tetraploids are clonal parthenogens, and we compare univariate and multivariate morphological variation among the tetraploids and their parental taxa. We also show that the extent of morphological variation in the clonal species can be equivalent to that of the gonochoristic A. gularis.
Figure 1.
Representatives of two species of Aspidoscelis and their F1 hybrid. Left: Adult male of A. gularis, SIMR 8868. Right: Adult female of A. exsanguis, MCZ R-198982 (SIMR 9049), snout–vent length approximately 80 mm. Middle: F1 hybrid female, MCZ R-199029 (SIMR 10140), an offspring of the other two individuals shown.
MATERIALS AND METHODS
Species, Populations Sampled, Specimens, and Museum Abbreviations for Specimens Examined
We examined specimens of the triploid parthenogen A. exsanguis (Chihuahuan spotted whiptail lizard); diploid gonochoristic A. gularis (common spotted whiptail lizard); laboratory F1 hybrids of A. exsanguis × A. gularis; and the new species, composed of clones of two F1 hybrids representing generations up to F7. Additionally, for comparisons of DNA quantification and gene sequences of adenosine deaminase (ADA), we examined specimens of diploid gonochoristic Aspidoscelis arizonae (no standard common name yet); diploid gonochoristic Aspidoscelis burti (canyon spotted whiptail lizard); tetraploid parthenogenetic Aspidoscelis neavesi (Neaves' whiptail lizard); diploid parthenogenetic Aspidoscelis neomexicanus (New Mexico whiptail lizard); and triploid parthenogenetic Aspidoscelis uniparens (desert grassland whiptail lizard). All specimens of A. gularis are from one locality in Texas that is inhabited by the subspecies A. gularis gularis. Because the name Aspidoscelis is treated as masculine (Tucker et al., 2016), we use specific epithets to match that gender.
The taxonomy of the genus Aspidoscelis, particularly with respect to species boundaries, often undergoes considerable revision and may be marked by disagreement among specialists. This reorganization occurs largely because of extensive hybridization, both present and in the past, that has resulted in exceptionally blurred species boundaries (Barley et al., 2019; 2021b). The absence of a uniformly applied species concept exacerbates the problem. Justification for the names we use here is found in Cole et al. (2014, for A. neavesi), de Quieroz et al. (2017, for A. exsanguis, A. gularis, Aspidoscelis marmoratus, A. neomexicanus, and A. uniparens), Barley et al. (2021b, for Aspidoscelis arizonae, A. burti stictogrammus, Aspidoscelis opatae, and Aspidoscelis preopatae), and Barley et al. (2021a, for A. gularis).
Methods for maintaining Aspidoscelis in a captive colony were those used by Lutes et al. (2010, 2011) and Jewell et al. (2015). For the current project, individuals of A. exsanguis females and A. gularis males were simply housed together in enclosures and allowed to mate naturally. All clutches of eggs were incubated at 28 °C. In the laboratory records, first generation hybrids are referred to as representing the H1 generation. Their cloned offspring are the H2 generation, the next is the H3 generation, and so on. However, in this publication, as before, we use a modified reference to the generations to emphasize both the origin of the new species and the fact that individuals of cloned generations are not hybrids (they had only one parent). Herein, first generation hybrids are referred to as F1 hybrids. The next generation, cloned by individual F1 females, is the P1 generation (parental) of the new species. The next generation, cloned by P1 females, is the F1 generation (filial) of the new species, their offspring are referred to as the F2 generation, and so on. Consequently, our reference here to the F1 generation of the new species is equivalent to the H3 generation as noted in the laboratory records, and so on. We do not consider F1 hybrids as representing the P1 generation of the new species because not all F1 hybrid females are known with certainty to have reproduced, and half of the F1 hybrids are males (chi-square test, P = 0.32).
All specimens examined and localities sampled are listed in Appendix 1. The specimens are in the following collections: American Museum of Natural History (AMNH specimens), New York, USA; Johannes Gutenberg University (D specimens), Mainz, Germany; Museum of Comparative Zoology (MCZ specimens), Harvard University, Cambridge, Massachusetts, USA; and Stowers Institute for Medical Research (SIMR specimens), Kansas City, Missouri, USA.
DNA Content in Whole Blood Cell Nuclei
Nuclear DNA content was quantified by propidium iodide staining and flow cytometry. Blood (approximately 50 µl) was isolated from the tail and immediately transferred to 1 ml of citrate-dextrose solution (ACD, C3821, Sigma). The BD Cycletest Plus DNA kit from BD Bioscience (cat. no. 340242) was used to isolate and stain the nuclei following the manufacturer's instructions. After counting on a hemacytometer and adjusting the cell density to 1 × 106 cells/ml in buffer solution (contains sodium citrate, sucrose, and dimethyl sulfoxide), the cell suspension (1 ml) was placed into a round-bottom polystyrene test tube (Falcon 352008). Cells were centrifuged for 5 min at 300 g at room temperature in an Eppendorf 5702R centrifuge with swing-out buckets. After decanting the supernatant, the cells were resuspended by vortexing at low speed in 1 ml of buffer solution. Cells were centrifuged again for 5 min at 300 g, and the washing step was repeated. Cells were counted again and the concentration was adjusted to 5 × 105 cells/ml in buffer solution. The cell suspension was now centrifuged at 400 g for 5 min, the supernatant was carefully decanted, 250 µl of Solution A (contains trypsin in a spermine tetrahydrochloride detergent buffer) was added, and the solution was gently mixed and incubated for 10 min at room temperature. Then, 200 µl of Solution B (contains trypsin inhibitor and ribonuclease A in citrate stabilizing buffer with spermine tetrahydrochloride) was added, and the solution was gently mixed and incubated for 10 min at room temperature. Lastly, 200 µl of ice-cold Solution C (contains propidium iodide and spermine tetrahydrochloride in citrate stabilizing buffer) was added, and the solution was gently mixed and incubated in the dark at 2–8 °C for at least 10 min.
A minimum number of 10,000 events were collected on the BD LSRFortessa™ cell analyzer with a YG561 (561 nm) laser and 610/20 bandpass filter. The BD FACS-Diva v9.0.1 software was used for acquisition with identical instrument settings for all relevant samples. A quality check for the laser was performed before all measurements. For data analysis FlowJo v10.8.0 software was used. The gating strategy includes forward scatter (FSC)-area/side scatter (SSC)-area and area versus height for doublet discrimination.
Microsatellite DNA
Microsatellite DNA (msDNA) analysis was performed as previously described with a sample volume of 1.5 µl per PCR reaction (Lutes et al., 2011). Primers used for each microsatellite are listed in Table 1.
TABLE 1.
Oligonucleotide primers used for amplification of microsatellite (MS) DNA sequences.
DNA Isolation, Cloning, and Sequencing
Genomic DNA was extracted from tail samples (160–200 mg) following the instructions of the Qiagen Blood & Cell Culture DNA Midi Kit (Qiagen, Cat. No 13343). PCR primers (281T: GTTTGAAAT CTACATGCCAG and 282B: GCTTTGG TTTCCACAAACTC) were designed based on ADA cDNA from Aspidoscelis arizonae llanuras to amplify parts of exons 9 and 10 and the intervening intron from genomic DNA. The ADA gene sequences of A. marmoratus (GenBank ON661265) and A. gularis (GenBank ON661263) were used to confirm sequence identity of the primer annealing sites among the three species. The region of interest was amplified by PCR with an Eppendorf Mastercycler X50s in 20-µl reactions containing 4 µl of 5× Phusion HF Buffer, 0.2 mM dNTPs, 0.5 µM of each primer, 20 ng of genomic DNA, and 0.2 µl of Phusion DNA polymerase. The initial denaturation was set to 30 sec at 98 °C. Thirty cycles of 10 sec at 98 °C followed by 15 sec at 59 °C and 40 sec at 72 °C were used in the amplification followed by a final extension at 72 °C for 10 min. PCR products were subjected to electrophoresis through 1% agarose gels in Trisborate EDTA and purified with the QIA-quick Gel Extraction Kit (Qiagen, Cat. No. 28706). Purified PCR products were cloned into the pCR 4Blunt-TOPO vector (Thermo Fisher 450159), transformants were picked, and plasmid DNA was isolated for Sanger sequencing.
Sequence Analysis
Sequence alignments were created by Muscle (v 3.8.1551) within SnapGene Version 5.2.4 with the program's default parameters (anchorspacing 32, clustering by upgmb, diaglength 24, diagmargin 5, gap open 1, objscore sp, maxiters 16). The evolutionary history was inferred by the Maximum Likelihood and Tamura-Nei model in MEGA version 10.1.8 (Kumar et al., 2018). The Tamura-Nei model was used with uniform rates among sites, using all sites (including gaps). The estimation of the statistical confidence in the branch points of the tree was obtained by bootstrap resampling (500 replicates). Initial trees for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated by the Tamura-Nei model and then selecting the topology with superior log likelihood values. The cutoff value for condensed tree was set to 55%.
Morphological Characters and Statistical Procedures
The univariate characters are listed in Table 2. Nomenclature for epidermal scales follows Smith (1946). The sex of individual specimens was determined by dissection and examination of primary sex characters, a history of oviposition, and examination of relevant external secondary characters.
TABLE 2.
Univariate meristic characters used in statistical analyses of tetraploid, triploid, and diploid species compared in this paper.
We used NCSS 12 Statistical Software (2018; NCSS, LLC, Kaysville, Utah, USA) for database construction, statistical procedures and analyses, and scatterplot construction. Our database for morphological characters comprised five samples: 1) 36 specimens of parthenogenetic A. exsanguis (all but one being progeny produced at SIMR); 2) 30 specimens of gonochoristic A. gularis (field captures and their progeny produced at SIMR); 3) 12 of the F1 hybrids from A. exsanguis × A. gularis (from a total of six females and 10 males designated as the H1 generation by SIMR), with males and females pooled into one sample because they lacked significant differences (t test probabilities from 0.13 to 0.52) for the five univariate characters used in this study; 4) 30 female specimens representing cloned progeny from the F1 hybrid females and representing the P1–F3 generations of the new parthenogenetic species described herein (designated as H2–H5 generations at SIMR); and 5) 52 specimens of A. neavesi, a tetraploid parthenogenetic species that also originated in a laboratory environment at SIMR (Lutes et al., 2011; Cole et al., 2014). All specimens examined are identified in Appendix 1.
Multivariate analyses such as canonical variate analysis (CVA) and principal component analysis (PCA) are sensitive to the presence of outlier specimens (Tabachnick and Fidell, 2013). For example, a CVA of three samples (as a priori groups) based on univariate meristic characters will generate a set of canonical variate (CV) scores (CV1 and CV2) for each specimen. A scatterplot of these scores generates a cluster of points for each sample. Each point represents a specific specimen, and its location in the scatterplot is the point of coincidence between its CV1 and CV2 scores projected from the CV1 and CV2 axes. Each of the three clusters has a centroid (i.e., a center of gravity for that set of points). However, an outlier specimen in a sample skews the centroid's location by moving it in the direction of the outlier. Because this changes the mean values of CV1 and CV2, outliers can lead to erroneous conclusions regarding similarities and differences among the samples. Consequently, outliers should be routinely excluded from the definitive CVA and PCA analyses (Tabachnick and Fidell, 2013). We screened each sample for multivariate outliers on the basis of a combination of granules around midbody (GAB), femoral pore (FP), circumorbital scale (COS), lateral supraocular granule (LSG), and subdigital lamella (SDL) univariate characters (Table 2). We used T2 to test the Mahalanobis distance (D2) of each specimen from the centroid of the remaining specimens in that sample and identified outliers as specimens with distances yielding probabilities < 0.05. We excluded these few outlier individuals from our definitive statistical analyses. We also used Grubbs' univariate outlier test to test each sample for CV1 and principle component (PC)1 outliers. There were no CV1 outliers and the only PC1 outlier was excluded from the statistical tests and scatterplots comparing A. gularis, A. exsanguis, and the new tetraploid species. Outlier specimens are mentioned with analyses discussed in the appropriate sections of Interspecific Comparisons.
In previous publications we used PCA to concentrate the variation in 10 univariate characters into a reduced number of PCs. The PCs were used in the descriptions of two parthenogenetic species (A. neavesi Cole et al., 2014, and Aspidoscelis priscillae Cole et al., 2017). We used the PCs as variables in CVAs to discriminate various groups (Jombart et al., 2010). Because we used only five univariate characters in the present study, we used preliminary CVAs to compare the discrimination effectiveness between a CVA based on PC characters versus a CVA based on raw univariate characters. The CVA model that used raw univariate characters provided better discrimination: 63.6% reduction in classification error for the CVA of univariate characters (27 misclassifications) compared with 59.6% reduction in classification error (30 misclassifications) for the CVA of PCs. Therefore, we used stepwise selection of the five univariate characters in our definitive CVAs. However, we used a modified Levene's test to test the null hypothesis of no significant differences in PC1 and PC2 variances between certain pairs of samples. This Levene option is based on differences between specimen scores and sample medians rather than specimen scores and means, a robust alternative that functions well despite various departures from normal distributions (Brown and Forsythe, 1974). One-way analyses of variance of the five univariate characters and CV1 and CV2 revealed the presence of at least one significantly different sample for each of these variables, and we used Tukey-Kramer multiple comparison tests to identify which samples were significantly different.
DESCRIPTION OF NEW SPECIES
Aspidoscelis townsendae
New Species
Townsend's Whiptail Lizard
Figure 2
Holotype. MCZ R-199041 (= SIMR 16662), an adult female of the F1 laboratory-reared generation that was naturally cloned by its mother and that also cloned itself at the SIMR, producing, for example, F2 offspring MCZ R-199054 (= SIMR 19476) and MCZ R-199057 (= SIMR 19837). This cloning occurs in the lizards without any laboratory manipulation and is the only way they are known to reproduce.
Paratypes. See Appendix 1, Specimens Examined. Most individuals of A. townsendae other than the holotype are paratypes. A few are not designated as such, as they were still alive at this writing, sampled only for DNA, or photographed alive but not examined thereafter.
Diagnosis. A species of the Aspidoscelis sexlineatus species group as reviewed by Lowe et al. (1970). The species is distinguished from all others in the genus except A. neavesi by the following combination of characters: abruptly enlarged mesoptychials; enlarged postantebrachials; two frontoparietals; usually three parietals; four supraoculars each side; unisexual (only females exist); body with six longitudinal pale stripes that fade and can disappear in large adults; hatchlings basically unspotted but adults with conspicuous pale dots and spots on body; maximum snout–vent length (SVL) about 92 mm; tetraploid number of chromosomes (inferred from DNA quantification [see below]). The species differs from the tetraploid A. neavesi in having large adults with rather large, bright spots on the dorsum and upper legs (spots fade in A. neavesi); distal tail tan in adults (bluish to greenish brown in A. neavesi); usually having a higher number of GABs, COSs, and LSGs, as shown in Table 8; and msDNA alleles inherited from A. gularis (in A. townsendae) versus those of A. arizonae (in A. neavesi).
Description of Holotype. The following is based on MCZ R-199041 (= SIMR 16662). Paired data presented in the form x-y are for scale counts on the left–right sides of the body. SVL, 84 mm; rostral large, visible from above, wider than high; nostril low, posterior to center of nasal; nasals with a long median suture behind the rostral; frontonasal hexagonal; pair of irregularly hexagonal prefrontals with a long median suture; frontal hexagonal, longer than wide, wider anteriorly than posteriorly; pair of irregularly pentagonal frontoparietals with a long median suture; 3 irregular-sized and irregular-shaped parietals in a transverse series, the medial one basically heptagonal; irregular-sized and -shaped occipitals, posterior to but in contact with parietals, much smaller than parietals, much larger than dorsal granules; posterior of head covered with small granules, smaller than middorsal ones on body; supraoculars 4; 1st and 4th smallest, 2nd and 3rd largest; separated from superciliaries by 1 or 2 rows of granules (LSGs, 16-20), except 1st supraocular broadly contacts 1st superciliary. Last supraocular separated from frontoparietals and parietals by 1 or 2 rows of small scales (COSs, 6-5); postnasal, 1 on each side; loreal, 1 on each side, large, more or less rectangular; preocular, 1 on each side, with distinct ridge (keel). A row of 3 suboculars, the 2 anterior ones with a suborbital ridge continuing anteriorly onto the preocular, the 2nd subocular longest. Postoculars irregular and varied in size and shape; superciliaries 8-6, the 2nd (left) or 3rd (right) longest, the anterior ones elongate, posterior ones basically quadrangular. A few somewhat enlarged supratemporals; somewhat enlarged scales anterior to ear opening; central region of temple with small, roundish granules. Ear opening large, surrounded by small scales forming a smooth edge; external auditory meatus short, tympanum clearly visible; large supralabials 6-6, followed by small ones; suture between 5th and 6th below center of eye; the 3rd and 4th largest. Lower eyelid with a semitransparent disc of 5 enlarged palpebrals; pupil shape basically round to somewhat oval horizontally, but difficult to see undistorted as eyelids tightly closed.
Mental trapezoid with convex anterior edge; postmental basically pentagonal; 6 pairs of chinshields (= sublabials of some authors) curving posteriorly and dorsally to the lower labials, only those of anterior pair (largest) in contact at midline; chinshields from 2nd or 3rd pair to 5th pair separated in part or completely from infralabials by the interlabial scales (6-4); 7 enlarged anterior infralabials, followed by small scales; 3rd infralabial each side largest.
Gulars small, flat, mostly rounded, juxtaposed to slightly imbricate, somewhat larger on the anterior part of the throat, smaller posteriorly, from level of ear openings. Mesoptychial scales (on anterior edge of distinct gular fold) abruptly enlarged, slightly imbricate, smooth. Scales on nape and side of neck similar to dorsal and lateral body scales but smaller.
Dorsal and lateral scales granular, in indistinct transverse and oblique rows, with some unusually small ones here and there; ventrals large, usually somewhat rhomboidal, usually wider than long, imbricate, smooth, mostly in 8 longitudinal and 28 transverse rows (axilla to groin), the anterior rows on chest interrupted by a small triangular area of smaller scales. Number of dorsal GABs 70; preanal area with 4 clearly enlarged, smooth, irregularly shaped, juxtaposed or slightly imbricate scales plus smaller ones (PAS pattern I).
Femoral pores 19-20; usually each pore surrounded by 3 small scales (medial one largest); midventrally, 3 scales separate the femoral pore series of each side.
Scales on dorsal and lateral aspects of tail basically rectangular, obliquely keeled, slightly mucronate, somewhat imbricate, in transverse rows, keels forming longitudinal ridges; scales under tail wider, smooth, more imbricate; tail round in cross section.
Scales on upper and anterior surfaces of upper arm, on upper and anterior surfaces of forearm, on anterior and lower surfaces of thighs, and on lower surfaces of lower legs large, smooth, imbricate. Scales on lower and posterior surfaces of upper arm, on posterior and ventral surfaces of forearm, on posterior and upper surfaces of thighs, and on upper, anterior, and posterior surfaces of lower legs small, granular, juxtaposed (but postantebrachials on forearms enlarged, angular, and irregular in shape). Lamellae on ventral surface of 4th finger 15-15. Lamellae on ventral surface of 4th toe 32-31, usually in single rows on fingers but often paired on toes; fingers and toes somewhat laterally compressed; palms and soles with small, irregular, juxtaposed, flat scales; upper surfaces of hands and feet with large, imbricate, smooth, flat scales.
Color and Pattern of Holotype in Preservative (70% Ethanol). Dorsal ground color on body, arms, legs, and top of head brown. The light dots and spots on those surfaces are distinct, numerous, pale gray or white. The six body stripes are gray and largely faded, more visible posteriorly than anteriorly. Dorsal surface of tail tan, but gray where superficial epidermis has shed.
Ontogenetic Development of Color Pattern in Life. Hatchings of A. townsendae have an SVL of about 30–35 mm. The immediate impression is of a conspicuously striped dorsum without light spots, although distinct spots are on the arms and legs. The youngest individual photographed was 1 day of age (SIMR 13085, Fig. 2A), a P1 generation offspring of AMNH R-179122 (= SIMR 11065). The dorsal ground color was dark brown with 6 conspicuous light stripes (paravertebrals and dorsolaterals yellow, lateral stripe basically white); paravertebral light stripes somewhat wavy or zigzag, others basically straight; body basically unspotted but with slight hint of small tan spots, especially posteriorly in the dark fields. Some individuals have slight traces of light tan spots or a partial stripe in the vertebral dark field. Top of head, arms, legs, dark brown; arms and legs with conspicuous tan spots and irregular tan markings; body color continues onto base of tail, but then tail becomes tan. Ventral surfaces are very light tan, but darker tan under the tail, hands, and feet.
At 8 months of age (Fig. 2B), lizards differ from their hatchling appearance as follows: brown dorsal ground color is becoming reddish brown, especially anteriorly; the light stripes may be somewhat less distinct and the paravertebral stripes may have a green tint anteriorly; dorsum covered with conspicuous light spots on and between light stripes. On young lizards around this age also, the color beneath the tail may take on an orange hue, being orangish tan.
At 18 months old (Fig. 2C), SIMR 13085 differed further from the juvenile state as follows: dorsum ground color more reddish brown, especially anteriorly; dorsal light stripes broken up and fading, especially posteriorly; more numerous conspicuous pale tan, white, or cream-hued spots on body, especially on sacral area; spots on arms fading. Distal tail is tan, but it may appear to be pale gray to pale blue if covered in dead skin that is about to be shed away. At about this age, the ventral color of the hands, feet, and tail is pale tan.
The largest adult condition is illustrated by F5 SIMR 25022 at an age of 4 years and 4 months and SVL of 80 mm (Fig. 2D). On this individual, the light stripes have nearly completely disappeared, the ground color is reddish brown, the light spots on the body, legs, and base of tail have become enlarged and most conspicuous, whereas spots on the arms have faded, and the distal tail remains tan. Some large spots on the dorsum are in linear series where juvenile stripes had been; some lizards have very pale blue on the ventral surfaces.
The ontogenetic changes can be summarized as follows: ground color changes from dark brown to reddish brown; paravertebral and dorsolateral light stripes change from conspicuous and yellow to tan and less conspicuous with a green tint anteriorly, then later essentially gone in large adults. The dorsal body surface of hatchlings is essentially unspotted, but later develops light spots, which in the largest adults become quite enlarged and conspicuous. The tail distal to the base remains tan, but on some young individuals appears gray to pale blue if unshed dead skin is adhering; light spots on arms fade while those on legs become large and conspicuous. Pale blue may be seen on the ventral surface of some individuals.
Tetraploidy. To assess directly that the animals are tetraploid and maintain this ploidy over several generations, we used flow cytometry to compare the DNA content in nucleated erythrocytes of F5 and F6 generation A. townsendae with samples from other species of Aspidoscelis of known ploidy used as controls (Fig. 3). As expected, nuclei from tetraploid A. neavesi contained twice as much DNA as those from diploid A. neomexicanus, with nuclei of triploid A. uniparens harboring an intermediate DNA content (Fig. 3). Nuclei from A. townsendae contained the same amount of DNA as those from A. neavesi, confirming that the new species is also a tetraploid.
Future research on the tetraploid karyotype should reveal that it is similar or identical to that of A. neavesi as illustrated by Lutes et al. (2011), because the female progenitor is the same species from the same locality as that of A. neavesi, and the male progenitor (A. gularis) has the same karyotype as that of A. arizonae, the male progenitor of A. neavesi. Characteristics of msDNA and other alleles sequenced are discussed immediately below.
Hybrid Origin. To determine the genetic fingerprint of the tetraploids and to assess whether they harbored a genetic contribution from A. gularis, seven highly polymorphic msDNA markers were analyzed. The putative A. gularis father was heterozygous at five of those loci, whereas the triploid A. exsanguis mother showed three different and distinct alleles for all seven markers, consistent with the high degree of heterozygosity that resulted from the two rounds of hybridization involved in the origin of A. exsanguis. All alleles present in the clonal maternal A. exsanguis were observed in her sibling offspring (MCZ R-199029 [SIMR 10140] and MCZ R-199030 [SIMR 10141]), indicating faithful inheritance of the entire maternal triploid genome, as commonly observed during reproduction of parthenogenetic whiptail lizards. Additionally, each offspring also harbored one of the two alleles present in the diploid gonochoristic A. gularis male at each of the seven loci, indicating that the A. gularis had indeed made a genetic contribution to the offspring (Fig. 4). For two of the five markers for which the A. gularis was heterozygous, the same allele was passed on to both hybrid offspring, and for the other three markers (MS10, MS12, and MS14), different alleles from the father were found in the two offspring. In summary, analysis of each microsatellite marker revealed four distinct alleles in the two offspring, further confirming a tetraploid genotype and hybrid origin.
Clonal Inheritance. We have observed individuals of A. townsendae produced in the laboratory up to the F7 generation, and all were females. To assess whether the paternal and maternal alleles of msDNA identified in first generation hybrids were faithfully inherited in subsequent generations, we subjected hybrid A. exsanguis×A. gularis AMNH R-179122 (SIMR 11065), two of her daughters (AMNH R-179128 [SIMR 15214] and AMNH R-179130 [SIMR 16111]), and four granddaughters to msDNA analysis (Fig. 5). Among the granddaughters, SIMR 20471 and SIMR 20472 are direct descendants of AMNH R-179128 (SIMR 15214), whereas SIMR 20610 and SIMR 20778 represent individuals from two separate clutches of eggs produced by MCZ R-179130 (SIMR 16111). The msDNA loci tested were the same seven loci used to study the hybrid origin of the species, discussed above. All alleles from the two parental species were faithfully inherited across these three generations. Apparent single-nucleotide shifts can result from the binning of signals in whole nucleotide increments and are therefore not necessarily indicative of genetic changes from generation to generation. Most importantly, four distinct alleles persisted at each of the seven msDNA loci across generations of A. townsendae, indicating maintenance of the tetraploidy and heterozygosity inherited from the original hybrids (compare Figs. 4, 5).
Etymology. The specific epithet, a noun in the genitive singular case, honors Carol R. Townsend, retired Senior Scientific Assistant at the AMNH. She pioneered development of methods to maintain laboratory colonies of whiptail lizards, starting in about 1971 (e.g., Cole and Townsend, 1977; Townsend, 1979; Townsend and Cole, 1985) and which, after year 2002, were modified and upscaled by one of us (D.P.B.) and colleagues at SIMR (Jewell et al., 2015). Work with lizards of known parentage and genealogy has allowed comparing offspring of known lineages with each other, which cannot be done with certainty with animals collected in nature. It also allowed obtaining critically important direct evidence on aspects of reproductive biology, patterns of inheritance, and evolution of clonal amniotes. Carol R. Townsend is the wife of C.J.C., but he did not participate in choosing the specific epithet. The majority of other coauthors simply announced their choice of the name and did not allow C.J.C. to discuss the matter.
Comments. The specimen selected as the holotype is an F1 individual that was cloned by a P1 female, which had been cloned by a hybrid A. exsanguis × A. gularis. Additionally, the holotype produced F2 offspring, all females. For this study, we have examined and photographed more than 109 individuals of A. townsendae representing specimens from the P1–F7 generations. Initially, four different individual A. exsanguis (MCZ R-198981 [= SIMR 7579], SIMR 8774, MCZ R-198982 [= SIMR 9049], and one of either SIMR 8284 or 8312) produced hybrids with four different males of A. gularis (MCZ R-199024 [= SIMR 8170], SIMR 8868, AMNH R-179105 [SIMR 8869], and AMNH R-179107 [SIMR 8996]). A total of 16 F1 hybrids were produced, six females and 10 males (which would have depended on whether the X-bearing or Y-bearing spermatozoan from the A. gularis was involved). Two of the F1 hybrid offspring (MCZ R-199029 [SIMR 10140] and AMNH R-179122 [SIMR 11065]) were known maternal parents of cloned offspring that represented the P1 generation of A. townsendae. See Cole et al. (2017) for a discussion of why we do not follow the suggestion of some authors (e.g., Frost and Wright, 1988) that each F1 hybrid should be considered the founder of a separate species if it reproduces successfully through parthenogenesis. We apply only one name for the tetraploid species described here. Aspidoscelis townsendae originated on 5 December 2013, when the first known F1 individual hatched.
Figure 2.
Color and pattern ontogeny of A. townsendae. A–C. SIMR 13085, P1 generation, at the following ages: A, 1 day after hatching; B, 8 months of age; C, 18 months of age. D, SIMR 25022, F5 generation, at 4 years and 4 months of age, snout–vent length, 80 mm.
TABLE 8.
Descriptive statistics and comparisons of means of morphological meristic characters among tetraploid Aspidoscelis townsendae and A. neavesi and triploid A. exsanguis, maternal progenitor of the two tetraploid species.
Figure 3.
DNA content in whole blood nuclei determined by propidium iodide staining and detection by flow cytometry. Samples from diploid A. neomexicanus (coral), triploid A. uniparens (dark green), and tetraploid A. neavesi (light green) were used as markers for the respective levels of ploidy. Samples from three individuals of A. townsendae represent the fifth generation of the lineage of MCZ R-199029 (SIMR 10140), SIMR 25022, and sixth generation of either the lineage of MCZ R-199029 (SIMR 10140) or AMNH R-179122 (SIMR 11065), D 500, and D 501. Number of events scored by flow cytometry were 10,986 (A. neomexicanus), 10,931 (A. uniparens), 11,633 (A. neavesi), 10,524 (D 500), 11,137 (D 501), and 10,548 (SIMR 25022).
Figure 4.
Microsatellite DNA analysis confirms hybridization. Animal identification numbers are shown to the left of each microsatellite (MS) data panel. Relationships among the four animals are shown in the lower right corner. The two hybrids originated from one clutch of eggs. For each marker, alleles present in the A. exsanguis and inherited by the hybrid offspring are shown in shades of pink; alleles from the A. gularis are underlaid in shades of blue. For MS1 and MS15 only one allele was detected in the A. gularis male, suggesting homozygosity for these loci. The numbers in the boxes correspond to the length of the analyzed DNA fragment rounded to the nearest integer. Single nucleotide deviations observed for MS7 and MS12 are probably rounding artefacts (see text).
Figure 5.
Clonal reproduction maintains ploidy over three generations. Upper. Microsatellite (MS) analysis of a first generation hybrid (AMNH R-179122 [SIMR11065]), one of its descendants AMNH R-178128 (SIMR15214), an A. townsendae founder, and two of its offspring (SIMR 20471 and 20472). For each marker, alleles originating from the A. exsanguis maternal ancestor are shown in shades of pink; alleles from the paternal A. gularis ancestor are underlaid in blue. Because of the inherent methodological noise affecting single-nucleotide differences, two MS15 alleles cannot be unequivocally identified and are shown in yellow. Lower. Microsatellite analysis of a first generation hybrid (AMNH R-179122 [SIMR 11065]), one of its descendants (AMNH R-179130 [SIMR 16111]), an A. townsendae founder, and two of its offspring from different clutches of eggs (SIMR 20610 and SIMR 20778). Colors are as in Upper, but allele 293 at MS12 would be from A. gularis.
Interspecific Comparisons among A. townsendae and Its Progenitors
Colors and Pattern. Aspidoscelis townsendae looks very similar to its maternal parent, A. exsanguis, throughout growth from hatchlings to young adults. However, the largest adults of A. townsendae have large and conspicuous light spots on the body and hind legs, whereas these are very much less conspicuous and smaller in A. exsanguis from Alamogordo (compare Figs. 1, 2). Adults of A. gularis have very different colors and pattern (Fig. 1). Additionally, A. gularis is significantly different from both A. exsanguis and A. townsendae in all five of the scalation characters compared herein (Table 3), illustrating that maternal inheritance overwhelms paternal inheritance in all these characters.
TABLE 3.
Descriptive statistics and comparisons of means of morphological meristic meristic characters among parthenogenetic Aspidoscelis exsanguis from alamogordo, New Mexico; diploid gonochoristic A. gularis from texas; their tetraploid F1 hybrids (pooled males and females); and tetraploid parthenogenetic A. townsendae cloned from F1 hybrids and representing the P1–F3 generations (pooled).
One additional scalation character involves sizes of the dorsal granules as seen on the side of the body. In some places in A. gularis, many small granules are much smaller than the more abundant larger dorsal granules. Most of these are in incomplete rows of small granules between the larger ones, but some of the small granules are scattered here and there among the others, which explains the higher and more variable GAB count in A. gularis (Table 3). Such small granules are seen on close inspection in A. exsanguis, but to a much smaller extent than in A. gularis, and A. townsendae has a condition intermediate to its parental species.
Multivariate Comparisons of Scalation. We used a CVA to reveal morphological relationships of A. townsendae, its progenitor species A. exsanguis and A. gularis, and the F1 hybrid generation from which A. townsendae was derived. Three outlier specimens were excluded: one A. exsanguis, two A. townsendae.
All five univariate characters met the stepwise selection criterion for inclusion in the CVA model, and the first two canonical variates (CV1 and CV2) summarized 99.1% of the intergroup variation (Table 4). The relative importance of each univariate character to the discrimination of these samples is indicated by the strength of its correlations with the canonical variates (CV1–CV3) generated by the CVA (Table 4). We used the correlation index of Comrey and Lee (1992) recommended by Tabachnick and Fidell (2013): a correlation larger than 0.71 (50% overlapping variance) is excellent; a correlation of 0.63 (40% overlapping variance) is very good; one of 0.55 (30% overlapping variance) is good; and a correlation of 0.45 (20% overlapping variance) is fair. None of the univariate characters was particularly distinctive for its contribution to the discrimination in CV1, but the contribution of CV2 was influenced primarily by SDL (Table 4).
TABLE 4.
Correlations between univariate characters (CVs and PCs) derived from canonical variate and principal component analysis of five and multivariate characters chasactersandongss derived from canonical variate and principal component Aspidoscelis exsanguis from New Mexico, gonochoristic diploid A. gularis from Texas, their tetraploid F1 hybrids (pooled males and females), and tetraploid parthenogenetic A. townsendae cloned from F1 hybrids and representing the P1–F3 generations (pooled). Eigenvalues and proportions of variation summarized by each multivariate character are also shown.
The absence of any misclassified specimens in the sample of A. gularis illustrates the meristic distinctiveness of this species (Table 5). Conversely, 31.9% of the A. townsendae sample was misclassified to the A. exsanguis and F1 hybrid groups, thereby reflecting the strong matriclinous resemblance of A. townsendae to its maternal parent and immediate ancestors (Table 5).
TABLE 5.
Classification results of a canonical variate analysis (CVA) of samples of triploid parthenogenetic Aspidoscelis exsanguis from New Mexico, gonochoristic diploid A. gularis from Texas, their tetraploid F1 hybrids (pooled males and females), and tetraploid parthenogenetic A. townsendae cloned from F1 hybrids and representing the P1–F3 generations (pooled). Rows show number of individuals assigned to each group by the CVA.
The pattern of morphological relationships is clearly shown by comparisons of CV1, the variable summarizing 96.3% of the intersample variation (Table 4). Canonical variate 1 shows a strong resemblance among A. townsendae, F1 precursors, and maternal progenitor A. exsanguis (Table 3). This pattern is reiterated by PC1 and GAB, FP, COS, and LSG univariate characters (Table 3). However, CV2, summarizing only 2.8% of the intersample variation, shows a stronger meristic resemblance of A. townsendae and its F1 precursors with paternal A. gularis (Table 3) because of the strong influence of SDL on CV2 (Table 4). Patterns of morphological relationships are depicted by scatterplots (Figs. 6, 7).
Figure 6.
Scatterplot of canonical variate scores of the tetraploid A. townsendae, its progenitor triploid species A. exsanguis and diploid A. gularis, and the tetraploid F1 hybrid generation from which A. townsendae was derived. The canonical variate analysis model is shown in Table 4. Axis percentages are proportions of among-group variance accounted for by CV1 and CV2.
Figure 7.
Scatterplot of principal component scores of the tetraploid A. townsendae, its progenitor triploid species A. exsanguis and diploid A. gularis, and the F1 hybrid generation from which A. townsendae was derived. The PCA model is shown in Table 4. Axis percentages are proportions of combined within-group and among-group variance accounted for by PC1 and PC2.
Univariate Comparisons of Scalation. Data for five meristic characters are presented in Table 3. For all univariate characters, diploid A. gularis differs significantly from triploid A. exsanguis, the tetraploid F1 hybrids and A. townsendae. However, there are no significant differences among A. exsanguis, the F1 hybrids, and A. townsendae in GAB, FP, COS, and LSG. Some significant differences are indicated for SDL (and SVL), but the differences are small and observed ranges of variation overlap widely, suggesting that the best way to distinguish individuals of A. townsendae and A. exsanguis in the field is by coloration. This preliminary identification should be followed up with molecular data in the laboratory.
Multivariate Comparisons of Scalation Among Tetraploid A. townsendae and A. neavesi and Their Maternal Progenitor, Triploid A. exsanguis
Thirteen outlier specimens were excluded from this CVA: two of A. townsendae and 11 of A. neavesi. Four of the five univariate characters met the stepwise selection criterion for inclusion in the CVA model, and the first two canonical variates (CV1 and CV2) summarized 100% of the intergroup variation; CV1 was influenced primarily by GAB and LSG, and CV2 by SDL (Table 6). Although both A. townsendae and A. neavesi were derived by hybridization from the same stock of A. exsanguis from Alamogordo, New Mexico, their morphological meristic relationship to A. exsanguis was quite different (Table 7; Figs. 6, 8). Aspidoscelis townsendae had 10 misclassifications to the A. exsanguis sample, whereas the sample of A. neavesi had none (Table 7), thereby providing a contrast to the matriclinous resemblance of A. townsendae to A. exsanguis, from which A. neavesi is more distinctive.
TABLE 6.
Correlations between univariate characters and canonical variates derived from a canonical variate analysis of four univariate characters among samples of tetraploid Aspidoscelis townsendae and A. neavesi and triploid A. exsanguis, maternal progenitor of the two tetraploid species. Eigenvalues and proportions of variation summarized by each canonical variate are also shown.
TABLE 7.
Classification results of a canonical variate analysis (CVA) of samples of tetraploid Aspidoscelis townsendae and a. neavesi and triploid A. exsanguis, maternal progenitor of the two tetraploid species. Rows show number of individuals assigned to each group by the CVA.
Figure 8.
Scatterplot of canonical variate scores of tetraploids A. townsendae and A. neavesi and their maternal progenitor, triploid A. exsanguis. The canonical variate analysis model is shown in Table 6. Axis percentages are proportions of among-group variance accounted for by CV1 and CV2.
Several alternative groups of samples are identified by various combinations of significant and nonsignificant differences for particular sets of univariate characters, thereby illustrating the value of CVA in revealing a comprehensive pattern of morphological variation that lends itself to interpretation (Table 8). Canonical variable 1, with the greatest discriminatory power, shows A. townsendae and A. exsanguis differing significantly from A. neavesi but not from each other, and CV2, which summarizes much less data, shows a pattern in which each of A. townsendae, A. exsanguis, and A. neavesi differs significantly from the other two species (Table 8; Fig. 8).
Patterns of Meristic Morphological Variability Among Clonal Tetraploid A. townsendae, Its Clonal Triploid Maternal Progenitor Species A. exsanguis, and Its Diploid Gonochoristic Paternal Progenitor A. gularis
The basic question addressed here is the null hypothesis that there are no differences among sample variances of clonal species versus gonochoristic species. One tends to think that the clonal species would be less variable, but previous comparisons have shown that variability in counts of scales and FPs can be similar in clonal and nonclonal taxa (Taylor et al., 2012; Cole et al., 2016, 2017), and this pattern is shown again here. All five univariate characters used in this PCA made significant contributions to the summary of meristic variation provided by PC1 (Table 9). The pattern of variability revealed was an absence of significant differences in PC1 variances among A. townsendae, A. exsanguis, and A. gularis (Table 10; Fig. 9), but there was a difference in variability between the two tetraploids (A. townsendae and A. neavesi; Table 10; Fig. 10), illustrating that two parthenogenetic species often have morphological variation equivalent to that expressed by sexually reproducing A. gularis. Of lesser importance (because less variation is summarized), A. exsanguis, A. gularis, and A. townsendae each differed significantly from the other two species in PC2 variances (Table 10).
TABLE 9.
Correlations between univariate characters and principal components from a principal components analysis (PCA) of tetraploid clonal Aspidoscelis townsendae, triploid clonal A. exsanguis, and diploid gonochoristic A. gularis. Eigenvalues and proportions of variation summarized by each principal component are also shown.
TABLE 10.
Comparative morphological variability assessed from two principal component analysesa of five meristic characters among 1) tetraploid clonal Aspidoscelis townsendae, triploid clonal A. exsanguis, and diploid gonochoristic A. gularis (total N = 88, Table 9) and 2) A. exsanguis, A. townsendae, and tetraploid clonal A. neavesi (total N = 169, Table 11). Relative variation summarized in principal components PC1 and PC2 is assessed by a modified Levene's test.b
Figure 9.
Scatterplot of principal component scores of the tetraploid A. townsendae and its progenitor species triploid A. exsanguis and diploid A. gularis. Axis percentages are proportions of combined within-group and among-group variance accounted for by PC1 and PC2. This PCA was based on the model shown in Table 9.
Figure 10.
Scatterplot of principal component scores of the tetraploids A. townsendae and A. neavesi and their maternal progenitor, triploid A. exsanguis. Ellipses define the 95% confidence limits for score distributions, and axis percentages are proportions of combined within-group and among-group variance accounted for by PC1 and PC2. This PCA was based on the model shown in Table 11.
In reptiles, temperature during egg incubation can have a significant effect on development of certain morphological characters (Fox, 1948; Fox et al., 1961). Consequently, one may question whether incubation temperature is a factor here and for the observations presented in the next paragraph. However, variation in temperature is not a factor in these comparisons because all clutches of eggs involved in this study were incubated at SIMR in incubation chambers maintained at 28 °C, whether of clonal or gonochoristic species.
Patterns of Meristic Morphological Variability Among Clonal Tetraploid A. townsendae and A. neavesi and Their Clonal Triploid Maternal Progenitor Species A. exsanguis
This PCA comparison contrasts with the one above in specific meristic characters being associated with one or the other principal component. PC1 summarized variation primarily in GAB, COS, and LSG characters, and PC2 summarized variation primarily in FP and SDL (Table 11). Both PC1 and PC2 revealed the same pattern: homogeneous variability between both A. townsendae and A. neavesi and their maternal progenitor A. exsanguis, but significant differences in the variances of A. townsendae and A. neavesi (Table 10). The explanation for this pattern may reside in the different haploid genomes contributed by different paternal gonochoristic progenitors in the hybridizations with A. exsanguis. Aspidoscelis townsendae arose from hybridization with A. gularis, and A. neavesi arose from hybridization with A. arizonae (sensu Barley et al., 2021b).
TABLE 11.
Correlations between univariate characters and principal components from a principal components analysis (PCA) of triploid clonal Aspidoscelis exsanguis, and tetraploid clonal A. townsendae and A. neavesi. Eigenvalues and proportions of variation summarized by each principal component are also shown.
IDENTIFYING POSSIBLE HYBRIDS FROM NATURE
To date, no specimens collected in nature have been reported as putative hybrids between A. exsanguis and A. gularis. However, these two species are in contact in an area in southeastern New Mexico, where either hybrids or specimens of A. townsendae may be found in the future. The Museum of Southwestern Biology (MSB), University of New Mexico, has numerous specimens of both parental species from in and around the city of Carlsbad, Eddy County, New Mexico, and within Carlsbad Caverns National Park. Because young representatives of A. exsanguis and A. townsendae are very similar to each other, the specimens on which to focus are the largest adults that resemble A. exsanguis but which differ from it as follows: light spots on dorsal body and hind legs are rather large and conspicuous (compare Figs. 1, 2). Although all scale counts are similar in A. exsanguis and A. townsendae, specimens with SDL greater than 32 could be A. townsendae (Table 3). Lizards that appear to be possible A. townsendae should have their ploidy checked by karyotyping or by DNA quantification if live individuals are available. Microsatellite DNA analysis and ADA intron sequencing might also be applied to preserved specimens of any age, which will reliably reveal differences in ploidy. J. Tom Giermakowski of MSB provided photographs of the six largest specimens from the Carlsbad area identified as A. exsanguis (MSB 61843, 61844, 61872, 61879, 66919, and 73615), which allowed C.J.C. to confirm the identification of each.
EVOLUTIONARY HISTORY OF A. TOWNSENDAE
The ancestry of A. townsendae involved three steps of hybridization over a period of at least hundreds, but probably thousands, of years as the distribution of habitats and lizards have shifted with environmental changes in what is now called northwestern Mexico and southwestern United States. In the last 10,000 years the major habitats of relevant woodland, chapparal, grassland, and deserts have significantly shifted in geographic distribution. For example, where deserts occur today, woodlands occurred in the recent past (Van Devender, 1977; Betancourt et al., 1990). Populations of whiptail lizards would have been shifting also, allowing for considerable chances for hybridization (Wright and Lowe, 1968). In the origin of A. townsendae, the initial two steps of hybridization in nature involved three species of Aspidoscelis, and the last step involved a fourth species in the laboratory at the SIMR.
Here we follow recent suggestions that best reflect the most reasonable taxonomy and nomenclature for the relevant species of whiptail lizards and the evolutionary origins of the relevant parthenogenetic species through hybridization events (Good and Wright, 1984; Dessauer and Cole, 1989; Moritz et al., 1989; Reeder et al., 2002; and Barley et al., 2021b). The taxonomic conclusions for treatment of the specific names A. gularis, and Aspidoscelis scalaris or Aspidoscelis septemvittatus follow de Quieroz et al. (2017), Barley et al. (2019, 2021a), and Barley (personal communication).
The oldest, first round of hybridization between two gonochoristic species in nature involved A. burti (female) × A. arizonae, probably in northern Sonora, Mexico, the result of which is the recently described diploid parthenogenetic species A. preopatae. The second round of hybridization in nature involved the parthenogenetic A. preopatae × A. scalaris, a diploid gonochoristic species, and the result of this hybridization, probably in northern Mexico, is the triploid parthenogenetic species A. exsanguis. The third round of hybridization is the round reported here of A. exsanguis × A. gularis. The DNA data that we have presented in this paper fully supports these conclusions.
All seven of the msDNA loci we analyzed show that A. townsendae inherited the three alleles present in the maternal A. exsanguis and one of the two alleles present in the paternal A. gularis (Fig. 4). Furthermore, the three alleles at each locus in A. exsanguis are consistent with the species having received them from three ancestral species, including data matching msDNA for A. arizonae and A. burti. We propose that the third msDNA allele was contributed by A. scalaris.
Additionally, we present new DNA sequence information and interpret this in the context of the evolutionary history. As A. exsanguis originated from two successive hybridization events involving three distinct species, we surmised that the three alleles of a single-copy gene should be distinguishable at the level of DNA sequencing of any genomic region where selection is weak, such as introns. We chose the well-characterized ADA locus and designed primers within exons 9 and 10 to amplify the intervening intron 9 from genomic DNA from an A. exsanguis (SIMR 22519). A 1.8-kb PCR product was isolated and cloned, and 40 clones were sequenced to obtain representations of all alleles present. As expected for a triploid species harboring the genomes of three distinct parental species, the dataset clustered into three groups corresponding to the three different alleles of the ADA locus (Fig. 11A). When the dataset was combined with sequences representing intron 9 derived from individuals representing two of the three parental species (A. gularis and A. burti; SIMR 007 and SIMR 12836, respectively), the group of 15 A. exsanguis sequences clustered with the 21 sequences from A. gularis (Fig. 11B). The group of seven A. exsanguis sequences clustered with two distinct A. burti groups containing 10 and nine sequences, respectively. The partition of the sequences derived from a single A. burti individual likely reflects a high degree of heterozygosity at the ADA locus in this animal; an equivalent split is not observed among 21 sequences derived from the A. gularis. Because we did not have tissue samples from A. scalaris, we hypothesize that this “green” group (Fig. 11B) was contributed by that species, with A. gularis representing its closest relative in this dataset. The remaining group of 18 A. exsanguis sequences clustered separately, and represents the allele derived from the A. arizonae ancestor (red).
Figure 11.
Analysis of ADA intronic sequences from A. townsendae and its ancestral species confirm the presence of four distinct genomes in A. townsendae. A. An alignment of 40 sequences from A. exsanguis (SIMR 22519) over 1,199 positions was used to generate a maximum likelihood tree. Sequences clustered in three groups, the numbers in parentheses denoting the number of sequences in each cluster. The color-coded boxes represent the ancestral species with the closest match: A. arizonae (red), A. burti (blue), and A. scalaris (green; inferred because no A. scalaris genomic sequence was available to us). Maximum likelihood values are shown next to the branches. B. Sequences of A. exsanguis (SIMR 22519; 40 sequences), A. burti (SIMR 12836; 19 sequences), and A. gularis (SIMR 007; 21 sequences; purple) aligned over 1,168 positions were used to generate a maximum likelihood tree. The A. burti sequences cluster in two groups reflecting a high degree of heterozygosity of the sequenced individual. Red and green boxes of A. exsanguis represent the ancestors of that species, with colors interpreted as in panel A. C. Sequences (118) of A. exsanguis (SIMR 22519), A. burti (SIMR 12836), A. gularis (SIMR 007), and A. townsendae (SIMR 25022) were aligned over 1,118 positions to generate a maximum likelihood tree. D. Sequences of two tetraploid individuals of A. townsendae, SIMR 25022 (38 sequences) and D 500 (41 sequences), were aligned over 1,176 positions to generate this maximum likelihood tree. The same clusters of four groups are visible in the two tetraploid individuals.
Next, we added 40 sequences of intron 9 of A. townsendae (SIMR 25022) to the dataset and aligned a total of 118 sequences over 1,118 positions (Fig. 11C). Two chimeric sequences were removed from the alignments at this stage. Although we cannot exclude the possibility that these sequences represent recombination events between homologous chromosomes within A. townsendae, the more parsimonious explanation is that the chimeras were generated during PCR amplification, a known artefact when amplifying multiple similar sequences within one reaction. Importantly, the remaining 38 A. townsendae sequences clustered as expected for a tetraploid harboring four distinct genomes, three representing the genomes present in the maternal A. exsanguis (green, red, and blue) and one genome derived from the A. gularis father (purple). Finally, we aligned the sequences derived from two individuals of A. townsendae (SIMR 25022 and D 500) representing generations 5 and 6 of the new species. Again, we found four distinct clusters, each containing a subset of the sequences from each of the two individuals (Fig. 11D). These observations are consistent with the conclusions that A. townsendae comprises genomes from four ancestral bisexual species and that the tetraploid condition and heterozygosity are maintained from generation to generation. The evolutionary history of A. townsendae is summarized in Figure 12.
Figure 12.
Phylogenetic origin of the tetraploid A. townsendae involving three hybridization events among four gonochoristic ancestral species and two parthenogenetic ancestors. Colors indicate transfer of specific ancestral genomes, with cloned genomes of parthenogenetic species in black. Each hybridization event resulted in a new parthenogenetic species, and the second and third events resulted in increases in ploidy. For clarity, only taxa involved in the origin of A. townsendae are included in the drawing (e.g., some taxa that could be plotted between A. arizonae and A. gularis are omitted). Origins of the other tetraploid species of Aspidoscelis (A. neavesi and A. priscillae) involved similar events but with a different mix of ancestral genomes contributed by the gonochoristic ancestors A. arizonae and A. burti. ARI, A. arizonae (red); PRE, A. preopatae (black); BUR, A. burti (blue); EXS, A. exsanguis (black); SCA, A. scalaris (green); TOW, A. townsendae (black); GUL, A. gularis (purple).
ACKNOWLEDGMENTS
We thank Rick Kupronis and the team of dedicated reptile technicians at SIMR—David Jewell, Alex Muensch, Christina Piraquive, and Kristy Winter—as well as Martin Fahr and his team at Johannes Gutenberg University for outstanding husbandry and herpetological skills. We also thank Troy Green and his colleagues in the Molecular Biology Core at SIMR, as well as Alex Muensch and Jillian Kupronis, for contributions to msDNA analysis. We thank Stefanie Möckel and the Flow Cytometry Core Facility at the Institute for Molecular Biology in Mainz for DNA content analysis and the Emergent AI Research Center for computational assistance. We also acknowledge computing time granted on the supercomputer MOGON 2 at Johannes Gutenberg University, Mainz (hpc.uni-mainz.de). The following individuals at museum repositories provided efficient curating and photography of specimens: David Dickey, David Kizirian, and Lauren Vonnahme (AMNH); James Hanken and Joe Martinez (MCZ); Tom Giermakowski and Howard Snell (MSB).
Charles Painter, Letitia Mee, and Samantha Ferguson, New Mexico Department of Game and Fish, helped with permit matters. The A. gularis were all collected at one locality in Texas (see Specimens Examined) with specific permission of the landowners, the family of W.B.N. As usual, Carol R. Townsend (AMNH) assisted in many ways.
This work was funded in part by the Howard Hughes Medical Institute, the SIMR, Johannes Gutenberg University, Mainz, and an Alexander von Humboldt Professorship awarded to P.B. We also gratefully acknowledge the support by the GenEvo RTG funded by the Deutsche Forschungsgemeinschaft (German Research Foundation) – GRK2526/1 – Projectnr. 40.7023052.
DATA AVAILABILITY
The ADA intronic sequences have been deposited in GenBank with the following accession numbers: ON661263, A. gularis; ON661264, A. arizonae; ON661265, A. marmoratus. The multiple sequence alignment FASTA files used to generate the maximum likelihood trees are available on GitHub (2022).
LITERATURE CITED
Appendices
APPENDIX 1. SPECIMENS EXAMINED
Aspidoscelis arizonae llanuras: SIMR 17121 was used for the DNA sequence analyses of ADA (adenosine deaminase). This specimen was collected at the following locality: New Mexico: Otero County; Alamogordo, at the coordinates 32°52′56.6″N, 105°57′45.9″W in WGS 84.
Aspidoscelis burti stictogrammus: SIMR 12836 was used for the DNA sequence analyses of ADA. This specimen was collected at the following locality: Arizona: Cochise County; vicinity of the headquarters of the Muleshoe Ranch Preserve of the Nature Conservancy, at the coordinates 32°20′14.8″N, 110°14′13.7″W in WGS 84.
Aspidoscelis exsanguis: Individuals with underline were used in the multivariate statistical analyses. MCZ R-100419; R-192263; R-192287–192292; R-198981–198990 (SIMR 7579; 9049; 9976; 10684; 10798; 11573; 11837; 12280; 12281; 12901, respectively); AMNH R-179076 (SIMR 13217); MCZ R-198991–198994 (SIMR 13550; 13621; 14235; 14315, respectively); AMNH R-179071 (SIMR 14316); AMNH R-179077 (SIMR 14317); AMNH R-179072–179073 (SIMR 14347; SIMR 14348, respectively); AMNH R-179074 (SIMR 14563); AMNH R-179078–179079 (SIMR 14588 and 14639, respectively); AMNH R-179075 (SIMR 14654); AMNH R-179081–179083 (SIMR 20133; 20335; 22262, respectively); AMNH R-179080 (SIMR 23029); AMNH R-179084 (SIMR 23793). Additionally, SIMR 22519 was used for the ADA DNA sequence analyses. All specimens were collected in the field in or within walking distance of the city park, Alamogordo, Otero County, New Mexico, or reared at SIMR from stock obtained at this locality.
Aspidoscelis gularis gularis: Individuals with underline were used in the multivariate statistical analyses. MCZ R-199013–199027 (SIMR 134; 186; 4191; 4225; 4859; 5471; 8002; 8003; 8019; 8032; 8098; 8170; 8213; 8809; 8810, respectively); AMNH R-179104–179118 (SIMR 8865; 8869; 8933; 8996; 9531; 9532; 9534; 9736; 9992; 10110; 10111; 10638; 17135; 17983; 17985, respectively). Additionally, SIMR 007 was used for the ADA DNA sequence analyses. All specimens were collected in the field 1 mi E of Spur, Dickens County, Texas, or reared at SIMR from stock obtained at this locality.
F1 laboratory hybrids of A. exsanguis × A. g. gularis: Individuals with underline were used in the multivariate statistical analyses. MCZ R-199028–199035 (SIMR 9478; 10140; 10141; 10152–10154; 10708; 10776, respectively); AMNH R-179119–179122 (SIMR 10777–10778; 10878; 11065, respectively).
Aspidoscelis marmoratus: SIMR 8450 was used for the DNA sequence analyses of ADA. This specimen was hatched and raised at SIMR from stock collected at the following locality: New Mexico: Sierra County; vicinity of Truth or Consequences.
Aspidoscelis neavesi: All these specimens were used in the multivariate statistical analyses, cataloged in the museums for use in a previous study, and hatched and raised in captivity from the same laboratory stock reported by Cole et al. (2014). AMNH R-176080–176085; 176087–176090; 176092–176096; 176099–176101; 176103–176105; 176110–176121; 176123–176125; 176126–176136; 176138–176146; 176148. MCZ R-192218–192219; 192224; 192226; 192233–192234; 192236–192237; 192241–192246; 192248; 192250–192252; 192255; 192258–192262; 192265–192267; 192269–192271; 192273–192275; 192277–192280; 192282; 192284–192286. Specimen D 459 was used for DNA quantification analyses.
Aspidoscelis neomexicanus: Specimen D 250 was used for DNA quantification analyses. It was hatched and raised in the laboratory from stock most likely collected in Socorro County, New Mexico.
Aspidoscelis townsendae: All these specimens are paratypes unless specified otherwise. Individuals with underline were used in the multivariate statistical analyses. MCZ R-199036–199039 (SIMR 12461; 12462; 13156; 13263, respectively); AMNH R-179123–179131 (SIMR 13541; 14551; 14867; 14943–14944; 15214; 15794; 16111–16112, respectively); MCZ R-199040–199044 (SIMR 16368; 16662–16663; 16705; 16836, respectively); AMNH R-179132; 179133 (SIMR 16837; 16988, respectively); MCZ R-199045–199049 (SIMR 17174–17175; 17869; 18077; 18441, respectively); AMNH R-179134–179140 (SIMR 18452–18453; 18499; 18617; 18671; 18718; 18800, respectively); MCZ R-199050–199059 (SIMR 19059; 19165; 19338; 19351; 19476; 19596; 19619; 19837; 21612; 21867, respectively); AMNH R-179141–179150 (SIMR 21868; 22029–22030; 22099; 22201–22202; 22290; 22393; 22700; 14942, respectively). Additionally, SIMR 25022 and D 500 were used for DNA quantification and ADA DNA sequence analyses, and D 501 was also used for ADA sequence analysis; these individuals are not paratypes.
Aspidoscelis uniparens: Specimen D 489 was used for DNA quantification analyses. It was hatched and raised in the laboratory from stock collected in Arizona: Cochise County; along Arizona Highway 186 SE Willcox, ca. 32°11′43.2″N, 109°45′5.1″W at 1,268 m elev., in WGS 84.
