Material and methods

Study species. Guaiacum sanctum (Zygophyllaceae) is a slow-growing tree distributed in tropical and subtropical forests from northern Costa Rica to southeast Mexico, the Florida Keys, and the Caribbean Islands (^{Holdridge et al., 1975}). In Mexico, G. sanctum is distributed in the states of Oaxaca, Chiapas, Quintana Roo, Yucatán and Campeche. It is an evergreen tree that grows up to 25-30 m in height and 60 cm diameter at breast height (dbh). Its compound leaves are paripinnate and opposite with 4-12 oblong to obovate leaflets. Perfect flowers are solitary or in terminal panicles with blue-violet petals 0.7-1.2 cm long and yellow stamens. Flowering occurs between March and July. Wasps and bees are the reported pollinators. Fruits are capsules 1.5-2.0 cm long, winged, that turn yellow-orange when ripe. Seeds are enclosed by an intensely red aril, which attract birds, the main dispersers (^{Chavarría et al., 2001}; ^{Grow and Schwartzman, 2001}).

Sampling. Within the Yucatán Peninsula 17 populations of Guaiacum sanctum were sampled as follows: six in Yucatán, three in Quintana Roo, and eight in Campeche (Figure 1). On average, from 7-10 individuals were sampled. In most of the cases, we collected the remaining reproductive individuals. In denser populations such as some populations of Campeche State, we randomly chosen ten reproductive trees with at least 10 m of separation among trees. Leaf samples were frozen until genetic analysis. Genomic DNA from each individual was extracted from 100 mg of frozen leaf following the protocol proposed by ^{Lefort and Douglas (1999)}. Seven 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 primers pairs for Gcoult_1, Gcoult_2 and Gcoult_5, whereas the second group included the primer pairs for Gcoult_10, Gcoult_12, Gcoult_14 and Gcoult_15 (^{McCauley et al., 2008}). PCR was performed using the QIAGEN Multiplex PCR kit (QIAGEN) in a volume of 5 µl containing 1X 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 ran in an ABI-PRISM 3100 Avant sequencer (Applied Biosystems). Fragments were analyzed and registered with the Peak Scanner program 1.0 (Applied Biosystems).

Figure 1 Potential geographic distribution of Guaiacum sanctum in Mexico. Red dots represent the localities obtained in Herbarium specimens and GBIF web site. 

Genetic diversity. To test for the presence of null alleles, upper alleles dropout or errors due to stutter in the microsatellite data obtained from all the individuals of 17 populations of Guaiacum sanctum, the MICRO-CHECKER 2.2.3 software was used (^{Van Oosterhout et al., 2004}) with 10² bootstrap simulations and a 95 % confidence interval. We also tested departures from Hardy-Weinberg equilibrium (e.g. heterozygosity excess and deficit F _IS ) with the GENEPOP 4.1 software (^{Raymond and Rousset, 1995}) using the Markov-chain approach (e.g. 10³ dememorization steps, 10² batches and 10³ iterations per batch). Additionally, we estimated the effective number of alleles (N _e ), observed heterozygosity (H _O ), and expected heterozygosity (H _E ) for each population were performed with 10³ iterations using GENETIX 4.05 program (^{Belkhir et al., 1996-2004}).

Genetic structure. A Bayesian cluster analysis was conducted using STRUCTURE version 2.3.3 (^{Pritchard et al., 2000}; ^{Falush et al., 2003}; ^{Hubisz et al., 2009}). In this approach, individuals were assigned probabilistically 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 to determine the maximum value of posterior likelihood [LnP (D)]. Each run was performed using 50⁴ burn-in periods and 10⁶ Markov Chain Monte Carlo (MCMC) repetitions after burn-in. We used a model allowing 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 vonHoldt, 2012}). A hierarchical test of population structure was estimated using both mutation models [i.e. stepwise mutation model (SMM), the infinite allele model (IAM)] performed with AMOVA in ARLEQUIN 3.5. (^{Excoffier et al., 2005}). The variance distribution between groups, among populations between groups, and within populations was compared. The statistical significance was tested using 10⁴ permutations utilizing the resulting groups of genotypes obtained by STRUCTURE.

