Introduction

The Mexican temperate forest is the ecosystem ranking second in cover loss due to anthropogenic activities (Guerra-de la Cruz & Galicia, 2017; Semarnat, 2016). Loss of habitat and connectivity seriously affect genetic diversity (Aavik et al., 2019; Monteiro et al., 2019). Changes in the landscape modify the quantity, quality, and configuration of the habitat and the matrix (i.e. the dominant land cover of a landscape, e.g. forest or farming), which influence the effective population size, mating, and migration (DiLeo & Wagner, 2016). All these factors affect evolutionary processes (gene flow and genetic drift) that, as a whole, determine genetic diversity and its spatial structure (Goncalves et al., 2019; Herrera-Arroyo et al., 2013; Moran & Clark, 2011). After a disturbance event, often original continuous populations decrease in size and are spatially isolated; the new conditions limit gene flow between populations and increase the effects of genetic drift and inbreeding, with consequences in the levels of genetic differentiation (Evanno et al., 2009; Goncalves et al., 2019; Manel & Holderegger, 2013).

In this scenario, the isolation by resistance (IBR) model proposes that gene flow is driven by the permeability of the matrix, or an environmental or landscape characteristic (McRae, 2006; McRae & Beier, 2007). The response of a species to landscape changes depends on the organism, its life-history traits, and its landscape context (Aavik et al., 2019). For example, in species with specific environmental requirements, alterations in connectivity associated with habitat fragmentation and loss are more severe due to their relative specialization (Aldrich & Hamrick 1998). Alternatively, the distribution of genetic diversity can also be explained by the classic isolation by distance (IBD) model, which assumes that geographic distance limits genetic dispersal and flow (Wright, 1943).

La Malinche National Park (LMNP) is a Protected Natural Area (PNA) that is part of the biological corridor for gene flow between volcanoes of the Trans-Mexican Volcanic Belt (TVB) (e.g., Popocatépetl, Iztaccíhuatl, Nevado de Toluca). This area is home to a temperate forest that covers 206.07 km² where a great diversity of species is conserved (Domínguez-Domínguez & Pérez-Ponce de León, 2009; Fernández & López-Domínguez, 2005). The first most important forest in the LMNP is dominated by Pinus montezumae, and the second by Abies religiosa (Kunth) Schltdl. & Cham. (oyamel fir), covering an area of 13.66 km² (6.6%) in the habitat of high mountain species (Anderson & Brower, 1996; Snook, 1993). This forest is characterized by a high carbon uptake capacity and soil formation, making it a forest of great importance for conservation in Mexico. Abies religiosa populations show a discontinuous distribution in the LMNP due to severe changes in land use, intentional fires, and illegal logging affecting the habitat of this species (i.e., high humidity, low temperature, and shade; Gallardo-Salazar et al. 2019) and having fragmented at least 56% of the original forest cover in the park (López-Domínguez & Acosta, 2005; Valdéz-Pérez et al., 2016).

Some genetic studies have been carried out in A. religiosa. For example, Aguirre-Planter et al. (2000) identified a relationship between gene flow and geographic distances in TVB populations based on isoenzymes. Heredia-Bobadilla et al. (2012) reported a high genetic variability between populations and low levels of gene flow in populations inhabiting the Nevado de Toluca National Park, using chloroplast DNA (cpDNA; cp-trnK, cp-TRNG-trnQ) and mitochondrial DNA (mtDNA; nad-5, coz3in). Also, Jaramillo-Correa et al. (2008), through mtDNA and cpDNA, revealed that the Mesoamerican Abies species share a recent common ancestor and that their divergence and speciation have been determined mainly by genetic drift and isolation during warm interglacial periods. Méndez-González et al. (2017) identified a high genetic diversity and a variable population genetic structure according to the genetic marker, i.e., four genetic groups with amplified fragment length polymorphism (AFLP) and two with chloroplast microsatellites (cpSSRs). Genetic diversity, genetic differentiation between populations, and the associated factors are unknown in LMNP. For this reason, the objective of the present study was to determine whether the remaining populations of A. religiosa are genetically structured and, if so, what factors best explain its gene flow.

