Biogeography and evolution of a widespread Central American lizard species complex: Norops humilis, (Squamata: Dactyloidae)

Background Caribbean anole lizards (Dactyloidae) have frequently been used as models to study questions regarding biogeography and adaptive radiations, but the evolutionary history of Central American anoles (particularly those of the genus Norops) has not been well studied. Previous work has hypothesized a north-to-south dispersal pattern of Central American Norops, but no studies have examined dispersal within any Norops lineages. Here we test two major hypotheses for the dispersal of the N. humilis/quaggulus complex (defined herein, forming a subset within Savage and Guyer’s N. humilis group). Results Specimens of the N. humilis group were collected in Central America, from eastern Mexico to the Canal Zone of Panama. Major nodes were dated for comparison to the geologic history of Central America, and ancestral ranges were estimated for the N. humilis/quaggulus complex to test hypothesized dispersal patterns. These lineages displayed a northward dispersal pattern. We also demonstrate that the N. humilis/quaggulus complex consists of a series of highly differentiated mitochondrial lineages, with more conserved nuclear evolution. The paraphyly of the N. humilis species group is confirmed. A spatial analysis of molecular variance suggests that current populations are genetically distinct from one another, with limited mitochondrial gene flow occurring among sites. Conclusions The observed south-to-north colonization route within the Norops humilis/quaggulus complex represents the first evidence of a Norops lineage colonizing in a south-to-north pattern, (opposite to the previously held hypothesis for mainland Norops). One previously described taxon (N. quaggulus) was nested within N. humilis, demonstrating the paraphyly of this species; while our analyses also reject the monophyly of the Norops humilis species group (sensu Savage and Guyer), with N. tropidonotus, N. uniformis, and N. marsupialis being distantly related to/highly divergent from the N. humilis/quaggulus complex. Our work sheds light on mainland anole biogeography and past dispersal events, providing a pattern to test against other groups of mainland anoles. Electronic supplementary material The online version of this article (doi:10.1186/s12862-015-0391-4) contains supplementary material, which is available to authorized users.

Anoles, present throughout Central America [10,30], are an ideal group to remedy this deficiency. There are several species groups in the genus Norops that are widespread, and present excellent opportunities to study the colonization of Mesoamerica. Among these taxa, we selected the Norops humilis species group as a model for refining biogeographic hypotheses for Central American lizards in general, and Central American anoles in particular. The N. humilis species group (described below) ranges from Mexico to Panama, although the timing of diversification within this group has not been determined. Our first objective was to estimate the date of origin for the group, which may in part provide support for one of two hypotheses. Hypothesis 1 (H1; see Figure 1) predicts a north-to-south dispersal pattern, while Hypothesis 2 (H2; see Figure 1) predicts a south-to-north dispersal pattern. Our second objective was to test biogeographic patterns within this group, as interpreted from the above hypotheses. The N. humilis species group [31] is an ideal group for examining these hypothesized patterns as well as their timing. As originally defined, the group includes N. compressicauda, N. humilis, N. notopholis, N. tropidonotus, and N. uniformis. To these we add the newly described N.
wampuensis [32] and N. quaggulus [33], plus N. marsupialis, a former subspecies of N. humilis that has recently been treated as a species [34]. Several studies suggest the N. humilis species group may be paraphyletic [35][36][37]. Therefore, our third objective was to confirm the polyphyly of the N. humilis species group. We also examined the N. humilis/quaggulus complex, which we define to include N. humilis plus N. quaggulus, a sister species described based on putatively unique hemipenial morphology and evidence of reciprocal monophyly based on mitochondrial genome sequences [33]. This complex also has a wide distribution (Panama to Honduras). Our analyses provide the opportunity to further evaluate the relationship between N. humilis and N.
quaggulus, as well as to test the specific status of N. marsupialis.

