1. Introduction

Loss of continuous natural areas and conversion of forests to smaller areas have drastically reduced the size of plant and animal populations, increasing the isolation between populations of the same species (Young and Clarke, 2000; Lowe et al., 2005). Plant reproduction may decrease in fragmented populations because plants may receive fewer flower visitors due to a decline in the richness and abundance of pollinators, modification in species composition and limitation in movement among landscapes (Young and Clarke, 2000; Cascante et al., 2002; Goverde et al., 2002; Frankham, 2010). Disturbance also changes the abundance, richness or behavior of animal seed dispersers affecting seed dispersal (Ghazoul, 2005; Lowe et al., 2005), recruitment and establishing of seedlings (Ewers and Didham, 2006), and the spatial genetic structure of plant populations (Hamrick, 2004). Shortly following a disturbance, levels of genetic homozygosity in affected plant populations are predicted to increase, and this may lead to the expression of recessive alleles deleterious to individual fitness (Hart and Clark, 1997; Reed and Frankham, 2003). In the long term, the loss of genetic variation may induce severe population bottlenecks due to genetic drift, high population divergence and low effective population size (Zhang et al., 2004; Moreira et al., 2009).

Dry ecosystems around the world are experiencing major environmental impacts due to climate change since water is the major limiting factor for plant productivity in such environments (Fischer and Turner, 1978; Anderson et al., 1994; Hamrick et al., 2002). Other factors, such as strong drying winds, extreme temperatures, limited nutrients, high light intensity, overexploitation of natural resources, and desertification have caused habitat degradation, fragmentation and loss of biodiversity (McNeely et al., 1990; Ayyad, 2003).

The Cactaceae family includes approximately 1600 species distributed in arid and semi-arid habitats (Gibson and Nobel, 1986; Dávila et al., 2002). Mexico is the most important center of diversification with almost 850 species and 54 genera (Bravo-Hollis, 1978; Arias, 1993; Anderson et al., 1994; Hernández and Godínez, 1994). The Balsas River Basin and the Tehuacán Valley in central Mexico have a great number of cacti with approximately 81 species, of which 25% are endemic, constituting one of the principal structural and floristic elements of the dry forests (Arias et al., 1997). The species of Mammillaria predominate in this area and have biological characteristics that make them more vulnerable to the effects of disturbance in their populations such as low growth rates and long life cycles (Gibson and Nobel, 1986). Particularly, Mammillaria pectinifera is an extremely narrowly distributed species with only 18 fragmented populations remaining in the Tehuacán-Cuicatlán Valley in the states of Puebla and Oaxaca in Central Mexico (Martorell and Peters, 2005; Valverde et al., 2009). This species grows on relatively deep alkaline calcareous soils and the plants are localized on small isolated patches in xerophyllous scrub and grasslands (Arias et al., 1997). The populations of M. pectinifera fae serious risk of extinction, mainly due to massive collection of plants for commercial purposes and the degradation of their habitats (Glass, 1998; Valverde and Zavala-Hurtado, 2006). Grazing also causes the mortality of plants of M. pectinifera in its natural habitats (Meyrán, 1973; Bravo-Hollis, 1978; Bravo-Hollis and Sánchez-Mejorada, 1991). Thus, M. pectinifera is included in the category of risk of extinction in international (IUCN; CITES, 2013) and national (NOM 059 ECO1 2001; SEMARNAT, 2010) regulations.

Most of the previous population genetic studies within the Cactaceae family have been mainly conducted on columnar and short globose cacti using allozymes (Nassar et al., 2001; Hamrick et al., 2002; Moraes et al., 2005; Parra et al., 2008), RAPDs (Clark-Tapia et al., 2005) and microsatellites (Figueredo et al., 2010; Castro-Félix et al., 2013; Solórzano-Lujano et al., 2014). Cornejo-Romero et al. (2014) found that populations of species of Mammillaria were genetically structured indicating that geographic isolation and limited dispersal were the primary causes of genetic population differentiation.

