Plants of Euphorbia tithymaloides (Euphorbiaceae) are evergreen woody succulents that inhabit dry environments, from tropical deciduous forests to inland and coastal xeric scrublands and are often associated with calcareous soils. This is the most variable and widely distributed species of the Pedilanthus clade of Euphorbia, a mainly Mexican clade that stands out due to its markedly bilateral inflorescences that suggest hummingbird pollination in an otherwise insect-pollinated genus (Dressler 1957; Cacho et al. 2010). In contrast to the rest of the species in the Pedilanthus clade that exhibit rather narrow geographic ranges (Olson et al. 2005), the range of E. tithymaloides spans the whole Caribbean Basin, including Florida, México, Central America, the northern portion of South America, and practically all Caribbean islands (Dressler 1957; Olson et al. 2005; Cacho and Baum 2012).

Euphorbia tithymaloides exhibits a remarkable morphological diversity for a single species, which includes variation in leaf shape as well as in other foliar attributes such as leaf base and apex, the presence and morphology of a keel, and features related to venation and indumentum (Fig. 1). In an attempt to make sense of morphological diversity in E. tithymaloides, Dressler (1957) proposed eight subspecies in this system, seven of which are located in the Caribbean. Subspecies assignation is complicated because of character intergradation among subspecies, to the point that geography plays a fundamental role in Dressler’s proposal (i.e. subspecies E. tithymaloides subsp. bahamensis and E. tithymaloides subsp. jamaicensis; Supplementary file 1). That is, the provenance of an individual is a key element for taxonomic determination.

Studies examining the historical biogeography of E. tithymaloides based on morphological (Dressler 1957), molecular phylogenetics (Cacho and Baum 2012), and landscape genetic approaches (Cacho et al. 2019) all point to E. tithymaloides originating in Mesoamerica (in particular, an area between Mexico and Guatemala) and expanding its range towards the Caribbean along two geographic fronts (Fig. 2A). One of these fronts colonized the Greater Antilles and continued east towards the Caribbean Basin, and the other travelled south through Central America, then east through northern South America, and north towards the Lesser Antilles. The two fronts appear to reunite near the Virgin Islands, where they remain morphologically, ecologically, and genetically distinct (Cacho and Baum 2012; Cacho et al. 2019; Fig. 2B).

Molecular phylogenetic and landscape genetic approaches support E. tithymaloides subsp. tithymaloides as the most variable and geographically extended subspecies, with a continental geographic range that spans from Mexico through Venezuela. Euphorbia tithymaloides subsp. jamaicensis, E. tithymaloides subsp. smallii, E. tithymaloides subsp. bahamensis, E. tithymaloides subsp. parasitica, and E. tithymaloides subsp. angustifolia are supported as part of the Greater Antillean front, and E. tithymaloides subsp. padifolia as part of the Lesser Antillean one. Subspecies angustifolia and padifolia are the most recent subspecies, and the extremes of either biogeographic front, thus representing lineages with independent evolutionary trajectories (Fig. 2).

In E. tithymaloides, clear geographic patterns in floral morphology have been documented: the floral involucre has shortened from the centre of origin towards the Anegada Passage, and this has happened in parallel along both biogeographic fronts (Cacho et al. 2019). On the other hand, there is no evidence supporting geographic patterns in leaf morphology in relation to the Anegada Passage (Cacho et al. 2019). Dressler (1957) presents some data on leaf length/width ratio in two subspecies (i.e. angustifolia and padifolia), also in the context of discussing floral morphology, but he does not discuss the leaf data. Apart from these two instances, there is no study focusing on variation in leaf shape in E. tithymaloides throughout its range.

We assembled a collection of images of E. tithymaloides leaves (with scale) representing as wide a geographic range for the species as possible. We included our own images, taken during various field trips, as well as pictures we took from herbarium specimens (COL, HUA, MEXU) and images obtained from websites of individual herbaria (AAH, FLAS, FSU, GH, NY, P) or portals like the Global Biodiversity Information Facility (GBIF).

We excluded specimens without leaves, specimens with leaves that were damaged to the point of limiting our ability to measure them, specimens with locality data that would not allow georeferencing, and specimens whose native status was questionable. Subspecies assignation followed information on the specimen label, or when absent, we followed the subspecies key and descriptions by Dressler (1957).

