Collecting Sites and DNA Amplification

Samples of Q. candicans (18 sites), Q. crassifolia (15), and Q. castanea (20) were collected in temperate forests across the main physiographic provinces in Oaxaca state, Mexico (Table 1, Figure 1), along an altitudinal gradient from 150 to 3,300 m. Oaks are mainly distributed in five of the biogeographic subprovinces in which Oaxaca has been divided: Sierra Madre of Oaxaca, Western Oaxacan Mountains and Valleys, Central Mountains and Valleys, Región de los Chimalapas, and Sierra Madre del Sur (Ortíz-Pérez et al., 2004; Figure 1). We collected leaves of 12 to 20 trees separated by at least 30 m in each location. The samples were frozen until DNA extraction.

Table 1.

Locality Name, Sampling information, N/W Geographical Coordinates (Degrees), Altitude, Mean Number of Alleles (N_a), Mean Observed Heterozygosity (H_O), F_IS the Mean Expected Heterozygosity (H_E) and Inbreeding coefficient (F_IS), for Populations of Q. candicans, Q. crassifolia, and Q. castanea in the Oaxaca state, Mexico.

Figure 1.

The physiographic provinces of the state of Oaxaca (drawn from Ortiz-Perez et al., 2004). The names of the numbered provinces and their abbreviations are as follows: 1, Western Oaxacan Mountains and Valleys (WOMV); 2, Tehuacan Depression (TD); 3, Sierra Madre of Oaxaca (SMOax); 4, Gulf Coastal Plain (GCP); 5, Central Valleys of Oaxaca (CVO); 6, Central Mountains and Valleys (CMV); 7, Tehuantepec Isthmic Depression (TID); 8, Sierra Madre del Sur (Sierra Madre del Sur); 9, Pacific Coastal Plain (PCP); 10, Tehuantepec Coastal Plain (TCP). Points represented the localities of Q. candicans (Yellow), Q. crassifolia (Red), and Q. castanea (Blue).

Genomic DNA was extracted from 100 mg of fresh leaf material using the protocol proposed by Lefort & Douglas (1996). The isolated DNA was diluted with deionized water to a final concentration of 20 ng/µL and stored at −20℃. For the polymerase chain reactions (PCRs), samples of Q. candicans and Q. crassifolia were amplified with nine nuclear microsatellite loci by multiplexing in two multiplex PCR reactions (QIAGEN Multiplex PCR Kit). The first multiplex included the following five loci: three (quru-GAOA01, quru-GA2M04, and quru-GAOM05) developed by Aldrich, Michler, Sun, and Romero-Severson (2002) and two (QpZAG110 and QrZAG4) by Steinkellner et al. (1997) and Kampfer, Lexer, Glössl, and Steinkellner (1998). The second multiplex also included four loci (quru-GAIC08, quru-GA2F05, quru-GAOM07, and quru-GAIC06) developed by Aldrich et al. (2002). For Q. castanea, nine nuclear microsatellite loci were amplified by multiplexing in two PCR reactions (QIAGEN Multiplex PCR Kit). The first multiplex included five loci: quru-GA2FO5, quru-GAOC19, quru-GAIF02, quru-GAOA01 and quru-GAIF07 (Aldrich et al., 2002). The second multiplex included four loci: quru-GAIC06, quru-GAIC08, quru-GAOC11, and quru-GAOE09 (Aldrich et al., 2002). The PCR was performed using the QIAGEN Multiplex PCR Kit (QIAGEN) in a 5 -µl volume containing 1 × Multiplex PCR Master Mix, 2 µM of each primer, dH₂O, and 20 ng of template DNA. The thermal cycling conditions consisted of 35 cycles of 94℃ for 1 min, annealing for 1 min, extension at 72℃ for 2 min, and final extension at 72℃ for 10 min. Multiplex PCR products were combined with a GeneScan-500 LIZ size standard, and the analyses were performed using an ABI-PRISM 3100 Avant sequencer (Applied Biosystems). Fragments were analyzed and recorded using the Peak Scanner program 1.0 (Applied Biosystems).

Genetic Diversity