To identify possible geographic and genetic discontinuities among populations of Guaiacum sanctum, we used the Monmonier’s maximum difference algorithm with BARRIER version 2.2 (^{Manni et al., 2004}). This program creates a map of the sampling locations from geographical coordinates, where barriers are represented on the map by identifying the maximum values within the population pairwise genetic distance. We used a matrix of average square distance (ASD) (^{Goldstein et al., 1995}; ^{Slatkin, 1995}), estimated for the 17 populations of G. sanctum. Resampling random subsets of individuals within populations provided 100 bootstrap replicate distances that were constructed utilizing MSA program (^{Dieringer and Schlötterer, 2003}) to achieve statistical significance for the predicted barriers.

The connectivity among Guaiacum sanctum populations was analyzed using the POPGRAPH software (http://dyerlab.bio.vcu.edu/software.html). This analysis creates a network of population linkages and describes the amount of genetic variation within populations (^{Dyer and Nason, 2004}). In the POPGRAPH framework, 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. Thus, this analysis examines specifically the connections among individual populations that maintain genetic connectivity among all populations (^{Dyer and Nason, 2004}). POPGRAPH identifies pairs of populations 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.

Population size change and bottleneck detection. The software BOTTLENECK 1.2 (^{Piry et al., 1999}) was used to detect recent population bottlenecks that could be defined as a population where the rare alleles are the first to be lost decreasing the mean number of alleles per locus. In contrast, heterozygosity is less affected, producing a transient excess in heterozygosity relative to that expected given the resulting number of alleles (^{Cornuet and Luikart, 1996}; ^{Luikart and Cornuet, 1998}). To test 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. Additionally, we also estimated the effective population size in populations of G. sanctum with the program LDNe (^{Waples and Do, 2008}), which implements the bias-correction method developed by ^{Waples (2006)} to obtain N _e from each sample of 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, 2009). Confidence intervals (CIs) for Ne were based in the chi-square approximation implemented by LDNe (^{Waples, 2006}).

Results

Genetic diversity. We did not find evidence for null alleles overall sample-loci combinations, and tests for error due to stutter and upper allele dropout resulted negatives in all cases. Populations of Guaiacum sanctum in the Yucatán Peninsula had, on average, 26.4 alleles for the seven loci. Average number of alleles per locus indicated that the loci Gcoult_1, 5, 10, 12 were the most polymorphic (e.g. 7.94, 7.88, 6.17 and 7.88, respectively), then followed moderate to low polymorphism shown by the loci Gcoult_2, 14, 15 (e.g. 3.29, 5.94 and 5.29, respectively). The average number of alleles per locus per population and the observed and expected heterozygosity (H _O and H _E ) were higher in populations from Yucatán (N _e = 7.56, H _O = 0.790, H _E = 0.80) followed by populations in Quintana Roo (e.g. N _e = 6.37, H _O = 0.600, H _E = 0.750), and Campeche (e.g. N _e = 5.42, H _O = 0.630, H _E = 0.740; Table 1). Significant deviation from Hardy-Weinberg equilibrium due to heterozygosity excess was observed in the loci Gcoult-2 and 15, whereas heterozygosity deficiency was more frequent in the locus Gcoult-1, 2, 5, 6, 12, 14 and 15 (Table 1). The highest positive values of the coefficient of endogamy (F _IS ) (i.e. heterozygote deficit) were observed in populations from Quintana Roo (F _IS = 0.182; P = 0.001), whereas in populations from Campeche (F _IS = 0.035; P = 0.001) and Yucatán (F _IS = 0.003; P = 0.001) significant negative values for F _IS were observed (i.e. heterozygote excess; Table 1).

Table 1 Nuclear microsatellite (i.e. nSSR) parameters of genetic diversity and the exact test of Hardy-Weinberg equilibrium estimated in 17 populations of G. sanctum in the Yucatán Peninsula Mexico. 

Genetic structure. The highest posterior probability obtained from Bayesian likelihood [LnP (D)] and ∆K with the ^{Evanno et al. (2005)} approach, revealed that K = 2 is the effective number of genetic clusters. Cluster 1 (in red) was more common in the populations from Yucatán and Quintana Roo whereas the Cluster 2 (in green) was well represented in populations from Campeche (Figure 2). Hierarchical analysis of molecular variance (AMOVA), performed for both mutation models (i.e. F _ST and R _ST ), indicated that most of the genetic variation resided within populations (Φ _ST = 88.2 %, P = 0.001; Φ _ST = 89.5 %, P = 0.001) followed by variation among populations within groups (Φ _SC = 8.7 %, P = 0.001; Φ _SC = 6.7 %, P = 0.001), while the differentiation among groups only accounted for the remaining variation (Φ _CT = 2.9 %, P = 0.001 and Φ _CT = 3.6 %, P = 0.001; Table 2).

