Morphological variation of leaf traits in the Ternstroemia lineata species complex (Ericales: Penthaphylacaceae) in response to geographic and climatic variation

Taxon sampling and morphological trait analysis

We collected leaves of 270 individuals (8–10 mature leaves per individual) from 32 populations located throughout the distribution range of the species (Fig. 1). A field permit of scientific collection (Official number SGPA/DGVS/12770/16) was issued by the Secretaria de Medio Ambiente y Recursos Naturales, of Mexico. Leaf samples were pressed and dried for further measurements in the lab and for herbarium specimens. For each specimen, we measured nine foliar traits: 1. total length including lamina and petiole (TL); 2. lamina length (LL); 3. maximum width of the lamina (MW); 4. petiole length (PL); 5. distance from the base to the maximum width of the lamina (DW); 6. petiole diameter (PD); 7. angle of the lamina apex (ALA); 8. ratio between MW and LL (WLR); and 9. ratio between DW and LL (DWLR) (Fig. 2A). All variables (except ALA) were measured using a Mitutoyo^® Vernier caliper (0.05 mm resolution) and are recorded in mm. For ALA, we used a Jeppesen PJ-1 Rotating Azimuth Plotter.

Variation of morphological traits

We performed a principal component analysis (PCA) with the nine traits using the software Statistica (Statsoft Inc, 2009) to group the variables according to their variability and to select those with the higher values per component to account for most of the morphological variation among individuals (Table 1). With the selected morphological traits (TL, WLR and DWLR), interpolated surface maps were generated in a geographic information system (GIS) (ESRI, 2011) using the Geostatistical Analyst extension (Johnston et al., 2001). The traits DWLR and WLR are descriptors of leaf shape. DWLR refers to the degree to which a leaf is spatulate or elliptical, and WLR is indirect evidence of the foliar area (Nautival et al., 1990; Cittadini & Peri, 2006; Singh, 2007), where high values indicate small foliar areas and low values indicate large areas.

Finally, the interpolated maps were clipped with a map of the selected hydrographic basins of Mexico (INEGI, INE & CONAGUA, 2007) based on the presence of T. lineata complex records and the altitudinal range obtained from the herbarium specimens and field collections to eliminate the interpolation areas outside the known distribution of the T. lineata complex.

We used an empirical Bayesian kriging method, which is based on a semivariogram estimated from the spatial structure of the data (Oliver & Webster, 1990; Kidd & Ritchie, 2006; Brito et al., 2008). This kriging method is a geostatistical tool that generates an estimated surface from a set of dispersed points with Z values. The method is based on statistical models that consider the spatial autocorrelation among data points and provide a measure of precision of the predicted values. We considered the mean error, the square root of the quadratic mean error, the mean standard error, and the square root of the quadratic mean error to evaluate the efficiency of the interpolation (Johnston et al., 2001).

Geographical areas of morphological homogeneity

The three interpolated maps were used to detect discrete areas delimited by zones of rapid change in the values of morphological traits because the variation inside the area is lower than that among the areas. For this method, we used a boundary delineation method with the BoundarySeer software (settings used: crisp boundary type, top 20% in boundary in threshold(s), thresholds 90 deg. vector-to-vector and 30 deg. vector-to-line gradient angle). Different boundaries are zones of rapid change, and for detection, wombling methods are used, which estimate the average amount of change in the variable across the space. The locations that have high values of change are referred to as boundaries (BioMedware, 2013).

Each map generated by the boundary analysis was vectorized, and all were intersected in a single map using ArcGIS Ver. 10.0 (ESRI, 2011). The resulting polygons represented areas of geographical coincidence of the distribution of the three morphological traits: those polygons where the presence of T. lineata complex was not corroborated were eliminated based on the records of biodiversity information systems (REMIB, GBIF, specimens of herbaria). To reduce the number of polygons and group them into homogeneous areas, a cluster analysis was carried out based on the mean values of each of the morphological traits present in each polygon. This analysis was performed using NTSYSpc ver 2.11 (Rohlf, 2000) with a taxonomic distance coefficient and the UPGMA algorithm, and group formation was taken at a value of 0.88. Finally, the resulting groups were evaluated using discriminant analysis to identify their differences (Table 2). Finally, the distributions of the values of each of the three variables within each group were represented using boxplots. Discriminant analysis and boxplot were carried out using SPSS Statistics ver. 19 software (IBM Corp. Released, 2010).

Relationship of morphological traits with climate and geography

