Landscape genetics reveals inbreeding and genetic bottlenecks in the extremely rare short-globose cacti Mammillaria pectinifera (Cactaceae) as a result of habitat fragmentation.

Mammillaria pectinifera is an endemic, short-globose cactus species, included in the IUCN list as a threatened species with only 18 remaining populations in the Tehuacán-Cuicatlán Valley in central Mexico. We evaluated the population genetic diversity and structure, connectivity, recent bottlenecks and population size, using nuclear microsatellites. M. pectinifera showed high genetic diversity but some evidence of heterozygote deficiency (FIS ), recent bottlenecks in some populations and reductions in population size. Also, we found low population genetic differentiation and high values of connectivity for M. pectinifera, as the result of historical events of gene flow through pollen and seed dispersal. M. pectinifera occurs in sites with some degree of disturbance leading to the isolation of its populations and decreasing the levels of gene flow among them. Excessive deforestation also changes the original vegetation damaging the natural habitats. This species will become extinct if it is not properly preserved. Furthermore, this species has some ecological features that make them more vulnerable to disturbance such as a very low growth rates and long life cycles. We suggest in situ conservation to prevent the decrease of population sizes and loss of genetic diversity in the natural protected areas such as the Tehuacán-Cuicatlán Biosphere Reserve. In addition, a long-term ex situ conservation program is need to construct seed banks, and optimize seed germination and plant establishment protocols that restore disturbed habitats. Furthermore, creating a supply of living plants for trade is critical to avoid further extraction of plants from nature.

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, Hernandez and Godinez, 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-Felix et al., 2013, Solórzano 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.

Materials and methods

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, 20–30 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 cm2 to get photosynthetic tissue, subsequently, the tissue was frozen at −70 °C 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 μl 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 °C for 1 min, annealing for 1 min (i.e. first and second primer groups 58 °C and 60 °C, respectively), extension for 2 min at 72 °C and a final extension at 72 °C 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).

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 102 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), mean effective number of alleles (N), mean observed heterozygosity (H), mean expected heterozygosity (H) and mean of fixation index (F) in 18 populations of M. pectinifera using the GENETIX 4 program (Belkhir et al., 2004).

Population structure and cluster analysis

Population genetic differentiation R was estimated by the stepwise mutation model (SMM) performed with 104 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). We compared the variance distribution between groups, among populations between groups and within populations. The statistical significance was tested using 104 permutations utilizing the resulting two groups of genotypes obtained previously by STRUCTURE.

A cluster 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 105 burn-in periods and 106 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 Schlötterer, 2003) to accomplish statistical significance 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 within-population 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). 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.

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 104 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).

Results

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 (Na) and 0.716 to 0.950 in mean observed heterozygosity (H) values (Table 1). Wright's inbreeding coefficients within populations (F) were positive in almost all cases, indicating a deficit of heterozygotes ranging from F = 0.071 to 0.621 and some populations with negative values ranging F = −0.010 to 0.446 (Table 1).

Table 1

Locality name, sample size, N/W geographical coordinates (degrees), mean number of effective alleles per locus (N), mean observed heterozygosity (H), mean expected heterozygosity (H), for 18 populations of M. pectinifera in the Tehuacán- Cuicatlán Biosphere Reserve, Mexico. Standard errors are included in parenthesis.

Locality	Genetic diversity
	Sample size	Coordinates	Ne	HO	HE	FIS
1. CO	11	18.413/−97.420	6.50 (2.07)	0.893 (0.10)	0.808 (0.08)	0.345 (0.05)
2. CH	7	18.619/−97.548	5.83 (2.78)	0.734 (0.24)	0.761 (0.17)	0.327 (0.07)
3. TE	10	18.403/−97.548	7.00 (2.53)	0.826 (0.11)	0.792 (0.14)	0.237 (0.01)
4. SM	13	18.703/−97.603	7.00 (3.03)	0.845 (0.12)	0.785 (0.10)	−0.010 (0.05)
5. RE	11	18.555/−97.635	6.50 (2.08)	0.932 (0.08)	0.808 (0.07)	0.017 (0.07)
6. TC	11	18.865/−97.702	5.83 (2.31)	0.842 (0.16)	0.758 (0.10)	0.128 (0.12)
7. SB	11	18.464/−97.568	6.50 (2.66)	0.950 (0.08)	0.795 (0.08)	−0.153 (0.06)
8. SI	7	18.409/−97.429	4.83 (0.98)	0.849 (0.16)	0.768 (0.10)	−0.162 (0.09)
9. ZA	10	18.358/−97.480	6.50 (2.07)	0.866 (0.10)	0.793 (0.07)	0.326 (0.11)
10. FR	11	18.476/−97.572	5.66 (1.21)	0.842 (0.10)	0.777 (0.03)	0.424 (0.01)
11. FT	8	18.302/−97.695	5.33 (2.06)	0.854 (0.16)	0.756 (0.12)	0.071 (0.09)
12. SU	10	18.619/−97.477	6.83 (2.40)	0.796 (0.12)	0.789 (0.07)	0.621 (0.08)
13. EC	11	18.605/−97.542	5.83 (1.94)	0.772 (0.21)	0.751 (0.10)	0.225 (0.04)
14. CC	11	18.636/−97.422	4.50 (1.37)	0.771 (0.14)	0.734 (0.03)	0.261 (0.14)
15. LA	11	18.587/−97.548	6.16 (2.31)	0.886 (0.10)	0.772 (0.07)	0.234 (0.06)
16. AZ	10	18.704/−97.414	5.00 (1.26)	0.948 (0.05)	0.753 (0.01)	−0.446 (0.19)
17. ER	10	18.444/−97.463	5.16 (1.94)	0.716 (0.21)	0.733 (0.12)	−0.116 (0.14)
18. CQ	11	18.638/−97.415	5.00 (1.67)	0.903 (0.14)	0.766 (0.07)	−0.254 (0.05)

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