We tested for the presence of null alleles, large-allele dropout, and errors due to stuttering in the microsatellite data using the MICRO-CHECKER 2.2.3 program (Van Oosterhout, Hutchinson, Wills, & Shipley, 2004) with 10² bootstrap simulations and a 95% confidence interval. For each of the three oak species, we estimated the following genetic diversity parameters: number of alleles per locus (A), observed heterozygosity (H_O), and expected heterozygosity (H_E), using the GENETIX 4.02 software (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 2004).

Genetic Structure and Bayesian Admixture Analysis

The genetic ancestry of each individual plant and the three species separately was analyzed with the STRUCTURE 2.3.4 software (Falush, Stephens, & Pritchard, 2003; Hubisz, Falush, Stephens, & Pritchard, 2009; Pritchard, Stephens, & Donnelly, 2000). STRUCTURE uses a Bayesian clustering model to determine the proportion of each individual’s ancestry originating from different populations (Evanno, Regnaut, & Goudet, 2005). The optimal number (K) of groups was determined by varying K from 1 to 10 and running the analysis 10 times for each K value to find the maximum posterior likelihood [LnP (D)]. Each run was performed using 10⁶ Markov chain Monte Carlo repetitions following a burn-in period of 50⁴ iterations. We used an admixture model that allows the correlation of allele frequencies without any a priori information. Following Evanno et al. (2005), we determined the most likely value of K based on the maximum value of ΔK as implemented in the Structure Harvester 0.6.9 software (Earl & von Holdt, 2012).

A hierarchical test of population structure was conducted using the stepwise mutation model with analysis of molecular variance (AMOVA) in ARLEQUIN 3.5 (Excoffier & Lischer, 2010). We grouped the populations of the three oak species according to their geographical distributions (e.g., biogeographic subprovinces). Statistical significance was tested using 10⁴ permutations.

To identify geographic barriers and genetic breaks among populations of Q. candicans, Q. crassifolia, and Q. castanea, we used the Monmonier’s maximum difference algorithm with the BARRIER version 2.2 software (Manni, Guerard, & Heyer, 2004). BARRIER creates a map of the sampling locations from their geographical coordinates. Barriers are then represented on the map by identifying the maximum values within the population-pairwise genetic distance matrix. We used a pairwise matrix of average square distances (Goldstein, Ruiz-Linares, Cavalli-Sforza, & Feldman, 1995; Slatkin, 1995) estimated for the populations of Q. candicans, Q. crassifolia, and Q. castanea. Resampling random subsets of individuals within populations with the MSA program (Dieringer & Schotterer, 2003) provided 100 bootstrap replicate distances that were used to achieve statistical significance for the predicted barriers.

Changes in Population Size

We estimated the effective size (Ne) of populations of Q. candicans, Q. crassifolia, and Q. castanea using the LDNe software (Waples & Do, 2008). This program implements the bias-correction method developed by Waples (2006) to obtain Ne from a sample of S individuals. We set Pcrit = 0.02 (i.e., alleles with a frequency <0.02 are excluded), which generally provides a good balance between accuracy and bias (Waples & Do, 2008). Confidence intervals for Ne were calculated with the chi-square approximation implemented in LDNe (Waples & Do, 2008).

Ecological Niche Modeling and Spatial Connectivity

A database was built with georeferenced occurrence data for Q. candicans, Q. castanea, and Q. crassifolia. The presence records were obtained from herbarium specimens (CHAP, ENCB, FCME, IEB, IZTA, LL, MEXU, MO, XAL, and UNL), the Global Biodiversity Information Facility ( www.gbif.org), and data from monographic and floristic studies reviewed by Torres-Miranda, Luna-Vega, and Oyama (2011). A total of 26 presence records were gathered for Q. candicans, 67 for Q. castanea, and 50 for Q. crassifolia; these records were used to build the ecological niche model for each species in Oaxaca state.

