Abstract

The general diversity pattern of the Caribbean anole radiation has been described in detail; however, the actual mechanisms at the origin of their diversification remain controversial. In particular, the role of ecological speciation, and the relative importance of divergence in allopatry and in parapatry, is debated. We describe the genetic structure of anole populations across lineage contact zones and ecotones to investigate the effect of allopatric divergence, natural selection, and the combination of both factors on population differentiation. Allopatric divergence had no significant impact on differentiation across the lineage boundary, while a clear bimodality in genetic and morphological characters was observed across an ecotone within a single lineage. Critically, the strongest differentiation was observed when allopatry and ecology act together, leading to a sharp reduction in gene flow between two lineages inhabiting different habitats. We suggest that, for Caribbean anoles to reach full speciation, a synergistic combination of several historical and ecological factors may be requisite.

1. Introduction

Speciation, the mechanism at the origin of species diversification, is one of the most studied subjects in evolutionary biology. Despite this enormous interest, very little is known about the factors needed for speciation to occur. For instance, the relative importance of ecological versus purely historical factors, as well as the geographic context of speciation, is still debated, and mechanisms that are sometimes invoked to explain lack of speciation, such as observation of a species-area relationship [1], are themselves not fully understood.

The most widely recognised speciation model is the allopatric model, where different populations that are geographically isolated develop genetic incompatibilities, purely by genetic drift or founder effects [24], by adapting to different habitats (by-product speciation [5]), by sexual selection [6], or by fixation of incompatible mutations through adaptation to a similar habitat (mutation-order speciation [7]). Because there is no gene flow to counter the differentiation of populations, this model is the easiest to explain speciation, especially if the isolated populations are exposed to distinct selective environments (reviewed in [8]).

The possibility of speciation in the presence of recurrent gene flow (sympatric or parapatric speciation) is much more debated. In this model, the divergence between populations exchanging migrants is driven by ecological differentiation (ecological speciation [9]), and/or sexual differentiation (speciation by sexual selection [10]). Because even a limited amount of gene flow is expected to counter differentiation, early theoretical studies rejected this speciation model [11]. However, other theoretical and empirical work suggested that nonallopatric speciation is possible under particular circumstances (reviewed in [12]).

A group in which both the geographic context of speciation, and the role of ecological speciation, is debated is the Caribbean anoles. The anole radiation is a highly diverse species group useful for understanding the mechanisms at the origin of diversification. Anolis is one of the most speciose vertebrate genera with more than 400 described species. Of these, ca. 150 are found in the Caribbean islands, and are thought to originate from just two independent colonizations from the mainland. This species diversity is accompanied by a very high morphological and ecological diversity. For these reasons, the anole radiation in the Caribbean has been the focus of intense studies in the last decades (reviewed in [13]).

One of the striking patterns in Anolis diversity is the contrast between the Greater and the Lesser Antilles. Whereas numerous species coexist within the large islands of the Greater Antilles, islands of the Lesser Antilles have at most two naturally occurring species with most of the islands having only one. Furthermore, if most of the diversification in the Greater Antilles occurred within island, most of the species pairs in the Lesser Antilles are not sister species, suggesting that these species pairs did not diverge within island but rather came together after dispersal from another island, with one possible exception in Saint Vincent [13].

The Caribbean anoles show a well-documented species-area relationship [14, 15]; however, the reason why so many speciation events happened in the Greater Antilles while almost none happened in the Lesser Antilles is still speculative. Losos and Schluter [14] proposed two explanations for this observation. First, larger islands could offer more opportunities of geographic fragmentation and hence lead to the formation of more species by allopatric speciation. Second, larger islands may have more habitat diversity, and hence populations submitted to divergent selective pressure in different habitats could lead to the formation of more species by ecological speciation. The observation that some of the largest islands in the Lesser Antilles, like Dominica, Guadeloupe, or Martinique, have a very high habitat diversity and hence should be able to support several species if ecological speciation was the driving force in anole diversification led Losos [13] to suggest that nonallopatric modes of speciation driven by ecological speciation are not supported in anoles.

Several studies suggest the opposite. For instance, because of its complex geologic history [16], the island of Martinique has offered plenty of opportunities for allopatric speciation to occur in its endemic anole, Anolis roquet. Present day Martinique was once formed of separate proto-islands where distinct mtDNA lineages of A. roquet evolved in allopatry for millions of years before coming back into secondary contact [17]. However, high gene flow is observed between previously allopatric lineages showing that this long geographic isolation did not lead to complete speciation [18]. Similarly, deep mtDNA lineages are observed in several species, both in the Lesser Antilles [1923], and in the Greater Antilles [2431]. In the cases where it has been studied, no restriction of gene flow has been detected between these previously allopatric lineages [18, 32, 33], with one possible exception in North-Eastern Martinique [18]. All this suggests that even if allopatric speciation probably occurred in some situations, geographic isolation alone is not sufficient to reach speciation in anoles.

