Oaxaca state is one of the main hotspots of biodiversity in Mexico, containing almost 40% of the Mexican vascular flora, due to its high variability in habitat and climatic conditions coupled with high elevations in mountains and low elevations in valleys. We studied the genetic diversity and population structure of Quercus candicans, Quercus crassifolia, and Quercus castanea across their geographical distribution in Oaxaca state to understand how the heterogeneous physiography had driven the genetic diversity and population differentiation in these three oak species. We found high levels of genetic diversity but ca. 40% of the populations had significant values of Wright’s inbreeding coefficient. The analysis of molecular variance indicated that most of the variation occurred within populations in the three oak species. Resistance analyses showed connectivity among almost all the populations but barrier analysis found genetic breaks that limited gene flow among some populations of the oak species. Even in a heterogeneous environment such as in Oaxaca state, the oak species still have high levels of genetic diversity and landscape connectivity. However, it is necessary to maintain the genetic connectivity through the preservation of natural corridors with forests in good condition, which is necessary to maintain the cohesiveness of the species in the long term. It is also important to protect the centers of species diversity in Oaxaca state located in the subprovinces of Western Oaxacan Mountains and Valleys, Sierra Madre de Oaxaca, and Sierra Madre del Sur because they harbor most of the population genetic diversity and oak species richness, as has been shown in previous studies.
Introduction
The population genetic diversity of temperate forest species is spatially structured at different scales, from local (Holderegger & Wagner, 2008) to landscape (Ashley, Abraham, Backs, & Koenig, 2015; Sork et al., 2010) and continental (Petit et al., 2002; Wang et al., 2014; Wang, Zhang, Liu, Xu, & Liu, 2015) scales. At each scale, the ecological and genetic factors that structure the genetic diversity operate through environmental heterogeneity, species life history traits, and interactions with other organisms (Bacles, Lowe, & Ennos, 2004; Galván-Hernández, Lozada-García, Flores-Estévez, Galindo-González, & Vázquez-Torres, 2015; Honnay & Jacquemyn, 2007) determining the contemporary geographic distribution of the species.
The genetic connectivity among populations of a single species can be explained as a result of historical processes and physical barriers (e.g., mountain chains, large watercourses, or valleys) that influenced the dispersal routes (Manel & Holderegger, 2013; Sork & Waits, 2010) and the disruption of gene flow between populations at a regional scale (Galván-Hernández et al., 2015; Wang et al., 2015). Ecological traits, such as population density, pollen movement, and seed dispersal, determine the capacity of individuals to move through different environments and affect the reproductive success and population size (Breed et al., 2012; Craft & Ashley, 2010; O'Connell, Mosseler, & Rajora, 2006; Sagnard, Pichot, Vendramin, & Fady, 2011). In addition, human-mediated habitat fragmentation interrupts the natural connectivity between natural populations decreasing the genetic variability through an increase in endogamy and genetic drift resulting in isolated fragmented populations (Breed et al., 2012; Farwig, Böhning-Gaese, & Bleher, 2006; Herrera-Arroyo et al., 2013; Jump & Penuelas, 2006; Lowe, Boshier, Ward, Bacles, & Navarro, 2005; O'Connell et al., 2006).
Particularly, in regions recognized with high species richness, habitat loss had a significant alteration on the population dynamics and the structure of natural communities (Farwig et al., 2006). Species that formerly were widespread but currently are restricted to limited number of forest patches had altered patterns of gene flow at a landscape scale (Harrison & Hastings, 1996) with significant impacts in genetic connectivity (Herrera-Arroyo et al., 2013; Oyama et al., 2017).
In Mexico, one of the main hotspots of biodiversity is located in the state of Oaxaca, which contains almost 40% of the vascular flora of Mexico (García-Mendoza, 2004; Ortiz-Pérez, Hernández-Santana, & Figueroa-Mah-Eng, 2004). This high species richness is distributed throughout the mountains, valleys, and basins in the confluence of three major Mexican physiographic systems, which creates a mosaic of heterogeneous environments (Ferrusquía-Villafranca, 1993; Ortiz-Pérez et al., 2004). In particular, temperate forests in Oaxaca state are dominated by oak species represented by 52 species of Quercus out of the 161 reported in all of Mexico (Valencia, 2004; Valencia, Gómez, & Becerra, 2002). This oak species richness occupies the great number of habitats distributed through the state of Oaxaca both altitudinally and latitudinally (Ramírez-Toro, Torres-Miranda, Ruiz-Sanchez, González-Rogríguez, & Oyama, 2017). However, land-use cover changes as a result of forest fragmentation and the increasing human-mediated pressures are factors that produce extremely divergence among populations (Oyama et al., 2017; Velazquez et al., 2003).
In general, oak species in America are characterized by having high genetic variation within populations but low genetic differentiation among populations (e.g., Albarrán-Lara, Mendoza-Cuenca, Valencia-Ávalos, González-Rodríguez, & Oyama, 2010; Aldrich & Cavender-Bares, 2011; Peñaloza-Ramírez et al., 2010) contrasting with the patterns found for European oak species that have low genetic variation and high population structure (e.g., Grivet, Deguilloux, Petit, & Sork, 2006). However, it is unknown if these genetic patterns are maintained at a regional scale. We can expect that in a heterogeneous landscape such as in Oaxaca state, the population genetic structure of oak species could be higher due to the natural isolation of oak trees through the physiographic systems but still maintaining high genetic diversity within populations. We also decided to include three oak species due to the lack of genetic studies that include more than one species but also in order to obtain comparative patterns of population genetic structure and connectivity of congeneric species.
The main aim of this study was to analyze the genetic diversity and structure of three red oak species (Quercus candicans, Quercus crassifolia, and Quercus castanea) to understand how the heterogeneous physiography determines the genetic connectivity among populations across the geographical distribution of these species. The specific objectives of this study were to (a) estimate the genetic diversity and differentiation within and among populations of the three red oak species; (b) determine the levels of inbreeding, genetic connectivity, and effective population size of the three oak species; and (c) to propose conservation practices using the population genetic connectivity of key biological species to preserve temperate forests at a regional scale.
To the best of our knowledge, this is one of the first studies to use genetic connectivity as a surrogate of corridors in conservation biology. Furthermore, this study complements the biogeographical analysis of oak species in the same region giving insights of both species richness and genetic diversity (Ramírez-Toro et al., 2017).
Methods
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 (Na), Mean Observed Heterozygosity (HO), FIS the Mean Expected Heterozygosity (HE) and Inbreeding coefficient (FIS), for Populations of Q. candicans, Q. crassifolia, and Q. castanea in the Oaxaca state, Mexico.
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, dH2O, 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 102 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 (HO), and expected heterozygosity (HE), 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 106 Markov chain Monte Carlo repetitions following a burn-in period of 504 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 104 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 km2), 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).
Results
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 (FIS) 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, RST), 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 RST.
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)).
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 (Ne) obtained using LDNe for Q. candicans populations showed that Sierra Madre of Oaxaca group had the highest value (Ne = 118 individuals), followed by Western Oaxacan Mountains and Valleys (Ne = 62) and Sierra Madre del Sur (Ne = 61). For Q. crassifolia, the populations showed that the estimates of Ne in the Western Oaxacan Mountains and Valleys group had the highest value (Ne = 167 individuals), followed by Sierra Madre of Oaxaca (Ne = 123) and Sierra Madre del Sur (Ne = 122). Finally, for the Q. castanea populations, the Sierra Madre of Oaxaca group had the highest value (Ne = 374 individuals), followed by the Western Oaxacan Mountains and Valleys (Ne = 296) and Sierra Madre del Sur (Ne = 83) groups. All Ne 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.
Discussion
Genetic variation is necessary for the future viability of populations of any species (Lowe et al., 2005). In general, we found high genetic variation in the populations of Q. candicans, Q. castanea, and Q. crassifolia distributed in three biogeographic subprovinces of Oaxaca State, Mexico. The genetic diversity observed in the populations could be related to the biological traits of long-lived species, such as wind pollination, outcrossing, large effective population size, and extensive pollen flow among populations (Sork et al., 2002; Sork & Smouse, 2006; Young & Pickup, 2010). In general, the values of genetic diversity found in several Mexican oak species were higher than those reported in other Quercus species in North America and Europe (Aldrich & Cavender-Bares, 2011; Ashley et al., 2015; Fernández & Sork, 2005; Vranckx et al., 2014).
Oak species in Oaxaca have been subject to a long process of adaptation to microclimatic heterogeneity coupled with a very complex geological history, generating a heterogeneous matrix of species richness in different regions through the state (Ramírez-Toro et al., 2017; Valencia, 2004). Interestingly, populations with high genetic diversity coincide with detected areas of high oak species richness. Such is the case of the Sierra Madre of Oaxaca that harbored populations with high genetic diversity (HO = 0.866) and high oak species richness (38 oak species), followed by the Western Oaxacan Mountains and Valleys (HO = 0.855 and 25 oak species) according to the results of recent biogeographic study of oak species in Oaxaca (Ramírez-Toro et al., 2017). In this sense, the study of Valencia-Cuevas, Piñero, Mussali-Galante, Valencia-Ávalos, and Tovar-Sánchez (2014) showed that genetic diversity of populations of Q. castanea appears to be strongly related to the number of oak species growing in sympatry. Vellend (2005) suggested that regions with higher species diversity and genetic diversity could be connected as a consequence of an evolutionary force, that imprinted divergent selection on populations and communities but also causing changes in genetic diversity (Papadopoulou et al., 2011; Vellend, 2005).
Natural forests in Oaxaca State have been dramatically diminishing over the past 50 years (Velazquez et al., 2003). Most of the deforested areas have occurred in forests that have produced 46% of the pine and 85% of the oak wood (Gómez-Mendoza, Vega-Pena, Ramırez, Palacio-Prieto, & Galicia, 2006; Palacio-Prieto et al., 2000; Velazquez et al., 2003). Populations of the three oak species have been threatened in Oaxaca state as a result of habitat destruction. However, despite the high rates of deforestation, there is still moderated evidence of inbreeding in the oak populations as has been observed in other species (Ashley et al., 2015; Sork & Smouse, 2006; Young & Pickup, 2010). In several cases, populations showed significant and highly positive FIS values for Q. candicans (e.g., 4 out of 18 populations), Q. crassifolia (e.g., 3 out of 15 populations), and Q. castanea (e.g., 11 out of 20 populations) and were detected mainly on subprovinces of Western Oaxacan Mountains and Valleys, Sierra Madre of Oaxaca, and in a less number in Sierra Madre del Sur. Significant effects were also observed in the reductions in population size in Western Oaxacan Mountains and Valleys (Ne = 62.4) and Sierra Madre del Sur (Ne = 61.5) for Q. candicans and (Ne = 83.9) for a Q. castanea. In the study of Ramírez-Toro et al. (2017), they detected that the same subprovinces had the higher loss and vulnerability of habitat during the period between 2000 and 2010. The same trend has been reported for mammal, reptile, and pine species (Ortiz-Martínez, Rico-Gray, & Martínez-Meyer, 2008). An extensive forest disturbance in Oaxaca, such as habitat loss, expansion of pasture for cattle raising, construction of new roads, and settlements have had a strong impact in biodiversity (Harrison & Hastings, 1996; Ramírez-Toro et al., 2017; Velazquez et al., 2003), limiting the connectivity among populations, with some initial genetic concerns such as the increase of inbreeding and low effective population size. Additionally, it has been suggested that even in wind-pollinated and mostly outcrossing species such as oaks, decreasing the local tree density in a community might alter mate availability such that the number of nearby pollen sources surrounding a mother tree decreases and the pollination distance between mates increases (Bacles et al., 2004; Breed et al., 2012; Eckert et al., 2010; Frankham, Brook, Bradshaw, Traill, & Spielman, 2013; Sork et al., 2002; Vranckx et al., 2014). Seed dispersal in oaks is limited because the acorns lack adaptations for long-distance dispersal (e.g., acorns germinate closely to the mother tree or are transported by birds to short distances; (Jump & Penuelas, 2006; Sork et al., 2002; Sork & Smouse, 2006; Vakkari, Blom, Rusanen, Raisio, & Toivonen, 2006; Vranckx et al., 2014).
We did not find population genetic structure in any of the three oak species; most of the genetic variation was harbored within populations. We also found a great amount of connectivity between the populations located in the Western Oaxacan Mountains and Valleys and the Sierra Madre of Oaxaca, but leaving unconnected populations of the three oak species in the Western Oaxacan Mountains and Valleys and Sierra Madre del Sur. Differences in land cover and elevation may create multiple climatic zones across the Oaxaca state, characterized by variation in climate, topography, and vegetation; therefore, it is possible that habitat discontinuity and physical geographical barriers both play a substantial role in genetic differentiation and genetic inbreeding.
Oaxaca state has an intricate landscape, with mountains with high peaks crossed by valleys and canyons. Heterogeneous environments promote the maintenance of high genetic diversity as well as deep processes of divergence between populations of species (García-Mendoza, 2004; Ortiz-Pérez et al., 2004). Additionally, we observed evidence of landscape differentiation in the barriers to gene flow between Q. candicans and Q. crassifolia that split the western and the eastern portions of Western Oaxacan Mountains and Valleys subprovince. This evidence was also supported by the results of the climatic spatial resistance surface analysis. The western and the eastern portions of Western Oaxacan Mountains and Valleys subprovinces are separated by intermontane valleys covered with tropical deciduous forests in the central portion of the province (Ortiz-Pérez et al., 2004) acting as a geographical barrier for the disjunct distribution of the oak species in this zone.
We also detected another geographic barrier to gene flow between oak populations that separate populations located in the northeastern part of the Sierra Madre of Oaxaca from those in the northwestern Western Oaxacan Mountains and Valleys. The spatial resistance surface showed less connectivity among Q. candicans, Q. crassifolia, and Q. castanea between Western Oaxacan Mountains and Valleys and Sierra Madre of Oaxaca. As the driest zone developed in the Tehuantepec depression in Oaxaca, forests contracted, fragmented or disappeared, which affected the species composition on both sides of Western Oaxacan Mountains and Valleys and Sierra Madre of Oaxaca (García-Mendoza, 2004; Ortiz-Pérez et al., 2004). This valley is located in the northern portion of the Sierra Madre of Oaxaca to the east and the northeastern portion of the Western Oaxacan Mountains and Valleys. Most of the valley lies at elevations below 1,000 m but rises to slightly more than 1,400 m toward the south elevations, where it is covered by xerophytic vegetation and many succulent plants (Ortiz-Pérez et al., 2004), that act as an additional barrier to gene flow for Q. candicans, Q. crassifolia, and Q. castanea in the Western Oaxacan Mountains and Valleys and Sierra Madre of Oaxaca. Therefore, we suggest that high variability in habitat and climatic conditions combined with high elevations and large geographic distances may explain the moderated genetic differentiation and the inbreeding of some isolated populations.
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).
Acknowledgments
This paper constitutes a partial fulfillment of the Graduate Program (Posgrado en Ciencias Biológicas) of the Universidad Nacional Autónoma de México for W. R. T.
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: W. R. T. was supported by a scholarship from CONACYT fellowship 229358 and by Posgrado en Ciencias Biológicas, Universidad Nacional Autónoma de México. This project was supported by DGAPA-PAPIIT UNAM IN209108, IN229803, and IV201015; SEMARNAT-CONACYT 2004C01-97, 2006-23728; and CONACYT 38550 -V grants.