Since chloroplast microsatellites (cpSSR) (a) are highly polymorphic, (b) are of a non-recombinant genome, and (c) reveal paternal inheritance (pollen) in conifers, they are considered suitable markers for studying the genetic diversity and structure associated with demographic changes and landscape factors such as habitat fragmentation (Mittal & Dubey, 2009; Navascues & Emerson, 2005; Wheeler et al., 2014). Therefore, we consider that cpSSRs are adequate markers to explore the consequences of landscape changes on gene flow in A. religiosa. Since LMNP maintains a high environmental heterogeneity given the altitudinal range and changes in land use (Semarnat-Conanp, 2013), we expected to find high genetic differentiation levels between populations, explained by isolation by resistance (IBR).

This study contributes to the knowledge of the impacts of anthropogenic activities on populations of native tree species distributed in areas with heavy forest cover loss from changes in land use. This information is essential for defining conservation strategies for vulnerable temperate forest species.

Methods

The LMNP is located in central Mexico, to the east of the TVB, in the states of Tlaxcala and Puebla, between coordinates 19°08’–19°20’ N and 98°08’– 97°55’. The park covers an area of 460.93 km² with an altitudinal range from 2,200 to 4,461 m a.s.l. (Figure 1). Depending on the altitude, the mean annual temperature varies between 2 °C and 17 °C, with rains in summer (Fernández & López-Domínguez, 2005). The local vegetation also varies according to altitude, slope, orography, and anthropogenic activities; the dominant species are Pinus montezumae, A. religiosa, Pinus hartwegii, Pinus leiophylla, and Alnus jorullensis (Rojas-García & Villers-Ruíz, 2008). The predominant land use is farming, followed by forest and secondary vegetation (Figure 1c).

Figure 1.

Location of the populations of Abies religiosa in La Malinche National Park used for collecting vegetative tissue and describing the spatial genetic structure.

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

Table 1.

Loci of Chloroplast Microsatellites (cpSSR) Designed for Pinus Thunbergii and P. Leucodermis (Vendramin et al., 1996) Used to Estimate the Genetic Diversity and Structure of Abies Religiosa in La Malinche National Park, Mexico. AT, Alignment Temperature; ✓, Successfully Amplified; ×, Not Amplified.

Results

Of the ten microsatellites tested, successful amplifications were obtained for only eight (Table 1). Of these, we identified 27 fragments or alleles that formed a total of 116 haplotypes; of these, only four haplotypes were shared among individuals.

The number of observed haplotypes (H_a) and the number of effective haplotypes (H_e) were higher in the populations inhabiting the East and North slopes (H_a = 20, H_e = 20) and lower in the populations of the South slope (H_a = 19, H_e = 18). The Shannon-Weiner index (I) ranged from 2.93 in the two populations of the South slope to 2.99 in the populations of the East and North slopes. Haplotype diversity (h) was high in all populations (h = 0.94–0.95) (Table 2).

Table 2.

Genetic Diversity of Six Populations (Pop) of Abies religiosa in La Malinche National Park, Based on Eight Chloroplast Microsatellites (Vendramin et al., 1996). H_a, Number of Observed Haplotypes; H_e, Number of Effective Haplotypes; I, Shannon-Weiner Index; h, Haplotype Diversity.

The AMOVA identified the highest percentage of variation within populations (97.11%), and a low but significant genetic differentiation (F_st = 0.03, P = <0.05). The DAPC grouped together the two populations of the South slope and the high-altitude population of the North slope (NP4); another group included the populations of the East slope (EP1, EP2) and the low-altitude population of the North slope (NP3) (Figure 2). Separately, the connectivity analysis identified a global spatial structure (max(t)= 0.007, P = 0.05), and the connectivity network identified two genetic groups: one including the populations of the South slope (SP5 and SP6) and the other the populations of the East and North slopes (Figure 3).

