Although our primary goal was to assess whether A. allisoni and the eastern and west-central populations of A. porcatus are phenotypically distinct, we included all described species from the A. carolinensis subgroup to evaluate distinctness between Cuban green anoles and all the other members of the subgroup. We obtained new phenotypic data from 80 specimens representing all eight described species of the Anolis carolinensis subgroup (Table S1). Our new dataset included 10 A. allisoni from across its range in central Cuba (8 ♂, 2 ♀), 10 A. brunneus from the eastern Bahamas (all ♂), 9 A. carolinensis from mainland North America (8 ♂, 1 ♀), 12 A. porcatus from across western Cuba, including representation from western and central Cuba (5 ♂, 7 ♀), 10 A. porcatus from across eastern Cuba, including representation from all provinces of this region (9 ♂, 1 ♀), 6 A. fairchildi from Cay Sal (4 ♂, 2 ♀), 10 A. longiceps from Navassa Island (8 ♂, 2 ♀), 3 A. maynardii from the Cayman Islands (2 ♂, 1 ♀), and 10 A. smaragdinus from the western Bahamas (all ♂; Table S1). Specimen availability influenced unbalanced sex sampling; for some species, only males were accessible in sufficient numbers. Although sexual dimorphism exists in green anoles, with males being larger, with larger snouts, and more robust overall bodies than females (Rodríguez Schettino 1999; Losos 2009), our primary focus is on population-level differences rather than examining dimorphic traits.

From each specimen, we recorded a total of 14 variables, most of which are classic meristic and mensural traits that are widely used for species diagnoses in anoles in general and in species of the A. carolinensis subgroup (Garman 1888; Ruibal and Williams 1961; Pérez-Beato 1996; Köhler 2014; Les and Powell 2014). Our variables included three mensural traits (overall body length measured from the tip of the snout to the vent [SVL], head length, and head width), nine meristic variables that quantified numbers of particular types of scales or scale rows (numbers of postmentals, postrostrals, supralabials, infralabials, loreals, loreal rows, scales between interparietal and interorbital [IP-IO], temporals, and lamellae of the fourth toe), and two categorical variables that we identified as showing some degree of variation within and among closely related species during the course of our work (sublabials smooth or keeled and caudal rings keeled or heavily keeled). We did not conduct analyses on coloration because our study was based on preserved specimens, but we did examine photos from 633 iNaturalist records identified as A. porcatus to assess color variation in western and eastern A. porcatus.

We conducted three separate analyses of these data that assessed different types of phenotypic distinctness. We first asked if we could identify non-overlapping diagnostic traits capable of reliably distinguishing the two populations, a clear phenotypic diagnosability. The identification of one or more such diagnostic phenotypic traits is considered highly desirable, if not strictly required, when describing new species. To identify non-overlapping traits, we simply examined the distribution of trait values scored across the A. carolinensis subgroup and asked if the values for any individual trait were non-overlapping.