Instead, ecological gradients seem to be driving population differentiation in Martinique anoles. This island is very heterogeneous, both topographically and ecologically. The mountains in the North are exposed to the trade winds and receive a very high amount of precipitation all year round. At the opposite, the northern Caribbean coast is in the “rain shadow” of these mountains and is much drier and seasonal. Hence, the habitat changes dramatically from a cool montane rainforest to a hot xeric scrubland in just a few kilometres. Previous studies [18] have shown that the divergent selective forces along this habitat gradient lead to significant morphological and genetic differentiation between coastal and mountain populations of A. roquet. On the neighbouring island of Dominica, a similar situation appears to occur with its endemic anole, Anolis oculatus [32], and indeed, in the single case to date where contact zones have been studied in the Greater Antilles, a significant reduction in gene flow between divergent lineages is only observed in an area with a steep ecological gradient [33].

In this paper we reanalysed data from a previous study [18] and added a new transect of A. roquet populations where the ecotone and the lineage boundary overlap. We compared the population structure of this new transect to the two previously published ones, one between two lineages within a single habitat, and the other between two habitats within a single lineage. We studied the population structure and admixture rates along these different transects to investigate the effects of geographic isolation, ecological isolation, and the combination of these two factors on anole population differentiation and speciation. We observe that population differentiation is at its highest when both factors act simultaneously and suggest that a possible explanation of the species-area relationship observed in Caribbean anoles is that the probability of both allopatric and ecological factors acting in synergy increases with island size.

2. Material and Methods

2.1. Sampling

The distribution of anoles in Martinique is more or less continuous. Hereafter we use the term “population” to refer to a discrete sampling site, not to genetically or ecologically differentiated entities. We sampled three distinct transects in northern Martinique (Figure 1). First, the “lineage transect” (transect II in [18]) consisted of eight populations sampled in similar habitat (transitional forest) across the lineage boundary between the North-West (NW) and the Central (C) lineages. Second, the “habitat transect” (populations 1 to 6 from transect IV in [18]) consisted of six populations sampled within the NW lineage and across an ecotone between coastal scrubland and montane rainforest. Finally, the “combined transect” (new data) consisted of eight populations sampled across the lineage boundary between the NW and the C lineages and across the ecotone between coastal scrubland and montane rainforest. For each population, tail tips from 48 individuals were sampled, and quantitative traits measurements (see below) were collected on ten adult males.

2.2. Genetic Analyses

For each individual, genotypes were scored at nine microsatellite loci (AAE-P2F9, ABO-P4A9, AEX-P1H11, ARO-HJ2 [34], ARO-035, ARO-062, ARO-065, ARO-120 [35], and ALU-MS06 [36]). Mitochondrial DNA lineage was inferred by amplifying the complete Cyt-b gene, digesting the PCR product with the restriction enzymes SspI and DraI (New England Biolabs) whose cutting pattern allows to distinguish the different lineages present in this species. More details on these molecular techniques can be found in [18]. Only individuals successfully genotyped at least at four loci were kept for analyses. Three sets of analyses were conducted.

(i) Admixture Analysis. We estimated the genetic structure within each transect using the Bayesian clustering method implemented in STRUCTURE v2.3.2 [37]. First, we estimated the most likely number of clusters (k) by running the analysis with k ranging from 1 to the real number of populations. For each k value, the analysis was run 10 times using the admixture model, with a burn-in of 100,000 steps for a total run length of 500,000 steps. The optimal number of clusters (K) was inferred following the method outlined in [38], choosing the number of clusters corresponding to the highest rate of change of the log probability of data between successive K values. We considered an individual as admixed if less than 80% of its genome was assigned to a single cluster. We also estimated the standardised pairwise Fst (Fst’ [39]) between adjacent populations along each transect using RecodeData v0.1 [39] and FSTAT v2.9.3 [40].