To characterize the environmental niches, we used the 19 bioclimatic variables obtained from the WorldClim Project ( www.worldclim.com; Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), with a spatial resolution of 0.0083° (∼1 km²), and the result of interpolating monthly averages from weather stations throughout the world from 1950 to 2000. Then, we defined a polygon that covers Oaxaca state to establish the area of accessibility in order to run the ecological niche model. We used the algorithm MAXENT v. 3.3.3 k (Phillips, Anderson, & Schapire, 2006) to model ecological niche model performance, which uses the maximum entropy principle to estimate, from the existing values in the records of climate layers, a probability distribution that ranges from 0 to 1 for each pixel, which can be interpreted as an index of habitat suitability for the population being modeled (Elith et al., 2011). The algorithm compensates for colinearity between variables using a regularization method that addresses feature selection; therefore, there is less need to remove correlated variables (Elith et al., 2011) because the contribution of each one is ranked through the analysis. For model performance, we used 75% of the occurrence records for each of the three oak species for training, and the remaining percentage was used for testing. We used 1,000 iterations to limit convergence and a 0.00001 convergence limit during modeling, as this is the default level used by background test models. We also used a regularization value of 1. The options for extrapolation and clamping were disabled.

For the performance of MAXENT models, we used the AUC values of the receiver operating characteristic curve (Phillips et al., 2006), which allows the evaluation of the correspondence of the climatic suitability generated by the model with the known occurrences, ranging from zero to one, where 1 indicates perfect discrimination between presences and absences.

To generate spatial resistance distances for dispersal among sites based on habitat suitability maps (conductance grid), we imported the potential current distribution map for the three Quercus species separately into Circuitscape 4.0 (McRae & Beier, 2007). This program calculates pairwise resistance to gene flow among populations based on all possible paths, not just the least-cost path, thus better explaining the movement of genes among widely separated regions over many generations (McRae & Beier, 2007; McRae, Dickson, Keitt, & Shah, 2008). The input for Circuitscape is a raster dataset (habitat map) in which each cell is assigned a conductance value corresponding to the probability of the study organism moving through the habitat type encoded by the cell. We chose a conductance grid in which higher cell values denote greater ease of movement and applied a connection scheme that allowed gene flow among the four nearest cells. We used the potential distribution map rasters as input maps to quantify pairwise resistance distances between Q. candicans, Q. castanea, and Q. crassifolia populations across their distributions in Oaxaca state. We evaluated the connectivity between the three species separately using the results of Circuitscape, where low-resistance distances between pairs of nodes indicate effective dispersal between them, suggesting major connectivity between the populations of each oak species (McRae et al., 2008).

Genetic Diversity

No evidence of null alleles was found in any of the sample-loci combinations. The tests for errors due to stuttering and large-allele dropout yielded negative results in all cases. The values of the genetic diversity parameters estimated and grouped into biogeographical regions (i.e., subprovinces) were high for all the species (Table 1). For the populations of Q. candicans, the mean number of effective alleles per locus (Ne) ranged from 6.0 to 12.11, the observed heterozygosity (Ho) ranged from 0.675 to 0.881, and the expected heterozygosity (He) ranged from 0.642 to 0.854. Quercus crassifolia had Ne values ranged from 5.0 to 12.11, Ho from 0.764 to 0.959, and He from 0.704 to 0.888. Finally, Q. castanea also had high values of genetic diversity; the Ne ranged from 5.88 to 8.22, Ho from 0.599 to 0.881, and He from 0.721 to 0.826. Wright’s inbreeding coefficient within populations (F_IS) showed positive values in the majority of cases, denoting heterozygote deficiency, with 22 (out of 55) populations with significant values using a jackknife test (Table 1).

Genetic Structure and Bayesian Admixture Analysis

The AMOVA performed for the stepwise mutation model (i.e., stepwise mutation model, R_ST), comparing the biogeographic subprovinces showed that most of the variation occurred within populations in the three oak species (Table 2).

Table 2.

Analysis of Molecular Variance Performed on the nSSR Data and Grouped According to Biogeographic Subprovinces for Populations of Q. candicans, Q. crassifolia, and Q. castanea and Using R_ST.