We next asked if Cuban green anoles could be distinguished by significant differences in individual trait values, even in cases where trait values overlapped between the two populations. While non-overlapping differences provide stronger evidence for divergence, statistically significant differences indicate that the populations are, on average, distinct in their trait values. This suggests a consistent, measurable pattern of differentiation that may reflect underlying genetic or ecological divergence, even if individual trait distributions overlap. We tested for differences in individual trait values across the entire A. carolinensis subgroup by conducting Analyses of Variance (ANOVA) with Tukey’s Honest Significant Difference (HSD) as a post-hoc test. Prior to these analyses, we tested for normality and homogeneity of variances and, in cases where these standard ANOVA assumptions were violated, we performed a Kruskal-Wallis test and Dunn’s test with Bonferroni correction as a post-hoc test instead. Prior to conducting all statistical analyses, we filtered our dataset to eliminate size outliers because body proportions can change over time with body size and sex, potentially affecting mensural variables as opposed to meristic variables, which remain constant (Butler and Losos 2002; Losos 2009; Sanger et al. 2012). We removed size outliers using the package “dplyr” (Wickham et al. 2023), which calculated the 15th and 85th percentiles of SVL to determine the interquartile range (IQR) and then identified lower outliers as values below a lower bound defined by the 15th percentile minus the IQR. This process resulted in the removal of four small males and 13 females, leaving us with a new dataset consisting of 64 male individuals representing the eight described species of the A. carolinensis subgroup (9 A. allisoni, 9 A. brunneus, 8 A. carolinensis, 4 A. fairchildi, 7 A. longiceps, 2 A. maynardii, 10 A. smaragdinus, 6 west-central A. porcatus, and 9 eastern A. porcatus). We conducted our univariate statistical analyses using the functions “aov,” “shapiro.test,” “bartlett.test,” and “TukeyHSD” for ANOVA, normality testing, testing for homogeneity of variances, and for post-hoc comparisons, respectively. For variables that had a non-normal distribution or different variances we conducted a Kruskal-Wallis test using the “kruskal.test” function from base R and the “dunn.test” function from the “dunn.test” R package (Dinno 2024) for the Dunn’s test with Bonferroni correction.

Third, we asked if green anoles populations could be distinguished based on their entire suites of phenotypic data by conducting multivariate statistical analyses. We specifically asked (1) if they form non-overlapping clusters in multivariate statistical space, and (2) whether they could be accurately classified based on multivariate data. As with statistical differences between individual trait values, such multivariate differentiation implies divergence and isolation of distinct species, but generally does not permit reliable species-level identification of individual specimens. For the clustering analyses, we conducted a Principal Components Analysis (PCA) using the R base function “prcomp.” PCA reduces dimensionality by linearly transforming the data into component axes that maximize the variance between samples (Lever et al. 2017). Prior to the PCA calculations, we centered and scaled the data to standardize the variables and ensure that each variable contributed equally to the analysis regardless of their original units or magnitudes. We then calculated PC scores and plotted these to visualize clusters in the multivariate phenotypic space and asked whether they overlapped or not.

For the classification analysis, we used Discriminant Analysis of Principal Components (DAPC) to ask if our multivariate phenotypic data could distinguish predefined groups corresponding with the previously diagnosed species and genetically distinct populations. DAPC works by identifying the principal components that best discriminate pre-defined groups (a priori species classification in this case). We assessed the ability of discriminant analysis to distinguish species using a cross-classification table generated with a jackknife resampling method (Jombart et al. 2010). This table shows the number of individuals accurately assigned to their species, as well as with which species each individual was misidentified. This pairwise comparison provided a quantifiable metric for assessing species similarity and the success of efforts to identify species based on multivariate data. If the two species were perfectly phenotypically distinguishable, we expected no misclassifications between them. The ability to distinguish one species could be assessed by dividing the number of correct classifications by the total number of observations. We retained the first 10 principal components to compute the first 5 discriminant functions. To conduct the DAPC we used the “dapc” function from the “Adegenet” R package (Jombart 2008).

Geographic isolation is often used as a proxy for species differentiation because it limits or prevents gene flow between populations, which can lead to genetic, ecological, and behavioral divergence over time (Liu et al. 2019). Molecular data suggest that west-central A. porcatus occupies the western two-thirds of the island, whereas the eastern population is restricted to the eastern third (Glor et al. 2004). Although A. porcatus has long been assumed to have an island-wide distribution (Rodríguez Schettino 1999; Rodríguez Schettino et al. 2010, 2013), superficial inspection of the available distributional records suggests that a real distributional gap might exist between the eastern and west-central populations of A. porcatus. To test for geographic discontinuity between eastern and west-central A. porcatus, we generated a new database of occurrence data and mapped the distribution of east and west-central A. porcatus. We also included the other Cuban species, A. allisoni, since it is also present in Cuba and is more closely related to eastern A. porcatus than to west-central A. porcatus (Glor et al. 2004). We gathered data from three sources: our own original field observations, the published literature, and the online Global Biodiversity Information Facility database (GBIF; Ruibal and Williams 1961; Pérez-Beato 1996; Glor et al. 2004, 2005; Rodrı́guez Schettino et al. 2013; Torres and Acosta 2014; Cajigas Gandia et al. 2018; Rodríguez-Cabrera et al. 2022; Rodríguez-Cabrera and Torres 2023). We downloaded the occurrences from GBIF using the “rgbif” package (Chamberlain and Boettiger 2017; Chamberlain et al. 2024). Many of the records in GBIF correspond with vouchered specimens identified by experts and permanently housed in biodiversity collections; thus, we trusted those records. Other GBIF records are from the iNaturalist database, which includes data on species observations from non-experts. We personally verified species-level IDs from photographs for all iNaturalist records included in the GBIF database using diagnostic traits (see Introduction).