(ii) We then estimated the influence of different potentially important factors on the genetic structure using the method implemented in GESTE v2.0 [41]. GESTE employs a hierarchical Bayesian framework to estimate population specific Fst (representing the differentiation of a given population relative to all other populations) and uses a generalized linear model in order to test the contribution of biotic or nonbiotic factors to genetic structuring. For our analyses we used the default settings (sample size of 10000, thinning interval 20), and in accordance with guidelines we allowed ten pilot runs to estimate means and variances for the required input parameters [41]. The analysis was repeated five times for each transect to ensure that the results were consistent. Three different factors were considered. First 19 bioclimatic variables were obtained for each population from the WorldClim database. These variables were annual mean temperature, mean diurnal range, isothermality, temperature seasonality, maximum temperature warmest month, minimum temperature coldest month, temperature annual range, mean temperature wettest quarter, mean temperature driest quarter, mean temperature warmest quarter, mean temperature coldest quarter, annual precipitation, precipitation wettest month, precipitation driest month, precipitation seasonality, precipitation wettest quarter, precipitation driest quarter, precipitation warmest quarter, and precipitation coldest quarter. A principal component analysis (PCA) of the log transformed variables was performed using the ade4 package [42] in the R environment [43]. The first component explained more than 75% of the variation and was used as a composite bioclimatic variable to describe the environment among each site. This first axis described the variation from hot and seasonally dry coastal sites to cool and wet montane sites. Second, the geographic connectivity of the populations was estimated as described in [44] to include a measure of geographic isolation for each individual population (mean geographical distance between a given population and all other populations). Third, each site was assigned to a lineage (lineage transect), a habitat (habitat transect), or both (combined transect) as follows. For the lineage transect, sites 1–3 were assigned to the NW lineage, and sites 4–8 were assigned to the C lineage. For the habitat transect, sites 1–3 were assigned to coastal habitat, and sites 4–6 were assigned to montane habitat. For the combined transect, sites 1–3 were assigned to NW/coastal group, while sites 4–8 were assigned to C/montane group.

(iii) Transects Comparison. To determine if the genetic structure was different among transects, we computed and plotted the pairwise Fst’ values between populations on each side of the contact zone within each transect. This allowed us to compare the level of across habitat/lineage differentiation among transects. Because the observations are nonindependent and hence violate the assumption of both parametric and nonparametric tests, these plots provide qualitative information that has not been tested statistically.

2.3. Combined Analysis of Genetic and Quantitative Data

We conducted a modified version of the Discriminant Analysis of Principal Components (DAPC, [45]) to describe the global structure of populations within each transect using both genetics and morphology.

First different morphological characters were recorded as previously described [22, 46]. These included body dimensions (jaw length, head length, head depth, head width, upper leg length, lower leg length, dewlap length), scalation (number of postmental scales, scales between supraorbitals, ventral scales, and dorsal scales), colour pattern (number of dorsal chevrons, chevron intensity, occipital “A” mark, black dorsal reticulation, white spots), and trunk background colour (percentage red, green and blue hue on the posterior trunk). These were combined to six independent hues (UV 330–380 nm, UV/violet 380–430 nm, blue 430–490 nm, green 520–590 nm, yellow/orange 590–640 nm, and red 640–710 nm) extracted from the spectrum of the anterior and posterior dewlap by a multiple-group eigenvector procedure [23]. This combined data set was subjected to a principal component analysis (PCA) using the ade4 package in R.

Then, the microsatellite data were also subjected to a PCA using the adegenet package [47] in R. Finally, the components from the genetic (14 components) and quantitative (4 components) datasets were combined and subjected to a linear discriminant analysis using the MASS package [48] in R, using the population of origin as the grouping factor.

3. Results

3.1. Lineage Transect

Two genetic clusters were identified in this transect, with individuals from populations 1 to 3 being assigned in majority to the first cluster and individuals from populations 4 to 8 to the second cluster (Figure 2(a)). This separation corresponds to what was observed with mitochondrial DNA, populations 1 to 3 being in majority from the NW lineage while populations 4 to 8 are in majority from the C lineage [18, 46]. However, the proportion of admixed individuals is very high and relatively constant all along this transect (between 35 and 56%), and there is no trend in the pairwise Fst’ that vary between 0.072 and 0.101 (Figure 2(a)). Morphological and genetic data analysed separately (Figure S1) show the same trend that is magnified in the combined analyses. We observe a slight differentiation between populations 1 to 3 and populations 4 to 8, but there is an overlap between these groups (Figure 3(a)). The population specific Fst estimated with GESTE (Table 1) that represents the level of differentiation of one population relative to all other [41] are very low (ranging from 0.008 to 0.013, Table 1) and do not correlate with any of the factors investigated (geographic connectivity, environment, lineage), underlining again the lack of genetic structure along this transect (Table 2).

3.2. Habitat Transect