Ecological studies on Mammillaria have assessed the risk of extinction of its species due to the high specificity of the habitats where they occur and the degree of disturbance of the natural environments (Contreras and Valverde, 2002; Martorell and Peters, 2005; Valverde and Zavala-Hurtado, 2006; Valverde et al., 2009; Peters et al., 2014). Many aspects about genetic structure of the populations of M. pectinifera are unknown. Therefore, it is important to conduct studies to assess the genetic diversity of their populations and provide some guidelines for conservation. The main objective of this study was to examine the levels of genetic diversity of populations of M. pectinifera and to characterize its population genetic structure, inbreeding, recent bottlenecks and fluctuations in the effective population size to identify the factors that are determining the evolution and maintenance of this endangered microendemic species, and to propose conservation strategies to preserve this species.

2. Materials and methods 2.1. Study species

M. pectinifera F.A.C.Weber in Bois (Cactaceae) is a small cactus that in habitat appears above the ground as a squat stem completely covered by flattened spines, and usually solitary. Flowers are medium-sized in a ring around the sides of the globe, white to pink in color with darker mid-strip, 20e30 mm long and stigmas are green (Arias et al., 1997). This species has a blooming season from December to March. Reproduction starts at an age of eight years. Fruits are small, red at maturity, barely emerging above the spines. Seeds, usually black, are retained among the plant tubercles and are released gradually (Arias et al., 1997).

We collected from 18 populations of M. pectinifera in the Tehuacan-Cuicatlan Biosphere Reserve in central Mexico; from 7 to 13 individuals were sampled randomly with at least 100 m of distance between each individual. In all the samples, we did cuts about 1 cm² to get photosynthetic tissue, subsequently, the tissue was frozen at-70 ℃ until used for DNA extraction. For the DNA extraction (184 individuals), we followed the protocols of Doyle and Dickson (1987) and Doyle and Doyle (1987). Six nuclear DNA (nSSR) microsatellite loci were selected and amplified in multiplex polymerase chain reactions (PCR). Two groups of primers were arranged according to allele size and fluorescent labels. The first group was formed by the primer pairs for Mam VTC2, Mam VTC8, Mam VTC9, Mam VTC11, previously designed for the species M. crucigera by Solórzano-Lujano et al. (2009), whereas the second group included the primer pairs for P44 and P50, previously designed for Polaskia chichipe (Otero-Arnaiz, 2005). PCR was performed using the QIAGEN Multiplex PCR kit (QIAGEN) in a volume of 5 ml containing 1× Multiplex PCR Master Mix, 2 μM each primer, dH2O, and 20 ng template DNA. The thermal cycling conditions consisted of 40 cycles, each at 95 ℃ for 1 min, annealing for 1 min (i.e. first and second primer groups 58 ℃ and 60 ℃, respectively), extension for 2 min at 72 ℃ and a final extension at 72 ℃ for 10 min. Multiplex PCR products were combined with a GeneScan-500 LIZ size standard and sequenced by an ABI-PRISM 3100 Avant sequencer (Applied Biosystems). Fragments were analyzed and registered with the Peak Scanner program 1.0 (Applied Biosystems).

2.2. Genetic diversity

To test whether there was evidence of null alleles, upper allele dropout, or small genotyping errors due to stutter in the SSR data that we accumulated in 18 populations of M. pectinifera, we used the Micro-Checker v 2.2.3 program with 10² bootstrap simulations and a 95% confidence interval (Van Oosterhout et al., 2005). We calculated the values and standard error of mean number of alleles per locus (N_a), mean effective number of alleles (N_e), mean observed heterozygosity (H_O), mean expected heterozygosity (H_E) and mean of fixation index (F_IS) in 18 populations of M. pectinifera using the GENETIX 4 program (Belkhir et al., 2004).

2.3. Population structure and cluster analysis

Population genetic differentiation R_ST was estimated by the stepwise mutation model (SMM) performed with 10⁴ permutations in the ARLEQUIN 3.5.1.2 software (Excoffier and Lischer, 2010). A hierarchical test of population structure was estimated using the stepwise mutation model (SMM) performed with AMOVA in ARLEQUIN 3.5. (Excoffier et al., 2005). We compared the variance distribution between groups, among populations between groups and within populations. The statistical significance was tested using 10⁴ permutations utilizing the resulting two groups of genotypes obtained previously by STRUCTURE.

