Materials and methods

Morphological data. We measured 207 individuals of Lycianthes moziniana from 21 populations. These included 102 individuals collected between 2017-2021. A voucher specimen for each population was deposited in the Luz María Villarreal de Puga Herbarium at the University of Guadalajara (IBUG, ^{Thiers 2023}). Additionally, we included 105 specimen images from ANSM, CIIDIR, DAV, GBH, HUAA, HUAP, MEXU, NY, QMEX, TEX, US and XAL (Thiers 2023). The taxa sampling included L. moziniana var. margaretiana (four populations and 35 individuals), L. moziniana var. moziniana (13 populations and 150 individuals) and L. moziniana var. oaxacana (four populations and 22 individuals).

Morphological analyses. Herbarium specimens were photographed with a Canon EOS Rebel T5 camera. Both, photos and herbaria specimen images were measured using ImageJ (imagej.net). Using various ImageJ functions, we obtained three leaf values (petiole length, leaf length and leaf width) from the second sympodium unit (^{Dean 2004}). Also, we included pedicel length, calyx diameter and length, calyx appendages length, corolla diameter, long stamen length and short stamens length, style length, calyx in fruit length and diameter, and fruit length and diameter. In addition, the ratio values of leaf length/width, calyx length/diameter, calyx in fruit length/diameter, and mature fruit length/diameter were calculated.

We used 11 characters to test the differences among varieties. Then, the characters were transformed with natural logarithms to fit normality and parametric assumptions. With the transformed values and the three varieties as the classificatory variable, we conducted a Discriminant Canonical Analysis (DCA) with the candics package in RStudio version 1.3.1093 (^{R Studio Team 2020}). Mahalanobis Square Distance (MD) between the centroids of the Canonical Classificatory Analysis (CCA) were calculated in SAS (^{SAS 2013}).

Climatic range, distribution models and niche differentiation. The geographic distribution was obtained from herbarium specimens, electronic databases, and literature. We examined data from the following 34 herbaria (either in person or using digital images): ANSM, CFNL, CICY, CIIDIR, CIMI, CIQRO, DAV, F, GBH, GUADA, HCIB, HGOM, HUAA, HUAP, HUAZ [not in Index Herbariorum (^{Thiers 2023}), Herbario de la Universidad Autónoma de Zacatecas], IBUG, IEB, INEGI, MEXU, MO, NY, OAX, P, QMEX, SERBO, TEX, UC, UCR, UNL, UNSIJ [not in Index Herbariorum (^{Thiers 2023}), Universidad de la Sierra Juárez, Oaxaca], US, WIS, XAL and ZEA (acronyms according to ^{Thiers 2023}). In addition, three electronic databases were consulted: Global Biodiversity Information Facility (GBIF, www.gbif.org), Missouri Botanical Garden (Tropicos www.tropicos.org) and Southwest Environmental Information Network (SEINet, www.swbiodiversity.org). Information contained in ^{Dean (2004)}, ^{Dean et al. (2020)}, ^{Martinez et al. (2020)} and ^{Nee (1986)} was added. We built a database of 420 records with 32 records of Lycianthes moziniana var. margaretiana, 349 of L. moziniana var. moziniana and 39 of L. moziniana var. oaxacana. The spatial autocorrelation of distribution records was eliminated manually in QGIS Las Palmas v. 2.18.3 (^{QGIS Development Team 2017}). In other words, from the projected records points and the climatic layers, a single point for each pixel was selected. The final database contained 348 records of L. moziniana var. margaretiana (25), L. moziniana var. moziniana (288) and L. moziniana var. oaxacana (35).

We obtained 19 climatic variables from WORLCLIM Database based on the CCSM4 model (^{Hijmans et al. 2005}, ^{Otto-Bliesner et al. 2006}) with a resolution of 30 arc seconds (www.wordclim.org). Using the presence points, we extracted the values for the climatic variables. Then, to select the variables, the correlation of the variables was tested for each variety with Person’s correlation using a 0.05 level of significance in the corrplot R package in RStudio. Later, we calculated the variance inflation factor for correlated variables and excluded the highly correlated variables using the threshold of 0.9 in the usdm R package in R Studio. Finally, the layers were trimmed out to the accessible area (M) (^{Soberón & Peterson 2005}, ^{Barve et al. 2011}). The M represents the MTZ of ^{Morrone et al. (2017)} plus a buffer of 50 km.

