Population expansion and genetic structure in Cephalocereus nizandensis (Cactaceae), a microendemic cactus of rocky outcrops of the Tehuantepec basin, Mexico

Background and aims – Cephalocereus nizandensis is a microendemic columnar cactus that grows isolated in xerophytic enclaves associated with rocky outcrops in the Isthmus of Tehuantepec, in the south of Mexico. Its demographic history and genetic structure were assessed to determine the main events that shaped its current restricted distribution. Material and methods – Chloroplast intergenic sequences of 40 individuals and inter simple sequence repeats (ISSRs) of 45 individuals from four isolated populations were used to estimate haplotypic and nucleotide diversity, using expected heterozygosity and the Shannon index. AMOVA, population pairwise FST, and Bayesian clustering analyses were performed to explore the genetic structure. Demographic history was estimated with neutrality tests, mismatch distribution analysis, and Bayesian skyline plots. Phylogenetic relationships and divergence times were determined using a median joining network and a Bayesian molecular clock. Key results – C. nizandensis has a high diversity and moderate genetic differentiation. The lowest elevation locality was found to be the most genetically distinct. The species has undergone a process of population expansion that began 150,000 years ago and has remained without evidence of a population contraction in the transition from the Pleistocene to the Holocene (11,700 years ago). Conclusions – C. nizandensis presents moderate but significant genetic differentiation, which may be due to an early divergence of its populations. Currently observed levels of genetic diversity are the result of historical maintenance of high population sizes and a population expansion approximately in the last 150,000 years, which was sustained independently of the climatic fluctuations of the Early Quaternary, due in part to the stability of the rocky habitat.