The Bayesian admixture analyses indicated that the maximum posterior likelihood [LnP (D)] and the maximum ΔK value was K = 3 for Q. candicans (Figure 2(a)) and Q. castanea (Figure 2(c)) and K = 2 for Q. crassifolia (Figure 2(b)). The pie charts in Figure 3(a) show the distribution of ancestry proportions for Q. candicans at each collection site. Cluster 1 (red) is widespread across the Oaxaca subprovinces but only present in one population of Western Oaxacan Mountains and Valleys. Cluster 2 (green) is consistently structured across the three subprovinces. Cluster 3 (blue) has a wider distribution (Figure 3(a)). For Q. crassifolia, both genetic clusters are widely distributed across the three subprovinces (Figure 3(b)). Finally, all genetic clusters for Q. castanea were widespread across the subprovinces (Figure 3(c)).

Figure 2.

Maximum probability to the most probable genetic group of Q. candicans (a), Q. crassifolia (b), and Q. castanea (c).

Figure 3.

Distribution of genetic ancestry proportion groups corresponding to Q. candicans (a), Q. crassifolia (b), and Q. castanea (c) populations across physiographic regions in Oaxaca. Black bars principal barriers detected with the BARRRIERS program.

The BARRIER analysis performed using 100 bootstrap replicates for Q. candicans showed five main breaks or genetic discontinuities (Figure 3(a)). The first break (B1), with 95% bootstrap support (BS), separated the populations of Sierra Madre of Oaxaca from those in Western Oaxacan Mountains and Valleys. B2 (BS = 89%) separated the Tlazoyaltepec population from the rest of the populations in Western Oaxacan Mountains and Valleys. B3 and B4 (BS = 79%) separated populations at the edge of the distribution from neighboring populations located in Sierra Madre del Sur and Sierra Madre of Oaxaca. B5 split populations in Sierra Madre del Sur and Central Mountains and Valleys (Figure 3(a)). For Q. crassifolia, six breaks were identified (Figure 3(b)). B1 and B4 (BS = 90%) separated the eastern from the western part of Western Oaxacan Mountains and Valleys and Sierra Madre of Oaxaca. B2 (BS = 80%) divided the Paxtlan population from the rest of the populations in the eastern portion of Sierra Madre del Sur. B3 (BS = 75%) separated the Coatlán population from the populations located in the southeastern portions of Sierra Madre del Sur. B5 (BS = 62%) separated the populations from the eastern portion of Western Oaxacan Mountains and Valleys from Tlaxiaco and Juxtlahuaca in the western portion of Western Oaxacan Mountains and Valleys (Figure 3(b)). For Q. castanea, four breaks were identified (Figure 3(c)). B1 and B2 (BS = 95 and 90%, respectively) separated the populations from the eastern portion of Western Oaxacan Mountains and Valleys from Tlaxiaco and Juxtlahuaca in the western portion of Western Oaxacan Mountains and Valleys. B3 (BS = 70%) divided the Colorada population from the rest of the populations in the southeastern portion of Sierra Madre del Sur. B4 (BS = 62%) separated the populations in northwest Sierra Madre of Oaxaca from those in the northeastern portion of Sierra Madre of Oaxaca (Figure 3(c)).

Changes in Population Size

Estimates of effective population size were performed using the subprovinces to get an idea about changes in population size across the landscape. The effective population size (N_e) obtained using LDNe for Q. candicans populations showed that Sierra Madre of Oaxaca group had the highest value (N_e = 118 individuals), followed by Western Oaxacan Mountains and Valleys (N_e = 62) and Sierra Madre del Sur (N_e = 61). For Q. crassifolia, the populations showed that the estimates of Ne in the Western Oaxacan Mountains and Valleys group had the highest value (N_e = 167 individuals), followed by Sierra Madre of Oaxaca (N_e = 123) and Sierra Madre del Sur (N_e = 122). Finally, for the Q. castanea populations, the Sierra Madre of Oaxaca group had the highest value (N_e = 374 individuals), followed by the Western Oaxacan Mountains and Valleys (N_e = 296) and Sierra Madre del Sur (N_e = 83) groups. All N_e estimates had high jackknife support and a good confidence interval (Table 3).

Table 3.

Analysis Results Obtained for the Estimation of the Population Effective Size, Tested for the Three Physiographic Groups, Values Obtained With the Program LDNe, for Q. candicans, Q. crassifolia, and Q. castanea Populations in the Oaxaca Mexico.