Methods
Samples from several of the species in the Norops humilis species group (N. humilis/quaggulus, N. marsupialis, N. tropidonotus and N. uniformis) were collected throughout Central America or acquired from tissue loans (Figure 2, S1). For the N. humilis/quaggulus complex, 147 specimens were sampled throughout Honduras, Nicaragua, Costa Rica, and Panama. DNA was extracted from liver or muscle tissues using Qiagen DNeasy kits (Qiagen, USA). PCR was conducted following protocol and using lizard-specific primers from Macey et al. [38] for mtDNA and Nicholson [35] for nucDNA. Purified PCR reactions were sent to Michigan State University's Research and Technology Support Facility for sequencing of the following gene regions: NADH-ubiquinone oxidoreductase chain 2 (ND2), tRNA Trp , tRNA Ala , tRNA Asn , tRNA Cys , tRNA Tyr , origin of light strand replication, and partial CO1 from mitochondrial DNA. A subset of these samples representing each major lineage detected by the mitochondrial analyses was selected for nuclear DNA analysis using a nuclear internal 5 transcribed spacer unit (ITS-1). A nuclear marker was included because many studies have shown the limitations of mtDNA to reflect levels of gene flow or the extent of reproductive isolation among populations (e.g. [39,40]; see [41] for a review). Analyses based solely on mtDNA can also provide results that are in conflict with the rest of the nuclear genome [42,43].
Within the N. humilis species group, a total of 1451 aligned bp of mitochondrial data were collected for 192 individuals and 1522 aligned bp of the nuclear gene region ITS-1 was collected for 48 individuals. All newly acquired data were combined with published sequences for 65 additional Norops species (S2) in order to investigate the monophyly of the N. humilis species group and of the N. humilis/quaggulus complex. Sequences were edited using Sequencher 4.9 (GeneCodes Corp., Ann Arbor, MI, USA) and aligned initially using MUSCLE in MEGA 5.2.2 [44], then adjusted manually.
The relatively continuous geographic distribution of the N. humilis/quaggulus complex combined with a lack of distinct phenotypic differences made a priori separation of sampling localities into populations somewhat arbitrary. In order to describe the genetic structure and identify the best maximally differentiated number of populations within the N.
humilis/quaggulus complex, a spatial analysis of molecular variance (SAMOVA 1.0) was used to group 15 localities selected in this study (each with n ≥ 4 specimens; S3) into a number of user-defined clusters (K). For each cluster, the proportion of total genetic variance (high F CT index) due to differences between populations [45] was estimated and evaluated to select the optimal number of genetic groups. A simulated annealing process for each cluster (K = 2 to 14) was repeated 1023 times for each of 100 sets of initial conditions to ensure that the final population groups were not affected by the initial configuration. Significance of the F SC index was used to obtain the suggested number of genetic groupings for the localities selected [46]. 6 This analysis was based on sequences for the mitochondrial region only, since ITS-1 data were only obtained for a limited number of specimens, and were not sufficient for population genetic analysis.
Phylogenetic estimations were conducted under Bayesian analytical methods. PARTITIONFINDER [47] was used to select models of evolution as well as to examine the suitability of partitioning each dataset (mitochondrial, nuclear, and concatenated). In all cases each gene (including tRNAs) was entered as a potential partition. Protein coding genes were further partitioned by codon position for the PartitionFinder analysis. Branch lengths were unlinked, all models of evolution available in MRBAYES 3.2.2 [48] were tested, and a BIC information criterion and greedy algorithm were used. The PARTITIONFINDER analysis recommended a HKY + I + model of evolution for the mtDNA segment, and GTR + for the nuclear data with no partitioning recommended within either region. For the combined dataset, two partitions were recommended (mt and nuc) with GTR + I + as the selected model of evolution for both partitions.
Bayesian analyses for each dataset (as above) were conducted using MRBAYES and BEAST 1.7.5 [49]. The phylogenetic hypotheses using MRBAYES 3.2.2 was developed using 20 million generations, sampling every 1000 generations for two independent runs using four Markov chains with node support evaluated via posterior probabilities (BAPP). We evaluated stationarity of variables by examining our output via TRACER v1.5 [50]. The first 20% of trees were discarded as burnin and a majority rule consensus tree was generated to summarize the post-burnin results. With the phylogeny constructed using BEAST, we estimated divergence dates, employing a lognormal relaxed clock and a calibration rate of 0.65% per million years for mtDNA [38] with a Yule Process speciation prior. This rate has been used for our selected mitochondrial gene region for reptile and amphibian groups [51,52] including anoles (e.g. [53,54]; mean of a prior distribution, SD = 0.0025 for ucld.mean parameter). The calibration rate was only applied to the mtDNA analysis, as there are no calibration rates available for the nuclear ITS-1 gene. In addition to Bayesian analyses, ML analyses were conducted on each dataset using MEGA 5.2.2 with node support evaluated via bootstrap analyses (MLBS) based on 2000 replications. In addition to the constructed topologies, pairwise genetic distance (uncorrected-p) was estimated between all individuals included in the phylogenetic analyses using both the mitochondrial and concatenated datasets. These distances were compared among all major lineages (as indicated by phylogenetic analysis) of N. humilis species group members, as well as within each species or lineage using MEGA 5.2.2.
Phylogeographic hypotheses were tested using likelihood-based inference in LAGRANGE [55,56]. The analysis was conducting using the phylogeny constructed with BEAST for mtDNA only, as divergence dating was not available for the nuclear (ITS) tree. Dispersal events between regions were examined and ancestral ranges were estimated within the N. humilis/quaggulus complex to evaluate the support for each of the hypotheses presented above ( Figure 1). The Nicaragua) with younger clades in the south, it would indicate a north-to-south dispersal, while a southern origin would support a south-to-north colonization route.

