Introduction

Dioon edule Lindl. (Zamiaceae) is the species with the widest geographical range in Mexico from the genus Dioon. It is distributed along the Sierra Madre Oriental (SMO), from southern Tamaulipas throughout San Luis Potosí, Hidalgo, Querétaro, and Veracruz, with 1,525 m asl (^{de Luca et al., 1982}; ^{Whitelock, 2004}; ^{Yáñez-Espinosa, 2006}). In the northernmost region of the Sierra Madre Oriental, in Tamaulipas and Nuevo León, there are small populations of D. edule var. angustifolium, disjunct and with restricted distribution areas (^{de Luca et al., 1982}; ^{González et al., 2008}; ^{Whitelock, 2004}), recognized recently with the species status as Dioon angustifolium Miq. (^{González-Astorga et al., 2003}). The populations of this species are confined to rocky canyons in transition zones between piedmont scrub and tropical dry forest or oak forest. The plants mainly develop in shallow soil (Lithosol), with a depth not exceeding 15 cm, clayey texture, pH of 7.54 and approximately 24% of organic matter (^{Chávez-Acuña, 2010}; ^{Rubio-Méndez, 2010}). This species has played a relevant role in the indigenous cultures of the Sierra Madre Oriental. It has been used as a ceremonial, ornamental, medicinal, and food when maize is scarce (^{Bonta et al., 2019}; ^{Chávez-Acuña, 2010}; ^{Tristán-Martínez, 2012}; ^{Yáñez-Espinosa, 2006}).

Due to the diversity of habitats in which Dioon edule grows and based on the observed morphological variation in populations (^{de Luca et al., 1982}; ^{Whitelock, 2004}) along a latitudinal gradient in the Sierra Madre Oriental, from the northernmost distribution area in Nuevo León and Tamaulipas to central Veracruz (Table 1), some studies have recognized that there is intraspecific variation in some morphological traits, proposing the classification of D. edule populations in 7 "forms". However, these studies warned readers about not having included more populations distributed in the Sierra Madre Oriental (^{Whitelock, 2004}). Particularly in the state of San Luis Potosí, ^{Whitelock (2004)} described the forms "rioverde" and "valles", which have been recognized in this study during the field trips carried out since 2008 to identify and characterize Dioon edule populations in the Sierra Madre Oriental of San Luis Potosí. However, there are more populations which show suggestive characteristics of other morphological forms not described by ^{Whitelock (2004)}. Some studies also highlight the importance of comparing the phenotypic variation among Dioon species and their correlation to habitats to understand the evolutionary response to the environment, particularly to aridity (^{Gutiérrez-Ortega, Yamamoto et al., 2018}).

Here, we analyze the main morphological traits of plants in the populations along the Sierra Madre Oriental of Dioon edule in San Luis Potosí and Querétaro and of D. angustifolium in Nuevo León and Tamaulipas, to recognize and classify populations into the 2 forms proposed by ^{Whitelock (2004)} in the state of San Luis Potosí, as well as possible new ones. Also, we related them to environmental and geographic variables that could explain this phenotypic variation.

Materials and methods

This study was conducted in the Sierra Madre Oriental in Mexico, considering 3 populations of D. angustifolium in Nuevo León and Tamaulipas, and 12 populations of Dioon edule located in the states of San Luis Potosí and Querétaro. A latitudinal range from 24°51' to 21°30' N, and an elevation range from 294 to 1,314 m asl were covered (Fig. 1). The study was conducted in the summer of 2007 to 2018. The populations were classified in the forms proposed by ^{Whitelock (2004)} according to their elevation range (Table 1). In each population, 2 sites of 100 m² each were randomly selected and mature plants, defined as those in the reproductive stage of the life cycle with a visible stem and leaves arranged in 1 or more well-developed crowns, were evaluated. The parameters that were measured in each plant were plant total height, stem diameter, leaves number and leaf length.

Figure 1 Study area in the Sierra Madre Oriental, México. 

A discriminant analysis (DA) was carried out to classify the 15 populations (n = 318) in the 3 forms defined by ^{Whitelock (2004)} for the Sierra Madre Oriental of San Luis Potosí, delimit the morphological traits of each form, and identify those that were not included in the range. A Kruskal-Wallis test was applied to assess for significant differences on elevation groups and populations for parameters because the measurement variables did not have a normal distribution not even with data transformation.

Populations were represented in Arcgis 10.3, to conduct a spatial analysis to identify the environmental features of the sites at which the populations grow. The environmental variables considered in the analysis were: annual mean temperature, and annual precipitation (1970-2010) obtained from the Worldclim database (^{Hijmans, 2012}). The geographical layers were geology, edaphology, altitude, slope and the use of the soil and vegetation of the series V, generated from the digital elevation model of ^{INEGI (2013)}. All layers were represented in raster format with a spatial resolution of 30 m. The maximum entropy of the Maxent 3.4.1 program algorithm was used to estimate the contributions of each variable, as well as their respective parameters. A Pearson correlation analysis among plant traits (leaf length and plant height), geography and environmental variables was performed to assess their association. All statistical analyses were performed with the software Xlstat v.2014.1.0.1 (Addinsoft, France).