To evaluate the relationship of the morphological traits with the climate across the range of the distribution of this species complex, we used the average values of each polygon, considering 19 bioclimatic variables (Hijmans et al., 2005) (Table 1) and three geographical variables. Values were obtained through a zonal analysis in the GIS. We used the redundancy analysis (RDA) implemented in the package ‘vegan’ (Oksanen et al., 2018) in R 3.5.1 software (R Core Team, 2018) to determine the combination of environmental and geographical variables that better explain the morphological variation. This analysis was carried out as follows: first, all the bioclimatic variables (BIO01 to BIO19) were included, and an iterative analysis was carried out so that the variables showing high collinearity were eliminated according to their values of variance inflation factor (VIF). Then, to evaluate the effects of climate and geography on leaf trait variation, we performed RDA analysis in three phases: (a) an analysis that includes both climatic and geographical variables (RDAfull); (b) a partial RDA analysis including only bioclimatic variables without collinearity (pRDA1); and finally, (c) a partial RDA analysis including only geographical variables (pRDA2).

We also conducted linear Pearson correlation analyses between the three morphological traits selected after the PCA and the bioclimatic variables without collinearity after an iterative analysis. We included geographic variables such as latitude, longitude and altitude in these linear correlation analyses to detect patterns of leaf morphological variation through latitudinal, longitudinal or altitudinal gradients.

Leaf morphological variation and interpolation maps

The principal component analysis of the morphological traits showed that the first three components explained 79.71% of the accumulated variance (41.83%, 24.32% and 13.56%). For the first component, the higher loading values were associated with leaf length (TL and LL) and width (DW and MW), and the second and third components were associated with leaf shape (WLR and DWLR, respectively) (Table 2). The relation among the components showed two sets of variables with maximum correlation associated with the leaf length and width (LL, TL and DW) and the lamina angle and leaf shape (ALA and WLR) (Figs. 2B and 2C).

For the three variables identified as important for each component (TL, WLR and DWLR), we obtained the interpolated surface maps (Figs. 3A–3C). We identified a coarse-grain geographical mosaic distribution in the variation of these foliar traits. The longest leaves occurred in southern Mexico (Sierra Madre del Sur), and the plants with the shortest leaves were mainly located in the southernmost region of the distribution (Sierra Norte y Los Altos de Chiapas and Sierra Madre de Chiapas and Cuchumatanes) and a small area in the western part of the Eje Neovolcánico (Fig. 3A). Plants with spatulate leaves (higher values of DWLR) were mainly located in the western part of the Eje Neovolcánico and in the southernmost part of the Sierra Madre de Chiapas and Cuchumatanes. Plants with elliptical leaves (lower values of DWLR) were located mainly in the Sierra Madre del Sur and Sierra del Norte y Altos de Chiapas (Fig. 3B). The regions with plants with the largest foliar area (low values of WLR) were mainly located in the Sierra Madre del Sur, Eje Neovolcánico and Sierra Madre Occidental. The highest WLR values (smaller foliar area) are located in Sierra Madre de Chiapas and Cuchumatanes (Fig. 3C).

Geographical areas of morphological homogeneity and the morphological delimitation of the T. lineata species complex

The boundary analyses and the vectorization of the boundary maps produced 228 polygons for the TL map, 395 polygons for the DWLR, and 270 for the WLR (Figs. 3D–3F). The intersection of these three maps and the subsequent elimination of polygons, using only the records of presence of the T. lineata complex, produced a single map with 108 polygons (Fig. 4A). The cluster analysis with the average values of the three morphological traits in these 108 polygons generated nine clusters (morphological groups) (Fig. 4B). The discriminant analysis showed significant differences among these nine clusters (Table 3) with the eigenvalues of the first two functions discriminated by 92.7% with small Wilks’ lambda values, which confirms the high variability among these groups. Finally, the percentage error among the groups was minimal (0.9%) (Fig. 5A). Figure 5B shows the differences of the three morphological traits among the nine morphological groups.

The nine morphological groups have some correspondence with the species and subspecies of the T. lineata species complex. T. lineata subsp. lineata were formed by the populations of group 1 but also by populations of groups 3, 7 and 9 and some populations of group 4 (in allopatry) and other populations of groups 4 and 5 (in parapatry). Populations of group 2 belong to T. dentisepala, which is distributed in Nayarit state (northwest Mexico), populations of group 6 belong to T. lineata subsp. chalicophila (Chiapas state in southern Mexico) and populations of group 8 belong to T. impressa, which is distributed in Guatemala (Fig. 5C). However, according to the phenogram of the relationships of the nine morphological groups based on morphological similarity, the putative species and subspecies did not show clear patterns of differentiation. For example, the morphological similarity within morphological group 2 (T. dentisepala) was higher and was also higher within the morphological groups of T. lineata subsp. lineata (groups 1, 3 and 4). Group 5 was an interesting case because the populations of this group were distributed parapatrically and differentiated from the main cluster of T. lineata subsp. lineata. Populations of the morphological groups 6, 7, 8 and 9 were geographically isolated and differentiated in different clusters, and plants of group 8 that were distributed in Guatemala were described as T. impressa (Fig. 5C).

