Conservation genetics of relict tropical species of Magnolia (section Macrophylla)

Special conservation efforts should be made for relict species, as they usually have small population sizes and restricted distributions, placing them in critical extinction risk. To achieve conservation, information about genetic diversity distribution is needed. Here, using nine nuclear microsatellites, we analyzed 23 populations of five recently described species of Magnolia distributed in Mexico, which were previously assigned to Magnolia dealbata. We aimed to determine the level of genetic diversity and the distribution of genetic variation and proposed conservation measures. Compared to other endemic and relict species, we found a moderate level of genetic diversity in most populations; however, we identified two populations with no genetic variation. Additionally, we found evidence of positive values of inbreeding likely due to geitonogamy. We found a strong population structure, low effective population size, and no evidence of bottlenecks. Patterns of genetic differentiation did not support the morphological distinction of five species, so we hypothesized that the gene pools may instead represent well-differentiated populations of a single species. We argue that the pattern of genetic differentiation is explained by the natural fragmentation of the cloud forests after glaciation events, and the effects of genetic drift in small populations poorly connected by gene flow. Despite the moderate levels of genetic diversity, special attention is needed to guarantee conservation, with emphasis on the populations in the central region of the country as well as the valuable populations identified in the southwestern region.


Introduction
The fact that relict species usually represent descendants of once widespread taxa (or populations) suggests that current populations are small in size and occupy disjunct and restricted distributions, placing them in an elevated extinction risk (Habel et al. 2010). Thus, in the face of global climate change and anthropogenic perturbations, special conservation efforts should be made for relict species. Knowledge of the genetic variation within and between populations is central to conservation and management strategies (Blouin et al. 2010;Gitzendanner et al. 2012). For natural selection to produce an adaptive evolutionary change, genetic diversity is needed; therefore, lower genetic variation may constrain the ability of populations to adapt to environmental changes (Frankham 2005(Frankham , 2012, and the maintenance of genetic diversity is essential to guarantee the long-term survival of species (Yamamoto et al. 2017;Lecocq et al. 2018). Moreover, spatial isolation, particularly habitat fragmentation, may result in population subdivision of many species, with low or no gene flow among populations, leading to an increase in inbreeding and genetic drift and directly impacting genetic diversity (de Ita et al. 2012;Nunziata et al. 2016).
The intriguing disjunct distribution of plants in the Northern Hemisphere have been of great interest for evolutionary studies. Such disjunction is thought to be relict of a more continuously distributed mesophytic forest during the Cenozoic (66 Ma), known as the Boreotropical flora, which became fragmented due to geologic and climatic changes (Azuma et al. 2001;Nie et al. 2008). One of the former Boreotropical element is the family Magnoliaceae, which subfamily divergences date back to the Cretaceous (ca. 100 my) (Azuma et al. 2001). It is thought that in response to climate changes, Magnolia migrated from North America to Europe and Asia in the early Eocene (Hebda and Irving 2004). Subsequently, severe cold events during this epoch caused tropical species to migrate to lower latitudes, resulting in the origin of disjunct distributions of tropical and temperate taxa (reviewed by Sánchez-Velásquez and Pineda-López et al. 2016). Although seed fossil records of Magnolia species suggested that they were abundant million years ago (Azuma et al. 2001), current populations remain in patchy and severely restricted distributional areas. Magnolia species conserve some primitive characteristics such as insect pollination interactions. Their flowers are highly specialized to be exclusively pollinated by beetles (Coleoptera); however, flies (Diptera), bees (Hymenoptera), and bumblebees (Bombus) are important pollinators to some species of Magnolia (Thien 1996(Thien ,1974. The flowers of Magnoliaceae are protogynous, although self-compatibility is common (Thien 1996). The fleshy and red-colored seeds of Magnolia are bird dispersed. Members of the genus produce secondary metabolites (e.g. alkaloids and flavonoids), which appear to be involved in defense against natural enemies playing an important role in the evolution of Magnoliaceae (Thien 1996).
Of the 352 extant species of magnolias worldwide, 40 species are distributed in Mexico (Vázquez-García et al. 2017, 2019, most of which are considered threatened or endangered because of human activities such as logging, agriculture and livestock production (Rivers et al. 2016). In addition, almost 80% of these species are important components of a highly threatened ecosystem, the cloud forest (Rodríguez-Ramírez and Luna-Vega 2020). Within the Magnolia genus, the section Macrophylla includes seven species characterized by large deciduous leaves with a glaucous underside. One representative species is Magnolia dealbata. For many years, researchers believed that M. dealbata has a wide distribution characterized by geographically isolated populations. However, some authors have argued that M. dealbata represents a case of postglacial latitudinal adaptive radiation from one ancestral species with a wide distribution in the United States and Mexico (Vázquez-García et al. 2015 and references therein). Consequently, Vázquez-García et al. (2013 and García-Morales et al. (2017) (Vázquez-García et al. 2013García-Morales et al. 2017;Gutiérrez-Lozano et al. 2020).
Population genetic studies of Magnolia, focused mainly on Asian and North American species, showed opposing genetic diversity patterns, with values ranging from high to low (Qiu and Parks 1994;Kikuchi and Isagi 2002;Setsuko et al. 2004Setsuko et al. , 2005Setsuko et al. , 2007Newton et al. 2008;Tamaki et al. 2008;He et al. 2009;Zhang et al. 2010;Azuma et al. 2011;Yu et al. 2011;Budd et al. 2015;von Kohn et al. 2018;Rico and Gutiérrez Becerril 2019;Veltjen et al. 2019). Because limited information is available regarding the genetic diversity, and distribution of genetic variation in Magnolia in Mexico (a priority taxon for conservation), this study focused on five of the seven species conforming to the Macrophylla section in Mexico: M. nuevoleonensis, M. alejandrae, M. rzedowskiana, M. vovidesii, and M. dealbata. These species are distributed across a latitudinal range and are found from 1400 to 2274 m in elevation and they have restricted distributions. Magnolia vovidesii, M. nuevoleonensis, and M. rzedowskiana are considered threatened species, and M. dealbata has been assessed as endangered (Rivers et al. 2016). Here, using microsatellite markers, we aimed to: (i) assess whether genetic patterns support the morphological distinctions of these five Magnolia species, (ii) characterize the level of intraspecific genetic diversity of these primary species, (iii) assess their genetic structure and the genetic differentiation between species and among populations, and in the light of these findings, and (iv) discuss and propose adequate conservation measures.

Sample collection and species identification
Samples consisted of foliar tissue of 230 adult trees from 23 localities (n = 10 per locality) sampled across the Macrophylla distribution in Mexico, including the Sierra Madre Oriental (SMO), Trans-Mexican Volcanic Belt (TVB), and Sierra Madre del Sur (SMS) ( Fig. 1a; Table 1). Sample size was a consequence of the scarce number of individuals across species distribution. To collect samples of Magnolia, we conducted extensive field surveys during 2014 and 2015 based on the geographical location of populations reported in previous works (Vega et al. 2000;Velazco-Macías et al. 2008;Vázquez-García et al. 2013, and herbarium records (Institute of Ecology herbarium (XAL), Xalapa, Veracruz, Mexico; and the National Herbarium of Mexico (MEXU), UNAM, Mexico City). We used the flower size, the number of carpels and stamens, the fruit shape and the geographical distribution to identify Magnolia species within the section Macrophylla (Table S1). Geographical coordinates were recorded for each locality, and the collected tissue was preserved in silica gel until DNA extraction.

DNA extraction and microsatellite genotyping
DNA was isolated following the CTAB extraction protocol (Doyle and Doyle 1987). DNA quality and concentration were determined using a NanoDrop spectrophotometer (Life Technologies). We used 16 nuclear microsatellites to genotype the samples (see Online Resource 1 for details; Tables S2, S3). For each primer pair, the forward primer was labeled with one of the following fluorescent labels: PET, NED, VIC, and 6-FAM (Applied Biosystems). Each PCR contained 2 μM of each primer, 2X of Multiplex PCR master mix buffer (Qiagen), and 5 ng/μL of genomic DNA. We determined five groups of primers for the multiplex PCRs (Table S3); the groups of primers were amplified with the following PCR program: an initial activation of 5 min at 94 °C, followed by 25 cycles of 94 °C for 1 min, annealing temperature (Table S3) for 1 min, and 72 °C for 1 min; and a final extension at 72 °C for 10 min. PCRs were performed using a thermal cycler from Applied Biosystems; we included a negative control to check for contamination. The microsatellite fragment lengths were determined using capillary electrophoresis with an ABI Prism 3100 Genetic Analyzer (Applied Biosystems) with a GeneScan-600 LIZ Size Standard (Applied Biosystems); we recorded fragment size with the program Peak Scanner V2.0 (Applied Biosystems). A positive control was used to correctly assign alleles, while genotypification was conducted according to Selkoe and Toonen (2006).

Data quality
We infer scoring errors and null alleles with Microchecker V2.2.3 (Van Oosterhout et al. 2004) using 10 4 randomizations and a confidence interval of 95%. We tested for deviations from the Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) with Arlequin V3.5.2.2 (Excoffier and Lischer 2010). We performed an exact test of HWE for each locus and each locality (hereafter populations), using 10 6 steps in a Markov chain and 10 5 dememorization steps. The LD between all pairs of loci in each population was evaluated with 10 4 permutations and five initial conditions in an expectation-maximation (EM) algorithm. We performed the analyses with a significance level of 0.05. The final data set was used for assessing the genetic diversity and genetic structure.

Genetic diversity
Parameters of genetic diversity were estimated for populations and individual species. The genetic diversity was assessed as the number of polymorphic loci (P), the total number of alleles (Na), the observed and expected heterozygosity (Ho and He, respectively), and the allelic richness (AR). For AR, the rarefaction method (n = 20 alleles) was applied for individual species. P, Ho, He, and inbreeding coefficients (F IS ) were estimated with Arlequin. Na and AR were calculated using the 'hierfstat' package (Goudet and Jombart 2015) for R V3.5.3 (R Core Team 2019).

Distribution of genetic diversity
The genetic structure was assessed using different methods. First, we evaluated if Magnolia populations adjust to an isolation by distance (IBD) pattern. To do this, we performed a Mantel test to estimate the correlation between Da genetic (Nei et al. 1983) and geographic (km) distances based on 10 4 permutations using the 'ade4' package for R environment. To infer the correspondence between species and populations, the program Structure V2.3.4 (Pritchard et al. 2000) and TESS3 (Caye et al. 2016) was employed to assign individuals to distinct genetic groups (K). For Structure, we used the admixture model with correlated allele frequencies, and all runs consisted of 10 6 Markov chain Monte Carlo simulations with a burn-in period of 10 5 . We assumed K values ranging from 1 to 15, with 10 replicates for each K value. We considered that maximum value of K = 15, can retrieve enough information about the different levels of organization present in the genetic structure. The number of clusters present in the data was assessed based on the mean posterior probability of the data (ln P(D)) (Pritchard et al. 2000) and ΔK (Evanno et al. 2005). The ΔK method generally identifies the uppermost hierarchical level of the genetic partition (Evanno et al. 2005). For ln P(D) if there is no clear maximum likelihood, the point at which the plot curvature plateaus or continue to increase slightly at larger values of K is suggestive of the optimal K (Pritchard et al. 2003).
Values of ln P(D) and ΔK were obtained from Structure Harvester V0.6.94 (Earl and VonHoldt 2012). To identify additional nested groups (fine scale population structure), we performed additional analyses on data subset obtained from the initial analysis: for cluster A (localities from TR to CU (excluding AH); n = 15 localities, 150 individuals) and cluster B (localities from XI to SJV3 plus AH; n = 8 localities, 80 individuals). Run conditions were identical to the initial analysis but testing with different ranges of K due to variation in number of localities; for cluster A, we tested K from 1 to 15 and for cluster B, K from 1 to 10. The individual Q-matrix was loaded in Structure Plot V2.0 (Ramasamy et al. 2014) to draw bar plots of the inferred K clusters.
To consider the geographic coordinates information of populations, we implemented a spatial clustering method. We used TESS3 by using the package 'tess3r' (Caye et al. 2016) for the R environment. For each value of K, which ranged from 1 to 15, we performed 10 independent runs retaining only those repetitions with the lowest root mean squared errors. We employed a cross-validation method to determine the optimum value of K. The best choice of K is when the curve exhibits a plateau or starts increasing (Caye et al. 2016).
To complement previous analyses, we performed a discriminant analysis of principal components (DAPC) (Jombart et al. 2010) using the 'adegenet' package (Jombart 2008) for the R environment. This multivariate method identifies and describes clusters of genetically related individuals, maximizing the genetic variation between groups (Jombart et al. 2010). After cross-validation, we retained 20 and the first three axes in the principal component and discriminant analyses, respectively.
Second, to better understand the patterns of the observed variation, we used a hierarchical analysis of molecular variance (AMOVA) implemented in Arlequin with 10 4 iterations to compute statistical significance. Different grouping hypotheses related to Structure, TESS K-groups, and the five distinct morphological species were assessed using an infinite allele model (IAM) and a stepwise mutation model (SMM).
Third, to infer the relationships between populations of different but closely related species, an unrooted neighborjoining tree based on Da genetic distance (Nei et al. 1983) was constructed using Poptree2 software (Takezaki et al. 2010). Bootstrapping with 10 4 replications was carried out to evaluate the strength of the evidence for the branching patterns. The neighbor-joining tree was visualized in The Interactive Tree of Life (Letunic and Bork 2019) (https :// itol.embl.de).
Finally, we used the program Arlequin to estimate differentiation between individual species and among populations through F ST and R ST . We performed 10 4 permutations to compute the significance level (P < 0.05).

Bottlenecks and effective population size
Populations that have experienced past reductions in size exhibit a decrease in the number of alleles and heterozygosity in polymorphic loci. We used Bottleneck V1.2.02 (Cornuet and Luikart 1996) with 10 4 iterations and default settings to examine whether changes in population size influence genetic patterns. To avoid bias related to population substructure, we performed a sign test and Wilcoxon range test to assess the excess of heterozygosity in each population cluster identified by TESS3 using the two-phase mutation (TPM), SMM and IAM models. In addition, a mode-shift test was performed when in at least one mutation model of both tests resulted in P < 0.05. Mode-shift tests determine whether the allele frequency distributions have an L-shape, which is expected under mutation-drift equilibrium (Luikart et al. 1998). Although we performed bottleneck test based on TESS K-groups, we focused on populations clusters with at least 20 individuals to achieve reasonable power in the inference (Cornuet and Luikart 1996).
A fundamental parameter in conservation biology is the effective population size (Ne). We estimated Ne by using linkage disequilibrium information (Waples and Do 2008) as implemented in NeEstimator V2.1 (Do et al. 2014). We chose a minimum allele frequency cutoff value of 0.05 to ensure that the results were not driven by the presence of rare alleles in the data. To avoid bias associated with population substructure or gene flow in Ne estimation, we estimated Ne for genetic clusters identified at a finer scale with the TESS3 program. The confidence intervals were obtained from parametric CI.

Data quality
No evidence for scoring errors due to stuttering or allele dropout in any of the loci was detected by Microchecker software. However, null alleles were present in seven of the 16 loci amplified, with frequencies between 0.17 and 0.33 (Table S4). We identified that one locus showed HWE deviations in most populations sampled (Table S5). LD pattern did not suggest a systematic bias across pairs of loci (Table S6). Finally, we removed six loci with allele frequencies higher than 0.17 and one additional locus that showed deviations from HWE across the populations. The genetic diversity and genetic structure analyses were assessed using a final dataset of nine loci (Table S3; see also Online Resource 2).

Genetic diversity
The estimates of the genetic diversity for five species in the section Macrophylla among 23 populations are shown in Table 2. At species level, we observed a notable discrepancy between the Ho and He and relatively high levels of inbreeding (F IS ) among all the species (P < 0.05). Magnolia rzedowskiana, M. vovidesii, and M. dealbata showed the highest values of allelic richness (5.6, 6, and 5.6, respectively).
Within the populations, we observed contrasting values of genetic diversity, ranging from very low (Ho = 0; He = 0) to moderately high (Ho = 0.55, He = 0.7); nonetheless, many of the populations had moderate (> 0.30) levels of heterozygosity (Table 2). We obtained a total of 177 different alleles, with a mean of 19.7 alleles per locus. Most of the populations showed polymorphic loci, although two populations of M. rzedowskiana (P1 and P2) showed no polymorphisms ( Table 2). The highest number of alleles

Distribution of genetic diversity
We detected strong positive correlation between genetic and geographic distances (r = 0.498; P < 0.001) for Magnolia populations (Fig. S1). For the initial analysis using Structure, we excluded a total of 10 replicates with substantial differences between replicates before run Structure Harvester (2 replicates for K = 5, 8 and 14; 1 replicate for K = 6, 7, 12, and 15). Whilst ΔK identified K = 2 as the best fit, we observed that ln P(D) increased gradually as increasing K, probably due to isolation by distance pattern (Fig.  S2a). Thus, we focus on the uppermost hierarchical level of genetic partition ( Fig. 1b). For the subsequent analysis of Cluster A, we excluded 11 replicates before run Structure Harvester (2 replicates for K = 13 and 15; 1 replicate for K = 4, 5, 6, 7, 9, 11, and 12). The inspection of ln P(D) and ΔK (Fig. S2b) suggested K = 5 as the optimum number of the genetic partition (Fig. 1c). Despite the higher value of ΔK (3.7) showed K = 2, ln P(D) increase slightly after K = 5 and also showed similar ΔK value (3.6) (Fig. S2b). Value of ΔK in subsequent analysis of the cluster B suggested K = 3 as the best fit, whereas ln P(D) showed a plateau at the same K value (Fig. S2c); K = 3 is shown in Fig. 1c. For this data subset, it was not necessary to exclude replicates since all of them showed similar variation between runs. For all Structure analyses, the estimated membership coefficient for most individuals was high for the inferred clusters (Q > 0.9). The cross-validation score from TESS3 analysis showed a continuum decrease until K = 13 (Fig.  S2d). However, higher number of groups than K = 11 showed no clear subdivision and increased variation; thus, we focus in a finer-scale population structure at K = 11 (Fig. 1d). Population groups from Structure and TESS3 were highly corresponding to the geographical distribution of the populations ( Fig. 1). However, for both approaches we observed different genetic groups despite their geographic proximity (Fig. 1a, c, d). Although the strong genetic differentiation, Structure detected some admixed individuals within the CH, TE, XI, SJV1, and SJV2 populations (Fig. 1c).
The AMOVA analyses indicated that under IAM model, most of the proportion of the genetic variation was found within individuals for most partitioning hypotheses, while for SMM model was among groups (Table 3). For both mutation models, we observed overall a trend in which the percentage of variation explained among groups increase as the number of genetic groups (from 14.74 to 30.63% for IAM, and from 0.59 to 78.96% for SMM). However, relative to other clustering hypotheses, the genetic divergence was maximized among groups (F CT = 0.789; P < 0.001) and minimized within groups (F SC = 0.210; P < 0.001) when TESS clustering (K = 11 under SMM model) was considered (Table 3).
Based on the neighbor-joining tree, we identified four major clades; however, only two of them showed statistical robustness higher than 60% (Fig. 2), and those clades included populations of M. dealbata, M. alejandrae and M. nuevoleonensis. The neighbor-joining tree corroborated the finer scale structure pattern found in a previous analysis (Fig. 1d), in which many populations, such as AX, CH, JCC, and AH, diverged into separate long branches (Fig. 2). The neighbor-joining tree also suggested that the northern populations, which consist of M. nuevoleonensis (TR) and M. alejandrae (SPE and PE), are closely related to the southwestern populations (M. dealbata). In turn, long branches of M. nuevoleonensis and M. alejandrae shared ancestry. We observed that one population, identified as M. rzedowskiana (P3), was more related to populations of M. vovidesii (Fig. 2), but there was low support for this relationship; only three populations of the same species (P1, P2, and P4) were nested within the same clade (Fig. 2).
The mean pairwise F ST (0.475) and R ST (0.582) showed high levels of differentiation among all the populations; the values ranged from 0.047 to 0.95 (P < 0.05) and from 0.064 to 0.99 (P < 0.05), respectively (Fig. 3). Although strong differentiation was found, we observed that the populations P2, P1, P4, P3, AX, and TR were more divergent based on the F ST values. In addition, the northern populations (TR, SPE, 1 3 and PE) and other populations geographically close to AX (CO, AT, and MO) were considerably more differentiated than the rest of the populations based on the R ST (Fig. 3). We determined that M. nuevoleonensis and M. alejandrae were the most divergent species based on both the F ST and R ST ; values ranged from 0.21 to 0.43 and from 0.27 to 0.96, respectively (Table S7).

Bottlenecks and effective population size
Overall, we did not find evidence of deviations from the mutation-drift equilibrium in Magnolia populations. We performed mode-shift tests in three population groups (cluster 3 (JCC), cluster 6 (AH), and cluster 9 (TE, MO, and CU)), which showed significant bottleneck mainly under the IAM mutation model (Table S8). However, the plot of the modeshift test confirmed normal L-shape of allele frequency distributions for cluster 9 but deviations for clusters 3 and 6 ( Fig. S5). None of the population cluster with at least 20 individuals showed deviations from the mutation-drift equilibrium (Table S8).
The estimates of the Ne revealed a very small population size across the 11 genetic groups of populations identified ( Table 4). The values of Ne ranged from 0.3 to 274.3; more than half of the populations had a Ne no more than 15 individuals, three had a Ne between 26.7 and 61.3, and only one group of population had 274.3 individuals of Ne (Table 4).

Discussion
For conservation biology and other disciplines, the delimitation of species as an operational taxonomic unit is fundamental (e.g., Duminil and Di Michele 2009). For plant species recognition, morphological characters have historically been used to group individuals into species. However, the morphological characters can fail to differentiate between cryptic species or incorrectly subdivide species because of a misinterpretation of natural morphological variation across their distribution range (Duminil and Di Michele 2009;Duminil et al. 2012). Since there is relative morphological homogeneity within the Magnoliaceae family as well as the highly conserved chloroplast DNA (Corneanu et al. 2004;Kim et al. 2004;Nie et al. 2008), higher mutation rates become microsatellites a useful and efficient tool for distinguishing between different related species (Duminil et al. 2012). High genetic differentiation between populations has been recorded in other Magnolia (Kikuchi and Isagi 2002), endangered and relict (Szczecińska et al. 2016;Vardareli et al. 2019;Wu et al. 2019) species. This result has been explained by the limited gene flow either currently or historically and the strong effect of genetic drift in small isolated populations. The effects of climatic changes in the recent past remain controversial; however, evidence for Mexico suggests that tropical mountains facilitated persistence in fragmented distributions and prompted differentiation because of different patterns of connectivity across distribution ranges Fig. 2 Neighbor-joining tree based on Nei's genetic distance was colored according to different species of Magnolia such as in Fig. 1. We only showed bootstrap values of 60% or higher; light red (bootstrap > 60%) and gray lines indicate major lineages identified 1 3 during the Pleistocene (Mastretta-Yanes et al. 2018). Our genetic data showed small effective population size of M. dealbata that were characterized by high levels of differentiation among all the populations, with several genetically distinct groups (K = 11) even over short distances (< 25 km). It is possible that during glacial periods in some regions of the species distribution, gene flow may occur more readily because of habitat continuity (e.g., in the YA, FA DU. SJV1, SJV2, and SJV3 areas), while others have remained isolated (e.g., in the TR, SPE, and PE areas). This scenario could be exacerbated with the continuous natural fragmentation since the last glacial maximum (Newton et al. 2008). Although repeated glaciation and bottlenecks were invoked to explain genetic structure in other tertiary relict species, (e.g. Glyptostrobus pensilis) (Wu et al. 2019), we failed to detect signs of recent declines in population size. However, it is possible that bottlenecks remain undetected due to time-lag, sample size or ancient high levels of genetic diversity previous to the bottlenecks (Perry et al. 2012). We think that the pattern of genetic structuring could be explained by the action of genetic drift on geographically restricted dispersal. The seeds of Magnolia are bird-dispersed (although dispersal by rodents and insects has also been recorded) and the flowers are pollinated mainly by beetles (Dieringer and Espinosa 1994;Thien et al. 1996;Hirayama et al. 2005;Zhang et al. 2010); thus, gene flow could be expected between the sampled populations. However, as shown by the genetic structure analyses, our results suggest that populations are poorly connected by gene flow; under this scenario, genetic drift can differentiate allele frequencies faster than dispersal can homogenize geographical distant populations. This result is not surprising since studies have found a tendency for the movement of beetles and birds to be sensitive to habitat discontinuity and their responses vary according to the degree of specialization (Devictor et al. 2008;Kobayashi and Sota 2019).
The interplay between evolutionary forces and life-history traits, such as mating system, longevity, and vectors of seed and pollen dispersal, plays an important role in determining the levels of genetic variation within populations (Hamrick et al. 1992;Hamrick and Godt 1996). The patterns of genetic   (Rico and Gutiérrez-Becerril 2019) and Canada (Budd et al. 2015). Moreover, the genetic relationship between southern and northern populations identified in the neighbor-joining tree could represent the signature of the once widespread ancestral population. If this scenario is correct, as was found in M. obovata, the natural and active interplant movement of beetles, which transport genetically diverse pollen, favors cross-pollination (Matsuki et al. 2008) and thus higher levels of genetic diversity in the ancestral population.
Although Magnolia species have flowers that encourage cross-pollination, there is evidence of inbreeding likely through geitonogamy for M. stellata and M. sieboldii (Kikuchi and Isagi 2002;Tamaki et al. 2008). Magnolia species produce protogynous flowers (Dieringer and Espinosa 1994;Gutiérrez and Vovides 1997); however, the asynchrony of flower anthesis results in the coexistence of female and male flowers that could result in geitonogamy (Hirayama et al. 2005). We found overall high levels of inbreeding in nine populations, which is also explained by their low population size (< 100 individuals), making it more likely that close related individuals will mate. On the other hand, in two populations (P1 and P2), the patterns of genetic variation reflected predictions of strong drift and inbreeding, and the frequencies of the heterozygotes are close to zero, causing a fixation of alleles across loci (Ellstrand and Elam 1993). Likewise, the combination of isolation, inbreeding, genetic drift, and limited gene flow between populations explains the substantially higher values of F ST and R ST found in some populations.

Conservation remarks
The Mexican cloud forest is particularly vulnerable to global climate change due to its fragmented nature, specific environmental conditions, and deforestation rate (Ponce-Reyes et al. 2013 and references therein). After mangroves, cloud forests represent the ecosystem with the lowest remaining surface area in Mexico (SEMARNAT, 2016). From the 17 274 km 2 of cloud forests currently found in Mexico, models predicted that only 32% will persist by 2080 due to climate change (Ponce-Reyes et al. 2012). Moreover, of the current total extent, only approximately 12% occurs in federally protected areas (Ponce-Reyes et al. 2012). Among the populations studied here, only seven are under federal protection: La Trinidad (M. nuevoleonensis) in Cumbres de Monterrey National Park and populations assigned as M. rzedowskiana (all but the Chilijapa location) in the Sierra Gorda Biosphere Reserve. However, according to predictions, by 2080, less than 1% of these protected areas are likely to persist (Ponce-Reyes et al. 2012). This highlights the need to increase the number of protected areas across the cloud forests in Mexico based on the existence of suitable habitat predicted by 2080 (Ponce-Reyes et al. 2012) and the importance of interconnecting habitats more effectively, because sites such as P1 and P2, which have heterozygotes frequencies close to zero and fixed alleles, without gene flow are highly susceptible to extinction.
Interestingly, the cloud forests in southwestern Mexico (Oaxaca state) is predicted to remain stable, without significant changes in its total area, but the distances between fragments are predicted to vary from 2.4 to 2.6 km (Ponce-Reyes et al. 2013), which is lower than the maximum pollen dispersal distance estimated for M. stellata (Setsuko et al. 2013). These populations are not under federal protection, but they are found in community reserves, in which the wood and flower extraction is prohibited; thus, this region is valuable in the face of global climate change, and federal protection is needed to maintain population viability. In contrast, in the eastern part of the country, the state of Veracruz has conserved less than 30% of its natural vegetation (SEMAR-NAT 2016), and forests are surrounded by coffee plantations or grazing areas. Although most populations studied in the region (consisting of M. vovidesii) have moderate genetic diversity, they are under special concern because inbreeding has been detected, and because of the susceptibility of further generations. Incorporating small patches of forest to allow the movement of pollen and seed vectors between fragments is necessary to conserve populations in this region as well as in nonprotected sites in the north (SPE and PE).
Conservation biology requires species lists as an accurate measure of different attributes of biodiversity, such as endemism and richness, comparable across taxa and localities (Isaac et al. 2004;Zachos et al. 2013). The increase of species number has important impacts on global targets for conservation because can mask the risk and rates of extinction, but also result in a higher proportion of threatened or extinct species (Isaac et al. 2004). Particularly, a large number of new species have been described in Latin America for Magnoliaceae, becoming the second center of species diversity for the family (Rivers et al. 2016). However, oversplitting of species suppose costly consequences for species conservation, because conservation resources are finite and resources allocated to one species do not necessarily benefit others (Stanton et al. 2019), and can lead to management actions at the wrong spatial scale and thus, incorrect management strategies.

3
In conclusion, we argue that the extant gene pools represent the legacy of postglacial fragmentation of an ancient population; nonetheless, other evolutionary processes can be invoked to better explain the observed genetic patterns. Nevertheless, a revision using morphological and other molecular markers (nuclear and chloroplast genetic sequences) from across the distribution range can help to elucidate the number of existing taxonomic units. Our results contributed to the understanding of the genetic structure of Mexican magnolias, highlighting a moderate genetic diversity with a strong genetic structure. We argue that despite the levels of genetic diversity recorded here, a greater effort should be made to guarantee conservation, with emphasis on those populations within highly fragmented habitats (e.g., Veracruz region), and for genetically eroded populations (e.g., P1 and P2). Similarly, populations identified in the Oaxaca region are valuable due to their predicted area stability in the near future.