Table 2 Analysis of molecular variance (AMOVA) performed on the nSSR data and using FST and RST for the two group genetic clusters obtained by means of STRUCTURE for populations of Guaiacum sanctum. Asterisks indicate statistically significant values (P < 0.01). Tests were based on 104 random permutations. 

Seven barriers with more than 50 % bootstrap support were detected (Figure 2). Single barriers were observed, such as Barrier 1 that divided Quintana Roo 3 from the rest of the populations, the Barriers 2 and 3 that divided Yucatán 1 from Yucatán 2, and this one from other populations of Yucatán, respectively, whereas Barrier 6 splitted Campeche 7 from Quintana Roo1. Barriers 4 and 7 divided red genotypes in northern populations (i.e. Yucatán and Quintana Roo) from green genotypes observed in southern populations from Campeche. Finally, the Barrier 5 is a complex barrier that divided population from Campeche 4 and Campeche 5 from the rest of the populations (Figure 2).

Figure 2 Each pie chart represents the proportions in each population of the two genetic groups as assigned by the program STRUCTURE. Green Genotype and red Genotype are representing the genetic ancestry groups corresponding to G. sanctum populations. Genetic discontinuities (bold lines B-1-7) obtained with Monmonier’s maximum difference algorithm on genetic distances derived from microsatellite allele frequencies. 

The POPGRAPH network of populations based on nSSR genotypes had 17 edges out of the 22 possible ones and indicates extensive historical connectivity among most populations (Figure 3). The most common colonization patterns occurs in a north- south axis along Yucatán and Quintana Roo throughout southern portions of Campeche to the Calakmul Reserve, but it is 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 central portion of Campeche to the north of Yucatán where the genotypes are significantly more similar than expected based on spatial distance. Our analysis also identifies a great amount of network connection between the Campeche and Quintana Roo, but not much connectivity with Yucatán that resulted more genetically dissimilar than predicted by spatial distance (Figure 3).

Figure 3 Map of inter-population genetic structure comput ed from seven microsatellite loci in G. sanctum, implemented in the PopGraph R package. Populations are represented by nodes connected by the minimal set of edges that sufficiently describes the total among-population covariance structure. Circles represent the differences in node size that reflects the within and between populations differences in their heterozygosities. Numbering of populations is as in Table 1. 

Bottleneck detection and population size change. The results of detection of a recent bottleneck in populations of Guaiacum sanctum in Yucatán and Quintana Roo were non-significant (P < 0.50) for any of the models (IAM, TPM and SMM) analyzed. Populations of Campeche exhibited evidence of bottleneck only with the IAM (P = 0.039) (i.e. particularly in the loci Gcoult_2 and 5), but not with the models TPM and SMM (Table 3).

The results of estimation of effective population size (N _e ) performed with the program LDNe in populations of Guaiacum sanctum showed that Yucatán populations had the highest values (N _e = 262 individuals), followed by populations of Campeche (N _e = 123) and Quintana Roo (N _e = 100); in all cases, estimates had high Jackknife support and a good confidence interval (CIs) (Table 3).

Table 3 Bottleneck analysis for G. sanctum populations in the Yucatán Peninsula 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. ** Indicate significant deviation from equilibrium as value less than 0.05. Also, we include the results obtained for the estimation of the population effective size for Yucatan, Campeche and Quintana Roo populations. Values obtained with the program LDNe. 

Discussion

Current situation of Guaiacum sanctum. ^{López-Toledo et al. (2011)} modeled the ecological niche of G. sanctum to know the historical changes in the extent of available habitat and to project future trends on the basis of contemporary rates of habitat loss. They found strong evidence of habitat loss in the states of Oaxaca, Chiapas, northern Yucatán and Quintana Roo. These results assessed the negative effects of dramatic changes in land use across southern and eastern Mexico during the last two decades (^{Miles et al., 2006}). In Mexico, G. sanctum has suffered a considerably habitat loss of about 28.2 % in the last decades and a great reduction in population size (^{López-Toledo et al., 2008}, ²⁰¹¹). Possible causes can be attributed to a rapid increase in deforestation rates, selective removal of trees in forest fragments, changes of forest use to agriculture or pasture for cattle raising, and other human activities (^{Oldfield, 2004}; ^{López-Toledo et al., 2008}).