INTRODUCTION
Rocky outcrops are sites where erosion has removed all the soil layers and the bedrock is exposed. They provide various microclimates that generate island habitats, increasing the heterogeneity in the landscape, levels of endemism and β diversity (Twidale 1982;Gibson et al. 2012). Plant species that occur on rocky outcrops are usually adapted to the edaphic conditions of their habitats, while being absent or rare in surrounding habitats (Gibson et al. 2012). In these habitats, the humidity is low, there is almost no water in the soil, and the exposure to the wind is higher compared to neighbouring areas. The concentration of phosphorus and nitrogen is very low as well (Porembski 2007).
Harsh edaphic conditions on outcrops hamper the establishment of seedlings that are not adapted to environmentally extreme conditions. Some adapted traits common in outcrop species are succulent water-storage trunks or caudex, and desiccation-tolerant vegetative tissues (Porembski 2007).
The inorganic and harsh nature of the outcrops provides them with long-lasting landscape characteristics and stable microclimatic conditions; thus, the rocky outcrops can function as ecological refugia for hundreds of thousands of years (Couper & Hoskin 2008;García et al. 2020). This long-term stability of rocky outcrop habitats has allowed the persistence of species and communities in latitudes where their presence cannot be explained by the regional climate. For example, in tropical regions of the Americas (Sarthou et al. 2003;Torres-Ribeiro et al. 2007;Silva et al. 2020), and Africa (Porembski 2007) where rainforests and cloud forests are the most common vegetation, rocky outcrops provide conditions for species adapted to xeric environments (Porembski 2007;Locosselli et al. 2016).
Rocky outcrops host communities of great floristic, ecological, and evolutionary interest. Plant communities of rocky outcrops have an island-like distribution and are thus good study systems to test different hypotheses related to their formation, distribution, and maintenance, particularly with regard to the xerophytic species in these communities (Porembski & Barthlott 2000;Byrne et al. 2019). Consequently, the vegetation of rocky outcrops has been highly studied around the world (Moraes et al. 2005;Bonatelli et al. 2014;de Aguiar-Campos et al. 2020). For example, in eastern Brazil's seasonally dry tropical forests, limestone outcrops have been found to have functioned as litho-refuges during the Quaternary (de Aguiar-Campos et al. 2020). In these ecosystems, population genetics and phylogeographic studies have been conducted in some rupicolous cacti, such as Praerocerus euchlorus (F.A.C.Weber ex K.Schum.) N.P.Taylor and the Pilosocereus auricetus (Werderm.) Byles & G.D.Rowley complex, which went through population expansion during the Quaternary glacial periods (Moraes et al. 2005;Bonatelli et al. 2014).
Some common characteristics that have been observed in studies of plant genetic variability in rocky outcrop species include moderate genetic diversity within populations, demographically stable populations, and little gene flow between locations (e.g. in Praerocerus machrissi Dawson and P. euchlorus, Moraes et al. 2005; Eucalyptus caesia Benth., Byrne & Hopper 2008; Acacia woodmaniorum Maslin & Buscumb, Millar et al. 2017; Pilosocereus auricetus, Bonatelli et al. 2014). These demographically stable populations have been maintained for long periods of time suggesting that genetic drift has not been intense in the rocky outcrop populations. In addition, the pattern of isolation by distance between populations is not common, which would suggest an old isolation pattern (Millar et al. 2017). The values of diversity and genetic structure show that the local populations of the species underwent expansion and contraction dynamics, which suggest these outcrops as possible vegetation refuges during the glacial and interglacial periods (Couper & Hoskin 2008;García et al. 2020).
In southern Mexico, on the Isthmus of Tehuantepec, limestone outcrops are found, with associated xerophilous vegetation with a high endemic richness (Pérez- García & Meave 2005;Pérez-García et al. 2009). There are similarities between the rocky outcrops of the Nizanda region in the Isthmus of Tehuantepec and the desert system of North America (for example, the Chihuahua and Sonora deserts; Pérez- . This similarity makes the Nizanda region an area of biogeographic and evolutionary importance because it is one of the narrowest low-elevation areas that connect the Nearctic and Neotropic biota of North and South America (Nizanda has an average elevation of 140 m a.s.l. and Sierra Madre de Oaxaca, which connects this area, is 3270 m a.s.l.; González-Medrano 1996). In addition, due to its high floristic richness, and its important endemic component it has been proposed that this region may have hosted refuges for Mexican xerophytic flora and limestone outcrops seem good candidate habitats for such refuges (González-Medrano 1996;Pérez-García & Meave 2005;Pérez-García et al. 2010). However, there are no studies on the genetic diversity and the structure or demographic history of any of the species inhabiting these xerophytic enclaves.
In this study, we assessed the population genetics of a dominant and endemic cactus species of these xerophytic enclaves, Cephalocereus nizandensis (Bravo & T.MacDoug.) Buxb., Cactaceae. This species grows in limestone outcrops on mountain tops and slopes. It only develops in cavities formed in calcareous rocky outcrops, where organic debris and water accumulate (Bárcenas-Argüello et al. 2010). We hypothesized that the species would display signatures of population expansion processes associated with glacial periods, as observed in similar systems in eastern Brazil (Moraes et al. 2005;Bonatelli et al. 2014), or of longlasting persistence with isolation, as has been found in other species dwelling on rocky outcrops (e.g. Byrne & Hopper 2008;Millar et al. 2017). In this study we evaluated (1) how much genetic variability this species contains, and how is it distributed, (2) when its populations arose historically, and (3) how its historical demographic dynamic in the middle Pleistocene and in the transition from the Pleistocene to the Holocene.

Study area
The study area is located in the southern state of Oaxaca, Mexico, in the southeast of the Tehuantepec Isthmus region ( fig. 1A, B). The climate is tropical, sub-humid, and highly seasonal. The elevation ranges between 90 and 700 m. The geomorphology is dominated by schist hills abruptly interrupted by exposed limestone outcrops. The geological history of the area shows that this part of Oaxaca was separated in the early Miocene, when the Tehuantepec Fracture caused a significant reduction in highlands (Barrier et al. 1998). The dominant vegetation is tropical dry forest, with considerable savanna areas; but as a result of the edaphic aridification processes facilitated by the limestone outcrops, there are patches of xerophytic shrubs where species such as Hechtia rosea, Comocladia engleriana, Pseudosmodingium multifolium, Mammillaria voburnensis, M. albilanata, and Agave ghiesbreghtii are common (Pérez- García & Meave 2005). In this region, four isolated patches of Cephalocereus nizandensis were chosen for sampling near the towns of Nizanda and La Mata and were named after the closest town ( fig. 1C). The coordinates, elevation, and sample sizes are indicated in table 1; each of these localities was considered a population.

Study species
Cephalocereus nizandensis (Bravo & T.MacDoug.) Buxb. (Buxbaum 1965) is a columnar growth cactus, rarely branched, and approximately 3 m in height and 15 cm in diameter. It is characterized by a woolly, differentiated reproductive zone in the apex, called pseudocephalium (Bravo-Hollis & Sánchez-Mejorada 1978;Anderson 2001). This species is micro-endemic to the Nizanda region of Oaxaca, Southern Mexico. It only develops in a few xerophytic enclaves, restricted to a few isolated populations on the peaks and shoulders of an exposed limestone corridor, narrower than 10 km, where the soil is barely a shallow accumulation of organic material inside some cracks (Pérez- García & Meave 2005;Bárcenas-Argüello et al. 2010). These cacti grow both in xerophytic shrubs and tropical dry forest on rocks. However, they are more abundant in xerophytic shrubs, where they are the densest species of the high stratum (Pérez- García & Meave 2005).

Sampling and production of genetic data
Approximately 100 g of photosynthetic tissue per individual was obtained. Only reproductive individuals were sampled, and 1.5 m were left between sampled individuals. Tissue was stored at -80°C until DNA extraction, which was performed using a DNeasy Plant MiniKit (Qiagen, Hilden, Germany). To explore variability in the chloroplast genome, two intergenic regions of chloroplast DNA (cpDNA) were amplified, rpl32-trnL (UAG) (Shaw et al. 2007) and trnT-trnL (Taberlet et al. 1991). The amplification program was as follows: 60 s of initial denaturation at 94°C; 15 s of denaturation at 94°C, 15 s of alignment with temperatures according to supplementary file 1; 15 s extension at 72°C; 30 cycles were used for rpl32-trnL region and 40 for trnT-trnL. A single sample of both Cephalocereus apicicephalium E.Y.Dawson and Cephalocereus totolapensis (Bravo & T.MacDoug.) Buxb. were collected and included in the amplification (samples were collected on highway 190 Juchitán-Oaxaca, near the town of Santa María Jalapa del Marqués). Both species are in the same clade as C. nizandensis (Bravo-Hollis & Sánchez-Mejorada 1978) and were included as outgroup. The cpDNA PCR products were purified and sequenced by Macrogen, Inc. (Seoul, South Korea). All amplified samples were sequenced in the forward direction, and two individuals per population were sequenced in the reverse direction to verify the sequences.
The sequences were aligned in MUSCLE software (Edgar 2004) and edited in MEGA (Kumar et al. 2016). Sequence ends, 5' and 3', were trimmed enough to avoid the introduction of spurious variability and the variable sites were contrasted against the electropherograms.
Regarding the nuclear genome, a set of nine ISSR primers was tested (in this case, samples of C. apicicephalium and C. totolapensis were not included). Of these primers, the five that showed variability and reproducibility were used (supplementary file 1). PCRs were carried out in a 25 µL volume for 30 cycles. DNA polymerase Taq Mastermix (RADIANT; Alkali Scientific Inc., Pompano Beach, FL, USA) was used according to the manufacturer's instructions. For ISSR amplification, a touchdown PCR was performed. The starting temperature varied from 54 °C to 56 °C, depending on the primer used (supplementary file 1), and the temperature was reduced by 1 °C in the first nine cycles until it reached 45 °C to 47 °C. ISSR products were separated using horizontal electrophoresis in a 1.5% agarose gel with Tris-Borate-EDTA buffer and visualized using GelRed (Biotium, CA, USA) intercalating agent in a UV light imaging system, Mini-Bis Pro (DNR Bio Imaging Systems, Neve Yamin, Israel), to create a presence/absence database where each band represented a locus. A 100 bp DNA ladder (Omega Bio-TEK M01-02, Atlanta, GA, USA) was used to determine the size of the bands.
The genetic structure of C. nizandensis was assessed for both marker types using analysis of molecular variance (AMOVA) and paired F ST in Arlequin v.3.5.5 (Excoffier & Lischer 2010). The AMOVA was initially used to test for differences between all populations, considering these as part of the same hierarchical level, and later, regionalization was tested according to the paired F ST results. In addition, Bayesian clustering was used to assess the genetic clustering of individuals without a priori consideration of geographic location. BAPS v.6.0 (Corander et al. 2003) andSTRUCTURE v.2.3 (Pritchard et al. 2000) were used for cpDNA and ISSR markers, respectively. BAPS was used to test for K = 2 to K = 6 cpDNA groups with 10,000 iterations, 10% burn-in period, and 10 replicates. STRUCTURE was used to test for 1 to 10 ISSR groups with 50,000 iterations, 10% burn-in period, and 20 replicates for each K. Default values were used for all other parameters. The most likely number of groups K assigned by BAPS was determined with the posterior probability distribution (Corander et al. 2003) and the most likely K for STRUCTURE analysis was estimated using the Evanno et al. (2005) method in STRUCTURE HARVESTER (Earl & von Holdt 2012) and plotted with DISTRUCT v.1.1 (Rosenberg 2004).

Phylogenetic inference and divergence time
Bayesian inference was used to estimate the phylogenetic relationships ) of the rpl32-trnL/trnT-trnL sequences. The ingroup comprised 14 sequences corresponding to the C. nizandensis haplotypes, the outgroup was the same as for the Network analysis.
Before doing the phylogenetic reconstruction and with the sequences already aligned, the nucleotide substitution model was defined in the program MEGA v.7.0 (Kumar et al. 2016). In this case, the model with the lowest value of the Bayesian information criterion (BIC) was the HKY (Hasegawa et al. 1985) with a gamma distribution and the shape parameter a = 0.05 (Kumar et al. 2016).
Phylogenetic reconstruction was performed in BEAST v.1.10.4 along with the associated software suite ) as described below. First, in the BEAUTi v.1.10.4 program ) the following information was defined: the evolutionary model already mentioned HKY + gamma (a = 0.05) with a substitution rate of 1.1 × 10 -9 to 1.6 × 10 -9 substitutions per site per year (Wolfe et al. 1987) with a strict molecular clock. The divergence time of the Pachycereinae subtribe (5.28 Mya (95% HPD 7.74-3.47); Hernández-Hernández et al. 2014), to which the genera Cephalocereus and Carnegiea belong, was used as a calibration point. The length of the Markov chain was defined as 20,000,000 steps sampled every 1,000 steps. The output file of BEAUTi was used to run BEAST program, and the results of the Markov chains generated in BEAST were visually verified in Tracer v.1.7.1 (Rambaut et al. 2018). Additionally, BEAST produces an output file with a set of phylogenetic trees that is summarized in the program TreeAnnotator v.1.10.4 , this program finds the best supported tree, with a burn-in of 10%. Using this tree, it was possible to generate information of the posterior probabilities of the nodes, and the HPD (highest posterior density) limits of the node heights. Lastly, the best supported tree was drawn and edited in FigTree v.1.4.4. (Rambaut 2018).