We measured the following five leaf traits: blade maximum length (Lmax) and maximum width (Wmax), distance from the base to the point of maximum width (DbWmax), leaf area (A), and leaf perimeter (P), and calculated the following four metrics, which are described below: leaf aspect ratio, leaf area-perimeter ratio, obovate index, and circularity.

Leaf aspect ratio (Lmax / Wmax) – This is a robust metric that reflects the ratio between the two axes of an ellipse that is related to how round or elongated leaves are (Li et al. 2018). For leaves of similar values of Wmax, the leaf aspect ratio will be determined by Lmax, so that it will be larger for longer leaves.

Area-Perimeter ratio (A / P) – A metric that captures aspects related to leaf shape. A circular shape maximizes the area for a given perimeter. Given a constant perimeter, leaves with larger areas will have larger A/P values and be the most circular, and those with smaller areas will be less circular, either through elongation or the presence of lobing or dissection, and this metric will have smaller values.

Circularity (4π * (area / perimeter²)) – Another metric related to foliar shape, also based on area and perimeter. This metric is more sensitive to lobing and captures aspects of the degree of how elongated leaves are (Li et al. 2018) and rather than being redundant with the metric AP, it can complement it in certain groups, as shown in Vitis (Chitwood et al. 2014).

Obovate index (DbWmax / Lmax) – This metric aims at capturing leaf shape along an ovate-obovate axis. Along this axis, for leaves that have the same values for Lmax, this index will be determined by DbWmax (the distance between the leaf base and the point of Wmax), so that ovate leaves will have smaller values and obovate leaves larger ones.

Whenever possible, we measured and calculated the metrics described above for as many as three leaves per individual. Averages per individual were calculated prior to subsequent analyses. Lmax, Wmax, and DbWmax were measured in Geogebra Classic v.5 (Supplementary file 2.1). Area and perimeter were measured in ImageJ v.2.0 (Schindelin et al. 2015), which required some image processing as illustrated in Supplementary file 2.2.

When we did not have geographic (latitude/longitude) information, we geo-referenced images based on information in the label of the specimens, using tools in Google Earth or Google Maps. Geo-referenced data were curated with QGIS v.3.16 (QGIS 2021) and geo-referencing errors (in cities, bodies of water, etc.) were corrected on a case by case basis.

We extracted climate data associated with curated occurrence points for 19 climatic variables and elevation from WorldClim v.2.0 (Fick and Hijmans 2017) at a scale of 30 s, which roughly corresponds to 1 km² at the equator.

To evaluate the extent to which leaf morphology alone predicts current taxonomy, as a first approximation we used a combination of Kruskal-Wallis and posthoc tests on leaf traits given that several variables did not conform to normality and homoscedasticity (per Shapiro and Levene tests, respectively). We then used linear models (of the type: morphology ~ subspecies) to estimate the proportion of variation in foliar morphology that is captured by subspecies assignation. To evaluate the correspondence between morphology and geography or climate, we used linear models and implemented Mantel tests using matrices of morphological, climatic, and geographic Euclidean distances. Variables were transformed as shown in Supplementary file 3.1 to improve the normality of the data, and we verified that the nature of leaves from where images were derived (fresh vs dry) would not introduce a systematic bias in our analyses (Supplementary file 4).

We used Principal Component Analyses (PCA) to account for collinearity among variables and reduce the dimensionality of our data. We implemented this approach for both morphologic and climatic data (which is multivariate by definition). PCA was based on the correlation matrix (function princomp, cor = TRUE) to ensure data would be all at the same scale.

Because we found significant differences among subspecies and because subspecies are not randomly distributed across the landscape, to evaluate the effects of geography or climate as possible predictors of foliar morphology, we eliminated the subspecies effect by using the residuals of linear models by subspecies (syntax of model: variable ~ subspecies) as the response variable in all our analyses.

Patterns of geographic variation in foliar shape in E. tithymaloides across the Caribbean were evaluated using three approaches. First, we used linear models to assess variation in leaf shape (using morphological PCs) in relation to latitude and longitude. Then, based on the results and the loadings of the morphology PCs, we assessed latitude and longitude as potential predictors of specific variables reflecting leaf shape variation. For this, we focused on the variables with the highest loadings for every PC of morphology: A/P ratio and Wmax for PC 1, L/W ratio and Lmax for PC 2, and obovate index for PC 3. Finally, to evaluate if individuals that co-occur would share foliar morphology, driven by geographical proximity alone, we used Mantel tests (method = Pearson, permutations = 999), implementing the model D_morphologic ~ D_geographic.