We executed an Ecological Niche Model (ENM) for each variety. The analysis included 25 records of Lycianthes moziniana var. margaretiana, 288 records of L. moziniana var. moziniana, and 35 records of L. moziniana var. oaxacana in MAXENT v. 3.3.3K (^{Phillips et al. 2006}). Each model was replicated 10 times. Samples collected between 2017-2021 were used to validate the models. These included four records for L. moziniana var. margaretiana, four presence points of L. moziniana var. oaxacana and 15 records for L. moziniana var. moziniana. The rest of the parameters had the default configuration. Finally, the models were evaluated using under the curve values.

The niche overlap among varieties was estimated with the values D (niche equivalency) and I (niche similarity) in ENMTools (^{Warren et al. 2008}, ²⁰¹⁰). These indices consisted of values between 0 and 1, indicating no overlap and full overlap among ENMs, respectively. This test was carried on among all binary combinations among the three taxa (Lycianthes moziniana var. moziniana and L. moziniana var. margaretiana, L. moziniana var. moziniana and L. moziniana var. oaxacana and L. moziniana var. margaretiana and L. moziniana var. oaxacana). Identity test was carried on with the average layer of ENMs constructed in MAXENT using four variables (BIO4, BIO7, BIO15 and BIO18) shared among the three varieties and the occurrences files mentioned above and 100 replicates. The test calculated the overlap between ENM of one taxon and the ENMs constructed using random occurrence point samples of the other taxon. To evaluate the differentiation among niches, we used the measured niche overlap as empirical estimation and compared with the distribution of overlaps from pseudo replicates.

Results

Morphological analyses. Discriminant Canonical Analysis recovered three groups. Two canonical functions explained 100 % of the morphological variation, 78.6 and 21.4 %, respectively (Figure 2H). Canonical structure coefficients showed that seven morphological characters (petiole length, leaf ratio, long stamen length, short stamens length, style length, calyx in fruit and fruit ratios) had the highest discriminatory power in the DCA (Table 1). On the first function (Can1), petiole length, leaf ratio and style length explained the variation. On the second (Can2), the short stamens length, the long stamen length, calyx in fruit and fruit ratios explained the variation (Table 1). The petiole length was important for Lycianthes moziniana var. margaretiana. Four characters (leaf ratio, style length, and the short and long stamens length) grouped the individuals of L. moziniana var. moziniana and two characters (calyx in fruit and fruit ratio) grouped the individuals of L. moziniana var. oaxacana (Figure 2).

Figure 2 Characters with the highest discriminatory power from Canonical Discriminant Analysis. Boxplot A. leaf ratio; B. petiole length; C. long stamen length. D. short stamens length; E. style length; F. calyx in fruit ratio; G. fruit ratio and H. Discriminant Canonical Analysis. 

The MD among all a priori groups were statistically significant (P < 0.05, Table 2). Between pairwise varieties comparisons, the centroids of Lycianthes moziniana var. margaretiana and L. moziniana var. oaxacana clusters were closer together. Followed by L. moziniana var. margaretiana and L. moziniana var. moziniana. In contrast, the furthest clusters were L. moziniana var. moziniana and L. moziniana var. oaxacana. The confusion matrix of CCA, including the seven characters with highest discriminatory scores, had 82 % recovery of groups corresponding to L. moziniana var. margaretiana, 80 % recovery of L. moziniana var. moziniana and 77 % of L. moziniana var. oaxacana (Table 3).

Climatic range, distribution models and niche differentiation. For each variety, the climatic variable selection was different. The analysis recovered five variables (BIO4, BIO7, BIO9, BIO15 and BIO18) for Lycianthes moziniana var. margaretiana, eight variables (BIO4, BIO7, BIO12, BIO14, BIO15, BIO17, BIO18 and BIO19) for L. moziniana var. moziniana and six variables (BIO4, BIO7, BIO9, BIO13, BIO15 and BIO18) for L. moziniana var. oaxacana. The ENMs of each variety showed high predictive power with AUC values 0.953, 0.912 and 0.984, respectively (Figure 3C). The jackknife test selected BIO15 and BIO4 as the most important variables of L. moziniana var. margaretiana. In contrast, L. moziniana var. moziniana was associated with BIO4, BIO14 and BIO7, and the presence of L. moziniana var. oaxacana was determined by BIO4 and BIO9.