Demographic history based on cpDNA
The demographic history was investigated using four indices.
(1) We used F s neutrality tests (Fu 1997), whose negative values indicate a recent increase in population size under the assumption of neutrality. (2) Using the mismatch distribution analysis (Rogers & Harpending 1992), we calculated the sum of squared differences (SSD) between the observed and the expected mismatch distributions, and (3) Harpending's Raggedness Index (R). These indices take large values for multimodal distributions commonly found in stationary populations, whereas unimodal and smoother distributions are typical of expanding populations (Rogers & Harpending 1992). These analyses were performed using Arlequin v.3.5.5 (Excoffier & Lischer 2010). (4) Historical demographic processes were also inferred by Bayesian Skyline Plot analysis (BSLP), implemented in BEAST , which indicates changes in effective population size over time (Drummond et al. 2005). A strict molecular clock was used, with a mutation rate of 1.1 × 10 -9 to 1.6 × 10 -9 substitutions per site per year (Wolfe et al. 1987) and an MCMC chain with 10,000,000 replicates sampled every 1,000 steps, with a burn-in of 10%.

Populations
Lat  from the M2 population and was found in 18 individuals (48% of the sample). The H04 haplotype was the second most common and was found in 8 individuals (20%), all from population M2, but it was also found in two individuals of C. apicicephalium and C. totolapensis tested as the outgroup. The H09 haplotype was found in 3 individuals (7.5%) and was private to the N2 population. The rest of the haplotypes were only found in a single individual ( fig. 1C).
In the median joining network, most haplotypes are closely related, with one mutational step separating them. Haplotypes from the same or a nearby population were located close to each other in the network, except for the H03 and H06 haplotypes ( fig. 1D). Abundance, distribution and connectivity of the H02 haplotype are characteristic of population expansion; this haplotype occupies a central position in the network, with multiple connections with haplotypes of the N1 (H09, H10, and H12) and M2 (H04) populations. Two groups outside the network centre are distinguished: one represents the H01, H05, and H06 haplotypes, found in the M1 and M2 populations, and the other group is formed by the haplotypes H03, H08, H11, Figure 2 -Cephalocereus nizandensis haplotype chronogram (A) and genetic clusters (B and C). A. Bayesian inference chronogram. The estimated divergence time for the rpl32-trnL/trnT-trnL sequences is illustrated. Posterior probability (on the branches); branch age, (bold text next to nodes), 95% highest posterior density interval (blue bars). Haplotype H04 includes the sequence of species C. apicicephalium and C. totolapensis. Ccol = C. columna-trajani, Cgig = Carnegiea gigantea. The boxes correspond to the genetic groups of shown in B in this figure, and the colours of the haplotypes correspond to fig. 1. B-C. Genetic groups. rpl32-trnL/trnT-trnL clusters obtained with BAPS (B) and ISSR clusters obtained with STRUCTURE (C), each bar represents an individual, and the colour indicates the probability of an individual belonging to a given group. H12, H13, and H14, all of which are from the N1 and N2 populations ( fig. 1D).
Cephalocereus nizandensis has moderately high levels of genetic diversity. Twelve polymorphic sites were found in the cpDNA sequences, of which six were parsimony informative. The total nucleotide diversity (π) was 0.0013; it varied from 0.0007 in M1 to 0.0015 in M2, and the total haplotype diversity was 0.75 (0.8 in M2 to 0.25 in M1). Seventy-seven ISSR loci were found in total: 76 (98.7%) were polymorphic in at least one population, and 19 loci (14.63%) were private. Both the expected heterozygosity and Shannon diversity index varied little between populations (H e = 0.296; I = 0.45; table 1).
The populations are genetically differentiated. Most of the variability lay within populations, although there was significant genetic differentiation among them, both for the chloroplast genome (F ST = 0.209; p < 0.05) and for the nuclear genome (F ST = 0.186, p < 0.005; supplementary file 3). Paired F ST indicated that for the cpDNA, populations M1, N1, and N2 showed no differentiation between each other, while the only comparisons significantly different from

