Origins and biogeography of the Anolis crassulus subgroup (Squamata: Dactyloidae) in the highlands of Nuclear Central America

Background Recent studies have begun to reveal the complex evolutionary and biogeographic histories of mainland anoles in Central America, but the origins and relationships of many taxa remain poorly understood. One such group is the Anolis (Norops) crassulus species subgroup, which contains ten morphologically similar highland taxa, the majority of which have restricted distributions. The nominal taxon A. crassulus has a disjunct distribution from Chiapas, Mexico, through Guatemala, in the highlands of El Salvador, and in the Chortís Highlands of Honduras. We test the relationships of these species using multiple mitochondrial and nuclear loci in concatenated and multispecies coalescent frameworks, in an effort to both resolve long-standing taxonomic confusion and present new insights into the evolution and biogeography of these taxa. Results Sequences of multiple mitochondrial and nuclear loci were generated for eight of the ten species of the Anolis crassulus species subgroup. We analyzed phylogenetic relationships and estimated divergence times and ancestral ranges of the subgroup, recovering a monophyletic subgroup within Anolis. Within the nominal taxon Anolis crassulus, we recovered multiple genetically distinct lineages corresponding to allopatric populations, and show that the Chortís Highland lineage split from the others over 13 MYA. Additionally, distinct mitochondrial lineages are present within the taxa A. heteropholidotus and A. morazani, and importantly, samples of A. crassulus and A. sminthus previously used in major anole phylogenetic analyses are not recovered as conspecific with those taxa. We infer a Chortís Highland origin for the ancestor of this subgroup, and estimate cladogenesis of this subgroup began approximately 22 MYA. Conclusions Our results provide new insights into the evolution, biogeography, and timing of diversification of the Anolis crassulus species subgroup. The disjunctly distributed Anolis crassulus sensu lato represents several morphologically conserved, molecularly distinct anoles, and several other species in the subgroup contain multiple isolated lineages. Electronic supplementary material The online version of this article (doi: 10.1186/s12862-017-1115-8) contains supplementary material, which is available to authorized users.