Linear models to evaluate association between foliar morphology and climate focused on the first three PCs for both climate and leaf shape. As with geography, we used the results of these models to guide variable selection for subsequent models focusing on the variables with the highest loadings on relevant PCs (for both morphology and climate), which due to the nature of PCA have a lower probability of being collinear. We selected one variable per climatic PC, and two for every morphological PC (except for PC 3 which is mostly correlated with only one variable, see loadings below). The three climatic variables selected with this approach were: Mean annual temperature (BIO 1, PC 1), mean annual precipitation (BIO 2, PC 2), and temperature seasonality (BIO 4, PC 3), and the five morphological variables were the same as above (A/P ratio and Wmax for PC 1, L/W ratio and Lmax for PC 2, and obovate index for PC 3). To test the hypothesis that similarity in leaf morphology could be explained by the occupation of similar climatic envelopes (i.e. proximity in climatic space), we implemented Mantel tests as above. Climate data are known to be spatially autocorrelated (Portier et al. 2018), and this is true for our data (model: D_climatic ~ D_geographic; adj. R² = 0.494, p < 0.0001). To account for such autocorrelation, we implemented the following partial Mantel test: D_morphologic ~ D_climatic + D_geographic.

All analyses were implemented using base functions in R (R Core Team 2018) except for distance matrix calculations which were derived with functions of the R package vegan v.2.5-6 (Oksanen et al. 2018).

Our image database consisted of a total of 578 leaf images of E. tithymaloides specimens spanning across this species’ range in the Caribbean (Fig. 3), including specimens from the following countries or territories:

From a taxonomic standpoint, all seven Caribbean subspecies described by Dressler (1957) were represented in our database. The three subspecies with the largest ranges and evolutionary significance, i.e. angustifolia, padifolia, and tithymaloides, also had a noticeable larger representation in our dataset (with 107, 98, and 351 leaves measured from 41, 42, and 146 specimens respectively) than the other four subspecies whose ranges are more limited (number of images: subsp. bahamensis, 7; subsp. jamaicensis, 3; subsp. parasitica, 3; subsp. smallii, 9). Thus, these three subspecies are the focus of most of our analyses, especially those related to geography and climate. On average, the number of leaves per specimen measured for the three main subspecies was 2.55 for subsp. angustifolia, 2.39 for subsp. padifolia, and 2.4 for subsp. tithymaloides.

Euphorbia tithymaloides subsp. tithymaloides was the most variable of all subspecies, which would be expected given its significantly larger geographic range. However, for certain metrics (e.g. circularity and obovate index), subsp. angustifolia was the most variable (Supplementary file 3.2). In contrast, there is not one variable for which subsp. padifolia was noticeably the most variable, and it is rather common that it is the least variable subspecies in the system.

Our PCAs on leaf form reveal that foliar variation in E. tithymaloides is well represented by three principal components (PCs) that are related to how elongated and how obovate a leaf is, as well as leaf dimensions (Lmax and Wmax). The three first PCs capture 96.7% of variation in leaf morphology in E. tithymaloides (Supplementary file 5). Based on the loadings of the variables on the PCs (Supplementary file 5), PC 1 (which captures 64% of total variance) is mostly positively correlated with leaf shape as summarized by A/P ratio (also with Wmax, area, and perimeter). PC 2 (22% of total variance) represents how elongated a leaf is due to its positive correlation with leaf aspect ratio (L/W) and Lmax, and negatively with circularity. PC 3 (10.9% of total variance) represents how ovate is a leaf (inverse relationship with obovate index).

When visualizing the foliar morphospace of Caribbean E. tithymaloides, either with a taxonomic perspective (subspecies) or a geographic one (focused on the three main regions: Greater Antilles, Lesser Antilles, and Mainland), it is possible to appreciate a substantial correspondence between three main subspecies and the main geographic areas (Fig. 4). Individuals of subsp. angustifolia and those that occur in the Greater Antilles occupy areas of morphospace that are highly coincident, but not identical. Similarly, individuals of subspecies padifolia occupy a portion of the foliar morphospace of E. tithymaloides that overlaps with the one occupied by individuals of the Lesser Antilles, and the morphospace occupied by individuals of subsp. tithymaloides corresponds broadly to that occupied by individuals form the mainland. This might be expected, in part because geography is a quite important component of current taxonomy.