Phylogenetic inference and divergence time
Bayesian analysis ( fig. 2A) showed a first divergence in the root of the tree where the genus Cephalocereus (with a posterior probability of PP = 0.91) is separated from Carnegiea gigantea 3.97 Mya (95% HPD 5.88-2.59). Next, in a clade with a posterior probability of PP = 1.0, C. nizandensis, C. totolapensis, and C. apicicephalium (sharing the H04 haplotype) were grouped and separated from C. column-trajani 2.15 Mya (95% HPD 3.47-1.11) in the Early Pleistocene. We can also see in the tree that haplotypes found in this work formed a clade with high statistical support (PP = 1.0). Subsequently, two clades were formed dated to 1.37 Mya (95% HPD 2.27-0.686), the first contains the haplotypes

Historic demography based on cpDNA
Our results suggest that C. nizandensis has undergone a process of population expansion, which was mainly registered in the N1 population, while other populations such as M1 showed demographic stability. The neutrality test, Fu's F s , indicated population expansion for the whole species (F s = -7.224, p < 0.001) and for the N1 population (F s = -2.594, p = 0.011; table 2). The mismatch distribution analysis for the entire species' dataset fits the spatial expansion, but there were significant differences from the demographic expansion model. At the population level, N1 was the only population that displayed a signature of spatial expansion (table 2). In the same sense, the Bayesian Skyline Plot (BSLP) showed population expansion from 150,000 years to the present, which increased the effective population size from N e = 100,000 to N e = 1,000,000 individuals, without evidence of contraction for the last 250,000 years ( fig. 3).