Results

The Dioon edule forms according to ^{Whitelock (2004)}. The DA showed that the vectors of the 4 forms were significantly different (p < 0.0001). The populations were classified into 4 groups according to their morphological traits. Three of the forms corresponded to those proposed by ^{Whitelock (2004)}: "edule" (3 populations), "rioverde" (4) and "angustifolium" (3). Five populations that were not classified because they did not fit the ranges established by Whitelock are now proposed to be a new form and designated "tamasopo" (Table 2, Fig. 2). It is remarkable that almost 40% of individuals of the "angustifolium" populations share morphological traits with the "tamasopo" form, but only 25% with the form "rioverde" and less than 4% with the "edule" form.

Figure 2 Discriminant Analysis biplot of the Dioon edule form "edule", "rioverde", "tamasopo" and D. angustifolium classification. 

Two discriminant functions explained 99% of the total variability, contributing significantly to the separation between forms for the 15 populations (p < 0.0001, n = 313). The first function (eigenvalue = 0.311) explained 79.55% of the total variation, the second (eigenvalue = 0.076) explained 19.55%, and the third (eigenvalue = 0.004) explained the remaining variation. Morphological traits that contributed significantly to the separation of the centroids were leaf length, stem diameter and leaves number (Table 3).

The species Dioon edule has the widest stem diameter and longer leaves than the other forms, while D. angustifolium presented the lowest values concerning stem diameter and leaf length. However, both species have the largest numbers of leaves, a trait which separates them from the other D. edule forms (Table 4, Fig. 3). Significant statistical differences were found in all parameters measured by population pairs. As for the base diameter, the edule form is different from the "angustifolium" (p < 0.005), "rioverde" (p < 0.0001), and "tamasopo" (p < 0.0001) forms. The "tamasopo" form is, in turn, different from "angustifolium" (p < 0.005) and "rioverde" (p < 0.005). In height, "edule" was different from "angustifolium" (p < 0.0001) and "rioverde" (p < 0.001) and "tamasopo" was different from "angustifolium" (p < 0.05) and "rioverde" (p < 0.05). Regarding the leaves, there was only a difference between the "edule" and "tamasopo" forms (p < 0.005). In the leaf number per flush, "angustifolium" is different from "rioverde" (p < 0.0001) and "tamasopo" (p < 0.0001), while "edule" is different from "rioverde" (p < 0.05) and "tamasopo" (p < 0.01).

Figure 3 Morphology of Dioon edule forms a) "edule"; b) "rioverde"; c) "tamasopo"; d) D. angustifolium (photo by María Magdalena Salinas Rodríguez), in their distribution range. 

The "tamasopo" form occurs in several locations at an altitude of 850-1114 m asl (Fig. 3c). Stems are about 0.1-1 m in height and 0.02-0.3 m in diameter. The leaves are flat, about 1.34 m in length. Leaflets are straight, revolute, spiny, of 7.29 cm (± 0.11) in length and 0.56 cm (± 0.01) in width, inserted in acute angle in the yellowish rachis, with 0.25 cm (± 0.005) between them. Female cones are ovoid, densely tomentose, light brown, 205-260 mm long and 122-175 cm in diameter, male cones are cylindrical with conical apex, whitish green. Seeds are ovoid, the sarcotesta is light yellow to orange.

In general, the largest diameter was found in population E2, as well as the largest number of leaves; on average, these were also the longest (Table 3). The tallest plants were found in T4 and the highest leaf production was found in E3 (Table 5). Significant differences were found in all parameters measured between population pairs (Table 6). For diameter, population T2 is different from RV4 (p < 0.05), RV2 (p < 0.0001), T5 (p < 0.001), A3 (p < 0.05) and T1 (p < 0.0001). In turn, T1 is different from T5 (p < 0.05) and E3 (p < 0.0001). This is different from RV2 (p < 0.05) and T5 (p < 0.05), which is also different from RV2 (p < 0.05) and A3 (p < 0.05).

For plant height, T5 was different from T1 (p < 0.05), RV2 (p < 0.05), and A3 (p < 0.05). T4 was different from T1 (p < 0.01) and A1 (p < 0.05), and T2 was different from A3 (p < 0.05). Regarding the number of leaves, T2 was different from A2 (p < 0.05), RV3 (p < 0.01) and E3 (p < 0.05). For leaves number per flush, T5 is different to A1 (p < 0.001) and E3 (p < 0.01). For leaf length, T5 is different from T2 (p < 0.001), RV4 (p < 0.001), A2 (p < 0.01), A3 (p < 0.05), A1 (p < 0.05) and T4 (p < 0.05). Population E3 is different from T2 (p < 0.05), RV4 (p < 0.001), A2 (p < 0.001) and A1 (p < 0.005), and population T4 is different from RV4 (p < 0.005) and A1 (p < 0.05).