To test if Cuban green anoles occupy distinct abiotic environments, we combined our compiled dataset of locality records with climatic data from GIS layers to evaluate niche overlap. We cleaned our dataset using an established pipeline for removing records with uncertain or missing coordinates (Cobos et al. 2018). We also conducted spatial thinning by randomly removing records that were less than 5 km away from other records of the same species to avoid model overfitting derived from the presence of multiple records in the same pixel in the raster variables (Anderson 2012). We implemented the spatial thinning in the R package “spthin” (Aiello-Lammens et al. 2015). All coordinates are in the World Geodetic System 1984 (WGS84).

We obtained current scenarios of 19 climatic variables at resolution of ~1 km from the WorldClim database, version 1.4 (https://www.worldclim.org). We excluded the four variables that could contain artifacts resulting from combining temperature and precipitation data (mean temperature of wettest quarter, mean temperature of driest quarter, precipitation of warmest quarter, and precipitation of coldest quarter) (Escobar et al. 2014). We then trimmed environmental layers to the Cuban Archipelago and conducted a Principal Components Analysis (PCA) on the fifteen remaining WorldClim variables to avoid model overfitting due to variable multicollinearity. PCA reduces the dimensionality of the data by creating a new dataset of uncorrelated principal components that reflect environmental variation, and which can then be used as inputs for downstream analyses (Zhang and Castelló 2017; Chan et al. 2022). We ultimately kept the first five PCs, which explained 98.7% of the variance, and we masked them to a hypothesis of accessible area for each species (see details below).

To characterize the suitable environments of A. allisoni, west-central A. porcatus, and eastern A. porcatus, we built niche models from our trimmed and cleaned locality dataset and the first five environmental PCs. We used ellipsoidal envelope ecological niche models (Farber and Kadmon 2003), constructed following Nuñez-Penichet et al. (2021) and implemented in the R package “ellipsenm” (Cobos et al. 2022). We excluded 5% of marginal data when fitting ellipsoids to avoid potential outliers. We produced 1000 replicate ellipsoid envelopes based on separate subsamples that included 75% of the occurrence data to include the variability contained in the data. We obtained the final models for each species by averaging the centroid and covariance matrices across all 1000 replicates.

To test if the Cuban green anoles occur in significantly different environments, we employed two complementary approaches. First, we conducted pairwise niche comparisons using a test based on ellipsoid envelopes, measuring the amount of overlap between ellipsoids corresponding to each species (points of an environmental background are checked to be inside or outside the ellipsoid of species A, B, or both) (Nuñez-Penichet et al. 2021). To measure overlap, we considered the background-overlap approach (Nuñez-Penichet et al. 2021). To this end, we defined an “accessible area” for each green anole lineage in Cuba using the R package “grinnell” (Soberón and Peterson 2005; Barve et al. 2011; Machado-Stredel et al. 2021). To define accessible area, we arbitrarily allowed 100 dispersal events with 25 individuals for all Cuban lineages. We then used the accessible areas to test for statistical significance of niche overlaps by creating 1000 ellipsoids with randomly sampled points from the background of each species and comparing these results with the observed overlap value. In this context, the null hypothesis is that the two ellipsoids fitted to actual observations overlap at least as much as random-data ellipsoids. To reject this hypothesis, the observed value of overlap must be as extreme or more extreme than the lower confidence limit (5%) of the overlap values derived from comparing random-data ellipsoids derived from each species background. We performed these analyses using the routines implemented in the “ellipsenm” package.

Second, in addition to comparing niches using principal components from all variables together, we conducted a univariate approach to separately evaluate each environmental variable between the three lineages of Cuban green anoles. We conducted ANOVAs to test for differences between lineages by variable, and Tukey’s HSD as post-hoc analyses to find which lineages were different. Statistical differences among the three Cuban green anoles would support an ecological distinctness hypothesis. We ran these analyses in a similar fashion to those of the phenotypic variables.

To test the hypothesis that A. allisoni, eastern, and west-central A. porcatus are distinct lineages, as previously supported (Glor et al. 2004), we reconstructed phylogenetic trees with molecular data for the first time representing all species of the A. carolinensis subgroup, including new samples of A. porcatus from across its range. We evaluated the expectations that eastern and west-central A. porcatus were reciprocally monophyletic and that they were not sister clades.

We sampled 105 individuals, three as outgroup taxa, and 102 representing all eight species from the A. carolinensis subgroup: A. allisoni (21 individuals), A. brunneus (2), A. carolinensis (14), A. fairchildi (3), A. longiceps (2), A. maynardii (14), A. porcatus (36), and A. smaragdinus (14). We also included samples from across the distributions of the four widespread species found on either mainland North America (A. carolinensis), Cuba (A. allisoni and A. porcatus), or the Great Bahama Bank (A. smaragdinus). Within the mainland North American species, A. carolinensis, we included sampling from each of five major mitochondrial lineages (Campbell-Staton et al. 2012). For Cuban A. porcatus, we sampled broadly from across the island, including representatives of both the eastern and west-central Cuban lineages previously identified as genetically distinct (14 and 22 samples, respectively) (Glor et al. 2004, 2005). For A. allisoni, we included 18 samples from across this species’ range in central Cuba as well as four specimens from the type locality on Roatan Island off the coast of Honduras. For A. smaragdinus from the Great Bahama Bank, we included 14 samples from six Bahamian islands (New Providence, Great Ragged, Staniel Cay, Long Island, and South Bimini). For A. fairchildi, which is only known from two small western Bahamian islands—Cay Sal Island and Cotton (= South Anguilla) Cay—we included three samples from Cay Sal (Reynolds et al. 2018). For the third Bahamian species from Crooked Bank in the eastern Bahamas, A. brunneus, we included two samples from Crooked Island. We included 10 samples of A. maynardii from the Cayman Islands—four from Little Cayman and six from Cayman Brac (introduced; Herrel et al. 2011). Finally, we included two samples of Anolis longiceps, which is endemic to Navassa, an uninhabited island 74 km west of Hispaniola.

We obtained genomic DNA from either a piece of liver, thigh muscle, or tail preserved in 95% ethanol and then frozen to between –20 and –80°C, or from flash-frozen liver. We first digested tissue subsamples overnight in 280 μl lysis buffer and 20 μl of proteinase K before extracting DNA using the Promega Maxwell® RSC Tissue DNA Kit (catalog number AS1610) with a Maxwell® RSC extraction robot. We obtained mtDNA sequences, including 18 from GenBank and 87 newly sequenced individuals (Table S2). We amplified a ~1200 base pair (bp) fragment of mtDNA that includes the complete sequence for NADH dehydrogenase 2 (ND2) and tRNA_Trp, as well as partial sequence for tRNA_Ala, using established forward (Macey et al. 1997) and reverse (H5730; Glor et al. 2004) primers. The PCR protocol consisted of a hot start for 180 s at 95°C followed by 34 cycles of 95°C for 30 s, 59°C for 30 s, and 72°C for 60 s, and ultimately 72°C for 300 s after the cycles were complete. We used PuReTaq Ready-To-Go (RTG) PCR Beads to prepare the reactions. To each tube with a bead, we added 1.25 µL of each primer, 2 µL of DNA extract, and 20.5 µl of water for a total volume of 25 µL per tube/sample. We checked whether we had amplified the desired sequence by inspecting fragment sizes of PCR products by 1% agarose gel electrophoresis before sequencing the entire 1200 bp fragment with the same primers used for amplification at Functional Biosciences. We assembled and aligned sequences using MAFFT (Katoh et al. 2009) in Geneious.

We reconstructed phylogenetic trees to assess the relationships across green anoles, with emphasis on eastern and west-central A. porcatus. We generated Maximum Likelihood (ML) phylograms using IQ-TREE v1.6 (Nguyen et al. 2015) from our mitochondrial concatenated dataset. We determined the best-fit substitution model for the dataset via ModelFinder, implemented within IQ-TREE (Kalyaanamoorthy et al. 2017) and performed a partitioned analysis by partitioning according to codon position within the ND2 marker and by marker (ND2, tRNA_Trp, and tRNA_Ala). We calculated branch support with 10,000 bootstrap replicates using the Ultrafast Bootstrapping algorithm (Hoang et al. 2018).

Overall, our phenotypic analyses largely confirm that the Cuban green anoles are phenotypically diagnosable, but eastern and west-central A. porcatus are cryptic with respect to traditional mensural and meristic traits. Our assessment of non-overlapping differences did not recover any diagnostic traits. However, one of our two binary traits—presence or absence of keels on the sublabial scales—did unambiguously diagnose the west-central and eastern A. porcatus, with the west-central lineage having keeled scales and the eastern lineage having either smooth scales or keels extending only from the first to the fourth scale (Fig. 1). More broadly, we found an overall absence of non-overlapping diagnostic mensural and meristic traits across the entire A. carolinensis subgroup. Although our second binary trait—keelation of caudal rings—was not useful for distinguishing the Cuban populations, it did unambiguously diagnose the Bahamian species A. smaragdinus from all other species in the A. carolinensis subgroup because A. smaragdinus has tall keels whereas the rest of the species have markedly lower keels.

Figure 1. 

Ventral view of heads of representatives of west-central and eastern A. porcatus showing diagnostic keelation of sublabial scales. Eastern A. porcatus has either smooth (KU 259008) or weakly keeled scales up to the third or fourth scale (KU 55635) (curly bracket). In west-central A. porcatus all sublabial scales are heavily keeled. Photo credit: Dexter Reilly.

Although extensive overlap among the two populations for each of our meristic and mensural traits limits their usefulness in species diagnoses, univariate traits did show significant variation across the A. carolinensis species group. For example, the central Cuban A. allisoni was the largest species in all mensural variables whereas A. carolinensis and A. smaragdinus were the smallest. Overall, our univariate analyses showed significant variation in all our mensural and meristic variables across the A. carolinensis subgroup. Post-hoc analyses sorted four meristic variables—infralabials, supralabials, lamellae, and loreals—into three levels, and the remaining six—IP.IO, loreal rows, postmentals, postrostrals, and temporals—into only two levels, suggesting high phenotypic conservancy across the subgroup. The west-central and eastern populations of A. porcatus were never assigned to different levels (Tables 1, S3; Figs S1–S12).

Our multivariate analyses further supported overall phenotypic conservatism across the A. carolinensis subgroup, and a lack of phenotypic differentiation between the eastern and west-central populations of A. porcatus. Plots of our PCA data recovered phenotypic overlap among most members of the A. carolinensis subgroup, and particularly extensive overlap between the eastern and west-central populations of A. porcatus (Fig. 2). Across the entire A. carolinensis subgroup, our multivariate DAPC analyses correctly classified species 71.6% of the time (Table 2). Two of the nine populations in the subgroup—A. carolinensis and A. longiceps—had perfect prediction accuracy, and just one individual of one other species (A. maynardii) was ever misclassified as A. carolinensis, whereas no other species were ever misclassified as A. longiceps. At the other end of the spectrum, we saw the highest rates of classification error for the two eastern and west-central A. porcatus, each of which was most frequently misclassified as belonging to the other geographic population of A. porcatus (Table 2).

Figure 2. 

Visualizations of multivariate analyses with meristic variables including all described species of the Anolis carolinensis subgroup: Principal Components Analysis (left) and Discriminant Analysis of Principal Components (right). Points delimited by polygons are eastern (E) A. porcatus (top in inset) and west-central (WC) A. porcatus (bottom in inset). Inset photo credit: Tomás M. Rodríguez-Cabrera.

West-central and eastern Anolis porcatus are geographically isolated

We gathered a total of 1,538 occurrences for the two Cuban species currently assigned to the A. carolinensis subgroup but separating A. porcatus in eastern and west-central Cuba. After removing duplicates, we kept a total of 777 localities, 224 for A. allisoni, 373 for west-central A. porcatus, and 180 for eastern A. porcatus. Of these, 100 records were from the field, 567 were from the literature, and 110 were from online databases.

Anolis porcatus has been recorded from across most of Cuba. We recovered numerous records from west-central A. porcatus ranging from Pinar del Río Province in the west to Camagüey Province in the east, whereas eastern A. porcatus was distributed across Gramma, Holguín, Santiago de Cuba and Guantánamo Provinces (Fig. 2). We only identified a single record for A. porcatus from Las Tunas Province, which could not be verified because it was based on an unvouchered personal observation from Rodríguez Schettino et al. (2013). Aside from this record, a gap of an entire Cuban province exists between west-central and eastern populations of A. porcatus, with a 50-km separation between the closest localities. Distributional overlap is evident among the other Cuban green anoles: west-central A. porcatus has a large overlap in central Cuba with A. allisoni, and the latter has a much more limited overlap with eastern A. porcatus in eastern Cuba (Fig. 3).

Figure 3. 

Occurrence records of the three lineages of Cuban green anoles. Numbers denote administrative regions: Pinar del Río (1), Artemisa (2), La Habana (3), Mayabeque (4), Matanzas (5), Cienfuegos (6), Villa Clara (7), Sancti Spíritus (8), Ciego de Ávila (9), Camagüey (10), Las Tunas (11), Granma (12), Santiago de Cuba (13), Guantánamo (14), Holguín (15), Isla de la Juventud (16). Color intensity in occurrences is directly proportional to the number of times a given species has been reported for a given locality.

Anolis allisoni, west-central and eastern Anolis porcatus are ecologically distinct

After removing duplicates and conducting spatial thinning on our full occurrence dataset, the occurrence dataset used for our analyses included a total of 483 records: 229 for west-central A. porcatus, 150 for A. allisoni, and 104 for eastern A. porcatus.

The suitable areas reconstructed by our niche modeling exercise largely reflected each species’ actual distribution across Cuba. Eastern A. porcatus had the lowest proportion of suitable areas (mean = 0.79), mostly restricted to eastern Cuba, with patches in the western side of the island. Anolis allisoni and west-central A. porcatus had higher proportions of suitable areas (mean = 0.82 and 0.89, respectively), with wider geographical connectivity across most of western and central Cuba (Fig. 4). All three lineages had suitable areas extending beyond their current known distributions (Figs 3, 4). For example, the centrally distributed A. allisoni had highly suitable areas in westernmost Cuba and portions of the eastern end of the Cuban main island. Similarly, the west-centrally distributed A. porcatus had areas of high suitability in portions of eastern Cuba, although a gap of lower suitability existed in transitional areas between central and eastern Cuba (Fig. 4). Anolis allisoni and west-central A. porcatus had more similar suitable areas than either of them with eastern A. porcatus. Despite having a much more geographically restricted distribution than the other two lineages (Fig. 3), eastern A. porcatus had the largest niche breadth, with a mean ellipsoid volume of 194.58 units, followed by west-central A. porcatus (59.03 units) and A. allisoni (47.34) (Fig. 4).

Figure 4. 

Predictions of ecological niche models for Anolis allisoni, eastern A. porcatus, and west-central A. porcatus. Predicted values of suitability are presented in geographic (left panel) and environmental space (right panel). Black color indicates the values of suitability below a 5% omission threshold.

In both of our approaches for comparing environmental differences, we found support for our hypothesis that west-central and eastern Anolis porcatus are ecologically distinct, both of which are also distinct from A. allisoni. First, we obtained non-significant niche overlap proportions between all Cuban lineages: 0.69 for eastern A. porcatus×west-central A. porcatus, 0.62 for eastern A. porcatus×A. allisoni, and 0.88 for west-central A. porcatus with A. allisoni (Fig. 5). Second, considering the environmental variables independently, all 15 variables were statistically different: nine were different between the three Cuban green anole lineages, two were not different for A. allisoni and eastern A. porcatus, two were not different for A. allisoni and west-central A. porcatus, one was not different for A. allisoni and either of the A. porcatus lineages but was different between them, and one was not different between A. porcatus lineages but it was for A. allisoni (Table 3).

Table 3.

Download as

CSV

XLSX

Univariate analyses of environmental variables. F is the ANOVA stadigraph. Different letters (a, b, and c) represent statistical differences between lineages (A = Anolis allisoni, E = eastern A. porcatus, WC = west-central A. porcatus).

Variable code	Variable meaning	F	P-value	A E WC
bio_1	Annual mean temperature	23.54	1.77×10-10	a b c
bio_2	Mean diurnal range	57.04	5.96×10-23	a b c
bio_3	Isothermality	20.00	4.54×10-09	a a b
bio_4	Temperature seasonality	523.84	2.14×10-121	a b c
bio_5	Max temperature of warmest month	72.58	2.89×10-28	a b c
bio_6	Min temperature of coldest month	10.74	2.73×10-05	a a b
bio_7	Temperature annual range	181.96	1.54×10-59	a b a
bio_10	Mean temperature of warmest quarter	40.42	5.96×10-17	a b c
bio_11	Mean temperature of coldest quarter	17.02	7.19×10-08	a b b
bio_12	Annual precipitation	20.56	2.72×10-09	a b c
bio_13	Precipitation of wettest month	13.89	1.36×10-06	a b b
bio_14	Precipitation of driest month	89.26	1.10×10-33	a b c
bio_15	Precipitation seasonality	124.75	2.35×10-44	a b c
bio_16	Precipitation of wettest quarter	5.27	5.46×10-3	a ab b
bio_17	Precipitation of driest quarter	78.82	2.51×10-30	a b c

Figure 5. 

Representation of niche overlap between Cuban green anoles: niche overlap between Anolis allisoni and eastern A. porcatus (top), A. allisoni and west-central A. porcatus (center), and west-central A. porcatus and eastern A. porcatus (bottom) in the environmental space (left panel), with a significance test evaluating whether overlap is significant or not (right panel). CL: Confidence limits.

Anolis allisoni, west-central and eastern Anolis porcatus are distinct lineages

Anolis allisoni, west-central and eastern A. porcatus met our predictions of reciprocal monophyly, and non-sister clades (Fig. 6). We recovered a well resolved phylogeny with high support at most nodes from our Maximum Likelihood analysis. We identified three deeply divergent clades within the A. carolinensis subgroup, each of which included both Cuban and non-Cuban populations (Fig. 6). We refer to these three clades based on their geographic distribution on Cuba. The first clade, the “western clade”, included western and central populations of Cuban A. porcatus, A. carolinensis from mainland North America, and A. fairchildi from Cay Sal. The second clade, which we refer to as the “central clade,” included A. allisoni from across central Cuba and Roatan Island, and A. smaragdinus from the Great Bahama Bank. The third clade, which we refer to as the “eastern clade,” included eastern Cuban A. porcatus, A. brunneus from Crooked Bank, A. longiceps from Navassa, and A. maynardii from the Cayman Islands.

Figure 6. 

Geographic sampling, clade representative and phylograms of all species in the Anolis carolinensis subgroup. Left panel: Maximum Likelihood tree consensus resulting from concatenated mitochondrial sequences. Three deep clades are coincident with the relative distribution of lineages in Cuba—western (circle), eastern (square), and central (triangle). Right panel, top: Map with the distribution of ingroup sampling localities (except some within the US), with symbols shaped and colored according to the tree. Light gray areas represent shallow waters. Discontinuous lines represent faults that are associated with paleoisland boundaries. Right panel, bottom: Sample size by species/population are also shape/colored coded (depicting males, no photos available for two species). Photo credit: Tomás M. Rodríguez-Cabrera (A. allisoni and A. porcatus), Robert W. Henderson (A. brunneus), Rachel E. Cohen (A. carolinensis), RGR (A. fairchildi), and Jake M. Scott (A. smaragdinus).

Six of the eight previously recognized species are monophyletic; however, A. allisoni is not monophyletic only with respect to A. smaragdinus. Anolis porcatus, meanwhile, is also not monophyletic. One reason for this paraphyly is that A. porcatus is divided between two of the three main clades in the A. carolinensis subgroup; individuals from eastern Cuba are in the eastern clade and individuals from across central and western Cuba are in the western clade. Eastern A. porcatus is monophyletic, but is more closely related to A. allisoni, A. smaragdinus, A. maynardii, A. brunneus, and A. longiceps than it is to central or western A. porcatus. Additionally, central and western A. porcatus are not monophyletic because they are further divided into two geographically circumscribed clades in west-central Cuba, each of which also includes non-Cuban taxa. The North American mainland A. carolinensis samples are reciprocally monophyletic with central A. porcatus, whereas A. fairchildi is nested within western A. porcatus, a pattern reported in previous studies (Glor et al. 2004, 2005; Reynolds et al. 2018; Wegener et al. 2019).

Our phenotypic analyses showed that A. porcatus and A. torresfundorai sp. nov. —and indeed, all members of the A. carolinensis subgroup—are generally very difficult to distinguish from one another using the type of mensural and meristic characters that have traditionally been relied upon to distinguish species of anoles. Beyond the color variation that distinguishes some species (e.g., the blue head and torso of male A. allisoni and the largely gray body coloration of A. brunneus), phenotypic traits are strongly conserved across the A. carolinensis subgroup. Overlap across species exists in all the traits we quantified, but significant variation in the distribution of trait values across the group revealed between two and four classes. We did find that A. porcatus has keeled sublabial scales, smooth in all or weakly keeled in the first three to four in A. torresfundorai sp. nov., but further sampling will be required to determine if this diagnosis is consistent. We also found that having very dark paramedial stripes seems to be diagnostic of A. torresfundorai sp. nov. but this trait is not always present, and a more detailed analysis will be necessary to test for differences in coloration between A. porcatus and A. torresfundorai sp. nov. The overall phenotypic conservatism we see across the A. carolinensis subgroup may persist due to limited interactions between species, which reduces the potential for phenotypic divergence driven by selection against competing phenotypes. All non-Cuban species are the only trunk-crown anoles naturally occurring in their respective ranges. In contrast, in Cuba, where two pairs of species co-occur, character displacement and hybridization with sporadic hybrids have been documented—but not between A. porcatus and A. torresfundorai sp. nov. (Schoener 1977; Glor et al. 2004), suggesting selective pressures operating in overlapping areas of A. allisoni with A. porcatus or A. torresfundorai sp. nov. Character displacement is an evolutionary pattern observed in similar coexisting species that evolve divergence in traits due to selective pressures that reduce resource competition and/or hybridization (Pfennig and Pfennig 2009; Stuart et al. 2017). This pattern was observed in areas with A. allisoni and A. porcatus where the former diverged to a larger size and the latter to a smaller size. Selection driving phenotypic divergence that reduces competition or hybridization with less-fit offspring could explain the similarity between A. porcatus and A. torresfundorai sp. nov. and the differentiation of A. allisoni.