Results
All phylogenetic hypotheses demonstrate the paraphyly of the N. humilis species complex (Figures 3-6). Each phylogenetic analysis is discussed in detail below in regard to the N.
humilis/quaggulus complex. Six main lineages were present when using the concatenated dataset (mt + nucDNA) for both Bayesian ( Nicaragua specimens was younger than any clades restricted to Costa Rican or Panamanian samples in both ML and Bayesian analyses. Coupled with our ancestral range estimation, the topology suggests that the N. humilis/quaggulus complex originated in the south before dispersing northwards.

Discussion
To address our first objective, we estimated the age of the Norops humilis/quaggulus species complex and evaluated two potential dispersal patterns for this group. The BEAST analysis yielded a mean crown group age of 17.2 Myr BP (range = 14.2-20.6) for the clade. We report this date cautiously because it was calculated from mitochondrial data only. To evaluate the second objective, we used the ancestral range estimation to evaluate dispersal patterns. Using the LAGRANGE output, we infer that the common ancestor of the Norops humilis/quaggulus complex originated in Panama before dispersing west to Costa Rica, and then north into  As in many other taxa ( [58] and sources within), LCA has served as a region of diversification for the Norops humilis/quaggulus complex (once again, distinct from the N. humilis species group). The northward distribution observed here is similar to that found in eleutherodactyline (genera Craugastor [19]; Pristimantis [59]) and hylid frogs (Dendropsophus [60]) that originated in South America and dispersed to Central America prior to the most recent completion of the Isthmus of Panama (3)(4). This is also similar to findings in other groups of squamates [2], although it remains to be seen if this pattern of dispersal is shared by other lizard species.
Our third objective was to examine polyphyly described by others [35,36] of the N. humilis species group. The lack of support for a monophyletic N. humilis species group is consistent with results from Nicholson [35] and Poe [36], further illustrating the extent of convergence in morphological characters of mainland anoles. All analyses agreed in assigning several species of the N. humilis species group into separate areas of the tree. Norops marsupialis, N. tropidonotus, and N. uniformis were placed in areas of the tree distant to a monophyletic N.
humilis/quaggulus complex. Initially N. marsupialis was included in this study as a member of the N. humilis/quaggulus complex. Taylor [61] described N. marsupialis as a subspecies of N. humilis in 1956, and it was not elevated to specific status until 2015 [34]. All of our analyses (Bayesian and ML for each gene region and combined) found support for the specific status of N. marsupialis, and agreed to its placement as sister to a clade containing N. aquaticus and N. woodi (approx. 30 Myr divergent from the N. humilis/quaggulus complex, Figure 6). Therefore, we confirm N. marsupialis as a distinct species morphologically similar to, but evolutionary distant from N. humilis/quaggulus, agreeing with Köhler et al. [34]. Morphological convergence within anoles has long been reported, particularly for Caribbean species [27,[62][63][64][65][66] and is 13 further demonstrated by our analysis of the N. humilis species group. Inclusion of the other three members of the group (N. compressicauda, N. notopholis and N. wampuensis) into the molecular phylogeny of Norops may yield further insight on their placement in the phylogeny, but their inclusion is not likely to significantly alter the results obtained here. Additional work within N. tropidonotus may also be necessary to examine the evolution that has occurred in that species, as it occupies a broad geographic range (Köhler, 2008). Our data suggest that N.
tropidonotus contains at least three deep mitochondrial divergences, and could potentially represent multiple cryptic lineages. While N. quaggulus is nested within N. humilis, we do not recommend synonymizing the two species until further work is done to investigate this complex. The deep divergences seen among the N. humilis clades may correspond to cryptic species, and if further work confirms the presence of two or more species, the name N.
quaggulus would have priority as the senior synonym.
The SAMOVA results indicate that limited mitochondrial gene flow is occurring among localities, suggesting 12 distinct genetic groupings among our sample localities. Therefore, we conclude that genetic differentiation within the N. humilis/quaggulus complex is significant enough to conclude that (1) population fragmentation has occurred and (2) the complex does not represent a single panmictic population. The isolation, coupled with the high genetic variance, supports all major lineages identified in the Bayesian analyses as being distinct from one another, as well as further subdivision within most lineages. While we do not suggest that the 12 groups correspond to separate species, the presence of at least six well supported, divergent clades within the N. humilis/quaggulus complex demonstrates that deep mitochondrial divergence has occurred within this group, although nuclear evolution appears to be more conserved. This result is similar to several studies on Caribbean anoles, which also display considerable mitochondrial differentiation, with much less diversity in the nuclear DNA [67][68][69][70]. This pronounced phylogeographic substructuring may be explained in part, by low vagility in lizards [9]. Deep mitochondrial divergences coupled with conserved nuclear evolution as seen here, may have at least two implications 1) ITS is more slowly evolving than the mitochondrial genome of anoles, which we consider to be a likely scenario given that nuclear DNA is generally regarded as experiencing slower rates of evolution [71][72][73] and 2). The mitochondrial-nuclear relationships observed here are characteristic of male-biased dispersal [74,75], indicating that female N. humilis are more philopatric than males as seen in Caribbean anoles [76,77].
The novel biogeographic pattern for Central American anoles revealed here illustrates a need for further work on mainland Norops. What remains to be tested is whether the south-tonorth dispersal route seen in the N. humilis/quaggulus complex is repeated in other Norops groups. In addition, it is important to clarify where the Central American Norops lineage originated, how it dispersed throughout the mainland, and when these events took place.
Investigating the phylogeography of other widespread anoles may be highly informative towards understanding other distribution patterns of the Central American herpetofauna. There are several widespread Norops species and species complexes that vary in their ecological roles with corresponding morphological features. These are grouped into designations called ecomodes, which are distinct from ecomorphs, to accommodate mainland anoles (see [35]).
Such an assortment of ecologically diverse anoles may provide good models for testing the biogeographic hypotheses discussed here. Further studies on mainland Norops species are needed to test if the cryptic diversity suggested here is present in other widespread species complexes within the genus.  humilis are separated out by clade as defined by the combined Bayesian analysis (Figure 3).
Distances for mtDNA data are located below the diagonal, while distances for the combined nuclear and mtDNA data are located above the diagonal. Cells are colored by species.       whereas in other figures are in regards to lineages. The tree used for this analysis was created in BEAST v1.7.5 and only contains mitochondrial data since divergence dating was not possible for ITS-1. All major nodes are strongly supported (see Figure 6). Appendix S1. Locality and sequence data for Norops humilis species group sensu Savage and Guyer [31]. ITS accession numbers are not applicable to all specimens, as many taxa were only sequenced for the mtDNA region.
Appendix S2. Genbank numbers for sequences used in this analysis from 65 Norops clade outgroups. mtDNA were from multiple studies [78][79][80] All ITS sequences are from Nicholson et al. [80]. Not all taxa used had available ITS sequences.  The authors declare that they have no competing interests.