Genetic diversity and structure
Cephalocereus nizandensis' genetic variability levels are within values registered for other cacti, but this cactus is distributed in a small area with few populations. The number of cpDNA haplotypes reported for cacti is variable (cpDNA of 11 species of the Cactoideae subtribe; table 3), ranging from 2 to 23. It is important to highlight that the number of haplotypes found in C. nizandensis (14) is close to the average, although almost all of these haplotypes have a restricted distribution to only one of the four sampled populations. The haplotype diversity in C. nizandensis (H d = 0.75) is considerably higher than that in other cacti (average of 0.415), and its nucleotide diversity (π = 0.0013) is lower than the average (π = 0.0022). However, its π values are similar to or higher than those reported in studies that evaluated the same cpDNA regions (Gitzendanner & Soltis 2000;Clark-Tapia & Molina Freaner 2003;Bonatelli et al. 2014;Ornelas & Rodríguez-Gómez 2015;Cornejo-Romero et al. 2017). It should be noted that C. nizandensis is as diverse as other Cephalocereus species. Nucleotide diversity in C. nizandensis was greater than that in Cephalocereus columna-trajani (π = 0.00023; Cornejo-Romero et al. 2017). The expected heterozygosity in C. nizandensis (H e = 0.296) was similar to that found for RAPD and ISSR markers of nine cactus species (average H e = 0.26) and slightly higher than that in C. totolapensis (H e = 0.233; Palleiro-Dutrenit 2008); the percentage of polymorphism (98.6%) was higher than the average (76% in RAPD), greater than that in C. totolapensis as well (67.8%; Palleiro-Dutrenit 2008) and similar to that reported in Neobuxbaumia mezcalaensis (Bravo) Backeb. (97.2%, table 3; Rivera-Montoya 2003; Tapia et al. 2017).