Figure 2.

Discriminant Analysis of Principal Components of six populations of Abies religiosa in La Malinche National, based on eight chloroplast microsatellites (Vendramin et al., 1996). DA, discriminant analysis; EP1, low-altitude population of the East slope; EP2, high-altitude population of the East slope; NP3, low-altitude population of the North slope; NP4, high-altitude population of the North slope; SP4, low-altitude population of the South slope; SP5, high-altitude population of the South slope.

Figure 3.

Spatial Principal Component Analysis (sPCA) using the Gabriel-type connectivity network of six populations of Abies religiosa in La Malinche National Park, based on eight chloroplast microsatellites (Vendramin et al., 1996). EP1, low-altitude population of the East slope; EP2, high-altitude population of the East slope; NP3, low-altitude population of the North slope; NP4, high-altitude population of the North slope; SP4, low-altitude population of the South slope; SP5, high-altitude population of the South slope. The color of boxes indicates genetic groups.

The partial Mantel tests showed a statistically significant association between genetic distances (GD) and elevation (IBR_e) (Table 3). The Reciprocal Causal Model detected altitude as the best explanatory variable, followed by distance (Table 4). Likewise, the MPLE determined that the model where genetic distances are explained by IBR_e is the best model (AICc = -60,088, logLik = 36.04, P = < 0.001; Table 5). The contribution of variables to explain GD was 75% for IBR_e and 20% for IBD in the MPLE. IBR_l-uc was not detected with either analysis (Tables 4 and 5).

Table 3.

Partial Mantel Tests Between Genetic Distances (GD) of Pairs of Populations of Abies religiosa in La Malinche National Park and Isolation by Landscape Resistance for Land-Use Changes (IBR_l-uc) and Elevation (IBR_e), and Isolation by Geographic Distance (IBD). y, Dependent Variable, Genetic Distances; x, Explanatory Variables for Genetic Distances; z, Controlled Variables for Partial Mantel Tests.

Table 4.

Partial Mantel Tests, with Reciprocal Causal Modeling (RCM) of Genetic Distances (GD) Between Populations of Abies religiosa in La Malinche National Park, Based on Eight Chloroplast Microsatellites, and Isolation by Resistance for Land-Use Changes (IBR_l-uc) and Elevation (IBR_e), and Isolation by Distance (IBD). The Most Supported Relationships for Each Pair of Tests According to the Relative Value of the Reciprocal Causal Modeling are Marked in Bold.

Table 5.

Maximum Likelihood Populations Effects Models to Assess Isolation by Resistance for Land-Use Changes (IBR_l-uc) and Elevation (IBR_e), and Isolation by Distance (IBD) of Six Populations of Abies religiosa in La Malinche National Park, Based on Eight Chloroplast Microsatellites. logLik, Model Fitted by Maximum Likelihood; AICc, Corrected Akaike Information Criteria (Hurvich & Tsai, 1989); Δ, Variance.

Discussion

We identified a spatial genetic structure in the populations of A. religiosa inhabiting LMNP, probably resulting from isolation by resistance to elevation. Gene flow is seemingly not prevented by changes in land use since high functional connectivity (gene flow performed) was observed between the populations of the East and North slopes; however, the evident isolation of the populations of the South slope is explained by the increase in genetic drift as a consequence of the reduction in population size.