Ecological Niche Modeling and Spatial Connectivity

The total number of occurrences allowed us to develop with good performance with high AUC values (0.862, 0.927, and 0.921) and shown the potential current distribution in Oaxaca for Q. candicans (Figure 4(a)), Q. castanea (Figure 4(b)), and Q. crassifolia (Figure 4(c)), respectively. The climatic spatial resistance surface for Q. candicans (Figure 4(d)) showed major connectivity among populations in Western Oaxacan Mountains and Valleys, Sierra Madre of Oaxaca, and Sierra Madre del Sur, leaving one unconnected population in Western Oaxacan Mountains and Valleys, with a resistance value ranging from 0 to 6.97 and an average resistance distance of 3.66 among populations. For Q. castanea, the climatic spatial resistance surface (Figure 4(e)) showed major connectivity between Western Oaxacan Mountains and Valleys and Sierra Madre of Oaxaca, leaving unconnected populations to the west of the Western Oaxacan Mountains and Valleys and Sierra Madre del Sur, with resistance values ranging from 0.48 to 15.23 and an average resistance distance of 3.30 among populations. Finally, the climatic spatial resistance surface for Q. crassifolia (Figure 4(f)) showed major connectivity among Western Oaxacan Mountains and Valleys, Sierra Madre of Oaxaca, and Sierra Madre del Sur, leaving two unconnected populations to the west of the Western Oaxacan Mountains and Valleys, with resistance values ranging from 0.42 to 6.99 and an average resistance distance of 3.36 among populations. The average resistance values among the three oak species were quite similar.

Figure 4.

Potential distribution of (a) Q. candicans, (b) Q. castanea, and (c) Q. crassifolia identified using ecological niche modeling, prediction of suitable habitat in the current environment. Black dots are the presence records used to build the ecological niche model. On the right side (d) connectivity map among populations of Q. candicans. (e) Connectivity map among populations of Q. castanea and (f) connectivity map among populations of Q. crassifolia. Black dots are the populations records used for each species in the genetic analyses to build the connectivity maps.

Implications for Conservation

We studied the genetic structure of populations of Q. candicans, Q. crassifolia, and Q. castanea in Oaxaca state, one of the areas with high biodiversity in Mexico. High levels of genetic diversity were observed in populations of these species across their distribution along some of the main physiographic subprovinces, indicating that sites with high oak species richness (Ramírez-Toro et al., 2017) also had high genetic diversity. Recent studies have suggested that species diversity and genetic diversity could be interconnected as a consequence of a process acting in parallel at these two levels (Papadopoulou et al., 2011; Vellend, 2005). However, this biological richness is threatened by the effects of land cover changes, such as deforestation; Oaxaca state has lost over half a million hectares of natural forests during the last 20 years (Gómez-Mendoza et al., 2006; Velazquez et al., 2003).

The conservation of forests in Oaxaca state is necessary because most of the natural areas are not currently protected by the government. Urgently, the implementation of real public policies to include more conservation areas is needed to preserve one the greatest hotspot in Mexico with the greatest amount of vegetation types and most of the forests harboring many endemic species. There are only five natural protected areas in this state. In Oaxaca, the initiative called Indigenous Conservation and Community Areas is another alternative to manage and preserve the remaining natural areas with different levels of disturbance (Anta-Fonseca & Sanchez, 2009; Martin et al., 2011; Robson, 2007). However, according to Martin et al. (2011), 126 indigenous conservation and community areas protect 375,457 ha, but only 22 include forests with oak species. Connectivity of forest fragments plays an essential role eluding the isolation of organisms in altered landscapes that allow their displacement (Ramírez-Toro et al., 2017; Torres-Miranda et al., 2011), and promoting the genetic exchange and maintaining of ecological processes (Oyama et al., 2017; Robson, 2007). Therefore, it is necessary to protect remnant natural areas with high species richness and high genetic diversity as well as regions that could serve as natural corridors to maintain natural biological processes, such as dispersal and pollination (Herrera-Arroyo et al., 2013; Oyama et al., 2017).