2002) or self-incompatibility
The study shows that C. nizandensis presents moderate but significant genetic differentiation, which may be due to an early divergence among its populations. cpDNA pairwise F ST analyses indicated that the most genetically different individuals were found in the M2 population ( fig. 2B), which is also found in the lowest elevation location. Environmental differences, such as humidity, temperature or places to germinate, could explain a pattern of local adaptation of plants in different sites Valiente-Banuet et al. 2002;Nassar et al. 2003;Gutiérrez-Flores et al. 2016), even giving rise to natural selection processes according to the microclimatic conditions of each population (Kim & Donohue 2013). Similarly, the genetic structure of Stenocereus stellatus Riccob. has been found to be determined by the evolutionary processes of the landscape, as this determines the recruitment dynamics and gene flow restriction in seeds (Cornejo-Romero 2004;Bustamante et al. 2016). To confirm the existence of local adaptation, it would be necessary to perform specific studies aimed at examining phenotypic differences, such as germination, and/or the genes subjected to natural selection, as well as the measurement of microclimatic variables, that may explain the differentiation (

Divergence time and demography history
Cephalocereus genus is divided into two subgenera: Neodawsonia (including three species; C. totolapensis, C. apicicephalium, and C. nizandensis) and Cephalocereus (including two species C. senilis Pfeiff. and C. columnatrajani) (Bravo-Hollis & Sánchez-Mejorada 1978), and it is probable that these subgenera diverged approximately 2.5 Mya (95% HPD 3.06-1.24; fig. 2A), during the transition from the Pliocene to the Pleistocene. During this period, the Late Pliocene climate optimum (3.3-2.4 Mya; Tripati et al. 2009) has been related to the most recent radiation of succulent plants (Good-Avila et al. 2006;Arakaki et al. 2011). Recent diversification is a common characteristic in cacti (Arakaki et al. 2011;Hernández-Hernández et al. 2014), as observed in Pilosocereus Byles & G.D.Rowley (Perez et al. 2016) and Cereus Mill. (Franco et al. 2017;Silva et al. 2020), both genera of South American deciduous tropical forests. The recent divergence time of C. nizandensis may explain the retention of ancestral haplotype H04 found in the M2 population of C. nizandensis, C. apicicephalium, and C. totolapensis (figs 1D, 2A). This haplotype could have been lost in the other populations of C. nizandensis due to founder events (Ellstrand & Elam 1993).
Other divergences within C. nizandensis were estimated to be approximately 0.968 Mya (95% HPD 1.7-0.363) and 0.833 Mya (95% HPD 1.43-0.362; fig. 2A). This period coincides with the most arid stage of the Middle Pleistocene (Russo-Ermolli & Cheddadi 1997). It must be considered that, in recent divergence processes, genetic divergence usually precedes population divergence (Arbogast et al. 2002), which suggests that the Quaternary climate cycles could have influenced the structure of the species (Nistelberger et al. 2015). These estimates precede the climatic events of the last glacial cycle (Last Glacial Maximum LGM, 0.021 Mya), which have been used to explain the population dynamics and genetic structure of plants in arid and semiarid environments at different latitudes Bonatelli et al. 2014;Cornejo-Romero et al. 2017;Quipildor et al. 2017).
On the other hand, Silva et al. (2020) suggested that the Early Quaternary climatic events were more important in the diversification and population divergence of plants from open vegetation, at least in South America. However, the region where C. nizandensis is distributed has presented changes in the composition of the vegetation during the Late Quaternary (Gámez et al. 2014), although the inferred demographic expansion for C. nizandensis appears to be independent of these climatic fluctuations ( fig. 3). An explanation for this may be that despite the altitudinal change in the elements of the tropical dry forest that surround the xerophytic enclaves, population expansion events of C. nizandensis were restricted to rocky outcrops. A similar case has been reported in Australia, where several species from the surrounding landscape show distribution changes, while the endemic elements of rocky outcrops have historically stable demographics (Byrne 2008;Broadhurst et al. 2017). The stability of rocky habitats, which reaches hundreds of thousands of years, has been used to explain the demographic stability and accumulation of genetic variability in other species of restricted distribution such as Stypandra glauca R.Br. (Tapper et al. 2017) and this is most likely what has happened in C. nizandensis.