Background
Anoline lizards are a well-known focal system for studying mechanisms and drivers of evolution [1,2]. Studies of the replicated adaptive radiations of anoles on Caribbean islands [3] have greatly expanded our understanding of evolution and biogeographic patterns (e.g. [4][5][6]). In contrast to their island counterparts, mainland anoles have been relatively understudied, despite hosting the majority of the diversity of these lizards [7]. Recent studies of South American anoles have provided a clearer understanding of the relationships and biogeography of those taxa (e.g. [8][9][10][11]). Central America has been a setting for biogeographic studies of a variety of taxa, including snakes (e.g. [12,13]), mammals (e.g. [14]), fish (e.g. [15]), and amphibians (e.g. [16]). Only recently, however, have studies begun to focus on the evolutionary history, diversity, and biogeography of Central American anoles, revealing a not wholly unexpected degree of complexity in patterns relating to these lizards (e.g. [17,18]). It is apparent that the diversity of mainland anoles has been greatly underestimated, and that comprehensive work remains across virtually all mainland taxa.
One particularly poorly resolved group of mainland anoles is the Anolis (Norops) crassulus subgroup, consisting of ten nominal species known from intermediate-tohigh elevations (1000-2500 m above sea level) in Central America: A. amplisquamosus, A. anisolepis, A. crassulus, A. haguei, A. heteropholidotus, A. morazani, A. muralla, A. rubribarbaris, A. sminthus, and A. wermuthi (sensu [19]). Anoles referred and bearing an affinity to A. crassulus have long been recognized as needing thorough resolution of their evolutionary relationships and taxonomy [19][20][21][22]. Seven of the ten species are endemic to the Chortís Block Highlands of Honduras and adjacent El Salvador and Nicaragua [19,23,24]. Anolis anisolepis and A. haguei are known only from Chiapas, Mexico, and Alta Verapaz, Guatemala, respectively [19,23], and are of questionable taxonomic validity [25]. Populations assigned to the nominal taxon A. crassulus can be found from southwestern Honduras to Chiapas, Mexico at elevations between 1300 and 3000 m [19,23]. However, A. crassulus has a disjunct distribution across Central America, with isolated populations ranging in Chiapas, Mexico, southern El Salvador, across central Guatemala, and western Honduras (Fig. 1) [23]. These species all inhabit a similar ecological niche largely restricted to Lower Montane Wet and Moist Forest formations, and are found active on the ground and low vegetation or tree trunks during the day and asleep on low vegetation at night [19]. This subgroup has long been a taxonomically confusing complex, and past studies using only morphological examinations have failed to resolve relationships of named taxa and populations. Meyer and Wilson ([20]: 108) stated "specimens of the crassulus group from Guatemala and Mexico have a bewildering array of admixtures of the distinctive characters observed in Honduras… The inter-relationships of [this subgroup] are exceedingly complex, and… we are unable to suggest a satisfactory arrangement." Subsequent authors, including Lieb [21], McCranie et al. [22], and McCranie and Köhler [19] have commented on the importance of a deeper investigation into this subgroup. Previously, these anoles been divided into 'crassulus-like' anoles (Anolis anisolepis, A. crassulus, A. haguei, A. morazani, and A. rubribarbaris) and 'sminthus-like' anoles (A. heteropholidotus, A. muralla, A. sminthus, and A. wermuthi) on the basis of morphological similarities [24,26]. As noted by Townsend and Wilson [24], the relationship of Anolis amplisquamosus to the A. crassulus subgroup remains poorly understood. The validity of A. anisolepis and A.
heteropholidotus have been questioned, and in some cases, these taxa are still not recognized. Köhler [27] demonstrated the validity of A. heteropholidotus, which has subsequently been recognized by Townsend and Wilson [24], Johnson et al. [28], McCranie and Köhler [19], and Poe et al. [29]. However, Lieb [25] stated that A. heteropholidotus was not distinct at the species level, Nicholson et al. [30] did not list A. heteropholidotus as a valid species, and both the IUCN [31] and Uetz et al. [32] continue to regard A. heteropholidotus as a synonym for A. sminthus. Anolis anisolepis was considered a valid species by Nicholson et al. [30], but not by Poe et al. [29].
Despite the poor understanding of the relationships within the Anolis crassulus species subgroup and recent increased sampling effort, the majority of previous anole phylogenies have included, at most, two samples of the Anolis crassulus subgroup: a sample considered to represent A. crassulus from Chiapas, Mexico, and another sample identified as A. sminthus from Olancho, Honduras. As this study was being finalized, two major anole phylogenies became available that expanded the sampling of this group: Nicholson et al. [18] and Poe et al. [29]. Nicholson et al. [18] used a continuous 1423 bp sequence of mitochondrial DNA to analyze the biogeography of mainland Anolis. They included a single sample of A. crassulus, A. amplisquamosus, A. wermuthi, A. heteropholidotus, A. morazani, and A. sminthus, and recovered them as a monophyletic subgroup within their larger "A. crassulus clade". The A. crassulus specimen, MZFC 6458 from Chiapas, Mexico, was the same sequence used in Nicholson et al. [7,18], and sampled from the same specimen sequenced by Gray et al. [33] under a different GenBank accession number. The A. sminthus sample, SMF 78830 from Olancho, Honduras, was the same sequence used in Nicholson et al. [7,18]. McCranie & Köhler [19] noted that that SMF 78830 is not actually a specimen of A. sminthus, and opined that it likely represents an undescribed species more closely related to A. crassulus. Poe et al. [29] included molecular samples (at least two mitochondrial and/or nuclear loci) of A. crassulus, A. amplisquamosus, and A. sminthus, and scored A. heteropholidotus, A. morazani, A. muralla, A. rubribarbaris, and A. wermuthi for morphological characters only (no molecular data). In their combined phylogeny, Poe et al. [29] did not recover a monophyletic A. crassulus subgroup, and the placement of the taxa scored only for morphology was not well supported.
The goals of this study were to resolve the evolutionary histories of taxa and populations assigned to the Anolis crassulus subgroup, and determine the subgroup's biogeographic origins and patterns of diversification in Central America. We assembled and analysed a multilocus phylogenetic dataset (three mitochondrial loci and three nuclear loci) of the majority of ingroup taxa (including multiple samples from across their distributions, when available), and estimated divergence times and ancestral distributions to provide a historical context for their diversification. We provide new insights into the evolutionary history of these taxa in an attempt to resolve some of the confusion that has plagued this group for decades.