Our analyses reveal differences in foliar variation among subspecies of Caribbean E. tithymaloides. Differences are subtle where our sampling is shallow (subspp. bahamensis, jamaicensis, parasitica, and smallii) possibly due to limited statistical power. However, for better-represented subspp. angustifolia, padifolia, and tithymaloides, our analyses reveal clear differences among subspecies (Fig. 5; Supplementary file 6.1). While subspecies differ for all three morphology PCs, posthoc tests reveal that not all three subspecies differ in all three PCs: there are significant three-way differences for PC 1 and 2, but not for PC 3, in which subsp. padifolia differs from subspp. angustifolia and tithymaloides, which are not statistically different for this axis (Fig. 5).

Statistically significant differences in PC 1 (correlated with Wmax, Lmax, A/P ratio, as well as A and P) among three main subspecies suggest that the leaves of subsp. angustifolia leaves are smaller in size (smaller values of Wmax, Lmax, A, and P) and more elongated (smaller A/P ratios) than those of subspp. padifolia and tithymaloides. On the other hand, leaves of subsp. padifolia are characterized by values of PC 1 in the opposite direction (larger values of Lmax, Wmax, A, P, and A/P ratio), suggesting that leaves are larger and rounder. Leaves of subsp. tithymaloides are intermediate for PC 1.

Values along leaf PC 2 (associated with leaf aspect ratio L/W, Lmax, obovate index, and P) reveal that leaves of subsp. angustifolia are more elongated (higher values of A/P ratio) than those of the other two subspecies. For PC 2, subsp. padifolia has intermediate values to those of subspp. angustifolia and tithymaloides.

Morphology PC 3 (how ovate a leaf is – negatively correlated with obovate index) has a lower power to discriminate among subspecies. Results reveal that on average, leaves of subsp. padifolia are statistically more obovate (less ovate) than those of subspp. angustifolia and tithymaloides.

In summary, our analyses provide statistical support for leaves of subsp. angustifolia being smaller and more elongate in shape than those of subspp. padifolia and tithymaloides, and less obovate than those of subsp. padifolia but not statistically different along an ovate-obovate axis from those of subsp. tithymaloides. On the other hand, the leaves of subsp. padifolia are statistically larger, rather round, and markedly obovate, and those of subsp. tithymaloides are quite variable and mostly elliptical in shape.

Results from a linear model approach are consistent with Kruskal-Wallis tests and provide evidence that a significant proportion of the variation in leaf morphology – ranging from 11% (PC 3) to 57% (PC 2) (Supplementary file 6.2) – is captured by taxonomy.

Mantel tests (D_morphologic ~ D_geographic) suggest a correlation between morphology and geography (adj. R² = 0.056, p = 0.0025). However, this relationship disappears when eliminating variation in morphology due to subspecies (p = 0.58; Supplementary file 7). Thus, there is no evidence supporting that similarity in leaf morphology could be explained by geographical proximity alone but rather that variation in leaf morphology is structured by processes that are captured by taxonomy (i.e. subspecies).

Analyses from linear models based on PCs and where variation from subspecies has also been accounted for (using as a response variable the residuals of models focused on subspecies: morphology PCs ~ subspecies) reveal that PC 1 and 2 are negatively related to latitude (p = 0.007 and p = 0.071, respectively) and positively to longitude (p = 0.031 and p = 0.002) (Fig. 6). PC 3 does not show a significant relationship to either latitude or longitude (Fig. 6C, F). These results suggest patterns by which leaves of E. tithymaloides tend to be larger and of somewhat roundish shape in the south and west portions of its range, of smaller sizes and more elongate with increasing latitude (towards the north), and of larger sizes and more rounded shapes with increasing longitude (towards the east) (Fig. 6). These patterns were confirmed by analyses of individual focal variables that showed the largest loadings on the PCs: Lmax, Wmax, and A/P ratio tend to decrease with increasing latitude and increase with decreasing longitude (Fig. 6G–L). Smaller and less round leaves (those that are like those of subsp. angustifolia) are found in larger latitudes and smaller longitudes (NW of the distribution of E. tithymaloides in the Caribbean) and leaves towards the SE of the species range tend to be of larger sizes and rounder shapes, as are those of E. tithymaloides subsp. padifolia.