Figure 3 Distribution of Lycianthes moziniana; A. Actual distribution of the varieties of L. moziniana; B. Potential distribution of L. moziniana; C. Potential distribution of each variety of L. moziniana. 

The niche identity test showed significant differences among ENMs. The empirical niche overlap index in all comparisons was significantly lower than the expected by chance. The empirical niche overlap between Lycianthes moziniana var. margaretiana and L. moziniana var. moziniana was I = 0.6098 and D = 0.3227 (Figure 4A). The same values between L. moziniana var. margaretiana and L. moziniana var. oaxacana were I = 0.2868 and D = 0.1186 (Figure 4B). Lastly, similarity and equivalency values between L. moziniana var. moziniana and L. moziniana var. oaxacana sumed I = 0.5279 and D = 0.2749 (Figure 4C).

Figure 4 Niche conservatism inference in Lycianthes moziniana; A. Niche identity test between L. moziniana var. moziniana and L. moziniana var. margaretiana; B. Niche identity test between L. moziniana var. moziniana and L. moziniana var. oaxacana; C. Niche identity test between L. moziniana var. margaretiana and L. moziniana var. oaxacana. The dotted line shows the empirical value of I index in blue and D index in red. The bars represent the null distribution of I index in blue and D index in red. 

Discussion

The morphological and geographical analyses recovered three groups within Lycianthes moziniana (Figures 2, 3). On the other hand, ENM comparisons within each variety showed low niche overlap (Figure 4). This evidence revealed the geographic and climatic isolation of the varieties. These results were supported also by the morphological and ITS sequences differences observed by ^{Dean (2004)} and the genetical structure found in populations distributed in different biogeographic provinces (^{Anguiano-Constante et al. 2021a}). The combined evidence suggests that each variety has been isolated long enough to generate genetic structure and morphological differences. We argue that all three are separately evolving metapopulations and in the process of speciation (^{de Queiroz 2007}).

Morphological analyses (DCA, MD and CCA) recovered three statistically supported groups that correspond to the three varieties (Figure 2 and Tables 2, 3). These analyses have been useful for delimiting species boundaries in other groups (^{Aquino et al. 2019}, ^{Gándara et al. 2021}, ^{Gutiérrez-Ortega et al. 2021}, ^{Sánchez et al. 2020}). Interestingly, the DCA explained the variation of Lycianthes moziniana in two canonical functions (Figure 2). Seven out of 11 examined diagnostic characters (petiole length, leaf ratio, short stamens length, long stamen length, style length, calyx in fruit and fruit ratios) were significantly different. Five of the seven were reproductive characters and could be closely related to pollinators and dispersal agents. The genetic variation along with geographic isolation has promoted morphological fixation and differentiation of populations.

Foliar characters changed among varieties. Petiole length in Lycianthes moziniana var. margaretiana was longer [0.1-10 (± 4.7 mm)] than L. moziniana var. moziniana and L. moziniana var. oaxacana. In the last two taxa, the leaves were sessile or the petiole length was less than 7 mm (± 2.8 and ± 3.4 mm). The leaves of L. moziniana var. moziniana tended to be lanceolate in shape and narrower [2.2-9.2 (± 4.7) × 0.7-4 (± 1.9 cm)], in comparison to the wider ovate to elliptic leaves of L. moziniana var. margaretiana [2.5-10 (± 6.2) × 1.7-6 (± 3.7 cm)] and L. moziniana var. oaxacana [3-7 (± 4.8) × 1-4.6 (± 2.3 cm)]. The leaf variation is indicative of adaptations to different climatic conditions. Temperature and precipitation have been recovered in other analyses as the main climatic factors modulating morphological variation (^{Rodríguez-Gómez et al. 2018}, ^{Maya-García et al. 2020}, ^{Martínez-García et al. 2022}).

