Introduction
The Dasylirion genus includes 22 species whose members are commonly called ‘sotoles’ in Mexico (Reyes-Valdés et al., Reference Reyes-Valdés, Palacios, Rivas-Martínez, Robledo-Olivo, Antonio-Bautista, Valdés-Dávila, Villarreal-Quintanilla, Benavides-Mendoza, Grumezescu and Holban2019). It is mainly distributed in the centre and north of Mexico and south of United States of America, places where it is endemic (Fig. 1). The Dasylirion genus has played an important role for the people inhabiting the arid lands of North America since prehistoric times. Paleofecal DNA analyses have demonstrated that sotol plants were used as food before the European colonization (Poinar et al., Reference Poinar, Kuch, Sobolik, Barnes, Stankiewicz, Kuder and Pääbo2001). Since the stems are rich in sugars, the cooked sotol hearts were an important caloric input for survival (Leach and Sobolik, Reference Leach and Sobolik2010). Currently, in several arid lands of Mexico, sotol flowers are used to prepare a soup or are combined with other foods for nourishment (López-Barbosa, Reference López-Barbosa, Contreras-Delgado and Orega-Ridaura2005). Sotol leaves have been used in handicrafts and are currently employed to weave baskets. Also, the basal parts of sotol leaves are widely used in Mexico during religious festivities to craft artificial flowers which adorn churches. The floral scapes are employed to construct fences and walking sticks, which are both light and strong (Reyes-Valdés et al., Reference Reyes-Valdés, Palacios, Rivas-Martínez, Robledo-Olivo, Antonio-Bautista, Valdés-Dávila, Villarreal-Quintanilla, Benavides-Mendoza, Grumezescu and Holban2019). However, the main current use of the genus is the production of an alcoholic beverage called ‘sotol’ as well, resulting from the fermentation and distillation of Dasylirion stems, with denomination of origin for three northern states of Mexico (IMPI, 2002), and it is commercialized both nationally and internationally (Cruz-Requena et al., Reference Cruz-Requena, De La Garza-Toledo, Aguilar, Aguilera-Carbo, Reyes-Valdés, Rutiaga and Rodríguez-Herrera2013). The three main species employed for the production of the sotol drink are D. cedrosanum, D. leiophylum and D. durangense (Reyes-Valdés et al., Reference Reyes-Valdés, Palacios, Rivas-Martínez, Robledo-Olivo, Antonio-Bautista, Valdés-Dávila, Villarreal-Quintanilla, Benavides-Mendoza, Grumezescu and Holban2019). Exploitation of sotol plants is mainly based on collection from natural populations, an activity that may erode genetic diversity and abundance (Reyes-Valdés et al., Reference Reyes-Valdés, Palacios, Rivas-Martínez, Robledo-Olivo, Antonio-Bautista, Valdés-Dávila, Villarreal-Quintanilla, Benavides-Mendoza, Grumezescu and Holban2019).
Female Dasylirion cedrosanum plant growing in a natural area of Northeastern Mexico.