The present study did not detect isolation by resistance associated with land-use changes (IBR_l-uc), despite the fact that a continuous and drastic change in land use occurs in the study area, likely because the effects of anthropogenic disturbance are delayed in conifer species due to their longevity (Aavik et al., 2019). Land-use changes disrupt the spatial continuity of the habitat, which inevitably decreases population size and isolation, both of which reduce gene flow and increase genetic drift, thus impacting genetic variation and differentiation (Eckert et al., 2010; Leroy et al., 2018). Conifers usually show high genetic diversity and low genetic differentiation with nuclear markers due to anemochoric dispersal. However, this is not the case with single-parent genomes such as chloroplast DNA, in which the effective population size is approximately half that of nuclear genes, reducing the time to fix an allele because of the greater influence of genetic drift (Clark et al., 2000). In the family Pinaceae, the chloroplast genome is subject to paternal inheritance and dispersed first by pollen and later by seeds. As a result, propagation by both vectors in each reproductive cycle means that the migration-drift balance is reached more rapidly in chloroplast than in other genomes (e.g., nuclear, mitochondrial) (Mogensen, 1996); consequently, markers efficiently reveal recent stochastic processes (Heuertz et al., 2003; Méndez-González et al., 2017; Petit et al., 1993, 2005), which could suppose the observation of the demographic dynamics associated to the resistance of the landscape.

In A. religiosa, each generation comprises about 27 years since sexual maturity is reached at 17 to 25 years, and the cycle from pollen release to seed takes another two years (Mantilla, 2006; Ortiz-Bibian et al., 2019). The detected genetic structure may not reveal the actual effect of land-use changes, so it would be important to include recent generations (i.e., seedlings) to determine potential changes in connectivity.

Despite we do not rule out the possible underestimation of genetic diversity due to genotyping technique used (capillary electrophoresis) which does not have a detection of variation < 10 bp (Qiagen, 2008), this study detected a high haplotype diversity (h = 0.94–0.95), at levels consistent with the results of previous studies using chloroplast microsatellites in Abies species. For example, Clark et al. (2000) reported an h value of 0.95 in A. balsamea, whereas Parducci et al. (2001) reported h values of 0.93–0.99 in A. alba, 0.97 in A. numidica, and 0.98–1.0 in A. cephalinica. For their part, Jaramillo-Correa et al. (2008) found haplotype diversity values of 0.93 in A. gualmensis, 0.94 in A. hickeli, and 0.91 in populations of A. religiosa in the TVB; the last data is similar to the one reported by Méndez-González et al. (2017) in populations of this same species at El Ajusco, in central Mexico (h = 0.94).

An important result of our research is the lower haplotype diversity identified in the populations of the South slope compared to h values in populations of the East and North slopes, which may be due to a more intense genetic drift associated with the significant reduction in population size (Leroy et al., 2018; Song et al., 2006); it is possible that the climatic conditions of the North and East slopes that favor the best development of A. religiosa populations (Allende et al., 2016) help to maintain genetic diversity in these populations. An intense change in land use has been observed in the LMNP since the decade of 1980 (Ríos, 2014). For example, Ern (1976) reported associations between A. religiosa and other conifers in areas where the former is currently absent (Cruz-Salazar, 2021, personal observation). The loss of tree cover, including A. religiosa forests, is even more serious on the South slope, where there has been uncontrolled illegal logging (López-Téllez et al., 2019; Rojas-García and Villers-Ruiz, 2008; Valdez-Pérez et al., 2016). This activity has dramatically reduced populations of tree species of economic interest, affecting the microhabitat conditions necessary for germination and recruitment of these species, particularly those that need shade and moisture, such as A. religiosa (e.g. shade, humidity) (Gallardo-Salazar et al., 2019; Rzedowski, 2006).

The structure genetic observed in this study is possibly due to the exchange of genes through pollen. In conifers, pollen dispersal can reach up to 10 km (Molina-Sánchez et al., 2019), while seeds reach only 31 m in A. alba, transported mainly by wind and some animals (Briggs et al., 2009; Ortiz-Bibian et al., 2019). At elevated sites such as volcanoes, anemochoric dispersal depends on regional wind patterns resulting from the global variation of wind regimes, which either facilitates or restrains the displacement of pollen and seeds (Kling & Ackerly, 2020). In the LMNP, the dominant winds flow from the southeast during autumn and winter and from the northeast in spring and summer (Fernández & López-Domínguez, 2005); this may result in constant gene flow and low genetic structure despite the changes in land use and the obvious geographic barriers (e.g. Barranca Grande).