Two genetic clusters were also identified in this transect (cluster 1: populations 1–3; cluster 2: populations 4–6, Figure 2(b)). This genetic division corresponds to the habitat division, the ecotone being situated between populations 3 and 4 [18]. The admixture rates are globally lower than on the lineage transect (range: 15–55%). It is relatively low at the two opposite sides of the transect (29% in population 1 and 15% in population 6) and gets higher in the centre of the transect, around the ecotone (55% in population 3 and 45% in population 4). Here again, there is no obvious trend in the pairwise Fst’ that vary between 0.059 and 0.112 (Figure 2(b)). The combined dataset shows a much higher level of variation than on the lineage transect (almost twice as high), and a marked differentiation between populations 1-2 and 4–6, with population 3 being intermediate (Figure 3(b)). The population specific Fst are higher than in any site of the lineage transect (ranging from 0.0148 to 0.0522, Table 1) confirming the higher genetic structuring along this transect. The best model to explain this genetic structure only incorporates a constant, but the second and third best models have a nonnegligible probability and incorporate, respectively, the habitat type and the bioclimatic variable (Table 1); these two factors have a combined probability of 0.175 and 0.134, respectively (Table 2).

3.3. Combined Transect

Here again two genetic clusters were identified (cluster 1: populations 1–3; cluster 2: populations 4–8, Figure 2(c)), but with a much higher differentiation. This genetic division corresponds both to the habitat and lineage boundaries, the ecotone being situated between populations 3 and 4, as well as the lineage boundary (Figure 4). On this transect, the admixture rate is much lower (range 4–27%) than on the other transects. It is somewhat higher in population 4, situated at the lineage boundary and at the ecotone, but even in this population, it is lower than in any of the lineage transect’s populations and than in most of the habitat transect’s populations (Figure 2(c)). There is also a marked increase of the pairwise Fst’ at the contact zone, with a value of 0.169 between populations 3 and 4, while the other values are much lower (range: 0.052–0.0760, Figure 2(b)), suggesting the existence of a barrier to gene flow at the lineage/habitat boundary. The combined dataset shows a similar level of variation than in the habitat transect, with a strong differentiation between coastal (1-2) and montane (5–8) populations, with the populations 3 and 4 being intermediate. Along this transect, the population specific Fst are similar or higher to what is observed in the habitat transect (Table 1). The genetic structure is significantly associated with the habitat type/lineage factor, and the second best model includes the bioclimatic variable with a nonnegligible probability (Table 2).

4. Discussion

In this paper, we describe the differentiation between two lineages of Anolis roquet living in contrasting habitats and coming into secondary contact at the ecotone between these habitats. As a comparison, we reanalysed two previously published transects, a “lineage transect” where two lineages meet within a same habitat, and a “habitat transect” where different populations of the same lineage live in contrasting environments and are in contact at the ecotone between these habitats. This design allowed to investigate the effects of allopatry, habitat, and of the combination of these two factors on anole population differentiation and speciation.

A marked difference could be observed in the population structure along these three transects (Figure 5), suggesting that they are at different stages of the speciation continuum [49, 50]. According to Hendry et. al [49], four stages can be distinguished along this continuum: “(1) continuous variation within panmictic populations, (2) partially discontinuous variation with minor reproductive isolation, (3) strongly discontinuous variation with strong but reversible reproductive isolation and (4) complete and irreversible reproductive isolation.” The weak structure observed in the lineage transect, and the high admixture level correspond the State 1 of the continuum. This pattern is similar to what was observed in a previous study [18] for all but one (transect I in [18]) transects sampled across lineages (transects II, III, IV, V, VI, VII, VIII in [18]). In the habitat transect, the variation is clearly bimodal, but the high level of admixture in the contact zone between the two habitats suggests that the reproductive isolation is still limited between the ecotypes, corresponding to a stage between State 2 and State 3 of the continuum. Here again, this pattern is similar to the other transect sampled across habitats (transect III in [18]). Finally, the combined transects present strongly discontinuous variation, but there is still a significant level of admixture at the contact zone suggesting that the reproductive isolation is not complete despite a strong reduction in gene flow. This would correspond to State 3 of the speciation continuum. The situation of this last transect, where the ecotone and the lineage boundary overlap, is unique on this island, and hence this result could not be replicated.