Taxon sampling
Most samples used in this study were collected across Honduras and northern Nicaragua by JHT and colleagues between 2005 and 2015 ( Fig. 1). Prior to formalin-fixing specimens, a tissue sample was removed and preserved in SED buffer (20% DMSO, 0.25 M EDTA, pH 7.5, NaCl saturated; [34][35][36] [37]. A list of all Fig. 1 Sampling localities of the Anolis crassulus species subgroup: A. crassulus sensu lato (top) and the remaining subgroup species (bottom). The type localities of two subgroup taxa, A. haguei and A. muralla, are also noted for completeness, though they were not able to be included in this study samples and their localities can be found in Additional file 1: Appendix 1.

DNA sequence generation
We generated sequences of three mtDNA loci and three nDNA loci for eight of the ten ingroup species (tissues of Anolis haguei and A. muralla were not available) and three outgroup taxa (A. cusuco, A. kreutzi, and A. laeviventris from Honduras; the A. laeviventris species subgroup sensu [19]). DNA extraction, amplification, and sequencing of mitochondrial genes 16S large subunit rRNA (16S) and cytochrome oxidase subunit I (COI) of some samples were carried out at the Smithsonian Institution Laboratory of Analytical Biology (LAB; Suitland, Maryland, USA) following standardized DNA Barcode of Life (BOLD) protocols. Template DNA was extracted by an AutoGen Geneprep 965 (Autogen, Holiston, MA) via a phenol-chloroform technique. The two mitochondrial genes were amplified via polymerase chain reaction (PCR). Two μl of ExoSAP-IT (Affymetrix, Santa Clara, CA) was used to remove unincorporated nucleotides from the PCR product, which was sequenced using a BigDye Terminator v3.1 Cycle Sequencing kit (ThermoFisher, Waltham, MA), cleaned using Sephadex spin column filtration, and electrophoresed on an ABI 3730xl DNA Analyzer.
All other sequences used herein, including additional loci generated for samples sequenced for only 16S and COI at LAB, were generated in the Townsend Lab at Indiana University of Pennsylvania (Indiana, PA, USA). Whole-genome DNA was extracted from tissue using PureLink Genomic DNA Kits (ThermoFisher). All available samples were first amplified for 16S and COI via PCR. A subsample of these was then amplified for an additional mitochondrial locus, NADH dehydrogenase subunit 2 (ND2), and three nuclear loci: prolactin receptor (PRLR), brain-derived neutrophic factor (BDNF), and protein tyrosine phosphatase non-receptor type 12 (PTPN12). PRLR and PTPN12 are considered relatively variable nuclear markers in squamate reptiles, while BDNF is considered relatively conserved [38].
Unincorporated nucleotides were removed from all PCR products using 2 μL of ExoSAP-IT (Affymetrix) per sample, and products were sent to Eurofins Scientific (Louisville, Kentucky, USA) for sequencing. They were electrophoresed on an ABI 3730xl DNA Analyzer. Chromatograms were checked manually and sequences were assembled using Geneious v.7.1.7 [39]. Sequence alignment was carried out in MEGA7 [40] using the Clus-talW algorithm [41], and redundant identical haplotypes were removed to streamline analyses. Transfer RNAs (tRNAs) at the ends of the ND2 sequences were removed due to inconsistencies in the sequencing of those regions. All sequences were deposited on GenBank (Additional file 1: Appendix 1). Additional ND2 sequences were pulled from GenBank (Additional file 1: Appendix 3) and analyzed with the ND2 sequences generated in this study.
Four datasets were then created from the full sequence data: ND2 only (75 samples), all mtDNA (16S + ND2 + COI; 43 samples), all nDNA (PRLR + BDNF + PTPN12; 37 samples), and all six loci concatenated (37 samples). The ND2 dataset included samples from GenBank and was analysed separately to determine if the Anolis crassulus subgroup is monophyletic with respect to other anoles. An additional ND2 dataset (19 samples) was used for divergence dating analysis (see below).

Phylogenetic analyses
In all analyses, unless otherwise stated, the datasets were partitioned by gene for 16S and codon position for the five protein-coding loci (ND2, COI, PRLR, BDNF, PTPN12). Best fit models of nucleotide substitution (Additional file 1: Appendix 4) were selected for each partition using PartitionFinder v.1.1.1 [42] via the Akaike information criterion (AIC), a greedy search scheme, and limiting the substitution models to those that could be implemented in MrBayes 3.5.2 [43,44]. The ND2only datasets (one including a larger sample of all Anolis and one including the A. chlorocyanus group; see below) were analysed separately to account for the additional sequences from GenBank.
To account for potential incongruences between mitochondrial-only and nuclear-only datasets due to incomplete lineage sorting, we inferred the species tree of our dataset under a Bayesian multispecies coalescent framework in StarBEAST2 [49], implemented in BEAST 2.4.7 [50]. Use of the multispecies coalescent framework have been shown to produce accurate estimates of species trees from independent gene trees, even when gene trees are not congruent [51,52]. This analysis requires a priori assignment of individuals to "species", which was done based on the results of our prior analyses. Two individuals per species were included when possible, though we included only samples that had at least two loci sequenced and some taxa had only single exemplars; in total, 35 individuals were used. We did not phase our nuclear data as many outgroup and ingroup taxa had only single representative individuals. We partitioned our dataset by locus, unlinking substitution models across each (based on PartitionFinder v1.1.1. ran as above, except limiting the models to those implementable in BEAST; Additional file 1: Appendix 4), and linking molecular clock and tree priors in the three mitochondrial loci while leaving them unlinked in the three nuclear genes. Four independent analyses were performed using a strict clock and Yule tree prior, for 100 million generations, sampling every 5000. We examined the resulting logfiles in Tracer v1.6 [53] to ensure effective sample sizes (ESS) were above 200, and combined trees and log files in LogCombiner v.2.4.7 [50], removing a burn-in of 10%. Finally, the combined tree files were annotated with TreeAnnotator v2.4.7 [50], using mean heights to calculate the maximum clade credibility (MCC) tree.

Divergence dating
We estimated the timing of diversification of the Anolis crassulus subgroup in BEAST 1.7.5 [54]. We analyzed the ND2 sequences from one representative of each ingroup lineage where possible (13 samples; Additional file 1: Appendix 1), with an uncorrelated lognormal relaxed clock model and a Yule tree prior. Initially, we chose only the calibration rate of 0.65% per lineage per million years (MY) [55]. This rate was calculated from approximately 1700 bp of ND1, ND2, and COI, including transfer RNAs [55], and has been used in several phylogenetic studies of anoles without fossil calibrations to determine the approximate timing of various taxa (e.g. [8,11,17,56]). However, as our ND2 sequences did not cover the entirety of Macey et al. [55]'s sequences, including the conserved trailing tRNAs, the 0.65%/MY rate alone greatly overestimated divergence dates in comparison to recent phylogenies [8,29,30,57]. Analyses using the published rate in addition to a fossil calibration (see below) continued to produce dates older than other recent phylogenies. Application of a calibration rate of 0.765% per lineage per MY (approximately 15% faster, accounting for the lack of conserved tRNAs in our sequences) estimated dates more comparable (though not identical) to those of previous studies [8,29,30,57], and was chosen for continued analysis together with a fossil calibration in order to provide a general temporal context for the diversification of this subgroup. As there are no known fossils within the A. crassulus subgroup, we added four ND2 sequences of the A. chlorocyanus species group (A. alinger, A. chlorocyanus, A. coelestinus, and A. singularis; Additional file 1: Appendix 1) from GenBank to our dataset, allowing one fossil calibration utilized in major anole phylogenies [29,30] to be incorporated: a Dominican amber anole assigned to this species group [58]. Following Poe et al. [29], we placed this calibration on the most recent common ancestor node of the A. chlorocyanus species group, and used a uniform prior distribution based on the stratigraphic information of the fossil (17-23 MYA).
We implemented this dating method in order to determine a general historical context for the diversification of the subgroup and report the resulting values cautiously, given they are estimated from a single mitochondrial locus with an approximated rate, and a single fossil calibration. Four independent runs were performed in BEAST with the same number of generations and frequency of sampling as the species tree analysis. As with the species tree analysis, the resulting log files were examined in Tracer v1.6 [53], log and tree files were combined in LogCombiner v.1.7.5 [54] and combined tree files were annotated with TreeAnnotator v1.7.5 [54].

Ancestral area reconstruction
To infer the ancestral distribution of the Anolis crassulus subgroup, we conducted ancestral area reconstruction using a Bayesian Binary Markov-chain Monte Carlo (BBM) [59] analysis implemented in RASP v3.2 [60] using the post burn-in trees obtained from BEAST to account for phylogenetic uncertainty. This method statistically infers an ancestral state at a node by averaging its frequency across all trees, using a full hierarchical Bayesian approach and hypothesizing a special distribution, the "null distribution", where an ancestral area does not contain any of the unit areas [60]. We assigned each ingroup taxon to one of four biogeographical regions: (A) Chiapas, Mexico; (B) Guatemala; (C) Salvadoran Cordillera; (D) Chortís Highlands. These boundaries isolate all of the ingroup taxa analyzed herein, and represent biogeographic breaks corresponding to speciation in reptiles (e.g. [12,13]). Outgroup taxa were removed from the analysis. The maximum potential areas occupied was set to four, the Fixed Jukes-Cantor model (with gamma shape parameter) was chosen, and the analysis was run for 1 million cycles using 10 chains, sampling every 100 generations and discarding the first 10% as burn-in.

Sampling success
In total, 107 samples of ingroup taxa were sequenced for 16S (444 bp) and 104 samples for COI (654 bp). This group was then subsampled to representatives of the recovered lineages, but due to problems with PCR amplification and/or gene sequencing, some samples failed to properly amplify for certain genes. Twenty-three samples of ingroup taxa were amplified for ND2 (1038 bp), 39 for PRLR (545 bp), 31 for BDNF (658 bp), and 25 for PTPN12 (832 bp), giving a final combined dataset of 4171 bp (2136 bp mtDNA, 2035 bp nDNA).

Phylogenetic relationships
All phylogenetic analyses recovered a monophyletic Anolis crassulus species subgroup, with multiple lineages assigned to some nominal taxa. The previous assignment of species to "crassulus-like" and "sminthus-like" series within the subgroup was not supported. Instead, in most analyses we recovered a clade containing A. anisolepis and A. crassulus sensu lato (hereafter the A. crassulus clade) as sister to a clade that contained most of the Chortís Highland endemics: A. sminthus, A. morazani sensu lato, A. rubribarbaris, and A. heteropholidotus sensu lato (hereafter the A. sminthus clade). Anolis aff. rubribarbaris and A. wermuthi were also recovered as part of this clade when included in mitochondrial only and species tree analyses. The placement of A. amplisquamosus was problematic and poorly supported across all analyses.
Analyses of the larger ND2 dataset (Fig. 2) recovered the Anolis crassulus subgroup as a well-resolved clade within the rest of Anolis. It is important to note that MZFC 6458, the sample used to represent A. crassulus in previously published anole phylogenies, was not recovered with any other samples of A. crassulus, including those from nearby localities in Chiapas, Mexico, and was instead recovered as a sister lineage to A. anisolepis. Additionally, SMF 78830, the sample of A. "sminthus" referred to as a potentially unrecognized species related to A. crassulus by McCranie & Köhler [19], was recovered as part of a lineage with other "crassulus"-like anoles from Olancho that is sister to A. morazani. Based on this information, we infer that the samples referred to herein as A. aff. morazani are conspecific with McCranie and Köhler [19]'s "crassulus-like" population from Olancho. Further investigation of this population is ongoing. Unfortunately, MZFC 6458 and SMF 78830 were not available to be sequenced for the remaining loci.
Analyses of mitochondrial-only and nuclear-only datasets were partially incongruent (Fig. 3). While both datasets recovered similar structure, the nuclear only dataset recovered a polytomy between the A. anisolepis, A. crassulus sensu lato, and the proposed A. sminthus clade. The placement of Anolis amplisquamosus, A. rubribarbaris, and A. sminthus were not congruent between analyses. Both datasets did recover multiple lineages assigned to A. crassulus, corresponding to geographic location, but the El Salvadoran population was nested within the Guatemalan samples in nuclear-only analyses. Interestingly, two samples previously identified as A. crassulus were recovered as a sister lineage to A. rubribarbaris (herein called A. aff. rubribarbaris). We were not able to include samples of A. aff. rubribarbaris or A. wermuthi in the remainder of the analyses due to a lack of viable tissue for sequencing. Four distinct lineages were recovered within A. heteropholidotus in mitochondrial-only analyses; this structure was not recovered in nuclear only analyses, which returned only a single lineage.
The species tree from the multispecies coalescent framework (Fig. 4) and the six-locus ML and BI topologies (Additional file 2: Figure S1) were largely congruent, supporting many of the relationships recovered from the other datasets analysed. Two clades within the Anolis crassulus subgroup were strongly supported (maximum likelihood bootstrap support (BS) = 100; Bayesian posterior probabilities (PP) = 1; species tree PP = 1). The substructure in A. crassulus was again recovered, with the Chortís population well-supported as basal to the other A. crassulus lineages (BS = 92; PP = 1/0.9), and the Mexican population strongly supported as distinct from the El Salvadoran/Guatemalan populations (BS = 100; PP = 1/1). As with earlier analyses, the placement of A. amplisquamosus differed between concatenated analyses (basal to the A. sminthus clade) and the species tree analyses (most closely related to A. sminthus and A. wermuthi), but was not well supported in either (BS = 43; PP = 0.52/0.50). Anolis aff. morazani was well-supported as a sister lineage to A. morazani (BS = 100; PP = 1/0.99). Four lineages were recovered within A. heteropholidotus, all well supported (BS = 100 for all splits except lineages 1 & 2, where BS = 83; PP > 0.99/0.87). The species tree analyses recovered A. wermuthi as sister to A. sminthus, and A. aff. rubribarbaris as sister to A. rubribarbaris. As only sequences of two mitochondrial loci were available for A. wermuthi and A. aff. rubribarbaris, they were not included in the concatenated, six-loci analyses. The addition of missing subgroup taxa (Anolis haguei, A. muralla, A. aff. rubribarbaris, and A. wermuthi) in the concatenated six-gene phylogeny may alleviate some of the incongruences apparent between the different analyses, and provide better resolution of the placement of some ambiguously placed taxa.
Within the A. crassulus lineages, average overall uncorrected pairwise (p) distances (of all six concatenated loci; Additional file 1: Appendix 5) ranged from 0.000-0.007, while between lineages, the average distances ranged from 0.032-0.071; between the Honduran population and the other populations of A. crassulus, these values ranged from 0.066-0.071. These inter-lineage distances are comparable to those between recognized taxa such as A. rubribarbaris and A. morazani   The divergence date estimates recovered here are similar to the larger squamate phylogeny of Zheng and Wiens [57] at two major nodes: the MRCA of the Anolis chlorocyanus outgroup and our ingroup taxa (43.1 compared to 43.6 MYA) and the A. crassulus and A. sminthus clades (22.3 MYA compared to 24.3 MYA). We recovered younger estimates within the A. chlorocyanus group: 20.9 MYA compared to 25.9 for the MRCA of A. coelestinus and the rest of the subgroup; 11.9 compared to 15.3 for the MRCA of A. chlorocyanus, A. aliniger, and A. singularis; and 6.5 compared to 8.8 for the MRCA of A. aliniger and A. singularis. Our estimates were comparable to Poe et al. [29], but not Nicholson et al. [18] (discussed below).

Ancestral area reconstruction
The BBM analysis (Fig. 5) inferred the ancestral distributions of ingroup ancestors, though these results were not entirely unambiguous. The most recent common ancestor (MRCA) of the Anolis crassulus subgroup was inferred to be in the Chortís Highlands (probability = 86.4%). As expected, the MRCAs of the Anolis sminthus clade (all Chortís Highland endemics) were reconstructed as the Chortís Highlands (all nodes >99%). The ancestral areas of the A. crassulus clade, on the other hand, were reconstructed with less certainty.
The MRCA of A. anisolepis and A. crassulus sensu lato was inferred to be in either the Chortís Highlands (54.2%) or in Chiapas, Mexico (28.9%), as was the MRCA of A. crassulus sensu lato (Chortís only: 52.5%; Mexico: 29.6%). The ancestor of the remaining A. crassulus lineages was likely found in Chiapas (82.3%), while the ancestor of El Salvadoran and Guatemalan populations was inferred to be in either Guatemala (37.4%) or the Salvadoran Cordillera (36.3%).

Discussion
Our results indicate that significant evolutionary diversity within the Anolis crassulus subgroup has been overlooked. Congruent results from mtDNA, nDNA, concatenated, and species tree analyses revealed at least three deep genetic lineages presently assigned to A. crassulus sensu lato, including one lineage in the Chortís Block, A. crassulus sensu stricto in Guatemala and El Salvador, and a third lineage in Chiapas, Mexico. Importantly, the sample of A. crassulus from Chiapas, Mexico used in many recent anole phylogenies [18,29,30,33] is not conspecific with A. crassulus sensu stricto , and instead represents a sister lineage to the Chiapan endemic A. anisolepis. Mitochondrial p-distances between Chortís Block A. crassulus and other populations of the species are notable in comparison to those of some recently described anoles, such as A. purpuronectes (ND2 p = 0.115 from its sister taxon A. barkeri; [33]), A. kathydayae (ND2 p = 0.125 from its congener A. brooksi; [61]), A. mccraniei and A. wilsoni (16S p = 0.032-0.044 from their congener A. tropidonotus; [62]), and A. elcopeensis (COI p ≥ 0.073 "from other included Anolis species"; [63]:4). Similar phylogeographic structure to Anolis crassulus sensu lato is apparent in the other taxa. Chortís Highland populations of the widespread pitviper Cerrophidion godmani, isolated from the Guatemalan populations by the same biogeographic boundaries that separate the lineages of A. crassulus sensu lato, were revealed to represent a distinct lineage and species (C. wilsoni; [64]). Similar structure is also apparent in montane mice of the Peromyscus mexicanus species complex, where numerous lineages are present throughout Chiapas, Guatemala, and western Honduras, corresponding to similar biogeographic breaks and highland localities [65].
Our BBM results suggest that the common ancestor of the Anolis crassulus subgroup was found in the Chortís Highlands, and the ancestor of non-Chortís A. crassulus lineages radiated north (Fig. 5). The BBM analysis suggested the possibility of a Mexico + Chortís Block distribution for the ancestor of A. crassulus and A. anisolepis lineages. While this seems unlikely given the modern placement of these land masses, it is plausible in light of currently accepted geologic hypotheses regarding the Chortís Block (as summarized in [66]). During the early Paleocene, the Chortís Block was located to the south of what is now mainland Mexico, disconnected from the rest of Central America [67]. Throughout the Paleocene, Eocene, and Oligocene epochs (between 65 and 20 MYA), the Chortís Block moved eastward along the southern margin of the North American Plate along the contemporary Pacific Coast of Mexico ( [68]). Given the divergence time estimates and BBM results, we hypothesize the A. crassulus ancestor was distributed in the Chortís Block and subsequently invaded Mexico before further diversifying. It is also possible it was distributed in Mexico and the Chortís Block simultaneously during a period of connectivity between the two bodies, before the continuous eastward shift of the Chortís Block isolated the two areas. During the Miocene, intense volcanic activity in the region that is now southwestern Honduras, El Salvador, and southeastern Guatemala, likely isolated the populations [69,70]. By the late Miocene, the Mexican lineage had subsequently invaded Guatemala and El Salvador, while the Chortís population remained isolated in the highlands of western Honduras. Rapid diversification also occurred in the Chortís Highlands within the A. sminthus clade, between 3.5-6.3 MYA.
Our results provide a clearer picture of the evolution of the Anolis crassulus subgroup than was available in either Nicholson et al. [18] or Poe et al. [29], due in large part to more complete molecular sampling of ingroup taxa. Our recovery of Nicholson et al. [18]'s "A. sminthus" (SMF 78830) as a divergent sister lineage to A. morazani sensu stricto makes sense in light of their phylogeny, which present the two as sister taxa. Our phylogenetic hypothesis is otherwise largely congruent with Nicholson et al. [18]. However, the timings of diversification presented by Nicholson et al. [18] are much older than those recovered here or by Poe et al. [29]. For example, Nicholson et al. [18] estimated the A. crassulus clade/A. sminthus clade divergence to have occurred at approximately 33 MYA, nearly 10 MY earlier than our estimate, and 17 MY earlier than the estimate of Poe et al. [29]. Their estimates for the A. morazani/A. aff. morazani and the A. heteropholidotus/A. morazani splits were between 7 and 10 MY earlier than our estimates. Our ancestral area reconstruction suggesting a Chortís origin and subsequent dispersal out for the A. crassulus subgroup is consistent with their results for the "A. crassulus clade" (consisting of more taxa than studied herein). However, they recovered a larger potential for a Mexican origin of A. crassulus, due to the fact their single sample of A. crassulus is from Chiapas, and not conspecific with A. crassulus.
Our phylogenetic hypothesis demonstrates numerous incongruences with those of Poe et al. [29], though our divergence time estimates were comparable. In an analysis that combined molecular-only, molecular + morphological, and morphological-only data into a single phylogeny, Poe et al. [29] presented a paraphyletic A. crassulus subgroup, recovering A. crassulus, A. sminthus, and A. amplisquamosus in a clade with A. petersii, distinct from any of the taxa they only scored for morphology. These were recovered in a different clade sister to A. ortonii and A. sulcifrons, along with some members of the A. laeviventris subgroup. These relationships were not well supported (posterior probabilities as low as 0.11), likely due to their use of morphology only for several taxa. Poe et al. [29] estimated the divergence between A. crassulus and A. sminthus clades to be 16 MYA, approximately 6.3 MY later than our estimate (22.3 MYA). Poe et al. [29] did not estimate other divergence times within our ingroup, but their estimates within the A. chlorocyanus group were similar to those recovered in this study, ranging from 2.5-3.9 MY later than our estimates. Given that we did not include the more conserved trailing tRNAs, and had no ingroup fossil calibrations, the dates we present for the A. crassulus subgroup may be slightly older than the actual divergence times in the group, and Poe et al. [29]'s estimate of the A. crassulus/A. sminthus clade split may be demonstrated to be more accurate. Until further calibrations within this subgroup can be implemented, however, our dates provide a general temporal context for the diversification of these taxa not previously sampled in any phylogeny.
While Anolis crassulus lineages are genetically distinct and have been isolated for millions of years, these allopatric lineages are highly conserved morphologically. Thorough morphological investigation into these lineages is underway, but initial examinations by the authors have failed to find consistent, discrete differences in morphological characters examined, including a lack of diagnosable differences in hemipenis or dewlap morphology. This is consistent with the results of past attempts to resolve these taxa (e.g. [20]) using morphology alone. Each lineage inhabits similar intermediate-to-high elevation Lower Montane Wet and Lower Montane Moist Forest formations [19]. This apparent lack of niche diversification has potentially limited selection that would lead to significant morphological differences, leaving these lineages with highly conserved morphologies. A thorough ecological study including niche modeling could find subtle but consistent differences between lineages, as much of the natural history of these animals remains unknown.
Based on the deep genetic differences between ingroup lineages, disjunct geographic distributions, the timing of diversification and potential distributions of the ancestor populations, and the apparent phenotypic and ecological stasis among populations, we hypothesize that the Anolis crassulus subgroup, and in particular A. crassulus sensu lato, represent a nonadaptive radiation of highland anoles. In contrast to adaptive radiations (such as those of the Caribbean anoles), where a single ancestor evolves to fill a variety of niche space, a nonadaptive radiation refers to "evolutionary diversification from a single clade, not accompanied by relevant niche differentiation" ( [75]: 264). Nonadaptive radiations have been not been well documented in reptiles, but some studies have noted similar levels of diversity among cryptic lineages. Barley et al. [76,77] found deep genetic diversification not accompanied by distinct morphological differentiation within cryptic lineages of the Asian skink genus Eutropis, which diverged between 1 and 10 MYA. A nonadaptive radiation has also been hypothesized as driving diversity in the South American lizard genus Phymaturus [78]. Within the Anolis crassulus subgroup, the majority of taxa and populations are allopatric, and have diverged within the last 15 MY. The subsequent genetic drift has left a series of genetically distinct, morphologically conserved, ecologically similar anoles across the Central American highlands.
It is clear that the Anolis crassulus subgroup is in need of a thorough taxonomic revision, as has been the case for some time [19][20][21][22]. Our results strongly support the validity of Anolis anisolepis and A. heteropholidotus, two taxa that previously have been questioned or inconsistently recognized. As for the rest of the subgroup, however, at this time we propose no other taxonomic changes until additional lines of evidence are fully examined. Continued morphological investigation into the Chortís populations of A. crassulus is underway, but given the distinct molecular substructure consistently recovered across all analyses, this population likely warrants species-level recognition. The unique Mexican lineages remain in need of further study. These populations are not as geographically isolated as Chortís populations relative to A. crassulus in Guatemala and El Salvador, suggesting additional modes of diversification within the group. Additional work is also underway to investigate populations assigned to A. aff. morazani and A. heteropholidotus.

Conclusions
Based on analyses of the first multilocus dataset of the majority of the Anolis crassulus subgroup, we provide a new hypothesis for the evolutionary relationships, timing of diversification, and ancestral ranges of these anoles. Our results show clear molecular structure for numerous cryptic lineages within this group, and lay the groundwork for future investigation and potential taxonomic revisions. Additionally, our results support the validity of A. anisolepis and A. heteropholidotus. Based on these results, we hypothesize that this species complex represents a nonadaptive radiation of highland anoles. However, further investigation and studies of morphology, behavior, ecology, and diet, need to be undertaken to test this hypothesis.
By building upon this work with better taxonomic sampling and additional molecular, morphological, and ecological data, new questions about the diversity and biogeographic histories of mainland subgroups can be explored. Taken with other recent studies [8,11,17,61], it is clear that the diversity in Mesoamerican anoles has been underestimated, and further investigation is needed to reveal the patterns of evolutionary diversification among these ubiquitous lizards.

Additional files
Additional file 1: Appendix 1. Localities, samples, and GenBank Accession numbers for sequences used in this study. Taxa are identified by the genetic lineage they were recovered as in the phylogenetic analyses. Appendix 2. Primers and cycling parameters used in this study.