In A. religiosa, sexual structures are produced in winter (December), while pollination takes place in spring (April-June), followed by seed dispersal (October-December). However, the limited period of pollen dispersal (from six days up to one month) and the lack of formation of female cone production due to limited resources may prevent effective cross-pollination and increase inbreeding levels (Mantilla, 2006). Likewise, it has been reported that low population density decreases seed viability and germination percentage with influence on inbreeding (Ortiz-Bibian et al., 2019); in this study low population density may be due to deforestation or habitat loss but we do not rule out the effect of sampling conducted at the extremes of the altitudinal distribution of the species, where population density might decrease.

Our results indicate a spatial genetic structure determined mainly by IBR to elevation and two genetic groups; these findings are consistent with those reported by Méndez-González et al. (2017) in populations of El Ajusco, Mexico, based on chloroplast microsatellites. Interestingly, the populations inhabiting the South slope of the LMNP were genetically connected with the higher-altitude populations of the other slopes, suggesting that dispersal seemingly occurs through the volcano peak. This assumption is supported by the presence of an important geographic barrier (Barranca Grande) between the East and South slopes (Semarnat-Conanp, 2013), which possibly limits gene flow.

Implications for Conservation

Studying landscape resistance to gene flow is fundamental for the management and conservation of viable populations in modified landscapes (Oyler-McCance et al., 2013). A management unit is a group of populations that maintain a constant gene flow (Barbosa et al., 2009; Coates et al., 2018). In the present study, we detected two management units, one formed by the populations of the South slope and the other by the populations inhabiting the East and North slopes. Although land-use changes do not seem to determine isolation, these are a threat that will eventually affect the population genetic structure. Therefore, starting assisted migration planning is advisable (Sáenz-Romero et al., 2012) to decrease inbreeding and promote the conservation of the genetic diversity of A. religiosa in the LMPN, as it has occurred in other conifers species (e.g. Mendoza-Maya et al 2015). Assisted migration may consist of two steps: (1) migration between populations with the least genetic differentiation, i.e., between the high-altitude populations of the East and North slopes and the populations of the South slope; and (2) migration between populations with the greatest genetic differentiation, that is, between low-altitude populations inhabiting the East and North slopes and populations of the South slope. In addition, it is essential to obtain information from other temperate forests in the region (e.g., Popocatépetl, Iztaccíhuatl, Mount Tlaloc) to set out conservation strategies for the oyamel fir at a regional level.

It is also essential to immediately halt land-use changes throughout the LMNP, especially on the South slope, because pollen scarcity in small, isolated populations increase self-pollination and inbreeding (Morales-Velázquez et al., 2010; Ortiz-Bibian et al., 2019). Otherwise, habitat loss and population decline of tree species in LMNP (Ern, 1976; López-Domínguez and Acosta, 2005; Valdez-Pérez et al., 2016) may adversely affect recruitment and genetic structure of recent populations because these conditions enhance genetic drift (Waples, 2016).

This study reports a spatial genetic structure derived from isolation by resistance to elevation and two genetic groups that show a clear separation of the populations inhabiting the South slope, which may indicate a greater effect of genetic drift on these. The insignificant effect of IBR to land-use changes on the genetic diversity of A. religiosa is possibly due to the late effect of fragmentation on the genetic diversity of long-live species or anemochoric dispersal that can preserve pollen movement despite high fragmentation. However, it is imperative to limit changes in land use to conserve the habitat required by this species and thus ensure the establishment of future populations. In addition, in order to define conservation strategies focused on the conservation and promotion of connectivity, it is necessary to include populations data from other temperate forests of the TVB. The information reported herein is highly relevant for the conservation of A. religiosa, the associated fauna, and the ecosystem services provided by the temperate forest of central Mexico.