All species in this group are dioecious, and there are no morphological, cytological or molecular criteria to date that can help to distinguish between male and female plants without reproductive organs (Cruz-Requena et al., Reference Cruz-Requena, De La Garza-Toledo, Aguilar, Aguilera-Carbo, Reyes-Valdés, Rutiaga and Rodríguez-Herrera2013; Hernández-Quintero et al., Reference Hernández-Quintero, Reyes-Valdés, Mendoza-Rodríguez, Gómez-Martínez and Rodríguez-Herrera2015; Reyes-Valdés et al., Reference Reyes-Valdés, Hernández-Quintero, Morales-Reyes, Mendoza-Rodríguez, González-Uribe, Ramírez-Godina and Villarreal-Quintanilla2017). Early studies of the genus reported that female flowers show non-functional staminoid vestiges, whereas male flowers show carpelloid structures (Bogler, Reference Bogler1995; Rojas-Piña et al., Reference Rojas-Piña, Olson, Alvarado-Cárdenas and Eguiarte2014).
Currently, phylogenetic relationships are reconstructed by taking advantage of polymorphisms existing in DNA sequences. This is a reasonably fast methodology for acquiring information on a large number of taxa. Furthermore, it allows avoiding errors of taxonomic characterization elaborated with morphological, cytological or biochemical characters (Liu et al., Reference Liu, Ci, Li, Li, Conran and Li2017). The use of chloroplast and nuclear DNA segments has made possible phylogenetic reconstructions. Multiple phylogenetic studies have been based on the chloroplast maturase-K gene (matK) and the large subunit of ribulose-1,5-bisphosphate carboxylase gene (rbcL), which have fairly conserved evolution rates and enough variation to provide useful data in intrafamily analysis (Tamura et al., Reference Tamura, Yamashita and Fuse2004; Selvaraj et al., Reference Selvaraj, Kumar-Sarma and Sathishkumar2008; Chetia et al., Reference Chetia, Phukan and Tamang2016). Molecular phylogenetic reconstruction and the progress of computational statistics tools have made possible to develop and test complex hypotheses concerning the origin and evolution of plants, their diversification and their patterns of phylogenetic relationships. For instance, molecular clock methods combine reconstructed phylogenies from genetic data and fossil records to infer evolutionary time scales (Chen et al., Reference Chen, Kim, Chase and Kim2013; Ho and Duchêne, Reference Ho and Duchêne2014). The relationship between the origin or diversification of endemic species and certain cataclysms has been suggested in the literature and has been made feasible hypotheses to propose their emergence (Scheinvar et al., Reference Scheinvar, Gámez, Castellanos-Morales, Aguirre-Planter and Eguiarte2017; Jiménez-Barron et al., Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020). Furthermore, studies have suggested that genotypes that co-inhabit in similar climatic conditions may influence joint characteristics such as dioecy (Henry et al., Reference Henry, Akagi, Tao and Comai2018). Phylogenetic estimations can provide major insights for further evolutionary biology analysis.
In this work we study the phylogenetic relationships, origin and diversification events through time in the Dasylirion genus, opening a path to understand the biological particularities of this group and to design conservation strategies like safeguarding seeds from rare, evolutionary distinct species.
Materials and methods
Plant materials
Field collection proceeded as follows: seven sotol species were sampled in well-defined areas, with consistent reports of their presence and confirmation of their identity by an expert taxonomist. The collected species comprised D. cedrosanum, D. berlandieri, D. leiophyllum, D. occidentalis, D. palaciosii, D. texanum and D. quadrangulatum. Tissue of at least five plants from each species was collected by cutting sections of three leaves from each randomly selected plant. To best represent possible variation, the sampled individuals typically had a separation of at least 50 m between them. Given the dioecy condition of this group, both sexual forms were included for each species. Each location was represented by an area of at least 10 km2. In all cases, the locations were selected by the consistency of previous reports for species prevalence, whereas marginal areas of distribution were avoided. All the samples were sequenced, but we found sequence homogeneity through samples within species. Thus, only one individual sequence represents each location and species. The specimens were labelled, georeferenced and registered with passport data. Voucher numbers were assigned to the reference plants (Table 1).
Reference location of the Dasylirion sampled plants and their voucher numbers