In terms of association between observed genetic structure and biotic and abiotic factors, no significant factors were detected on the lineage transect, all the models other than the one incorporating only a constant having a very low probability. To the contrary, for the combined transect both the habitat/lineage category and the bioclimatic composite variable are the best at explaining the genetic structure, suggesting the important driving force of the environment into population differentiation. The situation in the habitat transect is not so clear. The best model only incorporates a constant, but the combined probability of the habitat type and the environment factor are not negligible as they were on the lineage transect. The lack of significance of these factors could be due to the lack of power of this method when the number of populations analysed is low. Indeed, simulations showed that this method failed to identify the true model with five or seven populations [41]. However, the observed trend suggests that environment factors may play a role in the genetic structuring of the populations along this transect too.

As expected, we find the strongest signal of divergence in morphological data that are known to react quickly to selection [5154] (Figure S1). The genetic data is similar in pattern but smaller in magnitude. Recent simulation studies suggest that neutral markers, that is, markers that are not linked with the traits under selection, are not very sensitive to detect ecological speciation [55, 56]. The authors conclude that this would lead to false negatives (failure to detect ecological speciation) rather than false positives. Taking this into account, and the fact that we found a clear signal of gene flow reduction both in the habitat and the combined transect in accordance with an extreme environment gradient, the divergence we demonstrate in this study with these neutral markers is undoubtedly considerably less than the divergence of the traits and loci under selection.

Since the habitat along the lineage transect is very homogeneous we do not expect ecological speciation to currently play a role in this area. However, this does not mean it has always been the case. When the populations from the two lineages were isolated on different proto-islands, it is possible that they were submitted to different environmental pressures. In such a situation, several studies have demonstrated that isolated populations can rapidly evolve partial reproductive isolation as a byproduct of local adaptation [5760]. Such a mechanism could also explain the differences observed between the three transects. When populations from the C and NW lineages came back into secondary contact, they did so in two very different contexts: either along a very sharp environmental gradient (combined transect), where any preexisting reproductive isolation could be maintained or strengthened by current disruptive selection, or within an homogeneous habitat (lineage transect) where any preexisting reproductive isolation may have been lost. Breakdowns of reproductive barriers associated with ecological changes have been recently described in various species [61, 62].

The relative role of geography and ecology in speciation remains a subject of debate. The main discussion relates to whether or not ecological factors can drive speciation in the presence of gene flow. Several convincing empirical studies suggest that it is indeed possible to reach full speciation in sympatry by ecological speciation (e.g., [6367]), while others emphasize the combined role of historical and ecological factors in shaping species diversity [6870]. For instance, in Trinidadian guppies, ecological speciation played a role in premating isolation either in allopatry (byproduct speciation) or in parapatry (to avoid maladaptive matings) [69]. For Martinique anoles, it is clear that for the populations that reached secondary contact very little evidence of the role of geographic isolation exists, while ecological factors seem to play a more important role. However, without conducting mate-choice experiments, it is not possible to determine conclusively the role of ecological speciation in allopatry.

Despite the large number of studies on the Caribbean anole radiation, very little is known about the factors at the origin of their diversity. Several papers have described the diversity patterns, demonstrating a correlation between island size and species diversity in large islands, while no speciation events were recorded on islands below a threshold size [13]. Recent work demonstrated that within island diversity could be the result of ecological opportunities and that net speciation rate decreased with time as opportunities decreased to reach an equilibrium at the island carrying capacity [15]. However, these studies do not explain the mechanisms at the origin of these diversity patterns and only propose several hypotheses to explain these observations.

In line with previous studies, we demonstrate that indeed environmental factors have a strong effect on the genetic structure of Martinique anole populations, with a reduction of gene flow at the ecotone between coastal and montane habitats and that it does not seem sufficient to lead to full reproductive isolation. Furthermore, we refine our understanding of divergence in these anoles by demonstrating that when allopatric lineages come into secondary contact on an ecotone, the differentiation is much stronger, with a significant reduction in gene flow. The absence of replication of the combined transect does not allow the generalisation of these findings, but we hypothesize that to reach full speciation anoles need first to evolve in allopatry and then come into secondary contact in an area where divergent natural selection will allow them to stay separate and further reinforce their divergence. Because this combination of factors is more likely to be found on large islands, it could be the mechanism at the origin of the species-area relationship observed in Caribbean anoles.

Acknowledgments

The authors wish to thank the DIREN Martinique for providing permits to collect samples, and W. Grail for assistance with molecular work. This work was funded by a Marie Curie Intra-European Fellowship to Y. Surget-Groba, a BBSRC grant to R. S. Thorpe, and a NERC studentship to H. Johansson.

Supplementary Materials

Supplementary Figure 1 presents the separate morphology and genetic analyses for the three transects.

  1. Supplementary Material