A first visualization of the climatic space occupied by Caribbean E. tithymaloides suggests no marked differences with respect to portions of climate space occupied by different subspecies, nor that the main geographic areas in the system imply different climatic spaces (Fig. 7). Rather, subspecies of E. tithymaloides tend to occupy subsets of the climatic space already occupied by subspecies tithymaloides (Fig. 7A–C). From a geographic perspective, the Antillean individuals occupy only a portion of a climatic space that is otherwise already occupied by E. tithymaloides individuals from the mainland (Fig. 7D–F).

Because there are significant differences in leaf morphology across subspecies, to examine the proportion of variation in leaf shape that could be attributable to climatic factors alone, we implemented linear models in which we used as response variable the residuals of models where subspecies was used as a factor (X ~ subspecies). For these analyses, only data for the main three subspecies were included (subspp. angustifolia, padifolia, and tithymaloides).

Results reveal that there is an association between foliar morphology and climate. However, only a very small fraction of the variation in leaf morphology (between 1.4 and 2.8%) not accounted for by subspecies is explained by climatic variation (Fig. 8; Supplementary file 8). Larger and rounder leaves (Lmax, Wmax, A, P, A/P, all metrics related to morphology PC 1) are associated with sites that experience higher temperatures (mean annual, coldest month, driest quarter, and coldest quarter, BIO 1, BIO 6, BIO 9, and BIO 11, respectively, all variables related to climate PC 1; Fig. 8; Supplementary file 9) and less annual variation in temperature (BIO 7; p = 0.007; Fig. 8A). Also, PC 1 of leaf morphology is negatively related to climatic PC 3 (p = 0.039; Fig. 8C) suggesting that larger leaves are associated with larger isothermality (BIO 3) and precipitation seasonality (BIO 15) but with lower temperature seasonality (BIO 4).

Leaf shape in relation to aspect ratio (L/W) and degree of obovate shape as well as some leaf Lmax and perimeter (morphology PC 2) are related to higher annual temperatures (climatic PC 1; p = 0.033; Fig. 8D). Morphology PC 2 is also negatively related to climate PC 2 (p = 0.016; Fig. 8E), which is negatively correlated with precipitation (mean annual, of the driest month and quarter, and of the coldest quarter, BIO 12, BIO 14, BIO 17, BIO 19) and positively to maximum temperature of the warmest month (BIO 5).

The metric with the highest loading on morphology PC 3 is the obovate index (Supplementary file 5.1, negative relationship). And morphology PC 3 only is significantly related to climatic PC 2 (Fig. 8H, negative relationship), suggesting that leaves that are more markedly obovate are associated to sites that experience more precipitation (BIO 12, BIO 14, BIO 17, BIO 19) and lower temperatures in the warmest month (BIO 5).

Results of models focusing on the variables with the highest loadings on morphological PCs (Supplementary file 5.1) and one variable per climatic PCs (Supplementary file 9.1) reveal that while there are significant relationships between morphology and climate, the proportion of variation in morphology that is explained by the climatic component is quite low (~1–4.5%). Most significant correlations were those related to temperature seasonality (A/P ratio and Wmax, p = 0.0006 and p = 0.0005, respectively). Four of the seven models that were evaluated were marginally significant: A/P ratio, Lmax, and Wmax in relation to mean annual temperature (p = 0.038, p = 0.052, and p = 0.064; Fig. 9D–F), and Lmax in relation to annual mean precipitation (p = 0.072; Fig. 9G), and three were non-significant (Fig. 9). These models suggest that in E. tithymaloides, once the effect of subspecies has been accounted for, larger and rounder leaves (as per A/P ratio and Wmax) are associated to warmer and less seasonal climates, while leaves with rather elongate and narrow shapes are related to climates that are less warm but more seasonal (Fig. 9). Also, longer leaves (Lmax) are related to sites that receive more mean annual precipitation (Fig. 9). Analyses do not find a significant relationship between the degree to which leaves are obovate and climate.

Climatic distance is not a good predictor of leaf morphologic distance in this system (Supplementary file 10), and this does not change when considering the spatial autocorrelation in climatic data: the partial Mantel test ResidD_morphologic ~ D_climatic + D_geographic was not significant (p = 0.14).