Overall, 40% of the georeferenced sites are located at altitudes ranging from 497 to 925 m. In addition, 66% of field records indicate that the species develop on slopes lower than 14°, where soils tend to be deeper due to the accumulation of materials. Precipitation had the greatest contribution, followed by slope, vegetation, elevation and geology (Table 7). The remaining variables did not contribute to the model. The most relevant variables that influence the phenotypic differences among the forms are precipitation, physiography and vegetation. The correlation analysis showed a moderate negative linear relationship between plant height and latitude (r = -0.538, p < 0.05) and a moderate positive relationship between leaf length and annual precipitation (r = 0.538, p < 0.05). Likewise, there are moderate negative linear relationships between annual precipitation and elevation (r = -0.577, p < 0.05) and latitude (r = -0.520, p < 0.05), but a moderate positive relationship with longitude (r = 0.636, p < 0.05).

Discussion

The results of our study highlight the phenotypic variation of Dioon edule populations in the latitudinal gradient mentioned by ^{Whitelock (2004)}, as well as the altitudinal gradient. The inclusion of more populations in the sample allowed us to confirm the "rioverde" and "edule" forms, as well as to propose the new "tamasopo" form. The fact that D. angustifolium presents a large proportion of individuals that share morphological traits with D. edule in the "tamasopo", "rioverde" and "edule" forms, supports recent studies showing that D. angustifolium do not form a monophyletic group and individuals are nested within D. edule populations (^{Dorsey et al., 2018}; ^{Gutiérrez-Ortega, Salinas-Rodriguez et al., 2018}). However, the similarities of the morphological traits of some populations of the "rioverde" and "tamasopo" forms with D. angustifolium allow us to suggest the presence of the latter in San Luis Potosí, but this will only be resolved when corresponding genetic studies are carried out.

The morphological traits range of the sampled individuals are considerably smaller than those described by ^{Whitelock (2004)}, which may be due to a greater number of individuals evaluated in several populations, increasing the data variation. Likewise, D. edule stems bend or break when they fall to the ground due to the weight of the reproductive structures and are covered by the ground. The leaf length and plant height varied more across latitude than across elevation. Northern populations are subject to lower annual precipitation than southern populations (1:2), but higher interannual variation with 6 to 8 months of dry season, leading to greater phenotypic plasticity. Other studies about intraspecific phenotypic variation along gradients have shown a similar pattern of phenotypic plasticity associated to variation in precipitation across latitude than elevation (^{Jonas & Geber, 1999}).

The results confirmed that precipitation is the variable with the greatest influence on the phenotypic variation in Dioon edule forms in San Luis Potosí, as mentioned in a previous study for Dioon species (^{Gutiérrez-Ortega, Salinas-Rodriguez et al., 2018}). The annual precipitation is lower at higher latitude and elevation following an aridification pattern, suggesting that the aridification promoted the diversification of the Dioon edule complex. However, a previous study (^{Dorsey et al., 2018}) suggests that D. edule complex diversification is more consistent with fragmentation of a widespread ancestor or longdistance dispersal. The intraspecific phenotypic variability is associated with adaptations to environmental conditions derived from the isolation of populations, especially after the appearance of large mountain ranges, as well as to the reduction of genetic flow between populations (^{de Luca et al., 1982}; ^{Whitelock, 2004}). As found in other studies, the intraspecific traits variation among and within populations, is due to genetic variation and phenotypic plasticity (^{Souza et al., 2018}). Thus, the emergence of new abiotic conditions such as the effects of climate change are not only threats to the species, but also a trigger of biodiversity, causing natural selection pressures that can promote adaptations (^{Gutiérrez-Ortega, Salinas-Rodriguez et al., 2018}). Our study agrees that the variation of intraspecific traits is a relevant mechanism by which individuals can cope with environmental variations and could avoid local extinctions under climatic change (^{Souza et al., 2018}).

The populations of Dioon edule studied were classified into the 3 forms proposed by ^{Whitelock (2004)}: D. angustifolium form "angustifolium", and D. edule forms "edule" and "rioverde", and the new form "tamasopo". The forms "tamasopo" and "rioverde" present intermediate values between D. edule and D. angustifolium, with 63% of the individuals classified as "tamasopo". Our results suggest the presence of D. angustifolium in San Luis Potosí, but genetic studies must be carried out to corroborate it and find out if hybridization has taken place. The variables that define the distribution of the populations are vegetation type, slope inclination and annual precipitation.