Genetic diversity. Previous studies on a related endemic species, Guaiacum unijugum from the Cape region in Baja California, Mexico, showed low values of genetic diversity (N _e = 1.076, H _E = 0.070, H _O = 0.073) obtained with nuclear microsatellite loci (^{McCauley et al., 2010}). Populations of G. unijugum are also endangered due to human activities such as tourism and excessive logging. In contrast, in our study and using the same nuclear microsatellite loci, we found higher levels of genetic diversity in all populations of G. sanctum (Table 1) than G. unijugum, despite the fact that some populations at Yucatán Peninsula are strongly affected by forest fragmentation. Current levels of genetic diversity indicate that enough number of individuals of G. sanctum escape logging retaining moderate to high levels of genetic diversity, suggesting that only rare and low frequency alleles were lost. Populations of G. sanctum still had different private alleles in each region (Table 1). Therefore, we suggested that the effect of a recent habitat fragmentation processes, like in our case, may have left a moderate signature in the genetic diversity of this species as were also previously reported in other tree species (^{Aguilar et al., 2008}; ^{Moreira et al., 2009}; ^{Figueroa-Esquivel et al., 2010}).

Additionally, we observed the highest positive values of endogamy F _IS in populations from Quintana Roo (F _IS = 0.182, P = 0.001), whereas in populations from Campeche and Yucatán both significant positive values (i.e. heterozygote deficit) and negative (i.e. heterozygote excess) for F _IS index were observed (Table 1). The deficit of heterozygotes in populations out of Hardy-Weinberg equilibrium indicates the existence of biotic and abiotic factors causing change during the periods of pre-and post-disturbance populations (^{Andre et al., 2008}; ^{Schaberg et al., 2008}; ^{Jacquemyn et al., 2012}). Habitat disturbance seems to be one of the factors that caused the isolation of populations from one another, the progressive reduction in pollen flow into stands and the reduction in the population size (^{Young and Clarke, 2000}; ^{Vergeer et al., 2003}; ^{Young and Pickup, 2010}).

Genetic structure. At the landscape level, the distribution of genetic ancestry in populations of Guaiacum sanctum in southern Mexico showed two groups of populations. Both, STRUCTURE and BARRIERS inferences are coincident and may represent a case of habitat fragmentation or could indicate historic differences associated with the colonization of distinct geographic regions from southern populations. For instance, the Barriers 4 and 7 divides the genetic group 1 (red genotype) in northern populations of Yucatán and Quintana Roo from the genetic group 2 (green genotype) observed in southern populations in Campeche. Even if, we observed very different geographically distributed genotypes, the AMOVA indicated that most of the genetic variation resides within populations while the differentiation among groups only accounted for F _ST (Φ _CT = 2.9 %, P = 001) and R _ST (Φ _CT = 3.6 %, P = 0.001).

Given the general low genetic differentiation in Guaiacum sanctum, we suggest that the detected genetic connectivity is a consequence of the dispersion ability of this species. Previous studies indicated the potential of G. sanctum for long distance pollen and seed dispersal (^{Wendelken and Martin, 1987}; ^{Fuchs and Hamrick, 2010}). The known floral morphology with purple color of petals and yellow anthers suggest pollination by bees and wasps which are able to disperse pollen over quite large distances (^{Holdridge et al., 1975}). Seed dispersion in G. sanctum occurs mainly by mammals and birds (^{Fuchs and Hamrick, 2010}). The flying capabilities of seed dispersers range between 5 to10 km on average in populations of G. sanctum in Costa Rica (^{Wendelken and Martin, 1987}; ^{Fuchs and Hamrick, 2010}). However, the observed pattern of fragmented habitats may act to impede gene flow due to remaining heterogeneity of the landscape across the Yucatán Peninsula.