In the same way, flower characters differed among varieties. Lycianthes moziniana var. moziniana showed more variation in stamen length [long 6.5-12.8 (± 9.4 mm), short 5-10 (± 7.3 mm)] and style length [7.3-15 (± 10.6 mm)] than L. moziniana var. margaretiana [long 5.8-11.8 (± 7.8 mm), short 5.2-8.7 (± 6.2 mm) and style length 6.9-12.4 (± 8.7 mm)]. Lycianthes moziniana var. oaxacana [long 5.9-9.4 (± 7.9 mm), short 5-6.8 (± 6 mm), and style 7.1-10.6 (± 9.1 mm)] had the shortest stamens and style. In the case of the length of the calyx in fruit, the results did not overlap between L. moziniana var. margaretiana [7.5-10.9 (± 8.9) mm long × 18.3-19.8 (± 19 mm diameter)] and L. moziniana var. oaxacana [9.2-9.7 (± 9.4) × 16.9-18.4 (± 17.7 mm)]. Lycianthes moziniana var. moziniana [3.3-10.7 (± 6.5) × 7.8-33.8 (± 13.7 mm)] had the shortest length of the calyx in fruit (Figure 2A-F).

Finally, the fruit shape in Lycianthes moziniana var. moziniana was narrow [0.7-2.4 (± 1.6 cm)] and long [1.1-3.2 (± 2.1 cm)]. In contrast, L. moziniana var. margaretiana [1.9-2.2 (± 2) × 2.3-4.4 (± 2.7 cm)] and L. moziniana var. oaxacana [2-2.2 (± 2) × 3-3.4 (± 3.2 cm)] had ovoid fruits (Figure 2G). ^{Dean (2004)} found that the combination of fruit shape and appendages position in the fruiting calyx were important for delimiting the varieties. In fruit, the appendages of the calyx were appressed in L. moziniana var. moziniana and spreading in the other two varieties. The ovoid fruit in L. moziniana var. margaretiana and L. moziniana var. oaxacana was green but the former presented purple blotches. Moreover, the fruit in L. moziniana var. moziniana was green, round or ovoid. The morphological variation observed in the three varieties has been affected by gene flow, geographic distance, pollination ecology and climatic conditions (^{Rodríguez-Peña & Wolfe 2023}). Thus, these events have determined the independent evolutionary trajectories among lineages (^{Jacobo-Arteaga et al. 2022}).

The MD demonstrated significant differences among varieties (Table 2). Lycianthes moziniana var. moziniana was the most distant, relative to L. moziniana var. oaxacana, and L. moziniana var. margaretiana. This was congruent with the taxonomic treatments of ^{Dean (2004)} and ^{Dean et al. (2020)}. However, L. moziniana var. moziniana and L. moziniana var. oaxacana hybridized under greenhouse conditions (^{Dean 2004}) and L. moziniana var. moziniana and L. moziniana var. margaretiana were sympatric in Sierra de Álvarez, San Luis Potosí. The same analysis recovered L. moziniana var. margaretiana and L. moziniana var. oaxacana close to each other, but they were isolated by the TVB. None of the varieties show genetic flow among them (^{Anguiano-Constante et al. 2021a}).

The varieties grew in different biogeographic provinces, but in similar vegetation types. Lycianthes moziniana var. margaretiana is restricted to the SMOr, L. moziniana var. oaxacana to SMS and L. moziniana var. moziniana inhabit TVB and SMOc (^{Anguiano-Constante et al. 2021a}, Figure 3A). In addition, we obtained a different predicted set of climatic variables for each variety. Nevertheless, the temperature seasonality (BIO4) was important for the three varieties. This was congruent with the rest of the species included in the Lycianthes series Meizonodontae (^{Anguiano-Constante et al. 2018}). Other important variables were precipitation seasonality (BIO15) for L. moziniana var. margaretiana, temperature annual range (BIO7) for L. moziniana var. moziniana and mean temperature of driest quarter (BIO9) for L. moziniana var. oaxacana. Precipitation and temperature seasonality are important climatic variables for the evolution of geophytes (Sosa & ^{Loera 2017}, ^{Howard et al. 2019}). Thus, climatic variables in these geophytic taxa, could be an important cause of divergence.