Table 2

Pairwise population genetic differentiations among localities of M. pectinifera in the Tehuacán-Cuicatlán Biosphere Reserve using the matrix of Slatkin linearized R (SMM), over 104 replicates in the Arlequin 3.5 software (Schneider et al., 2000).

RST	CO	CH	TE	SM	RE	FR	TC	SB	SI	ZA	FT	SU	EC	CC	LA	AZ	ER	CQ
CO	–
CH	0.122	–
TE	0.008	0.171	–
SM	0.094	0.157	0.097	–
RE	0.223	0.172	0.050	0.011	–
FR	0.005	0.100	0.002	0.055	0.050	–
TC	0.047	0.305	0.044	0.094	0.281	0.047	–
SB	0.029	0.218	0.067	0.118	0.295	0.062	0.018	–
SI	0.058	0.172	0.053	0.006	0.086	0.044	0.111	0.186	–
ZA	0.190	0.256	0.047	0.071	0.024	0.024	0.197	0.197	0.104	–
FT	0.337	0.423	0.322	0.241	0.395	0.231	0.296	0.147	0.381	0.301	–
SU	0.153	0.317	0.006	0.017	0.144	0.068	0.047	0.067	0.088	0.001	0.325	–
EC	0.115	0.352	0.061	0.143	0.312	0.057	0.011	0.034	0.181	0.184	0.181	0.102	–
CC	0.161	0.366	0.007	0.025	0.219	0.106	0.041	0.087	0.075	0.075	0.212	0.068	0.011	–
LA	0.028	0.273	0.130	0.207	0.367	0.012	0.067	0.029	0.223	0.317	0.390	0.196	0.078	0.246	–
AZ	0.208	0.411	0.123	0.117	0.302	0.158	0.106	0.137	0.197	0.227	0.360	0.062	0.046	0.030	0.214	–
ER	0.276	0.445	0.199	0.187	0.362	0.184	0.149	0.051	0.294	0.215	0.093	0.179	0.028	0.021	0.278	0.146	–
CQ	0.225	0.435	0.145	0.191	0.372	0.118	0.109	0.078	0.271	0.249	0.171	0.100	0.025	0.078	0.181	0.010	0.024	–

The highest posterior probability was obtained from Bayesian likelihood [LnP (D)] and the ΔK approach (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.410–0.590 (Fig. 2). Table 3.

Fig. 1

Mean and standard deviation values of ΔK for ten independent runs of STRUCTURE plotted against the number of genetic groups K. The peak indicates the most probable number of genetic groups given the data using Structure Harvester (Earl and Von Holdt, 2011).

Fig. 2

Map shows sampling localities representing 18 populations of M. pectinifera in the Tehuacán-Cuicatlán Biosphere Reserve. Green genotype and red genotype represent the genetic ancestry groups corresponding to M. pectinifera populations obtained by the program STRUCTURE. Black bars show the barriers detected with the BARRRIERS program.

Table 3

Analysis of molecular variance (AMOVA) for populations of M. pectinifera. Analysis was performed on the nSSR data and R was used for the two genetic clusters obtained by means of STRUCTURE. Asterisks indicate statistically significant values (P < 0.01). Tests were based on 104 random permutations.

	Source of variation	Sum of squares	Variance components	Percentage of variation	Fixation index
RST
	Among groups	10231.3	47.211	9.04	ΦCT = 0.090***
	Among populations within groups	24313.4	53.881	10.32	ΦSC = 0.010***
	Within populations	147309.7	420.884	80.63	ΦST = 0.080***
	Total	181854.5	521.978

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 north–south axis along the Tehuacán-Cuicatlán Biosphere Reserve, throughout southern portions of Puebla to Zapotitlán, but it was not uncommon to see east–west gene exchanges, especially across the north–east to the south–west 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).

Fig. 3

Patterns of genetic connectivity among 18 populations of M. pectinifera in the Tehuacán-Cuicatlán Biosphere Reserve, Mexico performed with the program POPGRAPH.

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), performed with the program LDNe, in populations of M. pectinifera showed that both genetic groups had moderated values in the first group (N = 455 individuals), followed by the second group (N = 89.8); in all cases, estimates had high Jackknife support and a good confidence interval (CIs) (Table 4).

Table 4

Bottleneck analysis for M. pectinifera populations in the Tehuacán-Cuicatlán Biosphere Reserve using Wilcoxon rank test under infinite allele, stepwise mutation and two phase model. Parameters for TPM: variance = 10%, proportion of SMM = 90%, estimation based on 104 replications. P, probability. IAM, infinite allele model; TPM, two phase model; SMM, stepwise mutation model. Also, we include the results obtained for the estimation of the population effective size for the two genetic groups, values obtained with the program LDNe.

Models	Group1	Group2
IAM	0.007***	0.039**
TPM	0.781	0.656
SMM	0.046*	0.042*
LDNe	455.2**	89.8**

*, ** and *** Indicate significant deviation from equilibrium as value less than 0.05.

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-Arrastio, 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 Ellam, 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 Zavala-Hurtado, 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.