Herbarium ANSM Universidad Autónoma Agraria Antonio Narro.
DNA extraction and sequencing
Total genomic DNA was extracted from 2 mg of dried leaves by the Lopes et al. method (Reference Lopes, Takasaki, Bostwick, Helentjaris and Larkins1995). All PCR reactions were performed in 50 μl reaction mixtures containing 5 μl of 10× PCR buffer Taq standard (include Mg), 1 μl dNTPs 10 mM, 1 μl of 10 μM of each primer, 0.25 μl of Taq polymerase, 37.75 μl nuclease free H20 and 4 μl template DNA.
The primers used for the matK region were developed particularly for Dasylirion (Orozco-Sifuentes et al., Reference Orozco-Sifuentes, Mendoza-Rodríguez, Hernández-Godinez, Martínez, Villarreal-Quintanilla and Reyes-Valdés2019): matK 335-F 5′ TTTCATTCTCGTTGCGATTAG 3′ and matK 1327-R 5′ ACTTTATTCGATACAAACCCTGTT 3′. For the matK partial gene the reaction mixtures were heated to 94° C for 4 min; this was followed by 40 cycles of 94 °C for 30 s, 57.8 °C for 35 s, 72 °C for 35 s and a final extension at 72 °C for 10 min. The rbcL sequences were amplified with primers reported in literature (Levin et al., Reference Levin, Wagner, Hoch, Zimmer, Nepokroeff, Chris Pires and Sytsma2003; Kress and Erickson, Reference Kress and Erickson2007) as follows: rbcL-F 5′ ATGTCACCACAAACAGAGACTAAAGC 3′, rbcl-R 5′ GTAAAATCAAGTCCACCRCG 3′. For PCR of the rbcL partial gene, the mixes were heated to 94 °C for 4 min; this was followed by 40 cycles consisting of 94 °C for 30 s, 61 °C for 35 s, 72 °C for 35 s and a final extension at 72 °C for 10 min. Purified PCR products were sequenced in both directions with a 3730xl DNA Analyzer (Applied Biosystems Inc. Foster City, CA, USA) at the Advanced Genomic Unit – CINVESTAV, Irapuato, Mexico. Both sequences of five plants per species were concatenated. The clean fragment sequences ranged from 576 to 855 bp for matK and from 516 to 525 for rbcL. All sequences were deposited in the GenBank database (online Supplementary Table S1).
Phylogenetic, molecular clock and Yule analysis of molecular data
Additionally to our own sequencing results, sequences of the same genes for D. micropterum, D. miquihuanense, D. serratifolium and D. wheeleri were retrieved from the GenBank database of NCBI (http://www.ncbi.nlm.nih.gov; online Supplementary Table S1). The GenBank sequence data from D. micropterum and D. miquihuanense were previously uploaded by our research team. The size of the sequences ranged from 812 to 1560 bp for matK and from 511 to 750 for rbcL. GenBank data form several species of the Agave and Yucca genera (A. americana, A. sisalana, A. ghiesbreghtii, A. havardiana, Y. filamentosa, Y. schidigera, Y. necopina and Y. australis) and Acorus gramineus were used to establish external groups and to calibrate the molecular clock. Overall, the study included 20 matK and 20 rbcL sequences to carry out the phylogenetic analysis.
All statistical analyses were done with the language and environment for statistical computing R (R Core Team, 2021), mainly with the ape package v5.4-1 (Paradis, Reference Paradis2012) unless otherwise indicated. Multiple and complete sequence alignments were made from the set of matK + rbcL sequences using Clustal X (Chenna et al., Reference Chenna, Sugawara, Koike, Lopez, Gibson, Higgins and Thompson2003). The phylogenetic tree was constructed thought the maximum likelihood method for optimization of topology and parameters (Felsenstein, Reference Felsenstein1981) with the phangron package (Schliep, Reference Schliep2011). The GTR substitution model was the best fitted for the single and concatenated sequence justified by the lowest value of the Akaike index, based on likelihood test implemented in modelTest function. Cluster reliability in the maximum likelihood tree was evaluated through 1000 non-parametric bootstrap samples (Felsenstein, Reference Felsenstein1985), implemented in the bootstrap.pml function.
The molecular dating was estimated by penalized likelihood and maximum likelihood (Sanderson, Reference Sanderson2002), fitting a chronogram to a phylogenetic tree whose branch lengths represent the rate of substitutions per site. In order to choose an optimal value for the smoothing parameter lambda, we used a cross-validation technique testing values of 10−4, 10−3, …, 103, 104. The lambda value with the lowest cross-validated discrepancy between the predicted and the observed values was chosen. In this technique each terminal branch is removed from the data and then its length is predicted from the reduced data set (Sanderson, Reference Sanderson2002). The tree was dated by two calibration points. A. gramineus with its fossil record provided a crown age ranged in 130–134 Mya (Millions of years ago) (Bremer, Reference Bremer2000; Good-Ávila et al., Reference Good-Ávila, Souza, Gaut and Eguiarte2006). The other calibration point was based on the Yucca genus, with a stem age ranging from 12.86–18.34 Mya according to its fossil record (Tidwell and Parker, Reference Tidwell and Parker1990; Scheinvar et al., Reference Scheinvar, Gámez, Castellanos-Morales, Aguirre-Planter and Eguiarte2017).
Based on the chronogram we determined the diversification of the Dasylirion species by the Yule model (Paradis, Reference Paradis2005), which estimates the diversification rate through maximum likelihood assuming zero extinctions. The pattern of diversification was obtained with the lineage through time (Pybus and Harvey, Reference Pybus and Harvey2000) with a logarithmic scale with the phytools package (Revell, Reference Revell2012). The experimental distribution was tested by the gamma statistics, with the null hypothesis of the Dasylirion clade exhibiting constant rates of diversification, with a normal distribution with mean = 0 and standard deviation = 1 (Pybus and Harvey, Reference Pybus and Harvey2000).
Results
Due to sequence homogeneity among the individuals within species, one individual per species was used for the phylogenetic analysis. In fact, almost all species had identical rbcL and matK sequences among the sampled plants, with a single exception in D. texanum. In this case the most frequent sequence was chosen. This within species homogeneity, is a good indicator that the partial sequences of rbcL and matK are useful tools for intrafamily groups identification through barcode of life in the Dasylirion genus.
The mean sequence length (without gaps) of the aligned data set of the two partial plastid genes of Dasylirion was 1455 bp with 18 polymorphic sites. From the matK sequences with a mean of 894 bp, 16 sites were polymorphic, whereas those of the rbcL segment, with a mean of 561 bp, two sites were polymorphic. Despite the fact that roughly all variation is given by the matK segment, using the rbcL segment increases the chance to gain more robust hypotheses of the relationships among the Dasylirion species.
Among the Dasylirion species, the proportion of the number of sites differing between each pair of sequences ranged from 0 to 0.008, with an average of 0.003. D. palaciosii, D. quadrangulatum and D. miquihuanense did not show any site differences, being a group located in a well-defined geographic area near the boundaries of the states of San Luis Potosí, Nuevo León and Tamaulipas, close to the centre of Mexico in a semiarid region dominated by shrubs (Table 2). The highest proportion of site variation was found between D. berlandieri and D. leiophyllum. These two species are found in two contrasting ecological regions, with D. berlandieri growing in shady coniferous forests in the state of Nuevo León, whereas D. leiophyllum grows in a sunny, extremely arid region of the states of Coahuila and Chihuahua (Table 2).
Proportion of the number of sites differing among Dasylirion species for matK and rbcL partial sequences

D.ced= D.cedrosanum; D.ber, D.berlandieri; D.mic, D.micropterum; D.miq, D.miquihuanense; D.ser, D.serratifolium; D.whe, D.wheeleri; D.lei, D.leiophyllum; D.occ, D.occidentalis; D.pal, D.palaciosii; D.tex, D.texanum; D.qua, D.quadrangulatum.
Phylogeny and origin of Dasylirion genus
Figure 2 shows the phylogenetic tree of Dasylirion species with merged data of matK + rbcL sequences estimated by maximum likelihood and the GTR substitution model. The phylogenetic tree, constructed through maximum likelihood and linearized through time by penalized likelihood, is represented with a scale of millions of years. The topology of the tree shows that both A. gramineus and Agave and Yucca species appear outside of the group that comprises all Dasylirion species studied here with a 100% bootstrap support. The node corresponding to the putative point of divergence between the Dasylirion genus and the two reference genera, Agave and Yucca, shows an age of 78.1 Mya (BS = 100%). Our estimate for a lowest limit of the crown age of Dasylirion was 5.46 Mya (BS = 100%; Miocene 23 to 5 Mya), and the ages for two speciation events were placed in 3.14 and 2.31 Mya ago (Pliocene 5 to 2.5 Mya; Fig. 3).
Chronogram of Dasylirion genus and outgroups formed by Acorus gramineus, and Agave and Yucca species. Bootstrap values (blue) are shown above each branch and estimate ages (red) are shown at nodes. The scale bar is represented in million years. The tree was optimized with maximum likelihood and dated with penalized likelihood.

Chronogram of Dasylirion genus. Bootstrap values (blue) are shown above each branch and estimated ages (red) are shown at nodes. The scale bar is represented in million years. The tree was optimized with maximum likelihood and dated with penalized likelihood.

The D. cedrosanum and D. leiophyllum branch (BS = 92%) formed a sister clade to the rest of the species with a crown age of 2.32 Mya (Fig. 3). Both of them grow in desert scrub in northern Mexico. A second clade (BS = 61%) includes the remaining Dasylirion species. The crown age dated in this analysis for D. texanum and D. occidentalis (BS = 52%) was 1.64 Mya. Both species grow in rocky substrates in northern Mexico. An estimated divergence between D. wheeleri and D. serratifolium (BS = 95%) is estimated to have occurred 0.78 Mya; although they are within the same group, they belong to geographically remote regions. Since these two sequences were obtained from GenBank database of NCBI, geographic misplacement or identification issues cannot be discarded. In the same clade, D. micropterum and D. berlandieri (BS = 61%) appear with a crown age of 1.31 Mya and tend to grow in mountains, near or within coniferous-oak forests. The most recently divergence estimated in this study was for D. miquihuanense, D. palaciosii and D. quadrangulatum with a crown age of 0.004 Mya. Although, the clade of D. miquihuanense + D. palaciosii + D. quadrangulatum has a low bootstrap value (~13%) these species showed identical sequences. The three species grow in a well-defined semiarid region dominated by shrubs in Northeastern Mexico.
Substitution rates and diversification of Dasylirion genus
Accordingly with the time-linearized tree, we estimated an average substitution rate of 0.0005 nucleotides per million years for the concatenated sequence. The diversification rate of Dasylirion estimated from its origin to the present and under the Yule model was S = 0.0466 ± 0.011 species per million years, whereas the expected number of lineages was under a constant diversification rate as supported by the gamma statistic which scored 0.3512 (P value = 0.7254).
Discussion
Phylogenetic reconstruction based on molecular data relies directly on the polymorphic sites contained in the DNA fragments. Plastid genes as rbcL commonly present low percentages of polymorphic sites given their level of genetic conservation (Chetia et al., Reference Chetia, Phukan and Tamang2016). Reports of Agave sensu lato plastid genes show of 29 polymorphic sites from a sequence of 3670 bp (Jiménez-Barrón et al., Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020). In our study, 18 polymorphic sites for the concatenated plastid segments of 1456 bp were found and almost all variation was due to the matK partial gene (16 polymorphic sites). Information given by the DNA segments used in our study was high enough to distinguish intrafamily groups. Using more than one gene in phylogenetic analyses performed in previous works suggests that quality increases compared to single gene analyses (Hufford et al., Reference Hufford, Moody and Soltis2001; Tamura et al., Reference Tamura, Yamashita and Fuse2004). In fact, we observed a low substitution rate in the rbcL partial sequence, polytomy and low resolution trees constructed with the individual gene. The number of polymorphic sites and the low evolutionary rate shown by the matK and rbcL partial genes, allow to develop robust hypotheses of relationships among species. Moreover, data provides acceptable support values for most of the discovered clades, having support values above 50% (Liu et al., Reference Liu, Ci, Li, Li, Conran and Li2017).
The estimated topology by maximum likelihood achieved maximum bootstrap support values in the three for inferred clades of Dasylirion, Agave + Yucca and A. gramineus (Fig. 2). Dasylirion genera are separated into a monophyletic clade and Agave + Yucca resulted within the same paraphyletic clade, which are in accordance with previous studies (Rojas-Piña et al., Reference Rojas-Piña, Olson, Alvarado-Cárdenas and Eguiarte2014; Jiménez-Barron et al., Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020).
The molecular phylogenetic tree showed D. cedrosanum and D. leiophyllum as a separated sister group from the remaining Dasylirion species (Fig. 3). Interestingly, these two species are widely used in the sotol distillated drink preparation. They have the peculiarity of having stems with a high monosaccharide content, desirable in the distillate industry (Cruz-Requena et al., Reference Cruz-Requena, De La Garza-Toledo, Aguilar, Aguilera-Carbo, Reyes-Valdés, Rutiaga and Rodríguez-Herrera2013; Reyes-Valdés et al., Reference Reyes-Valdés, Palacios, Rivas-Martínez, Robledo-Olivo, Antonio-Bautista, Valdés-Dávila, Villarreal-Quintanilla, Benavides-Mendoza, Grumezescu and Holban2019). The clade group separation of D. cedrosanum and D. leiophyllum is not visible in the chloroplast phylogeny constructed by Bogler (Reference Bogler1995). They inferred a separation of the species into sister clades due to their geographic locations. For instance, three of the eleven Dasylirion species not analysed in this study as D. acrotrichum (Schiede ex Schultes) Zucc, D. parryanum Trel and D. lucidum Rose belong to a clade separated by their location in southern Mexico (Bogler, Reference Bogler1995). This three Dasylirion species are in the D. serratifolium clade which belong to the same region.
Lastly, the paraphyletic clade formed by D. berlandieri, D. palaciosii and D. quadrangulatum, was in accordance with Bogler (Reference Bogler1995).
Previous studies have estimated ages of origin for Asparagales that vary widely across literature. The highest suggested age is 125.3 Mya (Tank et al., Reference Tank, Eastman, Pennell, Soltis, Soltis, Hinchliff and Harmon2015) while the lowest is in the range of 69–60 Mya (Good-Avila et al., Reference Good-Ávila, Souza, Gaut and Eguiarte2006). Our estimated age of divergence between the Dasylirion and Agave genera of 78.14 Mya is placed within this range, and it is close to the previous estimated divergence between the two taxa in approximately 62 Mya (Jiménez-Barron et al., Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020). Compared with the estimated origin of Asparagales, Dasylirion is of recent emergence (minimum age of 5.46 Mya). A previous approximation to the origin of Dasylirion has been made in a recent work (1.2 ± 1 Mya) (Jiménez-Barron et al., Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020), where the genus was used as an external group in the reconstruction of Agave phylogenies. The Jiménez-Barron et al. (Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020) estimated is low, compared to ours. However, it should be considered that in this work 50% of the 22 known Dasylirion species are included, whereas Jiménez-Barron et al. (Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020) used only three species aimed to be used as a reference group.
The results indicate that the diversification rate is low in Dasylirion (S = 0.0466 ± 0.011 species per million years), compared to others estimates of species of the Asparagaceae family that suggest a remarkable diversification rate of 0.51 in Agave or 0.21 in Yucca (Good-Ávila et al., Reference Good-Ávila, Souza, Gaut and Eguiarte2006). Moreover, it is low compared to Flores-Abreu et al. (Reference Flores-Abreu, Trejo-Salazar, Sánchez-Reyes, Good, Magallón, García-Mendoza and Eguiarte2019) who estimated a diversification rate of 1.478 and 0.473 species/Mya for Agave and Yucca, close to those found by Jiménez-Barron et al. (Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020) of 1.50 for Agave sensu lato and 0.12 for Yucca. The number of sampled species has an important effect to obtain the diversification rate (Wertheim and Sanderson, Reference Wertheim and Sanderson2011). The Dasylirion species included in this study, belong to populations with fair to great abundance with various traditional uses in the arid lands of Mexico. Although, in phylogenetic analysis it is rarely possible to sampled all the species of a genus, it is important to use the largest number of species for more conclusive results.
The species not analysed in this study are: D. acrotrichum (Schiede ex Schultes) Zucc, D. durangense Trel, D. gentryi Bogler, D. glaucophyllum Hook, D. longissimum Lem, D. longistylum MacBride, D. lucidum Rose, D. parryanum Trel, D. sereke Bogler, D. simplex Trel and D. treleasei (Bogler) Hochstätter.
The ecological opportunity and speciation proposed for Dasylirion given by the inferred diversification rate, could have taken place during the late Miocene (5.3 Mya), which has been linked to fluctuations in the climate of North America implying the appearance of arid zones attributed to geological activity (Mastretta-Yanes et al., Reference Mastretta-Yanes, Moreno-Letelier, Piñero, Jorgensen and Emerson2015). The volcanic structures emerging in the Pleistocene (0.01 to 2.58 Mya) such as the Trans-Mexican Neovolcanic Belt, caused changes in humidity and wind flow, with the appearance of new types of plant communities adapted to these dry climates within the Mexican highlands (Jiménez-Barron et al., Reference Jiménez-Barron, García-Sandoval, Magallón, García-Mendoza, Nieto-Sotelo, Aguirre-Planter and Eguiarte2020; Mastretta-Yanes et al., Reference Mastretta-Yanes, Moreno-Letelier, Piñero, Jorgensen and Emerson2015).
Members of the Asparagaceae family such as Dasylirion, Calibanus, Beaucarnea and Nolina (sub-family Nolinoideae) occurring in the dry climate areas of the Mexican highlands, have the peculiarity to exhibit dioecious sexual systems (Rojas-Piña et al., Reference Rojas-Piña, Olson, Alvarado-Cárdenas and Eguiarte2014). If dioecy is expected to occur via diverse paths (Henry et al., Reference Henry, Akagi, Tao and Comai2018), it may be possible that the ecological opportunities that Dasylirion ancestors had in a new niche could have induced the appearance of separated sexes. This kind of assumption is not new. Related hypotheses have been outlined for shrubs that develop in adverse climatic regions with conditions of hydric stress (Aronne and Wilcock, Reference Aronne and Wilcock1996). It has been suggested that these adverse environmental factors influenced the evolution of dioecious sexual systems, because the separation of sexes allows a reduction in the total energy consumption of the plant (Aronne and Wilcock, Reference Aronne and Wilcock1996). This energy saving may be is necessary for the adaptation and survival probability of the species to new environmental conditions. Sexual differentiation might be genetically controlled nevertheless the environment is a big factor to trigger (Bawa, Reference Bawa1980; Renner, Reference Renner2014). However, data suggests that environment differences do not have any role in the sex determination in Dasylirion (Reyes-Valdés et al., Reference Reyes-Valdés, Hernández-Quintero, Morales-Reyes, Mendoza-Rodríguez, González-Uribe, Ramírez-Godina and Villarreal-Quintanilla2017).
Our results are consistent with phylogenetic dating from previous studies, albeit they indicate that the Dasylirion genus is a group with a low diversification rate. Nonetheless, this genus has been able to gain successful adaptation to several arid-land environments. Moreover, the close relationship between two species highly used for the drink production, is an indicative of the potential genetic basics for their productivity.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S1479262122000181.