Acluster analysis using a Bayesian approach was conducted using STRUCTURE version 2.3.3 (Pritchard et al., 2000; Falush et al., 2003; Hubisz et al., 2009). In this analysis individuals are probabilistically assigned to one of the predefined K populations (gene pools) to identify the optimal number of genetic groups (Evanno et al., 2005). The optimum number of groups (K) was determined by varying the value of K from 1 to 10 and running the analysis ten times per K value, in order to determine the maximum value of posterior likelihood [LnP (D)]. Each run was performed using 10⁵ burn-in periods and 10⁶ Markov Chain Monte Carlo (MCMC) repetitions after burn-in. We used a model allowing for admixture with correlated allelic frequencies without any prior information. Also, we determined the most probable value of K using the maximum value of △K according to Evanno et al. (2005) implemented in the program Structure Harvester 0.6.1 (Earl and Von Holdt, 2011).

To identify probable geographic and genetic breaks among populations of M. pectinifera, we used the Monmonier's maximum difference algorithm with the software BARRIER version 2.2 (Manni et al., 2004) that creates a map of sampling locations from their geographical coordinates. Barriers are then characterized on the map by identifying the maximum values within the population pairwise genetic distance matrix. We employed an average square distances (ASD) matrix of (Goldstein et al., 1995; Slatkin, 1995) estimated for 18 populations of M. pectinifera. Resampling random subsets of individuals within populations provided 100 bootstrap replicate distances were constructed with the MSA program (Dieringer and Schlotterer, 2003) to accomplish statistical signi € ficance for the predicted barriers.

