Adaptive signal coloration maintained in the face of gene flow in a Hispaniolan Anolis Lizard
© The Author(s). 2016
Received: 21 March 2016
Accepted: 7 September 2016
Published: 20 September 2016
Studies of geographic variation can provide insight into the evolutionary processes involved in the early stages of biological diversification. In particular, multiple, replicated cases of geographic trait divergence present a powerful approach to study how patterns of introgression and adaptive divergence can vary with geographic space and time. In this study, we conduct replicated, fine-scaled molecular genetic analyses of striking geographic dewlap color variation of a Hispaniolan Anolis lizard, Anolis distichus, to investigate whether adaptive trait divergence is consistently associated with speciation, whereby genetic divergence is observed with neutral markers, or whether locally adapted traits are maintained in the face of continued gene flow.
We find instances where shifts in adaptive dewlap coloration across short geographic distances are associated with reproductive isolation as well as maintained in the face of gene flow, suggesting the importance of both processes in maintaining geographic dewlap variation.
Our study suggests that adaptive dewlap color differences are maintained under strong divergent natural selection, but this divergence does not necessarily lead to anole speciation.
Geographic trait variation is widely found within species . The evolution of such variation can arise from divergent natural selection driving adaptation to local environments [2, 3], but whether adaptive trait variation persists through time depends on whether divergent selection is strong enough to counteract the homogenizing effects of gene flow [4–6]. Two alternatives may therefore result from ecologically-driven phenotypic divergence: speciation, whereby reproductive isolation has evolved, or selection may maintain a stable coexistence of two or more locally adapted phenotypes despite gene flow. Identifying the consequence of adaptive divergence is central to understanding the roles of local adaptation and speciation during the early stages of the formation of biological diversity.
Comparative studies of naturally replicated hybrid zones where adaptively divergent, closely related taxa come into contact and potentially interbreed can provide a powerful approach to study how patterns of introgression and adaptive divergence can vary with geographic space and time . In particular, fine-scaled molecular genetic analyses across such replicated contact zones can directly test whether divergent natural selection necessarily leads to divergence at neutral loci (e.g. [8–11]). For example, studies of mice and lizards that adaptively vary in coloration for crypsis show that not all geographic transitions in color are associated with reduced gene flow [8, 10]. This suggests that although phenotypically different populations have yet to evolve barriers to reproduction, strong selection is acting to maintain phenotypic divergence, most likely from visual predators [12–14].
A remarkably polymorphic lizard, Anolis distichus, provides an ideal system to investigate whether adaptive trait divergence consistently leads to the same genetic outcome. Across Hispaniola, A. distichus’ dewlap (an extensible throatfan) varies adaptively in response to heterogeneous environments, whereby dewlaps tend to be larger, orange and less bright in wetter habitats, while in drier habitats, dewlaps are smaller, yellow and relatively brighter . This geographic variation in dewlap color has led taxonomists to diagnose eight subspecies on mainland Hispaniola alone  that are characterized by deeply divergent mitochondrial DNA (mtDNA) haplotypes, and differentiation in allozyme and nuclear DNA loci [17–19]. We have also shown in previous work that in two areas of contact, subspecies with different colored dewlaps show a reduction in gene flow, consistent with the expected signature of speciation . However, it remains unclear if these reduced patterns of gene flow occur at other contact zones between divergent dewlap phenotypes, or whether other phenotypic disjunctions represent cases of strong natural selection on dewlap color but with ongoing gene flow. The many other naturally replicated shifts in dewlap color across the range of A. distichus provide an excellent opportunity to assess how selection on an adaptive trait may be associated with patterns of introgression across space and time.
Here, we significantly expand upon prior work by conducting replicated molecular genetic analyses of nine transects that span A. distichus populations. We contrast results from five geographically diverse transects that extend across phenotypically distinct populations that have been diagnosed using differences in dewlap color and pattern, and four control transects that span the same dewlap color. Specifically, we aim to test whether adaptive dewlap color divergence is associated with (i) speciation, whereby genetic differentiation is observed with neutral markers. Following the general lineage concept [21, 22], this result would suggest that populations differing in dewlap color represent independently evolving lineages. Alternatively, (ii) if genetic differentiation is not observed with neutral markers, this suggests that locally adapted traits are maintained in the face of continued gene flow.
For both transitional and control transects, we aimed to sample five or six evenly spaced localities along 10–20 km linear transects (Fig. 1, Additional file 1: Table S1). For the transitional transects, this sampling scheme captured the transition in dewlap color from one dewlap color extreme to the other. Inaccessibility or unsuitable habitat determined transect orientation and also led to some variation in distance between sites (1.74 km to 8.43 km; mean ± SD = 3.26 ± 1.55 km). We previously reported results from sampling along two of these transects, T1a and T2 . Here, however, we include data from an additional site between sites 2 and 3 along transect T2 that appears to be close to the point of contact between the two subspecies (T2-2.5; 0.89 and 1.55 km from sites 2 and 3, respectively). In addition, following initial sampling of T3 (T3-4 to T3-8), we subsequently extended the transect into the range of A. d. ignigularis with three additional sites sampled at larger intervals (T3-1 to T3-3) after preliminary analyses failed to recover evidence for geographic genetic variation. As such, these additional sites ranged from 5.28 - 19.23 km in distance from the next adjacent sites, resulting in T3 spanning a total distance of 42.20 km. All control transects spanned populations of the same subspecies and ranged in distance from 8.98 - 14.03 km, with four or five localities sampled along each transect (Fig. 1, Additional file 1: Table S1). We obtained tissue samples (either tail tips or livers) for molecular genetic analyses from a median of 20 individuals at each transect locality (sample sizes ranged from 10 to 38; Additional file 1: Table S1). In addition to transect sampling, we also sampled 20 individuals from three additional localities (herein called non-transect sites), each of which represent A. d. dominicensis, A. d. ignigularis and A. d. ravitergum elsewhere within their ranges to further investigate geographic genetic variation (Fig. 1).
Characterizing dewlap color variation
We ascertained phenotypic variation across all transects by visually categorizing patterns of dewlap color following previous studies [18, 20]: (i) primarily yellow (<10 % orange), (ii) yellow with a small orange spot (occupying 10–40 % of the dewlap), (iii) yellow with a large orange spot (occupying 50–90 % of the dewlap), or (iv) orange with a narrow yellow margin (>90 % orange). We further categorized any orange coloration as either diffuse (blush) or solid. Assignments of each dewlap were independently assessed and verified by JN and REG. Although we focus on visual assessments of dewlap color (as yellow or orange), this has previously been shown to be consistent with spectrometer measurements of hue, which is an axis of adaptive dewlap divergence [15, 20]. In addition, while some A. distichus dewlaps reflect ultraviolet (UV) wavelengths, UV reflectance is not correlated with different habitat types and is therefore not considered an adaptive trait in this species . Juveniles and females, which lack or have reduced dewlaps, were not scored for dewlap color.
Characterizing genetic variation
Mitochondrial DNA sequencing and phylogenetic analyses
We investigated patterns of mtDNA differentiation along transects by obtaining mtDNA sequence data for all individuals sampled from sites along each transect (N = 955). We focused on a 1147 basepair (bp) region of mtDNA extending from the beginning of ND2 (subunit two of NADH dehydrogenase) through to tRNAAla that has been widely used in previous phylogenetic and phylogeographic studies of anoles (e.g. [23–29]). We also obtained mtDNA sequence data from this same region for individuals at the three non-transect sites as well as 195 additional individuals from 50 other A. distichus populations across Hispaniola, including three Hispaniolan subspecies that did not occur in any of our transects: A. d. sejunctus, a subspecies found on Isla Saona (an island off the south-east coast of the Dominican Republic), and two Haitian subspecies, A. d. vinosus and A. d. aurifer. To represent outgroups, we obtained previously published mtDNA haplotypes from eight species: four species from the A. distichus species group (A. brevirostris, A. caudalis, A. websteri and A. marron), two representatives from clades closely related to the distichus series (A. cristatellus [cristatellus series] and A. bimaculatus [bimaculatus series])  and two more distantly related anoles (A. punctatus and A. occultus) . Our complete mtDNA dataset comprised haplotypes sampled from 1305 individuals, of which 675 were unique; 332 sequences were obtained from previous studies [15, 17, 20, 27] and 973 new sequences were generated for this study (GenBank accession numbers: KX854021-KX855205).
We generated new mtDNA sequences by first extracting genomic DNA from tissue samples using the Wizard genomic DNA purification kit (Promega). Using primers L4437  and H5730 , we conducted 25 μL PCR reactions with 1 μL genomic DNA, 2.5 μL of each 2 μM primer, 2.5 μL of New England Biolabs (NEB) 10× buffer (10 mM Tris–HCl, 50 mM KCl), 2.5 μL of 25 mM MgCl2, 2.5 μL of 0.5 mM dNTP solution, 0.125 μL of NEB Taq polymerase and 11.4 μL H2O. Our thermocycling protocol involved an initial denaturation at 95 °C for 4 min, followed by 30 cycles of 95 °C for 35 s, annealing at 52 °C for 35 s, and extension at 70 °C for 2 min 30 s. We sent successfully amplified PCR products to a commercial facility for purification and Sanger sequencing (Beckman Coulter Genomics, Massachusetts). We assessed each sequence chromatograph for quality in Geneious v4.6.1  and used the pairwise alignment tool in MacClade 4.0  to align sequences. With no indels in the protein coding sequences, alignment was straightforward. The tRNA genes were aligned using secondary structural models . As 43 bp of sequence data from tRNA genes were highly variable in length and a region of apparent ambiguous alignment, we excluded these regions from our dataset prior to analysis, leaving 1104 aligned characters without gaps.
We inferred phylogenetic relationships among A. distichus haplotypes using a maximum likelihood (ML) approach in RAxML v7.0.3  and Bayesian inference in MrBayes v3.2 [37, 38] (alignment and trees deposited in TreeBASE: study ID S19838). We used PartitionFinder v1.1.1  to determine the best fitting partitioning scheme with four data blocks: a separate partition for each codon position within ND2 and a fourth partition for the tRNA region. For RAxML, PartitionFinder identified three partitions to be the best fitting scheme, with GTR + I + G to be the best fitting substitution model for the first ND2 codon position and the tRNA region merged as one partition, as well as the second ND2 codon position, while GTR + G was identified as the most appropriate model of evolution for the third ND2 codon position region. For MrBayes, which allows analyses under more evolutionary models, PartitionFinder identified HKY + I + G to be the most appropriate model of evolution for the first ND2 codon position, GTR + I + G for the second ND2 codon position and tRNA region, and GTR + G for the third ND2 codon position. For both ML and Bayesian analyses, we did not include a parameter for the proportion of invariable sites (I) since the interactions between the I and G parameters can potentially lead to inaccurate estimates due to non-independent optimization . For our ML analysis, we conducted 1000 nonparametric bootstrap replicates and computed pairwise ML patristic distances from the best ML tree. For trees generated under Bayesian inference, we ran two independent analyses for 10 million generations and assessed convergence by checking the average standard deviation of split frequencies in MrBayes and using Tracer v1.6 . We considered there to be genetic differentiation in mitochondrial DNA along our sampled transects if populations from either ends of the transect belonged to distinct mitochondrial clades (node support as measured by bootstrap values and posterior probabilities ≥ 75 %).
For each transect site sampled, we calculated the number of haplotypes, number of polymorphic sites, nucleotide diversity and the mean number of pairwise differences between haplotypes using ARLEQUIN v3.11. Using the same program, we also assessed genetic differentiation between populations using ΦST (an analogue of FST) and tested for significant genetic differences between sites using 10,000 permutations. We also constructed a haplotype network that comprised haplotypes from all transect sites in PopART (http://popart.otago.ac.nz) using the TCS method .
Microsatellite genotyping and population structure analyses
To assess nuclear differentiation along transects, we genotyped all individuals using 7 highly polymorphic di- and trinucleotide microsatellite loci: DISTA2B12, DISTBB5, DISTBC4, DISTAG1, DISTCC7, DISTAH6 and BREV2E9 . We followed the same PCR protocols as described for mtDNA sequence data (above section) but, following Ng and Glor , we included a fluorescently-labelled CAG primer (6-FAM, VIC or NED) at a 1:20 ratio (0.01 μM CAG-tagged forward primer and 0.19 μM fluorescent CAG primer) to enable multiplexing. We followed thermocycling conditions detailed in Ng et al.  and visualized products on an ABI 3730 Genetic Analyzer (Applied Biosystems) located at the University of Rochester’s Functional Genomics Center. Genotypic data were analyzed with GENEMAPPER V3.7 software (Applied Biosystems). We tested all microsatellite loci for within population departures from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) using GENEPOP v4.0 , and accounted for multiple tests using sequential Bonferroni correction .
We tested for nuclear genetic differentiation along each transect using (i) genetic distance metrics and (ii) a Bayesian clustering approach. We did not use cline analyses (e.g. [46, 47]) because our microsatellite markers are not diagnostic; they do not show fixed, or almost fixed, differences between the two subspecies along transects (see Results). We calculated pairwise estimates of ΦST between adjacent sites along each transect using ARLEQUIN v3.11 and tested for significant genetic differences between sites using 10,000 permutations. In addition, we used a genetic distance metric based on the proportion of shared alleles between populations (DPS) , calculated using MSA v4.05 . DPS has an advantage over ΦST in that it does not assume population equilibrium [50, 51].
We conducted Bayesian clustering analyses using the program STRUCTURE v2.3.3, which uses a Markov chain Monte Carlo algorithm to probabilistically assign individuals to different genetic clusters (K) . We inferred genetic structure after combining data across all transects and the three non-transect sites. We used the admixture model, which allows individuals to have mixed ancestry, in addition to the correlated allele frequencies model, which accounts for potential linkage disequilibrium between markers that can arise due to admixture , and did not include a priori information about where individuals were sampled. We ran analyses for 105 iterations following a burn-in period of 104 iterations for each value of K ranging from 1–15. We chose a maximum K value of 15 because we expected that there would be at least 8 different genetic clusters across our transitional transects if genetic differentiation is associated with adaptive dewlap color transition (or at least 10 if individuals of the same subspecies were genetically different along T1a and T1b). We ran 10 independent runs of K to ensure consistent probability estimates and that particular runs were not trapped in different modes in the parameter space . The average of all ten runs were subsequently used to assess genetic structure. To identify the best estimate of K, we used log probabilities Pr(X|K) and calculated ΔK  using the web-based program STRUCTURE HARVESTER . We also investigated whether there was more subtle population structure by repeating the analysis on each genetic cluster identified. Following protocols outlined in Coulon et al. , we assigned individuals to a particular genetic cluster if their inferred ancestry to that group was higher than 0.6. With each genetic cluster, we sequentially ran STRUCTURE analyses using the same parameters as the initial run and continued this process until the optimal K was 1.
We interpreted results of distinct nuclear genotypic clusters at either ends of the transect to be evidence for neutral genetic divergence. We diagnosed how far along populations were in the speciation continuum by the frequency of admixed individuals. If geographic genetic structure was evident but the frequency of admixed individuals was high and widely spread across different populations, we considered this to reflect low levels of restricted gene flow, and that the populations were at an earlier stage in the speciation continuum. In contrast, a limited number or no evidence of admixed individuals indicated a strong reduction of gene flow, likely representing populations in the later stages of speciation.
Dewlap color variation
We found similar dewlap phenotypes across each control transect. All males sampled along C3 had dewlaps with greater than 60 % orange, while individuals sampled along all other control transects had very little orange in the dewlaps; the two individuals at C2-1 that exhibited greater than 50 % orange in the dewlap was diffuse in coloration (Figs. 1 and 2).
Mitochondrial phylogenetic relationships
Along the control transects, C1 and C3 comprised mtDNA haplotypes from the same clade (Figs. 1, 2 and 3). We found haplotypes from two closely related clades along C2 and C4 (4.7 % and 5.3 % average uncorrected sequence divergence, respectively) but these haplotypes did not show any geographic structure (Figs. 1 and 2).
Assessing overall phylogenetic relationships among A. distichus mtDNA haplotypes, we found A. d. favillarum to be the only monophyletic Dominican subspecies from mainland Hispaniola. Mitochondrial haplotypes from the other four Dominican A. distichus subspecies from mainland Hispaniola (A. d. dominicensis, A. d. ignigularis, A. d. properus and A. d. ravitergum) did not form monophyletic groups (Fig. 3). Instead, A. d. ravitergum and A. d. properus sampled along our transects were each rendered paraphyletic by A. d. ignigularis sampled along the same transects (T2 and T3 respectively). Our finding that most A. distichus subspecies are non-monophyletic is inconsistent with other phylogenetic studies of the group [17, 19], but is likely due to our inclusion of individuals from subspecific contact zones where mitochondrial introgression can confound phylogenetic placement. Our phylogeny, however, is consistent with Geneva et al.  in showing well-supported haplotype clades within A. d. dominicensis that appear to be geographically structured; haplotypes from transect C1 in central Hispaniola are most closely related to haplotypes sampled from A. d. favillarum and Haitian A. distichus subspecies, while A. d. dominicensis haplotypes sampled from northern Dominican Republic (T1a and b, non-transect A. d. dominicensis site) formed a separate group (Fig. 3).
Microsatellite loci tests of HWE and LD
The microsatellite loci exhibited significant deviations from HWE in 57 of 364 intrapopulation assessments (Additional file 1: Table S1). All loci but DISTAH6 deviated from HWE in only a few populations (3–10 populations), while DISTAH6 deviated from HWE in 23 of 57 sample sites. Thus, we ran analyses using all seven loci and then repeated analyses without genotypic data from DISTAH6. As analyses did not qualitatively differ, we here report results of analyses on all seven loci. We detected little evidence of LD among microsatellite loci; in 6 of 52 sampled sites, one pair of loci showed significant LD after Bonferroni correction (P < 0.007). These pairs differed between sites, except for DISTCC7 and DISTAH6, which showed significant LD in two sites along T4 (T4-1 and T4-3).
Patterns of nuclear genetic structure
The results of our STRUCTURE analysis identified K = 9 to be the best estimate of K. Of the transitional transects, three showed that individuals on either side of the transect belonged to different genetic clusters (T1a and b, T2; Figs. 1 and 2). Sampled sites in the middle of these transects included individuals from both genetic clusters as well as admixed individuals. While at least three localities included admixed individuals along T1a and T1b, we only found admixed individuals at one site along T2 (T2-2.5). The other two transitional transects (T3 and T4) did not show any genetic structure despite T3 being the longest transect; most individuals sampled along the transect were predominantly assigned to the same cluster. Additional runs on each genetic cluster revealed two further genetic clusters along T4 but there was no obvious geographic structuring (Fig. 2). Individuals along control transects were primarily assigned to the same genetic cluster as others sampled from the same transect (Fig. 2).
While STRUCTURE analyses showed T1a, T1b and T2 to all exhibit genetic differentiation, pairwise ΦST and DPS estimates showed that genetic differentiation along T1a and T1b were comparable to that found along the control transects and the transitional transects, T3 and T4, which did not exhibit genetic structure in the STRUCTURE analyses (Fig. 4, Additional file 1: Figure S2). Meanwhile, pairwise ΦST and DPS estimates along T2 showed moderately high genetic differentiation between T2-2 and T2-2.5 (ΦST = 0.12; DPS = 0.48) that, after correcting for differences in geographic distance, were at least 7 and 3 times greater respectively than that found along any other transect (geographic distance corrected: ΦST = 0.15; DPS =0.54; Fig. 4, Additional file 1: Figure S2).
Assessing general genetic structure within subspecies, we found that subspecies that have a widely distributed range exhibited nuclear genetic structure. The most widespread of the subspecies, A. d. dominicensis, comprised two clusters: one cluster comprised the control transect (C1) while the other included some T1a and b populations (Figs. 1 and 2). The non-transect A. d. dominicensis site comprised individuals, or proportions of individuals, assigned to either one of the two clusters found along T1a and b. We found three clusters within A. d. ignigularis: one cluster included some sites along T2, one included some sites sampled along T1a and b, while the other included the A. d. ignigularis control transect (C3) (Figs. 1 and 2). The non-transect A. d. ignigularis site looked similar to sites found along T3, whereby it comprised individuals, or proportions of individuals, assigned to either cluster found along C3 and C4. All A. d. ravitergum sampled were initially assigned to the same genetic cluster (T2, C2 and the non-transect A. d. ravitergum site) (Figs. 1 and 2). Additional STRUCTURE analyses showed further substructure into two clusters: C2 comprised one cluster, while T2 and the non-transect site were assigned to another cluster (results not shown). Further STRUCTURE analyses differentiated T2 A. d. ravitergum from the non-transect site.
Our study examining numerous transitions in adaptive dewlap coloration in Anolis distichus across the Dominican Republic suggests that adaptive trait divergence is associated with both speciation and local adaptation in the face of gene flow. We find evidence for local adaptation at two transects (T3 and T4), whereby there was no evidence of genetic differentiation despite bimodal phenotypic variation. Patterns of mtDNA variation suggest that introgression is occurring and nuclear variation at microsatellite loci does not show evidence of any genetic structure. This is further supported by pairwise estimates of ΦST and DPS that indicate similar levels of pairwise genetic differentiation to that found along the control transects. Despite the homogenizing effects of gene flow, the maintenance of dewlap color variation at such a fine spatial scale suggests that strong divergent natural selection is acting on loci underlying dewlap color differences between spatially close sites while neutral regions of the genome are homogenized with gene flow. This also suggests that the genes underlying dewlap color likely have different evolutionary trajectories that is not reflected by the results of our analyses of neutral markers. Together, these results support the importance of divergent natural selection, rather than neutral processes, in driving and maintaining dewlap color variation along these transects. These transects are similar to previous studies that have shown that adaptive trait divergence between populations can be maintained in the face of gene flow (e.g. [57–60]), including cases of adaptive color divergence [8–10, 61, 62]. It is also possible that while we did not find evidence of genetic differentiation with our markers, such ecologically-based trait divergence may represent the initial step towards speciation [63–66].
We find patterns consistent with the expected signature of speciation along three transects (T1a, T1b and T2) whereby adaptive dewlap color divergence is associated with genetic differentiation. The results of our control transects, showing no genetic differentiation where there was no change in dewlap color, further support the possibility of speciation occurring at these three transitional transects. The transects, however, exhibited differing levels of genetic differentiation and therefore likely represent different stages along the speciation continuum. T1a and b appear to be at an earlier stage along the speciation continuum relative to T2. Along T1a and b, phylogenetic and STRUCTURE analyses show distinct mtDNA clades and nuclear genetic clusters on either side of the transect, respectively. Both mtDNA and nuclear markers, however, suggest a low reduction in gene flow between the phenotypically divergent populations: mtDNA haplotypes prevalent on the western side of the transect were still found at a low frequency in the middle and other side of the transect, and nuclear markers showed evidence of admixed genotypes that decreased in frequency from the middle of the transect to the east. Furthermore, pairwise genetic distance estimates (ΦST and DPS) exhibit similar levels to that found along transects that show no genetic differentiation (T3, T4 and all control transects). Mitochondrial haplotype structure has been shown to reflect historic geological events while the more quickly evolving microsatellite loci are more likely to reflect contemporary gene flow [11, 67–69]. As such, the deeply divergent haplotypes we sampled from either side of the transects suggest hybridization upon recent secondary contact between A. d. dominicensis and A. d. ignigularis. While we know that these two subspecies along T1a were once isolated by a seawater channel that separated the Samaná Peninsula from mainland Hispaniola [20, 70], we are unaware of any historical barrier that isolated subspecies along T1b.
Closer to the other end of the speciation continuum, T2 exhibits evidence of a stronger reduction in nuclear gene flow; Bayesian clustering analyses suggest that hybridization is limited to one site (T2-2.5) whereby both parental genotypes and admixed individuals co-occur, and pairwise distance estimates between T2-2 and T2-2.5 are at least 3 times higher than that found along any other transect, despite closely related mtDNA haplotypes showing evidence of introgression. This pattern was also shown in Ng and Glor , but evidence of hybridization was only found with the additional sampling we include in this study. Our results showing high levels of genetic differentiation at nuclear loci along T2 is further supported by a previous phylogenetic reconstruction based on seven nuclear genes that showed that the same two subspecies collected at some of the same localities each formed monophyletic groups . Together, this suggests that the two subspecies along T2 represent independently evolving lineages. The maintenance of a narrow hybrid zone, with no evidence of hybrids at T2-2 or T2-3 (0.89 and 1.55 km away, respectively) suggests strong intrinsic or extrinsic selection against hybrids and almost complete reproductive isolation between A. d. ravitergum and A. d. ignigularis.
While three transects support the predictions of speciation, it remains unknown whether any of the transitional transects we investigated in this study will progress to speciation, or whether the current stages we observe represent stable states of migration-selection equilibria . The presence of admixed individuals at all transitional transects suggest that some level of gene flow is still occurring, and whether they progress to speciation will likely depend on the strength of divergent selection on dewlap coloration and/or whether divergent selection is also acting on other traits involved with fitness (‘multifarious selection’) (reviewed in ). Given that the degree to which dewlap color diverged across transects did not correlate with genetic divergence, it is likely that other traits are also involved in reproductive isolation and speciation in A. disitchus and that the dewlap does not serve as a ‘magic trait’ [66, 73]. For example, the abrupt reduction in gene flow across T2 was accompanied by a gradual shift in dewlap color, while T4, which showed the most abrupt shift in orange in the dewlap, did not show any genetic structure. Instead, our evidence of some genetic differentiation along T1a, b and T2, suggests that the reduction in gene flow may be more directly associated with another adaptive trait not quantified in this study. One such trait could be scale differences, which were also used to designate A. distichus subspecies . Differences in scale number or size has been shown to vary within and among other anole species, and may represent adaptive divergences to regulate water loss or body temperature in different habitats [62, 74, 75]. If weak selection on multiple traits can promote speciation , this may be one reason why the polymorphic A. d. favillarum along T4, which only differs in dewlap color, does not show genetic structure, despite exhibiting the most abrupt dewlap color shift of all transects.
Our analysis of multiple replicates of dewlap color divergence in A. distichus showed that geographic variation in adaptive dewlap color is associated with both local adaptation and speciation. At some transitional transects, there appears to be strong selection maintaining adaptive dewlap coloration in the face of gene flow. At other contact zones, we found signatures of speciation whereby adaptive dewlap divergence is associated with genetic differentiation, indicating some reduction in gene flow across the transects. Our results are consistent with previous work suggesting that within-island speciation in anoles is limited to the four large Greater Antillean islands . Despite striking adaptive variation in head, body and dewlap coloration in anoles on the smaller islands of the Lesser Antilles, there are no examples of phenotypically divergent populations on these islands that show as strong a reduction in gene flow as we have found between A. d. ravitergum and A. d. ignigularis across T2 [9, 11, 68].
We are grateful to Daniel Scantlebury, Anthony Geneva and Miguel Landestoy for their help collecting field samples, and Jorge Brocca for help obtaining permits. We thank members of the Glor Lab and Robert Laport for helpful comments on an early draft of this manuscript.
This work was supported by a National Science Foundation grant to R.E.G. (NSF-DEB 0920892).
Availability of data
All DNA sequence data that was newly generated for this study has been deposited into GenBank (accession numbers KX854021-KX855205). The DNA sequence alignment file and phylogenetic trees reported in this study have been deposited into TreeBASE (study ID S19838).
JN and REG designed and coordinated the study. JN, AGO and REG contributed to data collection and analyses. JN and REG drafted the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Tissue samples were collected with the authorization of the Ministerio de Medio Ambiente y Recursos Naturales in the Dominican Republic.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Endler JA. Geographic variation, speciation and clines. Princeton, New Jersey: Princeton University Press; 1977.Google Scholar
- Nosil P. Ecological speciation. New York: Oxford University Press; 2012.View ArticleGoogle Scholar
- Schluter D. The ecology of adaptive radiation. New York: Oxford University Press; 2000.Google Scholar
- Lenormand T. From local adaptation to speciation: specialization and reinforcement. Int J Ecol. 2012;2012(Article ID 508458):1–11.View ArticleGoogle Scholar
- Slatkin M. Gene flow and the geographic structure of natural populations. Science. 1987;236(4803):787–92.View ArticlePubMedGoogle Scholar
- Smadja CM, Butlin RK. A framework for comparing processes of speciation in the presence of gene flow. Mol Ecol. 2011;20(24):5123–40.View ArticlePubMedGoogle Scholar
- Harrison RG, Larson EL. Hybridization, introgression, and the nature of species boundaries. J Hered. 2014;105:795–809.View ArticlePubMedGoogle Scholar
- Hoekstra HE, Drumm KE, Nachman MW. Ecological genetics of adaptive color polymorphism in pocket mice: geographic variation in selected and neutral genes. Evolution. 2004;58(6):1329–41.View ArticlePubMedGoogle Scholar
- Muñoz MM, Crawford NG, McGreevy TJ, Messana NJ, Tarvin RD, Revell LJ, Zandvliet RM, Hopwood JM, Mock E, Schneider AL, et al. Divergence in coloration and ecological speciation in the Anolis marmoratus species complex. Mol Ecol. 2013;22(10):2668–82.View ArticlePubMedGoogle Scholar
- Rosenblum EB, Harmon LJ. “Same same but different”: replicated ecological speciation at White Sands. Evolution. 2011;65(4):946–60.View ArticlePubMedGoogle Scholar
- Thorpe RS, Surget-Groba Y, Johansson H. Genetic tests for ecological and allopatric speciation in anoles on an island archipelago. PLoS Genet. 2010;6(4):e1000929.View ArticlePubMedPubMed CentralGoogle Scholar
- Vignieri SN, Larson JG, Hoekstra HE. The selective advantage of crypsis in mice. Evolution. 2010;64(7):2153–8.PubMedGoogle Scholar
- Dice LR: Effectiveness of selection by owls of deer mice (Peromyscus maniculatus) which contrast in color with their background. Contributions of the Laboratory of Vertebrate Biology, University of Michigan 1947, 34:1–20Google Scholar
- Kaufman DW. Adaptive coloration in Peromyscus polionotus: experimental selection by owls. J Mammal. 1974;55(2):271–83.View ArticleGoogle Scholar
- Ng J, Landeen EL, Logsdon RM, Glor RE. Correlation between Anolis lizard dewlap phenotype and environmental variation indicates adaptive divergence of a signal important to sexual selection and species recognition. Evolution. 2013;67(2):573–82.View ArticlePubMedGoogle Scholar
- Schwartz A. Geographic variation in Anolis distichus Cope (Lacertilia, Iguanidae) in the Bahama Islands and Hispaniola. Bull Mus Comp Zool. 1968;137:255–309.Google Scholar
- Glor RE, Laport RG. Are subspecies of Anolis lizards that differ in dewlap color and pattern also genetically distinct? A mitochondrial analysis. Mol Phylogen Evol. 2012;64(2):255–60.View ArticleGoogle Scholar
- Williams EE, Case SM. Interactions among members of the Anolis distichus complex in and near the Sierra de Baoruco, Dominican Republic. J Herpetol. 1986;20(4):535–46.View ArticleGoogle Scholar
- Geneva AJ, Hilton J, Noll S, Glor RE. Multilocus phylogenetic analyses of Hispaniolan and Bahamian trunk anoles (distichus species group). Mol Phylogen Evol. 2015;87:105–17.View ArticleGoogle Scholar
- Ng J, Glor RE. Genetic differentiation among populations of a Hispaniolan trunk anole that exhibit geographical variation in dewlap colour. Mol Ecol. 2011;20(20):4302–17.View ArticlePubMedGoogle Scholar
- de Queiroz K. The general lineage concept of species, species criteria, and the process of speciation: a conceptual unification and terminological recommendations. In: Howard DJ, Berlocher SH, editors. Endless Forms: Species and Speciation. Oxford: Oxford University Press; 1998. p. 57–75.Google Scholar
- de Queiroz K. The general lineage concept of species and the defining properties of the species category. In: Wilson RA, editor. Species: New Interdisciplinary Essays. Cambridge, Massachusetts: MIT Press; 1999. p. 49–89.Google Scholar
- Castañeda MR, Kd Q. Phylogenetic relationships of the Dactyloa clade of Anolis lizards based on nuclear and mitochondrial DNA sequence data. Mol Phylogen Evol. 2011;61(3):784–800.View ArticleGoogle Scholar
- Glor RE, Losos JB, Larson A. Out of Cuba: overwater dispersal and speciation among lizards in the Anolis carolinensis subgroup. Mol Ecol. 2005;14(8):2419–32.View ArticlePubMedGoogle Scholar
- Jackman TR, Larson A, Queiroz KD, Losos JB. Phylogenetic relationships and tempo of early diversification in Anolis lizards. Syst Biol. 1999;48(2):254–85.View ArticleGoogle Scholar
- Kolbe JJ, Glor RE, Schettino LRG, Lara AC, Larson A, Losos JB. Genetic variation increases during biological invasion by a Cuban lizard. Nature. 2004;431(7005):177–81.View ArticlePubMedGoogle Scholar
- Nicholson KE, Glor RE, Kolbe JJ, Larson A, Hedges SB, Losos JB. Mainland colonization by island lizards. J Biogeogr. 2005;32(6):929–38.View ArticleGoogle Scholar
- Rabosky DL, Glor RE. Equilibrium speciation dynamics in a model adaptive radiation of island lizards. Proc Natl Acad Sci U S A. 2010;107(51):22178–83.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang IJ, Glor RE, Losos JB. Quantifying the roles of ecology and geography in spatial genetic divergence. Ecol Lett. 2013;16(2):175–82.View ArticlePubMedGoogle Scholar
- Brandley MC, De Queiroz K. Phylogeny, ecomorphological evolution, and historical biogeography of the Anolis cristatellus series. Herpetol Monogr. 2004;18:90–126.View ArticleGoogle Scholar
- Macey JR, Larson A, Ananjeva NB, Papenfuss TJ. Evolutionary shifts in three major structural features of the mitochondrial genome among iguanian lizards. J Mol Evol. 1997;44(6):660–74.View ArticlePubMedGoogle Scholar
- Glor RE, Gifford ME, Larson A, Losos JB, Schettino LR, Lara ARC, Jackman TR. Partial island submergence and speciation in an adaptive radiation: a multilocus analysis of the Cuban green anoles. Proc R Soc Lond Ser B Biol Sci. 2004;271(1554):2257–65.View ArticleGoogle Scholar
- Drummond AJ, Kearse M, Heled J, Moir R, Thierer T, Ashton B, Wilson A, Stones-Havas S. Geneious v4.6.1. 2006. Available from http://www.geneious.com/.
- Maddison DR, Maddison WP. MacClade 4: Analysis of phylogeny and character evolution. Version 4.0. Sunderland, Massachusetts: Sinauer Associates; 2000.Google Scholar
- Kumazawa Y, Nishida M. Sequence evolution of mitochondrial transfer-RNA genes and deep-branch animal phylogenetics. J Mol Evol. 1993;37(4):380–98.View ArticlePubMedGoogle Scholar
- Stamatakis A. RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22(21):2688–90.View ArticlePubMedGoogle Scholar
- Huelsenbeck JP, Ronquist F. MRBAYES: Bayesian inference of phylogeny. Bioinformatics. 2001;17:754–5.View ArticlePubMedGoogle Scholar
- Ronquist F, Huelsenbeck JP. MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572–4.View ArticlePubMedGoogle Scholar
- Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29(6):1695–701.View ArticlePubMedGoogle Scholar
- Sullivan J, Swofford D, Naylor G. The effect of taxon sampling on estimating rate heterogeneity parameters of maximum-likelihood models. Mol Biol Evol. 1999;16(10):1347–56.View ArticleGoogle Scholar
- Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6. Available from http://beast.bio.ed.ac.uk/Tracer. Accessed June 2016.
- Clement M, Snell Q, Walke P, Posada D, Crandall K. TCS: estimating gene genealogies, Proc 16th Int Parallel Distrib Process Symp, vol. 2. 2002. p. 184.Google Scholar
- Ng J, Perkins SL, Dussmann EJ, Glor RE. Eleven highly polymorphic microsatellite markers for the lizard Anolis distichus. Conserv Genet Resour. 2010;1(1):135–9.View ArticleGoogle Scholar
- Raymond M, Rousset F. GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism. J Hered. 1995;86:248–9.Google Scholar
- Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43(1):223–5.View ArticleGoogle Scholar
- Gompert Z, Buerkle CA. A powerful regression-based method for admixture mapping of isolation across the genome of hybrids. Mol Ecol. 2009;18(6):1207–24.View ArticlePubMedGoogle Scholar
- Szymura JM, Barton NH. Genetic analysis of a hybrid zone between the fire-bellied toads, Bombina bombina and B. variegata, near Cracow in southern Poland. Evolution. 1986;40(6):1141–59.View ArticleGoogle Scholar
- Bowcock AM, Ruizlinares A, Tomfohrde J, Minch E, Kidd JR, Cavallisforza LL. High resolution of human evolutionary trees with polymorphic microsatellites. Nature. 1994;368(6470):455–7.View ArticlePubMedGoogle Scholar
- Dieringer D, Schlotterer C. MICROSATELLITE ANALYSER (MSA): a platform independent analysis tool for large microsatellite data sets. Mol Ecol Notes. 2003;3(1):167–9.View ArticleGoogle Scholar
- Murphy MA, Evans JS, Storfer A. Quantifying Bufo boreas connectivity in Yellowstone National Park with landscape genetics. Ecology. 2010;91(1):252–61.View ArticlePubMedGoogle Scholar
- Storfer A, Murphy MA, Spear SF, Holderegger R, Waits LP. Landscape genetics: where are we now? Mol Ecol. 2010;19(17):3496–514.View ArticlePubMedGoogle Scholar
- Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.PubMedPubMed CentralGoogle Scholar
- Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies. Genetics. 2003;164(4):1567–87.PubMedPubMed CentralGoogle Scholar
- Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14(8):2611–20.View ArticlePubMedGoogle Scholar
- Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4(2):359–61.View ArticleGoogle Scholar
- Coulon A, Fitzpatrick JW, Bowman R, Stith BM, Makarewich CA, Stenzler LM, Lovette IJ. Congruent population structure inferred from dispersal behaviour and intensive genetic surveys of the threatened Florida scrub-jay (Aphelocoma coerulescens). Mol Ecol. 2008;17(7):1685–701.View ArticlePubMedGoogle Scholar
- Berner D, Grandchamp A-C, Hendry AP. Variable progress toward ecological speciation in parapatry: stickleback across eight lake-stream transitions. Evolution. 2009;63(7):1740–53.View ArticlePubMedGoogle Scholar
- Saint-Laurent R, Legault M, Bernatchez L. Divergent selection maintains adaptive differentiation despite high gene flow between sympatric rainbow smelt ecotypes (Osmerus mordax Mitchill). Mol Ecol. 2003;12(2):315–30.View ArticlePubMedGoogle Scholar
- Schneider CJ, Losos JB, de Queiroz K. Evolutionary relationships of the Anolis bimaculatus group from the northern Lesser Antilles. J Herpetol. 2001;35(1):1–12.View ArticleGoogle Scholar
- Smith TB, Schneider CJ, Holder K. Refugial isolation versus ecological gradients. Genetica. 2001;112–113(1):383–98.View ArticlePubMedGoogle Scholar
- Mullen LM, Hoekstra HE. Natural selection along an environmental gradient: A classic cline in mouse pigmentation. Evolution. 2008;62(7):1555–69.View ArticlePubMedGoogle Scholar
- Thorpe RS, Surget-Groba Y, Johansson H. Quantitative traits and mode of speciation in Martinique anoles. Mol Ecol. 2012;21(21):5299–308.View ArticlePubMedGoogle Scholar
- Funk DJ, Nosil P, Etges WJ. Ecological divergence exhibits consistently positive associations with reproductive isolation across disparate taxa. Proc Natl Acad Sci U S A. 2006;103(9):3209–13.View ArticlePubMedPubMed CentralGoogle Scholar
- Rundle HD, Nosil P. Ecological speciation. Ecol Lett. 2005;8(3):336–52.View ArticleGoogle Scholar
- Schemske DW. Adaptation and the origin of species. Am Nat. 2010;176:S4–S25.View ArticlePubMedGoogle Scholar
- Servedio MR, Doorn GSV, Kopp M, Frame AM, Nosil P. Magic traits in speciation: ‘magic’ but not rare? Trends Ecol Evol. 2011;26(8):389–97.View ArticlePubMedGoogle Scholar
- Malhotra A, Thorpe RS. The dynamics of natural selection and vicariance in the Dominican anole: Patterns of within-island molecular and morphological divergence. Evolution. 2000;54(1):245–58.View ArticlePubMedGoogle Scholar
- Stenson AG, Malhotra A, Thorpe RS. Population differentiation and nuclear gene flow in the Dominican anole (Anolis oculatus). Mol Ecol. 2002;11(9):1679–88.View ArticlePubMedGoogle Scholar
- Thorpe RS, Stenson AG. Phylogeny, paraphyly and ecological adaptation of the colour and pattern in the Anolis roquet complex on Martinique. Mol Ecol. 2003;12(1):117–32.View ArticlePubMedGoogle Scholar
- Grant C. Report on a collection of Hispaniolan reptiles. Herpetologica. 1956;12(2):85–90.Google Scholar
- Lenormand T. Gene flow and the limits to natural selection. Trends Ecol Evol. 2002;17(4):183–9.View ArticleGoogle Scholar
- Nosil P, Harmon LJ, Seehausen O. Ecological explanations for (incomplete) speciation. Trends Ecol Evol. 2009;24(3):145–56.View ArticlePubMedGoogle Scholar
- Gavrilets S. Fitness landscapes and the origin of species. Princeton, NJ: Princeton University Press; 2004.Google Scholar
- Calsbeek R, Knouft JH, Smith TB. Variation in scale numbers is consistent with ecologically based natural selection acting within and between lizard species. Evol Ecol. 2006;20(4):377–94.View ArticleGoogle Scholar
- Malhotra A, Thorpe RS. Microgeographic variation in scalation of Anolis oculatus (Dominica, West Indies): a multivariate analysis. Herpetologica. 1997;53(1):49–62.Google Scholar
- Nosil P, Funk DJ, Ortiz-Barrientos D. Divergent selection and heterogeneous genomic divergence. Mol Ecol. 2009;18(3):375–402.View ArticlePubMedGoogle Scholar
- Losos JB, Schluter D. Analysis of an evolutionary species-area relationship. Nature. 2000;408(6814):847–50.View ArticlePubMedGoogle Scholar