Bottleneck detection and population size change. We observed evidence of recent bottleneck only for Campeche populations with the IAM model but not with TPM and SMM models (Table 3). This discrepancy in results between the IAM and TPM and SMM models is the consequence of different heterozygosity expectations at mutation equilibrium (^{Shaffer, 1981}; ^{Luikart and Cornuet, 1998}). Given that microsatellite mutation is thought to occur largely through the stepwise process, a combination of the SMM and IAM is expected to provide the best estimate of equilibrium heterozygosity for the bottleneck analysis. The absence of heterozygosity excess using both the strict SMM and the mixed TPM suggest that the contemporary population is at mutation-drift equilibrium (^{Luikart and Cornuet, 1998}). On the other hand, several authors have pointed out that these tests often failed to detect bottlenecks in populations known to have experienced population reductions because they have low statistical power as a result of limited sample size (^{Piry et al., 1999}; ^{Peery et al., 2012}). Power to detect significant bottlenecks based on the levels of heterozygosity is generally limited (^{Williamson-Natesan, 2005}; ^{Bouzat, 2010}). Therefore, results on bottleneck estimations should be taken with caution.

With respect to the effective population size for Guaiacum sanctum, ecological studies are in agreement with genetic estimates. Both estimations showed a tendency to low population size; populations from Yucatán have the highest N _e = 262 individuals, followed by Campeche with a N _e = 123 individuals and by Quintana Roo populations with a N _e = 100 individuals. Ecological studies registered historically higher densities of large trees of G. sanctum (up to 70 cm diameter at breast height - dbh) in Campeche, Yucatán, Quintana Roo, Oaxaca, and Chiapas, but remaining trees currently available in Campeche do not reach more than 55 cm in dbh (^{Oldfield, 2004}; ^{López-Toledo et al., 2008}). Population densities currently registered in Campeche seems to be acceptable with up to 1,200 potentially reproductive individuals per hectare (^{López-Toledo et al., 2008}). However, towards the edges of its distribution the situation drastically change in Oaxaca, Yucatán, and Quintana Roo with densities of only 150-470 individuals per hectare (^{Vester et al., 2007}; ^{López-Toledo et al., 2008}). Also, the reproductive system of G. sanctum seems to be impacted by fragmentation because fruits of trees that occurred at lower densities tend to produce less seeds than in non-fragmented habitats (^{Vester et al., 2007}). According to the rules used by conservation practitioner to estimate the population numbers, 50 individuals are needed to prevent an unacceptable rate of inbreeding, 500 to retain evolutionary potential and 5,000 to ensure overall genetic variability (^{Franklin, 1980}; ^{Frankham, 1996}; ^{Leimu et al., 2006}; ^{Frakham et al., 2014}). Several authors have suggested that population size, plant fitness and genetic diversity are generally associated, which probably due to the negative effects of small population size on genetic variation and plant fitness (^{Vergeer et al., 2003}; ^{Traill et al., 2007}). In G. sanctum, despite the severity of habitat fragmentation we observed a moderate effect on population genetic diversity and population bottleneck.

Concluding remarks. Greater population densities of Guaiacum sanctum occurred in porous karst terrain soils which are prevalent in northern Yucatán, the eastern portion of Quintana Roo and the southern portion of Campeche (^{Orellana et al., 1999}). In contrast, very low population densities occurred at the seasonally flooded terrains that prevail in the western portion of Campeche and northern and central Tabasco (^{Martínez and Galindo-Leal, 2002}; ^{López-Toledo et al., 2008}; ²⁰¹¹). The current geographic distribution of G. sanctum is the result of combined effects of environmental pressures (i.e., soil types) and extensive habitat loss (^{White and Hood, 2004}; ^{Vester et al., 2007}; ^{López-Toledo et al., 2008}; ²⁰¹¹) that in turn, may explain the observed population genetic diversity, connectivity and differentiation.

Guaiacum sanctum is facing a high risk of extinction and it has been classified as an endangered species by the International Union for the Conservation of Nature (IUCN) (^{CITES 2000}, ^{Vivero et al., 2006}). Successful strategies to conserve this species depend on detailed knowledge of the levels and distribution of genetic diversity within and between populations but most importantly we need to protect large fragments of remnant forests such as the population of G. sanctum at the natural biosphere reserve of Calakmul, Campeche to maintain the actual levels of genetic diversity of this species. Smaller fragments are also important to protect in order to maintain the genetic connectivity among populations of different regions in the Yucatán Peninsula.