A genomic assessment of species boundaries and hybridization in a group of highly polymorphic anoles (distichus species complex)

Abstract Delimiting young species is one of the great challenges of systematic biology, particularly when the species in question exhibit little morphological divergence. Anolis distichus, a trunk anole with more than a dozen subspecies that are defined primarily by dewlap color, may actually represent several independent evolutionary lineages. To test this, we utilized amplified fragment length polymorphisms (AFLP) genome scans and genetic clustering analyses in conjunction with a coalescent‐based species delimitation method. We examined a geographically widespread set of samples and two heavily sampled hybrid zones. We find that genetic divergence is associated with a major biogeographic barrier, the Hispaniolan paleo‐island boundary, but not with dewlap color. Additionally, we find support for hypotheses regarding colonization of two Hispaniolan satellite islands and the Bahamas from mainland Hispaniola. Our results show that A. distichus is composed of seven distinct evolutionary lineages still experiencing a limited degree of gene flow. We suggest that A. distichus merits taxonomic revision, but that dewlap color cannot be relied upon as the primary diagnostic character.


| INTRODUCTION
The formation of new species is typically a gradual process that occurs over thousands or even millions of generations. As this makes speciation difficult to observe experimentally, investigating how and why speciation occurs tends to rely heavily on observations of species at varying stages of the speciation process-the snapshot approach to studying speciation. Relatively young species are particularly important but also the hardest to identify because they often fail to meet one or more of the criteria expected of deeply divergent species. They may, for example, exhibit incomplete reproductive isolation, readily hybridize with other species, or be difficult to distinguish morphologically or genetically (Coyne & Orr, 2004;Knowles & Carstens, 2007;Maddison & Knowles, 2006;Shaffer & Thomson, 2007). Speciation is a continuum under the general lineage concept, and criteria are expected to accumulate gradually (de Queiroz, 2007).
Our goal here is to use genomic data to identify candidate species within a polytypic lizard species that may include a number of young lineages at varying stages of divergence (Geneva, Hilton, Noll & Glor, 2015;Glor & Laport, 2012;Ng & Glor, 2011). The Hispanolian bark anole (Anolis distichus) is a trunk-dwelling lizard species that currently includes more than a dozen subspecies distributed across Hispaniola and the Bahamas (Schwartz, 1968). These subspecies are primarily delimited by differences in the color and pattern of their dewlaps, throatfans that are extended by males during behavioral displays (Schwartz, 1968). Dewlaps in A. distichus range from pale yellow to wine red, with many variants in between; most dewlap color and pattern variation occurs among geographically circumscribed populations, but considerable variation can also exist within some populations (Lambert, Geneva, Mahler & Glor, 2013;Schwartz, 1968). Because the dewlap is thought to play a critical role in species recognition and sexual selection, dewlap divergence has been used as a proxy for reproductive isolation and is often used to delimit species boundaries (e.g. Lotzkat, Bienentreu, Hertz & Köhler, 2011;Poe & Yañez-Miranda, 2008;Velasco & Hurtado-Gómez, 2014). However, in the case of A. distichus and a few other polymorphic anole species, populations with strikingly different dewlaps have been recognized as subspecies or unnamed geographic populations rather than distinct species because they appear to hybridize where they come into contact (Heatwole, 1976;Schwartz, 1968;Schwartz, 1974;Underwood & Williams, 1959). This decision to recognize dewlap color variation at the subspecific level (or not at all) is supported by more recent evidence that dewlap color and pattern variation in A. d. distichus may represent an adaptive response to local signaling conditions rather than an indicator of reproductive isolation (Ng, Kelly, MacGuigan & Glor, 2013;Ng, Landeen, Logsdon & Glor, 2012;Webster, 1977).
Prior molecular genetic studies of A. distichus provided mixed support for the evolutionary independence of the subspecies diagnosed by differences in dewlap color and pattern. Early allozyme studies revealed molecular differentiation and reduced gene flow at the contact zone between some subspecies (Case & Williams, 1984) but not others (Case, 1990;Case & Williams, 1984;Williams, 1977;Williams & Case, 1986). Meanwhile, mitochondrial DNA (mtDNA) sequence data suggested that each of the subspecies found in the Dominican Republic form distinct and deeply divergent clades, with the exception of the widespread subspecies A. d. dominicensis, which is associated with multiple mtDNA clades (Glor & Laport, 2012). Fine-scale studies of contact zones between pairs of subspecies involving phenotypic, mitochondrial, and microsatellite data have uncovered evidence for abrupt phenotypic and genetic divergence along narrow hybrid zones, but also evidence for extensive introgression and relatively shallow genetic differentiation (Ng & Glor, 2011;Ng, Ossip-Klein & Glor, 2016).
Multilocus phylogenetic analyses have found that while most subspecies of A. distichus are genetically distinct, these differences were mostly restricted to mtDNA, and several subspecies were not monophyletic (Geneva et al., 2015).
The multilocus phylogeny of Geneva et al. (2015) also suggested for the first time that geography may be more important than dewlap color variation for delimitation of Anolis species.
Modern Hispaniola formed when a North and a South paleo-island merged approximately 15 mya (Graham, 2003;Iturralde-Vinent & MacPhee, 1999). The boundary between these paleo-islands, also known as Mertens' Line, has long been recognized as one of the most important biogeographic boundaries on Hispaniola (Schwartz, 1980). The current boundary between the paleo-islands (Figure 1, black dashed line) has likely remained a biogeographic barrier since the merger because it coincides with a lowlying xeric valley that is periodically inundated with seawater, and is relatively inhospitable to lizards adapted to the more mesic environments flanking the valley (Glor & Warren, 2011;Townsend, Rimmer, Latta & Lovette, 2007). Anolis distichus populations appear to have diverged across Mertens' line, with the deepest phylogenetic split dividing clades of subspecies found primarily on either the North or the South paleo-island (Geneva et al., 2015).
In spite of this prior work, no study of A. distichus has involved range-wide assessment of genomic variation with the goal of identifying candidate species. Utilizing amplified fragment length polymorphism (AFLP) genome scans, we apply a two-step process of candidate species discovery and validation (Carstens, Pelletier, Reid & Satler, 2013 (Ng & Glor, 2011;Ng et al., 2016; (Ng & Glor, 2011;Ng et al., 2016;  We extracted DNA from tail tips or liver samples stored in 95% ethanol at −80°C using either a Wizard SV Genomic DNA Purification System kit (Promega Corp.) or via a phenol chloroform extraction protocol modified from Laird et al. (1991). For phenol chloroform extractions, we combined up to 20 ng of tissue with 250 μl of TENSII (base solution 4 ml 5 mol/L NaCl, 50ml 1 mol/L Tris pH 8, 2 ml 0.5 mol/L EDTA pH 8, 844 ml H 2 O, 100 ml 10% SDS), 20 μl of proteinase K (20 μg/μl), and 5μl RNase A solution before incubating for 16-18 hr at 55°C. Following incubation, we transferred this solution to a prespun (15,000 g for 1-2 min) Phase Lock Gel TM (PLG) 2 ml heavy tube (5 Prime, Inc), added 0.5 ml of phenol:chloroform:isoamyl alcohol (PCI, 25:24:1), and mixed via repeated inversion. We then centrifuged at 14,000 g in an Eppendorf model 5,424 microcentrifuge for 5 min before transferring the resultant aqueous phase to a fresh prespun PLG 2 ml tube heavy tube. We next added 0.5 ml of chloroform:isoamyl alcohol (CI, 24:1) to the sample in the PLG2 ml tube, mixed by repeated inversion, and centrifuged the tube at 14,000 g for 5 min before transferring the resultant aqueous phase to a fresh microcentrifuge tube. We then added 30 ml of sodium acetate 3 mol/L, pH 5.2) and 1.25 ml of 95% ethanol before incubating at −20°C overnight. Finally, we centrifuged this mixture at 14,000 g for 20 min, rinsed with 1 ml of 95% ethanol, centrifuged again at 14,000 g for 10 min, and ultimately re-suspended the resulting DNA pellet in 200 μl H 2 O.
F I G U R E 1 Distributions of A. distichus subspecies, sampling for Set 1, and results from genotypic clustering analyses conducted in STRUCTURE. Each column on the bar plots represents an individual sample. Different colors correspond to different genetic clusters. Shading of each column represents the proportion of the genome for that individual assigned to one of the genetic clusters identified by STRUCTURE. Each point on the map is a locality included in Set 1, labeled with corresponding locality numbers. The color of each locality reflects the genetic cluster to which the majority of the genomes at that locality were assigned. Locality 150 is colored gray because the genome of the specimen at that locality was not assigned to a single genetic cluster. Localities 686, 687, 818, and 819 are colored light orange to represent their admixed status. The colored regions on the map represent approximate subspecies ranges, with white space where no subspecies is present. Subspecies ranges are based on the maps of Ng et al. (2013) and Schwartz (1968

| AFLP genotyping
AFLPs can provide a large amount of genomic data to address questions about population structure and hybridization (Bensch & Akesson, 2005;Mueller & Wolfenbarger, 1999) and are more costeffective than other methods for acquiring genomic data (e.g., GBS or similar SNP-based approaches). Because AFLPs are dominant markers, they generally suffer from reduced information content relative to sequencing or SNP-based approaches (Elshire et al., 2011). AFLPs may also suffer greatly from genotyping error (Crawford, Koscinski & Keyghobadi, 2012), but these errors are accounted for by the methods outlined below.
AFLP genotyping involved four steps: (1) digestion of genomic DNA and ligation of adaptors, (2) preselective PCR amplification of genomic DNA, (3) selective PCR amplification with fluorescently labeled primers, and (4) scoring of fragments resulting from selective amplification.
We followed Lambert et al. (2013) in using AFLP protocols modified from Vos et al. (1995). 500 bp. Samples that failed the preselective amplification were rerun until successful or excluded from the final dataset.
Following the preselective amplification, we performed selective PCR amplification on products of the preselective amplification using primers identical to the preselective primers, but with the addition of two nucleotides at the 3' end. We used a total of six primers for selective In the selective amplification step, we utilized two EcoR1 florescent primers (one labeled with VIC and one labeled with 6-FAM) and three Mse1 primers for a total of total of six unique primer pair combinations.
For our first set of samples, we used a fourth Mse1 primer for a total of eight unique primer pairs to increase the number of loci (Table 1).

All fragment analyses were performed by the Functional Genomics
Center at the University of Rochester Medical Center using an Applied Biosystems 3730 Genetic Analyzer with a LIZ500 size standard.

| AFLP scoring and error analysis
We individually analyzed every primer pair for each set of samples, as well as for a concatenated dataset containing all samples. For the combined dataset, only the first six primer pairs were used. We first visually inspected and analyzed AFLP electropherograms using PeakScanner v1.0 (Applied Biosystems) with light peak smoothing and default settings. We analyzed results from PeakScanner using a modified version of the R package RawGeno (Arrigo, Tuszynski, Ehrich, Gerdes & Alvarez, 2009). For analyses of individual sets, we set the maximum bin width to two base pairs (bp), minimum bin width to one bp, minimum fragment size to 50 bp, and maximum fragment size to the observed maximum fragment length. To remove low-intensity peaks, we set the minimum peak height threshold to 100 relative florescence units (RFU). For analyses of the concatenated dataset, we used the same settings except for maximum fragment size, which we set equal to the smallest of the observed maximum fragment sizes from the individually analyzed sets.
We used the "visualize samples" tool in RawGeno (Arrigo et al., 2009) to identify samples that had fewer AFLP peaks than expected, potentially indicating a methodological or analytical failure. The "visualize samples" tool produces a binary matrix of AFLP loci, with an AFLP peak either present or absent at a particular fragment size ( Figure S1).
If a sample failed in any of the marker preparation steps, it would have very few strong AFLP peaks and the "visualize samples" tool would call most loci as "absent." For a given sample, we identified primer pairs as problematic if they had few "present" loci compared to the rest of the samples. We then removed all problematic primer pairs for that sample from the final dataset. If this procedure resulted in diagnosis of three or more problematic primer pairs for a particular sample, that sample was completely removed from the dataset prior to downstream genetic clustering and species delimitation analyses.
We exported raw peak height output data from RawGeno and used a custom R script to convert this output into a format accepted by the R package AFLPScore (Whitlock, Hipperson, Mannarelli, Butlin & Burke, 2008). We used ALFPScore to assess the mismatch genotyping error rate for duplicate samples (Whitlock et al., 2008). We duplicated 16 samples, representing 12.2% of our dataset and exceeding the recommended 5-10% (Bonin et al., 2004). To generate duplicate samples, we repeated selective amplification for eight samples from Set 1, 14 samples from Set 2, and eight samples from Set 3 (30 total samples). All duplicate samples were randomly selected. We scored AFLPs for these duplicate samples following the same protocol discussed above. We analyzed error rates for each primer pair at phenotype threshold values of 100, 250, 500, 750, and 1,000 (the minimum RFU required to call the phenotype at a specific locus as present) and at locus threshold values of 100, 250, 500, and 750 (the minimum RFU required to call a peak as a locus). We selected and applied the least strict threshold values that produced an error rate <0.05 as suggested by Zhang and Hare (2012). Our acceptable error rate is consistent with the 6-18% error rates reported for optimized phylogenetic resolution (Holland, Clarke & Meudt, 2008) and only slightly higher than the 2-5% error rates reported in several other studies using semi-automated AFLP scoring (Bonin et al., 2004).

| Species delimitation & species tree inference
We used two methods to infer boundaries between candidate species from the AFLP data acquired for Set 1, which included broad geographic and taxonomic sampling. The first method was largely exploratory and relied on the clustering algorithms implemented in the program STRUCTURE (Pritchard, Stephens & Donnelly, 2000) to ask whether some populations or sets of populations correspond with genotypic clusters that may represent distinct species. The second method used Bayes factors and a coalescent-based framework to quantitatively evaluate and compare a set of alternative species delimitation scenarios derived a priori from taxonomy, biogeography, or the genotypic clusters identified by STRUCTURE (Leaché, Fujita, Minin & Bouckaert, 2014). In addition to identifying species boundaries, we also used our AFLP data to infer a species tree that reflects evolutionary relationships between candidate species.

| Genotypic clustering analyses
To determine whether genotypically distinct populations exist in the A. distichus species group, we conducted genotypic clustering analyses with our Set 1 AFLP data using the program STRUCTURE (Pritchard et al., 2000). STRUCTURE uses a Bayesian Markov chain Monte Carlo (MCMC) algorithm to probabilistically assign individuals to genetic clusters (Pritchard et al., 2000). We also performed the same STRUCTURE analyses independently on the AFLP data from Set 2 and Set 3 to examine the degree of admixture along our two hybrid zone transects.
Strong confounding effects prevented us from combining all three sample sets for one STRUCTURE analysis (see Section 3). For each STRUCTURE analysis, we utilized the admixture and correlated allele frequency models. Because AFLP markers are dominant, we also employed the recessive alleles model. We ran each analysis for 1,000,000 generations, excluding 100,000 generation as burn-in. We then used the ΔK method in Structure Harvester to determine the optimal number of genetic clusters (K) (Earl & vonHoldt, 2012;Evanno, Regnaut & Goudet, 2005). The ΔK method examines the rate of change in lnP(D), where lnP(D) is an estimate of the posterior probability for K genetic clusters. The optimal number of genetic clusters is identified as the breakpoint where the slope of lnP(D) vs. K begins to plateau (Evanno et al., 2005).
The ΔK method alone can underestimate the actual number of genetic clusters in a dataset (Evanno et al., 2005), so we employed a hierarchical version of the ΔK method (Coulon et al., 2008). We assigned individuals to genetic clusters based on majority (>0.5) inferred ancestry. Individuals that did not have the majority of their genotypes assigned to one cluster were excluded from subsequent hierarchical analyses, but were included in the final nonhierarchical analyses at a fixed value of K. Individuals were divided into subsets based on their majority cluster assignment and subjected to another round of STRUCTURE analyses using the ΔK method. For each round of analyses, we performed 10 independent replicate STRUCTURE runs at each value of K, with K values ranging from 1 to the maximum number of subspecies included in the dataset plus two. We set the maximum K to a value larger than the number of subspecies represented in each analysis to avoid forcing individuals into inappropriately few clusters (Kalinowski, 2011).
Adapting the general guidelines outlined by Coulon et al. (2008), we employed this hierarchical ΔK method until K = 1 had the highest posterior probability, <5 individuals were assigned to a genetic cluster, or no individuals had the majority of their genotypes assigned to any cluster. When no further population subdivision was possible, we calculated the total number of clusters identified during the hierarchical analyses. To generate our final clustering results, we performed 100 replicate STRUCTURE runs on the complete Set 1 dataset with K fixed at the total number of clusters identified by the hierarchical analyses.

| Bayes factor delimitation
We tested alternative species delimitation scenarios using the coalescent-based model comparison framework outlined by Leaché et al. (2014). We generated eight species delimitation scenarios that consisted of between two and thirteen species based on (1) current taxonomy, (2) biogeography, and (3) genotypic clusters identified by STRUCTURE (Figure 2). The three models based on the current taxonomy were (I) two species corresponding to the two currently  (Geneva et al., 2015). The two biogeographic models were (IV) three species, corresponding with Hispaniola's   (Kass & Raftery, 1995). Thus, our BF model selection statistics indicated the degree of support for the best fitting model relative to each alternative model. A BF model selection statistic between 0 and 2 reflects weak support, between 2 and 6 reflects positive support, between 6 and 10 reflects strong support, and greater than 10 reflects decisive support (Kass & Raftery, 1995).

| Species tree inference
To estimate a species tree for our optimal species delimitation model, we used SNAPP with the Set 1 AFLP data. We employed default priors for mutation rates and ancestral population sizes, estimating both the forward and reverse mutation rates from the data. We used a gamma prior on the ancestral population sizes, with a shape parameter (α) of 11.75 and a scale parameter (β) of 109.73. We used a Yule prior for the species tree height and branch lengths with a lambda parameter of 0.00765. We ran two independent MCMC chains for 1 × 10 6 generations, with parameters and trees sampled every 1,000 generations. To assess MCMC mixing and convergence, we visualized the output using Tracer v.1.6 (Rambaut, Suchard, Xie & Drummond, 2014

| Interactions between candidate species at areas of contact
Our delimitation of candidate A. distichus species was restricted to analyses of Set 1, which included broad geographic and taxonomic sampling. To test whether these candidate species are currently experiencing gene flow, we conducted separate analyses of Set 2 and Set 3, each of which contain samples from across contact zones between candidate species that have traditionally been recognized as subspecies.
For Set 2 and Set 3, we tested for the presence of distinct genotypic clusters corresponding with candidate species by conducting the same type of hierarchical STRUCTURE analyses used for the Set 1 analyses. By analyzing genotypic assignment proportions across the hybrid zone, we determined whether hybridization is ongoing, as well as the extent to which admixture is evident outside of the contact zone.

| Genetic diversity and pairwise F ST calculation
We calculated genetic diversity (H e ) within and pairwise F ST among each of the species identified in the optimal species delimitation model for Set 1 (see below) and for the genotypic clusters identified by independent STRUCTURE runs for Sets 2 and 3 using AFLP-Surv v.1.0 (Vekemans, Beauwens, Lemaire & Roldan-Ruiz, 2002). We ran 5,000 permutations of the Bayesian method with a nonuniform prior for allele frequencies (Zhivotovsky, 1999) to estimate F ST under the assumption of Hardy-Weinberg equilibrium. For all three datasets, we calculated these statistics both with and without individuals from localities at known hybrid zones.

| Error rates and quality control
We determined AFLP scoring error rates for 20 locus/phenotype threshold combinations for each primer pair. Ultimately, we only used phenotype thresholds of >500 or greater because lower thresholds tended to result in noninterpretable results in downstream analyses, indicative of low-quality data. Results were considered noninterpretable when STRUCTURE failed to assign more than 50% of the genomes of most individuals to any cluster. We were unable to determine primer pair specific error rates for sample Set 1 for two primer pairs (M53/E1 and M53/E2) due to a technical change at the core facility conducting our AFLP fragment analyses. As a result of this change, duplicate samples run with these two primer pairs had significantly higher AFLP peaks relative to the original samples, making comparison impossible. Thus, for all M53 primer pairs, we applied the filtering threshold most commonly used for all other primer pairs: a phenotype threshold of 500 and a locus threshold of 100.
Presence of an unusually low number of AFLP peaks resulted in complete exclusion of 10, 15, and 8 individuals from our three sample sets, resulting in 82, 77, and 51 retained individuals in Sets 1, 2, and 3, respectively. Due to low AFLP peaks, we excluded one primer pair from 7, 3, and 1 individual(s) and two primer pairs from 8, 1, and 0 individuals in sets 1, 2, and 3, respectively. Following exclusion of these primer pairs, the three sets were 97.6%, 99.2%, and 99.7% complete, respectively. The total number of loci retained was as follows: Set 1 included 534 loci, Set 2 included 552 loci, and Set 3 included 836 loci. North paleo-island.

| Genotypic clustering analyses
The first hierarchical STRUCTURE analysis of the South paleo-island cluster did not result in any further subdivision (lnP(D) greatest for K = 1, Figure S2a). The first hierarchical analyses of the North paleo-island cluster suggested additional subdivision, with the optimal K = 2.
However, the two genotypic clusters identified by this analysis were not easily interpretable, with many individuals from the same subspecies and locality assigned to different clusters. Additionally, most individuals' genotypes were not strongly assigned to any one cluster (average of maximum genotype assignment proportions = 58.5%).
Closer inspection of other values of K from the hierarchical analysis of North paleo-island cluster revealed a more readily interpretable pattern at K = 4, the value of K with the highest overall lnP(D). For K = 4, most individuals had the majority of their genotype assigned to a single cluster (average of maximum genotype assignment proportions = 83.8%, Figure S2a). After recovering a total of five genotypic clusters (one in the South paleo-island and four in the North paleo-island) via hierarchical analyses, we ran 100 replicate STRUCTURE runs on Set 1 with K fixed at 5. However, this final analysis produced slightly different clustering than was suggested by the hierarchical analyses (Figures 1 and S2b).
These differences included (1)  We evaluated the fit of the these two clustering schemes separately (models VI and VII) and in combination (model VIII) using BFD*.

| Bayes factor delimitation
We compared eight different species delimitation scenarios under a coalescent-based framework (Figure 2 All of the populations assigned to this candidate species tend to have relatively pale dewlaps. One of the Bahamian subspecies included in this group (A. d. ocior) has the most green body coloration of any distichoid population (Schwartz, 1968).

| Species tree inference
We used the BEAST package SNAPP to infer phylogenetic relationships among the candidate species identified by model VIII (Figure 4).

| Population structure statistics for candidate species
Pairwise F ST and H e values for Set 1 are reported in Table 2

| Interactions between candidate species at areas of contact
Our second set of samples consisted of 552 AFLP loci for 77 A. d. ravitergum and A. d. ignigularis from a transect that spans a zone of contact between the two subspecies. The first round of ΔK analyses with this dataset identified an optimal K = 2, with clusters corresponding largely with subspecies ( Figure 3). We also calculated pairwise F ST and H e values for the clusters identified by independent STRUCTURE analyses of sets 2 and 3. Pairwise

| DISCUSSION
Using genotypic clustering and species delimitation methods, we recover strong support for the hypothesis that A. distichus is comprised   (Lambert et al., 2013). This observation supports the hypothesis that divergence within populations currently recognized as A. distichus is younger than that observed between the four distinct species previously recognized as A. brevirostris (Arnold, 1980).

| Biogeography
Although pure biogeographic scenarios were among the worst performing delimitation models (Figure 2), our results support prior hypotheses (Geneva et al., 2015;Glor & Laport, 2012) that suggest divergence of populations on Hispaniola's North and South paleoislands has contributed to diversification in bark anoles (Figures 1   and S2). The first division in our hierarchical STRUCTURE analyses distinguishes populations found primarily on Hispaniola's North and South paleo-islands. Our analyses are unable to determine whether this divergence across the paleo-island boundary occurred prior to the paleo-island merger or from restricted gene flow since the merger due periodic inundation of the paleo-island boundary or the inhospitable environmental conditions of this region (Glor & Warren, 2011

| Dewlap color in species delimitation
The historic use of dewlap color as the primary taxonomic character in the A. distichus complex has led to recognition of many subspecies that may not reflect true evolutionary lineages. We identified several candidate species that contain a broad array of dewlap colors.
For instance, A. d. favillarum appears to be a single genetic population with impressive dewlap color polymorphism, consistent with prior phylogenetic (Geneva et al., 2015) and allozyme studies (Case, 1990;Williams & Case, 1986 ation, yet make up a single genetic cluster (Figure 1). This "genetic continuity" of the three Tiburon subspecies was previously hinted at by the unfinished allozyme work of Webster in the 1970s (Williams, 1977). On the other hand, at least one previously delimited subspecies with similar dewlap color across its range appears to represent multiple independent evolutionary lineages; populations of A. d. dominicensis were split across four separate candidate species, in agreement with prior phylogenetic results (Geneva et al., 2015). Together these results suggest that dewlap color is not by itself a reliable diagnostic trait in the A. distichus complex, and perhaps in anoles more broadly.
Other polymorphic anoles may also be composed of multiple genetically divergent species, which implies that the biodiversity of anoles is currently underestimated. Future studies should explicitly quantify both dewlap color variation and genetic variation to determine whether other anole species exhibit a similar disassociation between dewlap color and population structure (e.g. Ng et al., 2016).

| Hybridization and introgression
We examined two A. distichus subspecies pairs for evidence of hybridization at contact zones. A. d. ignigularis (Figure 4, candidate species B) and A. d. ravitergum (Figure 4, part of candidate species A) come into contact in the southern Dominican Republic along the Baní River. This contact zone, first described by Williams (1977)  the south to mesic habitat in the north (Ng et al., 2016). Our genotypic clustering analyses reveal a strong signal of admixture in the middle of this transect with very little admixture at either end ( Figure 3). This pattern is indicative of hybridization without substantial gene flow into the home range of either subspecies (Ng et al., 2016). Despite low pairwise F ST estimates between the two subspecies, we conclude that there is a strong genetic break between A. d. ignigularis and A. d. ravitergum, with admixture at the hybrid zone but limited gene flow between the subspecies.
The second subspecies pair we examined was A. d. ignigularis ( Figure 3, candidate species B) and A. d. dominicensis (divided among four candidate species). Our transect for these subspecies runs eastwest, spanning a recently recessed marine channel that separated the Samaná Peninsula from mainland Hispaniola (Grant, 1956;Ng et al., 2016). Unlike the transect between A. d. ravitergum and A. d. ignigularis, this transect does not encompass any obvious environmental gradient. While there is signal of admixture in the middle of this transect, the hybrid zone is not as well defined as the hybrid zone between A. d. ravitergum and A. d. ignigularis (Ng et al., 2016). There appears to be significant admixture well into the range of A. d. dominicensis at the western edge of the transect. Thus, we conclude that A. d. ignigularis genetic material has effectively introgressed into A. d. dominicensis beyond the contact zone. However, without further geographic sampling of A. d. dominicensis populations in eastern Hispaniola, it is difficult to determine the extent of this gene flow.
Previous phylogenetic analyses found that A. d. dominicensis consists of three or four geographically distinct and deeply divergent polyphyletic lineages (Geneva et al., 2015). Their species tree analyses recovered a clade of northern Haitian/central Dominican A. d. dominicensis and a separate clade of northern Dominican A. d. dominicensis whose most recent common ancestor was that of all A. distichus (Geneva et al., 2015). This deep divergence within A. d. dominicensis is reflected in our own analyses, with a distinct genetic break between populations from northern Haiti/central Dominican Republic and populations from the northeastern Dominican Republic ( Figure   3). Comparatively, A. d. ignigularis located on mainland Hispaniola shows little genetic differentiation across its entire range, from the Samaná Peninsula to the southeastern Dominican Republic. Thus, despite the relative uniformity of dewlap color, there is genetic evidence for multiple independent lineages within the widespread A. d. dominicensis.

| TAXONOMIC RECOMMENDATIONS
Our results together with prior work strongly suggest that formal taxonomic revision of populations previously recognized as A. distichus is needed because this species is comprised of numerous distinct populations that likely warrant recognition as distinct species. We have not undertaken such a taxonomic revision here because we are unable to provide diagnostic phenotypic traits to distinguish the candidate species identified on the basis of genomic differentiation. Additionally, delimiting the geographic boundaries between these putative species requires more extensive geographic sampling of genomic variation.
The fact that A. altavalensis is genetically indistinguishable from A. d. ravitergum could be used to argue in favor of no longer recognizing the Alto Velo populations of bark anoles as a distinct species.
However, we agree with the suggestion by Geneva et al. (2015) that A. altavalensis warrants continued recognition because it is clearly geographically isolated and phenotypical distinct from A. d. ravitergum.

| CONCLUSIONS
Our study provides a geographically broad first-take genomic perspective on a young species complex of anoles with remarkable dewlap color polymorphism. Consistent with results from mitochondrial DNA (Geneva et al., 2015;Glor & Laport, 2012) and several nuclear genes (Geneva et al., 2015), we find strong evidence for genetic differentiation despite some gene flow between the lineages of Anolis distichus. We identify six new candidate species with our molecular species delimitation and suggest that A. altavalensis should be maintained as a seventh species. The genetic breaks and candidate species we recovered are largely unassociated with shifts in dewlap coloration. We conclude that dewlap color is a highly labile trait that may be misleading if used as the primary diagnostic character for species delimitation. Thus, there is likely substantial unrecognized biodiversity within other polymorphic anole species.
In contrast to the lack of genetic divergence between populations differing in dewlap coloration, we find support for several biogeographic hypotheses. First, we find evidence for a genetic break between populations of A. distichus on the North and South paleo-islands of Hispaniola. We also observe that the Hispaniola satellite island endemic A. d. sejunctus appears to be the result of colonization by the nearest mainland subspecies, A. d. properus, suggesting reconsideration of satellite island endemics as distinct subspecies. In an example of long-distance dispersal to a satellite island, A. altavalensis was likely founded by A. d. ravitergum traveling at least 100 km over-water.
We also posit that the Bahamian distichoids are the result of colonization by A. d. dominicensis from northern Hispaniola. Our insight into such biogeographic patterns will only grow clearer as future studies increase in genomic, taxonomic, and geographic scope.

ACKNOWLEDGMENTS
We would like to thank Julienne Ng, Daniel Scantlebury, Miguel Landestoy, and D. Luke Mahler for assistance with field collections.
We would also like to thank two anonymous reviewers whose comments helped improve the quality of this article. This work was supported by the University of Kansas (REG), NSF DEB #0920892 (REG), NSF DEB #1457774 (REG), NSF DEB #1500761 (REG and AJG), and a Sproull University Fellowship (AJG).

CONFLICT OF INTEREST
None declared.