Population connectivity among M. pectinifera was scrutinized with the POPGRAPH library (http://dyerlab.bio.vcu.edu/software.html), this approaches generate a network of population linkages and defines the degree of genetic variation within populations (Dyer and Nason, 2004). In the POPGRAPH context, the set of nodes represents sampled populations, and the edges represent the multivariate measures of genetic covariance among populations. The difference in node size reflects differences in withinpopulation genetic variability, while the edge length represents the among-population component of genetic variation due to the connecting nodes (Dyer and Nason, 2004). POPGRAPH, recognize population pairs where long distance migration may have occurred by indicating population pairs with significantly greater inter-site distances (Dyer et al., 2010). It also identifies pairs of populations that are located significantly closer than predicted by inter-site separation, which suggests that a direct barrier might exist between the linear distances.

2.4. Population bottleneck and effective population size

To detect evidence of recent bottlenecks in M. pectinifera populations, we used the BOTTLENECK 1.2 software (Piry et al., 1999) to examine the two genetic groups previously obtained by STRUCTURE. Due to small sample sizes, our analysis was limited to. Also, England et al. (2006) showed that the disequilibrium method leads to severe downward bias in the estimate of Ne when sample size (S) is small relative to effective size. Therefore, the estimation by population will result in bias due to the small sample size of the population. Recent bottlenecks could be defined as a population where the rare alleles are the first to be lost, diminishing the mean number of alleles per locus. In contrast, heterozygosity is less affected, generating a transient excess in heterozygosity compared to that expected given the resulting number of alleles (Cornuet and Luikart, 1996; Luikart and Cornuet, 1998). To analyze the data set, we used 90% stepwise and 10% multistep mutations 10⁴ iterations with the Wilcoxon signed-rank test, and the stepwise mutation (SMM), the infinite allele (IAM) and two-phase mutation (TPM) models. Moreover, we also assessed the effective population size in populations of M. pectinifera with the program LDNe (Waples and Do, 2008), which implements the bias-correction method developed by Waples (2006) to obtain Ne from each sample of S individuals. For LDNe, we used the criterion Pcrit=0.02 (alleles with frequency ^˂ 0.02 are excluded), which generally provides a good balance among accuracy and bias (Waples and Do, 2008). Confidence intervals (CIs) for Ne were based in the chi-square approximation implemented by LDNe (Waples, 2006).

3. Results 3.1. Genetic diversity

In general, we did not find evidence for null alleles over all sample-loci combinations, and tests for error due to stutter and upper allele dropout were negatives in all cases. The values of the genetic diversity parameters estimated in 18 populations of M. pectinifera are reported in Table 1. We observed a wide range of values from 4.50 to 7.00 in mean number of effective alleles per locus (N_a) and 0.716 to 0.950 in mean observed heterozygosity (H_O) values (Table 1). Wright's inbreeding coefficients within populations (F_IS) were positive in almost all cases, indicating a deficit of heterozygotes ranging from F_IS=0.071 to 0.621 and some populations with negative values ranging F_IS=-0.010 to 0.446 (Table 1).

3.2. Population structure and cluster analysis

Pairwise genetic differentiation among populations indicated low to moderate genetic differentiation based on R_ST (Table 2). Hierarchical analysis of molecular variance (AMOVA), performed for the SMM mutation models (i.e. R_ST), indicated that most of the genetic variation resided within populations (Φ_ST=80.63%, P=0.001) followed by variation among populations within groups (Φ_SC=10.32%, P=0.001), while the differentiation among groups only accounted for the remaining variation (Φ_CT=9.04%, P=0.001) (Table 2).

The highest posterior probability was obtained from Bayesian likelihood [LnP (D)] and the △Kapproach (Evanno et al., 2005) implemented in the program Structure Harvester (Earl and Von Holdt, 2011). This statistic determined that K=2 is the optimum value for the number of genetic clusters included in the analysis (Fig. 1). Both genetic groups showed a very widespread distribution and also some of them exhibited a great amount of exchange among groups. For instance, the first genetic group (i.e. Cluster 1, green) of populations were consistently structured and very widespread over all the distribution of M. pectinifera. In addition, some of the populations had a greater ancestry coefficient Q from 0.669 to 0.847 such as the populations CO, CH, TE, SM, RE, SI, ZA and LA (Fig. 2). Populations from the second genetic group (i.e. Cluster 2, red), such as FR, SU, EC, CC, AC, ER and CQ had a widespread distribution and were very well structured with ancestry coefficient (Q) going from 0.710 to 0.904. The rest of populations TC SB and NO showed a greater amount of exchange among genetic groups with an ancestry ranged from Q 0.410e0.590 (Fig. 2). Table 3.

In total, four barriers were identified among the 18 populations of M. pectinifera using the Barriers analysis with 100 bootstrap replicates of ASD genetic distance matrixes (Fig. 2). The most significant barrier, with a bootstrap support of 98%, separated TC from SM in the northern part of the distribution. The second barrier, with 90% bootstrap value, indicated that there is a complex barrier that divides SM from the rest of the populations in the center and the south-western region. The third barrier, with 87% bootstrap support, separated populations in the center of the distribution, including, CH, EC, LA and SU. The fourth barrier, with a 79% bootstrap support, separated the populations from RE and NO located in the center-west of the distribution of M. pectinifera (Fig. 2).

The POPGRAPH network of populations based on nSSR genotypes had 17 edges out of the possible 22, which indicates extensive historical connectivity among most populations (Fig. 3). The most common colonization patterns occurred in a northesouth axis along the Tehuácan-Cuicatlán Biosphere Reserve, throughout southern portions of Puebla to Zapotitlán, but it was not uncommon to see eastewest gene exchanges, especially across the northeeast to the southewest part of the species range. Our data indicate occasional long-distance gene exchange across the southern portion of Tehuacán, FR, to the north toward TC and CH, where the genotypes are significantly more similar than expected based on spatial distance. Our analysis also identified a great amount of network connection between the center of the distribution towards the northwest, south, but not much connectivity with the southwest portion, which were more genetically dissimilar than predicted by spatial distance (Fig. 3).

3.3. Population bottleneck and effective population size

The results of the analysis to test the evidence of recent bottlenecks (excess of heterozygosity) using the infinite allele model (IAM), two-phase model (TPM) and stepwise mutation models (SMM) are presented in Table 4. Results were significant (P < 0.050) for the IAM P=0.0003 and the SMM P=0.006 models observed in the loci VTCMAM11, P50 and P44. Our analysis revealed that some populations of M. pectinifera departed from neutrality, exhibiting excess heterozygosity. This suggests that the populations have recently undergone a bottleneck. The results of estimation of effective population size (N_e), performed with the program LDNe, in populations of M. pectinifera showed that both genetic groups had moderated values in the first group (N_e=455 individuals), followed by the second group (N_e=89.8); in all cases, estimates had high Jackknife support and a good confidence interval (CIs) (Table 4).

4. Discussion

Anthropogenic activities, such as massive extraction of plants and disturbance of natural habitats, may lead to a decrease in population size and changes in the genetic structure of species populations (Lowe et al., 2005). M. pectinifera has been subject to these pressures, particularly in reproductive stages, which has led to a significant reduction in their populations. Unexpectedly, the populations of M. pectinifera have been shown to have similar (López-Ortiz et al., 2013; Macías-Arrastrio, 2013) or higher (Ibarra-Suárez, 2009; Tapia-Salcido, 2011) values of genetic diversity compared to other short-globose cacti species of the same genus. However, we found significant evidence of recent bottlenecks in some populations of M. pectinifera. Lowe et al. (2005) suggested that after an episode of habitat disruption, populations that remained small for many generations suffer loss of allelic diversity and random genetic drift (Hamrick, 2004; Ellstrand and Elam, 1993). The destruction and removal of reproductive adult trees might result in fecundity variance, reduced effective population size and limited gene flow among the remnant populations (Ewers and Didham, 2006). Furthermore, M. pectinifera has life history traits that make it more vulnerable to the effects of disturbance over the course of generations such as low growth rates, long life cycles and serotiny (Gibson and Nobel, 1986; Bravo-Hollis and Sánchez-Mejorada, 1991; Valverde et al., 2009).

High inbreeding values, moderate to high heterozygosity and low to moderate effective population size are positive indicators that predict short-term population viability in M. pectinifera. While only a relatively small number of individuals may be needed in the short term to prevent inbreeding, greater genetic variation is thought to increase the likelihood of survival over much longer time scales (Frankham, 2010; Jamieson and Allendorf, 2012). Adaptation or resilience to novel or changing habitats is the most important feature (Williams and Hoffman, 2009). Thus, management strategies should be implemented that maintain, as far as possible, the full complement of genetic diversity within the Tehuacán-Cuicatlán Biosphere Reserve.

Hierarchical analysis of molecular variance revealed that most of the genetic variation is within populations, and there is low genetic differentiation among populations. This observation was confirmed by the POPGRAPH analysis that also indicated extensive historical connectivity. The connectivity values obtained for M. pectinifera were relatively high and it is likely the result of historical events of gene flow through pollen and seed dispersion by biotic agents. However, extant populations of M. pectinifera are distributed in isolated patches with few opportunities for seed dispersal across valleys and foothills (Valverde and ZavalaHurtado, 2006; Valverde et al., 2009; and Peters et al., 2014). Furthermore, populations of Mammillaria in the Tehuacán Valley present a restricted spatial distribution, high habitat specificity in restricted altitudinal ranges, soil and vegetation type (Hernandez and Godinez, 1994; Hershkovitz and Zimmer, 1997). In particular, M. pectinifera plants occurred only at altitudinal gradients from 1778 to 2100 m in mild to moderate slopes, on relatively deep alkaline calcareous soils, with high water retention capacity on stony hills (Valverde and Zavala-Hurtado, 2006; Valverde et al., 2009; Peters et al., 2014).

The massive extraction of plants, deforestation, changes in land use, and over-exploitation of areas by human activities have damaged the natural habitats and the populations of M. pectinifera (Zavala-Hurtado and Valverde, 2003; Martorell and Peters, 2005). Despite this, we observed that populations of M. pectinifera still show moderate to high levels of genetic diversity due to their strong capacity to adapt and survive in localities with high water deficits (Hernandez and Godinez, 1994; Hershkovitz and Zimmer, 1997). However, if not properly preserved, this species will become extinct in the future. For instance, according to Peters et al. (2014), 45% of all known populations of M. pectinifera are located outside of natural protected areas. For that reason, we suggest in situ conservation to prevent the decrease in population sizes and loss of genetic diversity. Natural reserves such as the Tehuacán-Cuicatlán Biosphere Reserve must play an essential role in protecting and restoring the habitat of M. pectinifera (Contreras and Valverde, 2002; Peters et al., 2014). Also, a long-term ex situ conservation program which includes collecting seed from different sources and optimizing seed germination and plant establishment protocols is needed to restore disturbed habitats. Also, a supply of living plants for trade would prevent further extraction of plants from nature.

Acknowledgements

We thank to D. Lugo-Aquino, N. Pérez-Nasser and Ana Luisa Albarrán-Lara for technical assistance and suggested improvements to the manuscript. This project was supported by SDI-UNAM (Secretaria de Desarrollo Institucional, Universidad Nacional Autónoma de México, Octavo piso Torre de Rectoría, Ciudad Universitaria).