Conservation considerations
Maintaining genetic diversity is essential for the conservation of species, as this enables the long-term viability of populations and survival in response to environmental changes, disease resistance, and chances of survival in general (Toro & Caballero 2005). Isolated populations are much more prone to extinction; however, a high dispersal capacity or colonization processes from different sources can reduce the probability of extinction (Frankham 1997). The restricted and fragmented distribution of C. nizandensis does not seem to be a limiting factor for its conservation, as population densities and genetic connectivity, together with its life history characteristics, have allowed it to sustain a high effective population size that mitigates the effects of genetic drift. Therefore, it maintains moderate to high levels of genetic diversity, and its population has increased in the last thousands of years, as shown in this work.
Despite the fact that populations of C. nizandensis naturally show a population expansion, their restricted distribution makes these cacti susceptible to anthropogenic effects due to habitat destruction. Even though the soil where they grow is not functional for agriculture, their populations have suffered from the burning of the nearby forest for the creation of paddocks. Other risks observed were rock extraction and pruned plants (Aldo Isaac Juárez-Miranda pers. obs.). Another danger for this microendemic species is the Transisthmic railway expansion, since there are populations described on both sides of the railway (Bravo-Hollis & Sánchez-Mejorada 1978).
Thus, one possible conservation strategy is to name the highest area where the rocky enclaves are located, a reserve of diversity. This area is not in Mexico's priority protected area, so it is important to include it (Arriaga et al. 2000). It is worth highlighting the tourist attraction and the unique landscape beauty of these "succulent rock gardens" and how this can help preserve these ecosystems similar to the N1 population that shows the highest degree of conservation and least damage, being located in proximity of the tourist trails signposted from the town of Nizanda. Additionally, the M2 population could represent a unique and important conservation unit due to its genetic difference, in terms of some haplotypes. However, given that this species has very few localities, it is recommended to preserve the largest number of populations to achieve effective and efficient conservation strategies (Pezoa 2001).

CONCLUSION
The species that inhabit rocky outcrops are of great interest either because of their restricted distribution or because of their adaptations to these arid environments (Crowther 1982;Porembski 2007). In some cases, rocky outcrops have been shown to be buffers of short-and long-term environmental changes that also act as filters for species unadapted to these areas (Porembski 2007). We can think the outcrop habitats have served as ecological refuges for certain species, as can be the case of C. nizandensis that has expanded its population in the past 150,000 years (Porembski & Barthlott 2000;Pérez-García & Meave 2005;Porembski 2007;Couper & Hoskin 2008). These sites still provide the opportunity to ask unexplored questions, from the point of view of population genetics (Torres-Ribeiro et al. 2007).
In Mexico, the investigation of rocky outcrops is still incipient, and in reality, there is no research on population genetics or phylogeography in species that inhabit outcrops. As far as we know, this is the first work to address this question. Therefore, it would be an important approach that can be used as a model for other species living in these enclaves, such as Agave nizandensis Cutak (Agavaceae) which is also associated with rocky substrates in the Nizanda region (Pérez-García 2002;Pérez-García & Meave 2005