Genetic structure of a population can be defined by the resistance of the landscape, the distance between populations, or geographic barriers. We analyzed the population genetic structure of Abies religiosa on a fine spatial scale and examined isolation models by resistance, distance, and barrier. We collected vegetative tissue from populations located at the altitudinal extremes of the distribution range of the species on three slopes of La Malinche National Park (LMNP) (South, North, and East) in central Mexico. Genomic DNA was obtained using the CTAB 2X method, and eight microsatellite chloroplast loci were amplified. The genetic structure was identified based on a Discriminant Analysis of Principal Components with cross-validation and a spatial Principal Component Analysis using the Gabriel-type connectivity network. The isolation hypotheses were evaluated by constructing partial Mantel tests using Reciprocal Causal Modeling and Maximum Likelihood Population Effects models. A genetic structure of isolation by barrier was identified, and two genetic groups were recognized: one including populations of the South slope and the other comprising populations of the North and East slopes. The greatest genetic exchange between groups was recorded in populations located at higher altitudes. It is suggested to promote the connectivity between slopes through assisted migration and immediately halt land-use changes, as part of the actions to preserve genetic diversity at LMPN. This study contributes to the knowledge of the spatial genetic structure of species at risk that are components of the Mexican temperate forest. 

Vegetative Tissue Collection

We toured three slopes of the LMNP that are relatively easy to access (North, East, and South) because of the orography and the type of ownership (social) of agricultural plots (López-Téllez et al., 2019), aiming to identify the relative abundance of the species within its altitudinal range (2400–3600 m a.s.l.; Sáenz-Romero et al., 2012). Two collection sites were selected on each slope, which were assumed to be distinct populations located at the upper and lower limits of its altitudinal range, with a total of six populations throughout the LMNP (2 populations × 3 slopes) (Fig. 1). We established a central point in each population, from which two transects were drawn in opposite orientations. Along each transect, vegetative material was collected from 20 individuals separated by at least 30 m (Fig. 1). The geographic location of each individual was recorded; then, young undamaged needles were collected, transported, and stored in 1.5 mL Eppendorf tubes with Tris-EDTA pH 8.0 (Sigma-Aldrich) buffer. 

Laboratory

Tissue samples were first ground with liquid nitrogen (Doyle, 1990); afterward, genomic DNA was extracted using the CTAB 2X method. The extractions were visualized in 1% agarose electrophoresis to determine their viability. Unsuccessful samples were purified using Wizard ® SV Minicolums (Promega) columns and Wash Buffer 2 (Qiagen) following the vendor-standardized procedure, although modifying the concentrations to improve the result. 

Ten chloroplast microsatellites (cpSSR) designed for Pinus thunbergii and P. leucodermis were amplified (Vendramin et al., 1996). Amplifications were run with a final volume of 14 µl using the Master Mix (Taq DNA Polymerase; Qiagen) solution in a T100 Touch^TM Thermal Cycler (Bio-Rad). The PCR conditions proposed by Vendramin et al. (1996) were modified: initial denaturation at 95 ºC for 5 min, then 32 cycles with denaturation at 94 ºC for 1 min, alignment at 50 ºC–58 ºC for 1 min (see Table 1), followed by extension at 72 ºC for 1 min. The final extension step was at 72 ºC for 8 min. The PCR product was read by capillary electrophoresis (QIAxcel, Qiagen) using the method OM500, with a 10-bp resolution for fragments of 100–500 bp. Amplicons were considered different when they were >10 bp (Qiagen, 2008); amplicon size was determined with the ScreenFel software (Qiagen v 1.0.2.0; Ambion Inc., Austin TX) provided by the QIAxcel system, using the 15 bp/500 bp QX Alignment Marker and the 25–500 bp QX DNA Marker. The binning of the fragments obtained was carried out with the program Allelogram v. 2.2 (Morin et al., 2009).

Genetic Diversity

Since the chloroplast has a haploid genome, it has paternal inheritance in conifers and is not subject to recombination (Neale and Sederoff, 1989; Watano et al., 1996); hence, it can be considered as a locus (Rai and Ginwal, 2018). It can also be perceived as a circular haploid chromosome in which sequence variation generates different alleles within individual non-recombinant loci (Echt et al., 1998). Therefore, this study analyzed genetic diversity based on haplotypes formed by combining fragments from the eight successfully amplified microsatellite loci (cpSSR) (Rai and Ginwal, 2018). Genetic diversity was described through the number of observed alleles (N_a), number of effective alleles (N_e), Shannon-Weiner index (I), haplotype diversity (h), percentage of polymorphic loci (%P), and genetic distances (Nei, 1987), using the program GenAlEx v. 6.4 (Peakall and Smouse, 2006). 

Genetic Structure

The pattern of spatial genetic variation was evaluated using a Discriminant Analysis of Principal Components (DAPC) with cross-validation (Jombart et al., 2010). Connectivity was determined by constructing a spatial Principal Components Analysis (sPCA) (Jombart et al., 2008) through a Gabriel-type network. This network was selected because it has moderate saturation, which makes networks more informative than saturated ones (Dyer and Nason, 2004; Naujokaitis-Lewis et al., 2013). All analyses were carried out with the software R v. 4.0.3 (packages adegenet, akima, poppr, ResistanceGA, and tess3r; R Core Team, 2020). 

To evaluate the hypotheses of isolation, we constructed matrices of resistance distances (IBR) and of barrier distances (IBB) between pairs of populations, with the program Circuitscape v. 4 (McRae and Shah, 2009). To this end, for IBR we elaborated a raster from the visual categorization of land use (farming, forest, secondary vegetation, urban, bare ground) with the program ArcGIS v. 10.4 (ESRI, 2018), from Landsat 5 images (30 m × 30 m resolution) downloaded from GloVis (https://glovis.usgs.gov). And for IBB, we obtained an elevation raster of the study area using the Digital Elevation Model downloaded from ASF Data Search Vertex (https://search.asf.alaska.edu/#/; Alos Pasar images, 12.5 m resolution). In the case of the IBD, we used linear geographic distances between population pairs calculated with the program GenAlEx v. 6.4 (Peakall and Smouse, 2006). 

Once the matrices of each isolation hypothesis were constructed, partial Mantel tests were performed with the genetic distances (Nei, 1987), using Reciprocal Causal Modeling (RCM) to avoid erroneous conclusions due to the high correlation of alternative models (Cushman et al., 2006; Cushman and Landguth, 2010). In addition, to identify the factors that contribute to gene flow, Maximum Likelihood Population Effects (MLPE) models were performed, these linear random effects models account for the lack of independence between pairwise comparisons (Row et al. 2017; Burgess and Garrick 2020). The analyses were carried out using the software R v. 4.0.3 (packages MuMIn, vegan, and ComMLPE; R Core Team, 2020).