The three varieties shared four climatic variables (BIO4, BIO7, BIO15 and BIO18). Except for the precipitation of the warmest quarter (BIO18), all of them were important predictors for the ENMs. However, the niche overlap analysis among varieties showed significant differences. Lycianthes moziniana var. moziniana overlapped the most with L. moziniana var. margaretiana (I = 0.6098 and D = 0.3227), followed by L. moziniana var. oaxacana (I = 0.5279 and D = 0.2749, Figure 4A, B). This result makes sense, because L. moziniana var. moziniana has a wide range relative to the other two varieties. A similar case was observed in Pinus montezumae Lamb. and P. pseudostrobus Lindl. (Pinaceae) which were distributed in the SMOr, SMS and TVB, and had intermediate niche overlap and the ability to hybridize (^{Manzanilla-Quiñones et al. 2019}).

Lycianthes moziniana var. margaretiana and L. moziniana var. oaxacana had the lowest niche overlap (I = 0.2868 and D = 0.1186, Figure 4C). The morphological evidence suggested significant differences (Figure 2 and Table 2). The isolation between the taxa could be responsible for this divergence. Several Mexican plants show a similar pattern when their populations are distributed in the SMOr and SMS (^{Suárez-Mota et al. 2015}, ^{Alvarado-Sizzo et al. 2018}, ^{Rosas-Reinhold et al. 2022}). This pattern suggests that both taxa diverged from an ancestral population and isolation generated and maintained the genetic and morphological differences.

Most likely, differences observed in Lycianthes moziniana resulted from geological, topographic, and climatic variation. These conditions isolated the populations, and each evolved separately, forming, and fixing the variation. The Mexican biogeographic provinces have different origins and ages (^{Mastretta-Yanes et al. 2015}). Also, the climatic conditions have changed drastically in the last 100 Myr (^{Cevallos-Ferriz et al. 2012}). All together contributed to the current biodiversity of Mexico and were directly related to genetic and morphological variation of the taxa. Important patterns have been documented. First, the geographic barriers along the MTZ promoted the speciation in several Mexican plant groups (^{Alvarado-Sizzo et al. 2018}, ^{Gutiérrez-Ortega et al. 2020a}, ^{Romero-Soler et al. 2022}). Second, the climate refugia created genetic variation in Ephedra compacta Rose (^{Loera et al. 2017}), Psittacanthus sonorae (S.Watson) Kuijt (^{Ornelas et al. 2018}), Podocarpus spp. (^{Ornelas et al. 2019}) and L. moziniana (^{Anguiano-Constante et al. 2021a}). Lastly, our results suggest that the geographic barriers play an important role for the differentiation among Lycianthes moziniana varieties.

Current climatic and soil conditions of the SMOr, SMS and TVB are different. Even though, they have similar vegetation types, the taxonomic composition is different. These facts might promote the morphological differences observed in Lycianthes moziniana varieties. The TVB soils are volcanic in origin, with annual temperature ranging from 12 to 18 ºC and 200 to 4,500 mm of annual precipitation (^{Hernández-Cerda & Carrasco-Anaya 2007}). In contrast, the SMOr has limestone soils derived from sedimentary rocks, with an annual temperature of 12 to 25.5 ºC and 1,000 to 4,000 mm of annual precipitation (^{Hernández-Cerda & Carrasco-Anaya 2004}). The SMS and the TVB have similar soil types, temperature values and precipitation ranges. In the SMS, the temperature varies between 18 and 20 ºC and precipitation ranges from 400 to 4,000 mm (^{Hernández-Cerda et al. 2016}). These may have facilitated ecological responses and contributed to the morphological variation among varieties.

In conclusion, the three varieties displayed an incipient process of speciation. The morphological and geographic differences among the varieties support this conclusion. The morphological evidence was statistically significant and showed that Lycianthes moziniana var. moziniana was more different relative to the others. Lycianthes moziniana var. margaretiana and L. moziniana var. oaxacana were morphologically similar. Furthermore, the climatic and geographical data showed low niche overlap among varieties. Lycianthes moziniana var. moziniana had a widespread range and had intermediate overlap with the other two varieties. Lycianthes moziniana var. margaretiana and L. moziniana var. oaxacana were isolated in different biogeographical provinces and showed the lowest overlap. Previous studies showed DNA sequence variation (^{Dean 2004}, ^{Anguiano-Constante et al. 2021a}). For these reasons, we consider that each variety is a separated evolving metapopulation. Further plastome analyses, could help to decide if each variety merits species status.