Relationship of geography and climate to morphological traits

We selected six bioclimatic variables (BIO01, BIO04, BIO13, BIO14, BIO15 and BIO19) based on an iterative analysis. The results of the full redundancy analysis (F = 10.037, p < 0.001, N = 999 permutations), including geographic and bioclimatic variables, explained 47.79% of leaf variation (p < 0.001) with two significant axes RD1 and RD2; the first axis explained 83.72% of the variance, and the second axis explained 13.61%. The most significant variables for the first axis were precipitation of the driest month (BIO14, p < 0.001), precipitation seasonality (BIO15, p < 0.043), precipitation of the coldest quarter (BIO19, p < 0.007) and longitude (p < 0.001); for the second axis, the most significant variables were latitude (p < 0.001), precipitation of the wettest month (BIO13, p < 0.024), temperature seasonality (BIO04, p < 0.014) and annual mean temperature (BIO01, p < 0.001). The most important morphological variables were total leaf length (TL) and width-length ratio (WLR). The first axis of the full RDA separated those populations mainly located in the south in Sierra del Norte y Los Altos de Chiapas and Sierra Madre de Chiapas and the populations of the Sierra Madre del Sur (in the state of Oaxaca) and the most western area of the Eje Neovólcanico (state of Nayarit) (Fig. 6).

The analysis of leaf trait variation indicated that climate had a higher impact than geography. The percentage of the variance explained only by climate was 31.08%, the percentage explained by geography was 13.38% and the join effect, geography—climate, explained 2.38%.

The partial redundancy analysis using climate as control of geography (pRDA1) was significant only in its first axis (p < 0.001), explaining 87.26% of the variation, with precipitation of the coldest quarter (BIO19, p < 0.009) as the most important variable. The partial redundancy analysis using geography as a control climate (pRDA2) was only significant in its first axis (p < 0.001), which explains 99.45% of the variation, with latitude being the variable of greater importance (p < 0.001) (Fig. 6).

Relationship between morphological traits and environmental and geographic variables

After the principal component analysis defined TL, WRL and DWLR as the morphological traits that explain most of the variance in the T. lineata species complex and the iterative analysis selected the climatic variables without collinearity, we conducted linear correlations between both groups of variables. Total leaf length (TL) was positively correlated with annual mean temperature (BIO01) (r = 0.30; p < 0.01) and precipitation seasonality (BIO15) (r = 0.48; p < 0.01) and negatively correlated with the precipitation of the driest month (BIO14) (r =  − 0.53; p < 0.01), precipitation of the wettest month (BIO13) (r =  − 0.13; p < 0.01) and precipitation of the coldest quarter (BIO19) (r = 0.45; p < 0.01) (Fig. 7).

The ratio between the maximum width of the lamina (MW) and leaf length (LL) (WRL) as a descriptor of leaf shape was positively correlated with the precipitation of the wettest month (BIO13) (r = 0.39; p < 0.01) and precipitation of the driest month (BIO14) (r = 0.31; p < 0.01) and was negatively correlated with temperature seasonality (BIO04) (r =  − 0.27; p < 0.01) and precipitation seasonality (BIO15) (r =  − 0.37; p < 0.01) (Fig. 8).

The ratio between the distance from the base to the maximum width of the lamina (DW) and leaf length (LL) (DWLR) was positively correlated with temperature seasonality (BIO04) (r = 0.40; p < 0.01) and precipitation seasonality (BIO15) (r = 0.25; p < 0.05) and negatively correlated with the precipitation of the wettest month (BIO13) (r =  − 0.21; p < 0.05), precipitation of the driest month (BIO14) (r = 0.32; p < 0.01) and precipitation of the driest quarter (r =  − 0.41; p < 0.01) (Fig. 9).

Additionally, TL and DWLR were positively correlated with latitude (r = 0.20; p < 0.05; r = 0.22; p < 0.05, respectively) and negatively correlated with longitude (r =  − 0.30; p < 0.05; r =  − 0.24; p < 0.05, respectively), but WLR was negatively correlated with latitude (r =  − 0.53; P < 0.01) and positively correlated with longitude (r = 0.49; p < 0.01). The three morphological traits were not correlated with altitude (TL: r =  − 0.12; p > 0.05; WLR: r =  − 0.08; p > 0.05; DWLR: r =  − 0.14; p > 0.05) (Fig. 10).