Magnoliaceae Juss. is an old family in the evolutionary history of flowering plants and displays an important Holarctic fossil record of over 100 million years (APG IV 2016, Romanov & Dilcher 2013), estimated to have originated between 113.15 - 104.71 Mya (Magallón et al. 2015). It includes ca. 350 species of trees and shrubs within two subfamilies: Liriodendroideae (containing only Liriodendron L. with two species) and Magnolioideae (comprising only Magnolia L.) (Figlar & Nooteboom 2004, Figlar 2006, Kim & Suh 2013, Vázquez-García et al. 2016). Magnolias are currently found in temperate and tropical areas of Asia and the Americas (APG IV 2016). About half of the species occur in the New World, from temperate eastern North America through Mexico, Central America, the Caribbean islands and northern South America to Bolivia and Brazil (Lozano-Contreras 1994, Vázquez-García 1994, Frodin & Govaerts 1996, Vázquez-García et al. 2016). Mexico and Colombia, with nearly three dozens of species each, are the two magnolia richest countries in the continent (Lozano-Contreras 1994, Vázquez-García 1994, Vázquez-García et al. 2017). Their showy and fragrant flowers make magnolias of great horticultural and ethnobotanical value throughout the world.
The Magnolia pacifica species complex, subg. Magnolia, sect. Magnolia (Figlar & Nooteboom 2004) consists of four morphologically close and currently recognized endemic species of Mexico (Vázquez-García 1994, Vázquez-García et al. 2002, 2012, 2013): Magnolia tarahumara (A.Vázquez) A.Vázquez, Magnolia pacifica A.Vázquez, Magnolia pugana (Iltis & A.Vázquez) A.Vázquez & Carvajal, and Magnolia vallartensis A.Vázquez & Muñiz-Castro. M. tarahumara inhabits northwestern Mexico (Vázquez-García et al. 2013) whereas M. pacifica, M. pugana, and M. vallartensis belong to the “Pacific southwest” group sensuVázquez-García (1994). Magnolia species in the Pacific southwest region are distributed along a continentality and moisture gradient because of their different distances to the Pacific Ocean, and only survive in permanent humid microclimates such as very moist ravines (on the coastal mountains), or adjacent to springs and permanent streams (in the interior seasonally-dry valleys and ravines). The lands of the interior of the continents have a different climate than the lands near the coasts because they do not receive (or receive in a smaller proportion) the moderating effect of the humidity coming from the sea. The interior, unlike the coasts, presents a drier and a more variable and extreme climate, with a greater range of variation of diurnal and annual temperatures (Snow 2005).
Most reports within the M. pacifica complex are limited to the description of species based on qualitative morphological and geographical information as well as ecological differences. In particular, M. pugana was first treated as a subspecies of M. pacifica (Vázquez-García 1994); afterward, it was described as a species because of its distinctive characters such as the numbers of carpels per fruit and stamens per flower, glabrous peduncles, and tepals shape, which allowed to clearly separate it from M. pacifica (Vázquez-García et al. 2002). Afterward, M. vallartensis was also described as a new species based on the form and size of its leaves, number of carpels per fruit, number and shape of tepals, and height of the trees (Vázquez-García et al. 2012). Since phenotypic plasticity is one of the major traits by which plants can cope with the variability of environmental factors (Gratani 2014), morphological differences between species of the M. pacifica complex might be the result of ecological differences. Thus, in order to get a more reliable species delimitation within this complex, morphological boundaries must be validated with molecular data (Schlick-Steiner et al. 2010).
Population genetic analyses are considered complementary and useful tools for species delimitation in complex groups (Sites & Marshall 2004). The population approach for species delimitation is based on the unified, general species concept proposed by de Queiroz (2005, 2007). This concept describes species as population-level evolutionary lineages. Under this notion, species generally show a higher genetic divergence from one another than populations within species (Drummond & Hamilton 2007). Nevertheless, the levels of differentiation among species depend on the time of separation and amount of gene exchange (Hey & Pinho 2012). In trees, factors such as long generation time, large effective population sizes (Rosenberg 2003), and frequent introgression, increase the opportunity of sharing polymorphisms (Zhou et al. 2017), which makes the definition of species even more problematic. Thus, delimitation of recently diverged tree species based on DNA markers is particularly challenging, and sequence data (nDNA, mtDNA, cpDNA) used in traditional phylogenetic analysis may not be informative enough. In particular, mitochondrial DNA is generally insufficiently polymorphic in plants, and incongruences between phylogenies obtained from nDNA and cpDNA seem to be frequent (Wang et al. 2014).
The population genetic approach is widely used for species identification and evolutionary analysis of closely related plant species (Medrano et al. 2014, Stiehl-Alves et al. 2017, Esfandani-Bozchaloyi et al. 2018, Sheidai et al. 2018, Pinangé et al. 2019). For this purpose, it is recommended to use highly polymorphic markers with large genome coverage (Duminil & Di Michele 2009) which are not influenced by environment nor affected by natural selection (Holderegger et al. 2006), but experience high rate of intraspecific gene flow (Petit & Excoffier 2009). Dominant markers and microsatellites (a codominant marker) have demonstrated to be useful for species delimitation, however, dominant markers have a higher genome coverage compared with microsatellites, and have been efficient to provide information on species boundaries (Hausdorf & Hennig 2010). Inter simple sequence repeats (ISSR) are dominant and neutral multi-locus markers amplified by polymerase chain reaction (PCR) in the presence of single primers complementary to simple sequence repeats (SSR) or microsatellites (Zietkiewicz et al. 1994, Reddy et al. 2002). The high evolutionary rate, nuclear origin and a large number of loci scored, make ISSR markers useful to define inter-specific boundaries between closely related plant species (Fernandes-Lima et al. 2015, Rodrigues et al. 2015, Stiehl-Alves et al. 2017). The use of Bayesian model-based clustering methods for the analysis of multi-locus markers has proven to be a useful and recommended approach to identify species boundaries (Petit & Excoffier 2009, Hausdorf & Hennig 2010). These tools can be particularly powerful when combined with other methods, such as phylogenetic analyses (Noble et al. 2010).
Due to the large number of threatened Magnolia species worldwide, the evaluation of its genetic diversity to propose and implement conservation measures has become a major international task (Cires et al. 2013). Out of the 306 Magnolia species registered in the IUCN 2019 Red List of Magnoliaceae, 121 are listed as critically endangered or endangered; in particular, the Neotropical region holds the highest proportion (75 %) of threatened magnolias (Rivers et al. 2016). Based on its restricted geographical distribution (4,732 km2) and a continued decline in the extent of occurrence and quality of its habitat, M. pacifica has been classified as endangered according to Cicuzza et al. (2007) and Rivers et al. (2016). M. pugana has also been categorized as endangered because of its small extent of occurrence (2,460 km2) and habitat degradation (Rivers et al. 2016). M. vallartensis has been considered as critically endangered as its extent of occurrence was estimated to be less than 100 km2 (Rivers et al. 2016). The conservation of genetic diversity is critical for the long-term survival of a species (Schemske et al. 1994), however, information about the levels and distribution of genetic variation of these three taxa is not available. Here we evaluate the genetic structure, differentiation and diversity of the M. pacifica species complex using ISSR markers in order to 1) analyze species boundaries and 2) contribute to the conservation status assessment of each taxon.
Materials and methods
Study species and distribution. The studied taxa of the Pacific southwest M. pacifica complex exhibit differences in morphology, phenology, and habitat along a continentality and moisture gradient (Table 1) (Figures 1 - 3). Plant voucher specimens were collected in 10 localities for reference and morphological revision, and were deposited at the IBUG Herbarium of the Universidad de Guadalajara (Supplementary data 1).
M. pugana is endemic to the center of the state of Jalisco and southern Zacatecas, in the confluence of the Trans-Mexican Volcanic Belt (TMVB) and the southern end of the Sierra Madre Occidental (SMO) physiographic provinces (INEGI 2001) (Figure 3). This taxon, being nearly 170 - 215 km away from the sea, thrives in a seasonally-dry continental climate (AAR ca. 833 - 980 mm), in gallery riparian forests, on the side of permanent streams and springs; surrounded by tropical seasonal dry forests and pine-oak forests, usually exposed to high water stress (Vázquez-García et al. 2002).
Traits | M. pugana | M. pacifica s.s. | M. vallartensis |
---|---|---|---|
Mature tree height (m) | 10-25 (30) | 15-28 (35) | 8-15 (23) |
Mature tree dbh (cm) | 30-80 (160) | 40-80 (160) | 20-35 (47) |
Leaf shape | Narrowly elliptical to elliptical, to lanceolate | Narrowly elliptical to elliptical, to obovate | Elliptic to wide elliptic or oblong |
Leaf apex | Frequently acute | Frequently acute | Frequently obtuse |
Leaf size (cm) | 12-22 × 2.5-8 | 8-17 (18) × (3) 5-6 (8) | 13.5-27.8 × 6-14.8 |
Peduncle | Essentially glabrous to pubescent at the base | Pubescent throughout or at least at the nodes | Internodes glabrous, nodes pubescent |
Sepals shape | Widely obovate | Oblong | Oblong |
Petals (number) | 6-7 | 6 | 6-8 |
Stamens (number) | 92-100 | 100-110 | 75-82 (118) |
Carpels | 25-35 | (18) 23-28 (37) | 10-19 |
Fruit (polyfollicle) shape | Oblongoid | Narrowly oblongoid | Obovoid to narrowly ellipsoid |
Locule shape | Oval to almost round, slightly longer than wide | Elliptic, much longer than wide | Elliptic, much longer than wide |
M. pacifica s.s. is distributed in montane cloud forests of western Jalisco and Nayarit, in the confluence of three mountain physiographic provinces, the TMVB, SMO and Sierra Madre del Sur (SMS). Because of its proximity to the Pacific Ocean (40-70 km from the sea) it is exposed to higher levels of precipitation (AAR ca. 1,100 - 1,364 mm), both vertical (rain) and horizontal (clouds, mist, and fog) (Vázquez-García 1994).
M. vallartensis is found only in the west of Jalisco state, in the municipalities of Puerto Vallarta and Cabo Corrientes, in the northwestern limit of SMS; this species inhabits montane cloud forests and its ecotone with tropical sub-evergreen forests, close to the Pacific Ocean (0.5-20 km from the sea) and receives moisture throughout the year (AAR ca. 1,348-1,591 mm) (Vázquez-García et al. 2012).
Plant material and DNA extraction. This study includes 278 samples from 10 localities: four belong to M. pugana, three to M. pacifica s.s., and three to M. vallartensis (Figure 3). These selected locations cover the geographical distribution range of each taxon. Ten field trips were conducted from February 2012 to August 2015. In each locality 21 to 33 georeferenced (GPS Garmin, 60 CS) adult magnolia trees were randomly sampled; individuals were separated by a minimum distance of 100 m to reduce the probability of sampling from family clusters. From each tree, fresh mature leaves were collected, stored in plastic bags filled with silica gel and taken to the laboratory. Genomic DNA was extracted using the CTAB method reported by Cota-Sánchez et al. (2006) with slight modifications. The DNA concentration was quantified at A260 with a Spectro UV-VIS RS digital spectrophotometer (LaboMed Inc., Los Angeles, CA, USA).
ISSR-PCR. In an initial step, twenty-three anchored and one un-anchored ISSR primers were screened using DNA extracted from five randomly selected samples. Out of the screened primers, six were selected based on their amplification pattern. In order to test for reproducibility of the PCR with the selected primers, duplicate amplifications were performed with fifteen DNA samples chosen at random. Amplifications were carried out in 20 µL reaction volumes containing 75 ng of genomic DNA, 1X PCR buffer, 2.5 mM MgCl2, 0.2 mM of each dNTP (Invitrogen), 1 μM of the primer, and 0.8 U Taq polymerase (Invitrogen) in a PTC-100 thermal cycler (MJ Research, Inc.) programmed for an initial denaturing step of 3 min at 95 °C and 40 cycles of the following temperature profile: 45 s at 95 °C; 45 s at 50 °C (primers 810, 814 and 857) or 56 °C (primer 834) and 60 °C (primers 836 and 855); and 2 min at 72 °C. Cycling was concluded with a final extension at 72 °C for 10 min. Amplification products were separated by electrophoresis on 2 % agarose gels (standard and low melting point agarose at 2:1 ratio) with 1X TBE buffer, under 100 V for approximately 2 h. The gels were stained with ethidium bromide (1 µg/mL) and photographed under UV light with a Kodak photo documentation system (Kodak ID Image Analysis Software, version 3.5). Lengths of the DNA fragments were estimated using a 100 bp DNA ladder (Invitrogen).
Data analysis. The presence or absence of consistently reproducible DNA bands in the stained gels was scored manually. All smeared and weak bands were excluded. A binary data matrix was constructed for each sampling site, with fragments of the same size considered as the same locus. The species boundaries in the M. pacifica complex were examined through a multi-analytic strategy. First, Bayesian clustering analyses were performed with the software STRUCTURE 2.3.4 to infer the probability to assign the genotype data set to a given number of clusters (K) (Pritchard et al. 2000). Admixture proportions (Q) of all localities and individuals were estimated. Fifteen independent runs were conducted for K values between one (panmixia) and 11, using 150,000 MCMC (Markov Chain Monte Carlo) iterations and a burn-in period of 20,000, admitting the admixture model of ancestry and correlated allele frequencies. This analysis was performed in a hierarchical framework to identify the complete substructure in the data set. In a first run, STRUCTURE was executed with the full dataset, and then independent runs were executed for each genetic group revealed by this analysis. The results of all STRUCTURE analyses were clustered and averaged using CLUMPACK (Cluster Markov Packager Across K) (Kopelman et al. 2015). Evannos’s index (Evanno et al. 2005) was used to detect the best K in the dataset. Second, a dendrogram based on Nei’s genetic distance (Nei 1972) was constructed by using TFPGA 1.3 (Miller 1997), implementing an unweighted pair group method of cluster analysis that uses arithmetic averages (UPGMA). A bootstrap analysis using 1,000 replicates was performed to assess support for the inferred groups. Third, in order to identify genetic barriers, Barrier 2.2 software (Manni et al. 2004) was used to correlate geographical and genetic distances between populations using Monmonier's maximum difference algorithm (Monmonier 1973), with significance tested by means of 1,000 bootstrap matrices of genetic distance. Geographical coordinates of populations were used to obtain a Voronoï tessellation where barriers were delineated. Finally, analyses of molecular variance (AMOVA) were performed using GenAlEx 6.5. The distribution of the genetic variance was assessed in 1) morphology-based taxonomic groups and, 2) Bayesian analysis-based groups. AMOVAs were also run for each putative studied species.
To assess the genetic diversity of genetic groups, each taxon, and each locality, the percentage of polymorphic bands (P) (0.95 level) and Shannon’s diversity index (I) were estimated with POPGENE v. 1.31 (Yeh et al. 1999). Nei’s gene diversity (HE), intrapopulation genetic diversity (Hs), and total genetic diversity (HT) were estimated using GenoDive Beta 2.0 (Meirmans & Van Tienderen 2004). GST (Nei 1972) and D (Jost 2008) differentiation indexes were also calculated. In order to test for differentiation between localities, Exact Tests (Raymond & Rousset 1995) for locality pairs were performed with TFPGA 1.3 (Miller 1997), with 1,000 dememorization steps, followed by 20 batches of 2,000 permutations per batch; significance values were obtained following the Fisher’s method that combines probabilities of exact tests (Sokal & Rohlf 1995). Finally, pairwise genetic distances among localities were estimated to test for significant correlation with the corresponding geographical distances, a Mantel test was performed with TFPGA 1.3 (Miller 1997) using 999 random permutations. Because the three studied taxa belong to a complex of closely related species, this analysis was done for all localities grouped as one species complex and for each main genetic cluster.
Conservation status assessment. Extinction risk evaluation of the three studied taxa was performed using the IUCN Red List Criteria (criteria B1 + B2, IUCN 2019) and the GeoCAT software tool (Kew-IUCN-VIBRANT 2019). The records used for the delimitation of the Extent of Occurrence (minimum convex polygon, EOO) and the Area of Occupancy (grid cell area with occupancy, AOO) of each taxon were obtained from IBUG herbarium specimens, the GBIF, Tropicos.org, and REMIB-CONABIO databases. AOO calculation was based on the IUCN default cell width of 2 km. Additionally, the criterion C-2 (Genetics) of the Method of Extinction Risk Evaluation of Plants in Mexico (MER) from the Mexican Official Norm 059 (SEMARNAT 2010) was used to contribute partially to the MER assessment. Criterion C-2 proposes that (1) if species has heterozygosity (He) < 10 - 20 % (depending on the molecular marker used) and (2) a genetic differentiation (Gst or Fst) > 20 %, it has a higher threatened status or extinction risk.
Results
ISSR diversity. Primers 810, 814, 834, 836, 855 and 857 produced the best amplification patterns. Each primer amplified between eight and 16 DNA fragments, ranging in size from 300 to 1,650 bp (Table 2). A total of 76 clear and reproducible amplified DNA fragments were analyzed, of which 41 were polymorphic among all samples tested, with each primer displaying from three to 10 polymorphic fragments. The 900 bp (primer 834) and 850 bp (primer 855) ISSR loci did not amplify M. pugana samples.
Primer code | Sequence (5´-3’) | Annealing temperature (°C) | Fragments | |||
---|---|---|---|---|---|---|
Number | Polymorphic | Exclusive | Size (pb) | |||
UBC810 | (GA)8T | 50 | 16 | 5 | 380-1350 | |
UBC814 | (CT)8A | 56 | 8 | 3 | 630-1200 | |
UBC834 | (AG)8YT | 56 | 13 | 10 | 1 | 410-1650 |
UBC836 | (AG)8YA | 60 | 14 | 9 | 290-1150 | |
UBC855 | (AC)8YG | 60 | 10 | 4 | 1 | 320-1600 |
UBC857 | (AC)8CTC | 50 | 15 | 10 | 320- 870 | |
Total | 76 | 41 | 2 |
Species boundaries and genetic clustering. The Bayesian analysis with the full data set identified a value of K = 2 as the most likely number of genetic clusters (Figure 4). Alpha values were monitored and convergence was reached following the 20,000 iterations in the burn-in period. The analysis revealed one genetic cluster for the four M. pugana localities and another cluster including the six localities of M. pacifica s.s. and M. vallartensis. Thus, the localities of M. pugana can be assigned to one population, and the localities of M. pacifica s.s. and M. vallartensis to a second population. The M. pugana cluster was highly homogeneous internally with almost no uncertain assignments at the localities (≥ 93 %) or individual level (> 97 %); an individual membership coefficient threshold of Q > 0.8 was used to distinguish between pure individuals and individuals of admixed ancestry (Q < 0.8). The M. pacifica s.s.-M. vallartensis group was less homogeneous than the M. pugana group; the localities of BA and CSJ exhibited 78 and 75 % of assignment at this cluster, respectively; the rest of the localities were exclusively made up of this genetic cluster (≥ 93 %); admixed individuals were also more common in this group, as only 88 % of individuals were assigned to this cluster.
In order to reveal further genetic substructure within each of the two main clusters, a second Bayesian analysis was run, and a value of K = 2 was the most likely number of genetic subgroups within each cluster (Figure 4). The M. pugana main group was divided into the ALV-ASL and RV-APV subgroups, whereas the M. pacifica s.s.-M. vallartensis group was subdivided into M. pacifica s.s. and M. vallartensis subgroups. The proportions of individuals assigned to each the M. pugana and the M. pacifica-M. vallartensis subgroups were 70 and 74 %, respectively. Similar to STRUCTURE results, two major clusters were identified in the UPGMA dendrogram (Figure 5), the first one included the four localities of M. pugana, and the second one the six localities of M. pacifica s.s. and M. vallartensis. The bootstrap analysis showed a value of 100 % at this node, indicating strong support for these groups. Within the cluster of M. pugana, the localities ALV and ASL form one subgroup, and RV and APV form another one, however, the node had a very weak bootstrap support value (< 50 %). Within the second main cluster, the localities of M. pacifica s.s. form one subgroup and the localities of M. vallartensis form a second one, with a bootstrap support value of 74 %.
Three geographical boundaries were detected, using Monmonier's algorithm (Figure 3). One barrier separates the two subpopulations of M. pugana. A second geographical barrier was revealed between M. pugana and the genetic group M. pacifica s.s.-M. vallartensis and finally another barrier separates M. pacifica s.s. from M. vallartesis. All geographical boundaries had 100 % bootstrap support.
The AMOVA on the morphology-based taxonomic groups showed that 84 % of the genetic variation was found within localities, and that only 8 % was explained by differences among species, and a similar proportion was due to differences among localities within each taxon. Similar results were revealed by the AMOVA on the Bayesian genetic analysis-based groups (Table 3). All components of molecular variance were significant (p value 0.001).
Percentage of variation | |||||
---|---|---|---|---|---|
Groups | Taxa | ||||
Variation source | Bayesian analysis (2) | Taxonomic (3) | M. pugana | M. pacifica s.s. | M. vallartensis |
Among | |||||
Taxa / Genetic groups | 9 (0.001) | 8 (0.001) | |||
Localities | 9 (0.001) | 8 (0.001) | 12 (0.001) | 7 (0.001) | 7 (0.001) |
Within | |||||
Localities | 82 (0.001) | 84 (0.001) | 88 (0.001) | 93 (0.001) | 93 (0.001) |
(2) two Bayesian groups, (3) three taxa, p value is given in parentheses.
Genetic diversity and differentiation. The percentage of polymorphic loci (P), Shannon’s index (I) and heterozygosity (HT) were P = 51 %, I = 0.268, HT = 0.158 ± 0.023 and P = 53 %, I = 0.283, HT = 0.174 ± 0.024 for the M. pugana and M. pacifica s.s.-M. vallartensis genetic groups, respectively. The P value per taxon varied from 46 % to 51 %, with an overall mean of 64 % (Table 4). Shannon’s index (I) ranged from 0.268 to 0.275, with a total I of 0.309. Heterozygosity (HT) per taxon ranged from 0.158 ± 0.023 to 0.175, with an overall mean of 0.178 ± 0.023. Genetic diversity within populations (HS) ranged from 0.134 (M. pugana) to 0.156 (M. pacifica s.s.-M. vallartensis) with an average of 0.153 ± 0.022. Genetic diversity also varied among localities; M. pugana localities had the lowest genetic diversity parameters, whereas BA locality of M. pacifica s.s. showed the highest genetic diversity.
Genetic groups | Taxa | Locality | P (%) | I | HE | HT | HS | GST | D |
---|---|---|---|---|---|---|---|---|---|
APV | 38 | 0.227 | 0.136 | ||||||
RV | 38 | 0.218 | 0.121 | ||||||
ASL | 39 | 0.230 | 0.140 | ||||||
ALV | 38 | 0.208 | 0.138 | ||||||
M. pugana | M. pugana | 51 | 0.268 | 0.158 (0.023) | 0.134 (0.020) | 0.120 (0.021) | 0.028 (0.007) | ||
SS | 38 | 0.209 | 0.140 | ||||||
BA | 47 | 0.280 | 0.178 | ||||||
CSJ | 41 | 0.233 | 0.159 | ||||||
M. pacifica s.s. | 51 | 0.272 | 0.175 (0.025) | 0.159 (0.023) | 0.060 (0.016) | 0.018 (0.005) | |||
PV | 39 | 0.244 | 0.159 | ||||||
LL | 39 | 0.241 | 0.147 | ||||||
APM | 39 | 0.236 | 0.152 | ||||||
M. vallartensis | 46 | 0.275 | 0.171 (0.024) | 0.153 (0.022) | 0.073 (0.016) | 0.021 (0.006) | |||
M. pacifica s.s - M. vallartensis | 53 | 0.289 | 0.174 (0.024) | 0.156 (0.022) | 0.106 (0.016) | 0.026 (0.006) | |||
Total | 64 | 0.309 | 0.178 (0.023) | 0.147 (0.020) | 0.222 (0.039) | 0.040 (0.009) |
P proportion of polymorphic loci, I Shannon index, HE expected heterozygosity, HT total heterozygosity, HS intrapopulation heterozygosity, GST fixation index, D Jost’s differentiation index, standard deviation in parentheses.
Differentiation indices (GST and D) (Table 4) confirmed the results of the genetic structure analyses; moderate differentiation among localities was observed within each genetic group. Differentiation was higher in M. pugana (GST = 0.120 ± 0.021, D = 0.028 ± 0.007) than in M. pacifica s.s.-M. vallartensis (GST = 0.106 ± 0.016, D = 0.026 ± 0.006). Pairwise matrix analyses among localities using the Exact Test for population differentiation showed significant genetic differences (p < 0.05) among the localities that constitute the two subgroups of M. pugana (Table 5). This same result was observed when localities of M. pugana were compared to all of M. pacifica s.s.-M. vallartensis studied localities. In contrast, pairwise matrix analysis within each M. pugana subgroup, and within the M. pacifica s.s.-M. vallartensis group showed no significant genetic differences. Mantel tests revealed a significant correlation between geographical and genetic distances across the M. pacifica species complex r = 0.80, p < 0.001, showing a significant effect of isolation by distance (IBD) (Figure 6). However, when Mantel tests were applied to each of the two genetic groups separately, no significant correlation between geographical and genetic distances was detected.
M. pugana | M. pacifica s.s. | M. vallartensis | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
APV | RV | ASL | ALV | SS | BA | CSJ | PV | LL | APM | |
APV | - | 0.039 | 0.051 | 0.050 | 0.069 | 0.057 | 0.060 | 0.098 | 0.097 | 0.088 |
RV | 0.079 | - | 0.054 | 0.052 | 0.074 | 0.070 | 0.073 | 0.110 | 0.112 | 0.088 |
ASL | 0.006 | 0.001 | - | 0.032 | 0.060 | 0.065 | 0.050 | 0.077 | 0.079 | 0.072 |
ALV | 0.013 | 0 | 0.144 | - | 0.070 | 0.069 | 0.066 | 0.099 | 0.097 | 0.080 |
SS | 0 | 0 | 0.002 | 0 | - | 0.031 | 0.030 | 0.058 | 0.042 | 0.049 |
BA | 0.002 | 0 | 0 | 0 | 0.520 | - | 0.032 | 0.065 | 0.042 | 0.045 |
CSJ | 0 | 0 | 0 | 0 | 0.312 | 0.703 | - | 0.056 | 0.042 | 0.040 |
PV | 0 | 0 | 0 | 0 | 0.060 | 0.052 | 0.050 | - | 0.039 | 0.043 |
LL | 0 | 0 | 0 | 0 | 0.060 | 0.475 | 0.119 | 0.976 | - | 0.038 |
APM | 0 | 0 | 0 | 0 | 0.360 | 0.210 | 0.162 | 0.550 | 0.360 | - |
Conservation status assessment. The estimated EOO for M. pugana, M. pacifica, and M. vallartensis were 2178, 3196 and 124 km2, respectively, and AOO were 100, 104, and 44 km2, respectively.
Discussion
Species delimitation and assessment of species diversity have broad implications to biological conservation and evolutionary analysis (Coates et al. 2018). The population genetics approach with ISSR markers used in this work was useful to delineate species boundaries and relationships within the M. pacifica species complex. The genetic groups identified through a hierarchical Bayesian clustering analysis partially supported the current taxonomic delimitation among the studied species. In a first structure analysis, M. pugana was recognized as a genetically distinct group, whereas M. pacifica s.s. and M. vallartensis were clustered together. Barrier analysis results revealed that the TMVB acts as geographical barrier to gene flow between these two groups. The M. pugana group was internally more homogeneous (assignment proportion > 97 %) than the M. pacifica s.s.-M. vallartensis group; reduction in effective gene flow facilitates genetic and morphological divergences between these groups (Rieseberg & Willis 2007).
Interestingly, the second genetic cluster analyses showed substructures within each main groups; M. pugana was divided in the ASL-ALV and RV-APV subgroups, whilst in the second main cluster M. pacifica s.s. was separated from M. vallartensis. These results were confirmed by the dendrogram resulting from the UPGMA analysis. However, AMOVA results were not fully conclusive for the identification of species boundaries in the complex, and although a significant proportion of genetic variation either among taxa or genetic groups was shown, it was similar to the proportion of variation found within taxa or genetic groups. Most of the variation (≥ 82 %) was explained by differences within localities, which suggests high levels of cross-pollination for the studied taxa (Loveless & Hamrick 1984, Hamrick & Godt 1989). This result coincides with a recent molecular study indicating a predominance of outcrossing reproductive mode in some Neotropical magnolias (Veltjen et al. 2018). In agreement with our results, population genetic approaches with arbitrary dominant markers such as ISSR and AFLP have proved to be useful in the delimitation of diverse plant species. Bayesian analysis of genetic population structure using STRUCTURE and dominant markers has been conducted by several authors. In particular, Medrano et al. (2014) recognized two main genetic lineages in an endemic group of Narcissus species; Rodrigues et al. (2015) identified three main evolutionary lineages in Cattleya species; Esfandani-Bozchaloyi et al. (2018) delimitated 10 Geranium species; Sheidai et al. (2018) applied this strategy in Crocus species delimitation; and Pinangé et al. (2019) identified two population lineages in closely Dyckia species. Similar to our findings, AMOVA results on some of the studies enlisted above were less informative than Bayesian model-based clustering.
The present study revealed that genetic differentiation among allopatric populations is occurring in M. pugana. Gene barrier analysis suggests that the Santiago river canyon (184 km long, 3-15 km wide and 0.5-0.7 km deep) may act as a barrier to gene flow by limiting pollination and seed dispersal between the two subpopulations of M. pugana. Localities ASL (type locality) and ALV found southwest of the canyon (Barranca del Río Santiago) form a subpopulation, hereafter referred to as M. pugana, while the RV and APV found northeast of the canyon, represent another subpopulation, hereafter referred to as the Barrancae group. The Santiago river canyon has also been reported as an important barrier that limits the gene flow between Nolina parviflora populations (Ruiz-Sanchez & Specht 2013). According to Vázquez-García et al. (2016), allopatric speciation seems to be a major driver of Magnolia diversification in Neotropical Magnoliaceae. On the other hand, genetic differentiation between M. pacifica s.s. and M. vallartensis is also occurring, but this is lower than expected according to the morphological and ecological differentiation reported by Vázquez-García et al. (2012) and Dahua-Machoa (2018). Gene flow between these subgroups (assignment proportion 74 %) is higher than between the M. pugana subgroups (assignment proportion 70 %). The barrier detected between M. pacifica s.s. and M. vallartensis may correspond to several river canyons that stand among the sampled localities. Particularly CSJ is separated from the M. vallartensis localities by canyons of the rivers Ameca, Mascota, Cuale, and Pitillal; and SS is separated from M. vallartensis by the last three canyons. The climate in the lower zones of these canyons is too warm for these species. In contrast, there is no clear geographical or geomorphological barrier between some localities of these two taxa in the upper zones (i.e., BA and all the M. vallartensis localities). Therefore, genetic differences between these localities might be rather explained by a process of parapatric ecological differentiation (Nosil 2012). Sister species divergence in response to ecological factors is more frequent than formerly thought (Chapman et al. 2013, Anacker & Strauss 2014). M. vallartensis and M. pacifica s.s. occupy different habitats (ecotone of tropical montane cloud forest and tropical sub-evergreen forest at 100-1,120 m elevation vs. montane cloud forest at 750-2,250 m elevation), they also differ in flowering seasonality (absent vs. high), seasonality of seed dehiscence (high vs. moderate), and production of flowers, fruits and seeds (low vs. regular) (Dahua-Machoa 2018). Furthermore, M. vallartensis grows closer to the Pacific Ocean (0.5-20 vs. 40-70 km), with higher relative humidity and warmer climate (tropical maritime vs. temperate montane), higher mean annual temperature (20.9-25.6 °C vs. 16.6-19.9 °C) and higher annual rainfall (1,348-1,591 vs. 1,100-1,364 mm) (Fick & Hijmans 2017).
Molecular analysis of closely related taxa that are still exchanging genes is a recent approach in speciation research (Gao et al. 2019, Chapman et al. 2013, Feder et al. 2012). Gene exchange among closely related taxa happens in at least 25 % of plant species (Mallet 2005), therefore, further studies are necessary to assess the importance of parapatric and sympatric speciation in the Neotropical Magnolia species group. Complementary phylogenetic, biogeographic and phylogeographic studies are required to confirm species boundaries and relationships among the M. pacifica species complex revealed by ISSR dominant markers.
ISSR markers were found to be very effective in revealing the genetic variation within the M. pacifica complex. The genetic diversity estimates for the studied taxa (HT = 0.158-0.175, I = 0.268-0.275), were lower than the average values usually reported for plant genetic diversity based on ISSR (H = 0.22) (Nybom 2004), and much lower than those reported for two threatened eastern Mexican Magnolia species: M. sharpii Miranda (I = 0.56) and M. schiedeana Schlecht. (I = 0.50) (Newton et al. 2008). This is despite the fact that M. sharpii has an EOO (2,228 km2) similar to that of M. pugana (2,178 km2) and larger than that of M. vallartensis (124 km2). The three western Magnolia species of this study have much smaller EOO than that of the eastern M. schiedeana (17,411 km2) (Rivers et al. 2016); this difference in range size is in accordance with their lower genetic diversity. Genetic diversity of the taxa studied here is similar to that estimated for the rare and endangered Chinese species Magnolia cathcartii (HT = 0.162, I = 0.27) using AFLP data (Zhang et al. 2010). Some other studies based on ISSR have also shown low levels of genetic diversity on endangered and endemic plant species (Luan et al. 2006, Li & Jin 2007). M. pugana displayed the lowest genetic diversity (HT = 0.157, HS = 0.134) within the M. pacifica complex, this result is in accordance with expectations from population genetics theory for taxa with small populations in small geographic ranges (Godt & Hamrick 1999, Cole 2003). The highly fragmented small subpopulations of M. pugana are the result of both natural and anthropogenic causes; because of its higher continentality, M. pugana is exposed to greater water stress (less rain and fog), and more forest fires than M. pacifica s.s. and M. vallartensis. M. pugana is not a cloud forest species, but an element of riparian gallery forest, usually surrounded by a matrix of drier habitats such as Quercus-Pinus and tropical seasonal dry forests, and anthropogenic pasturelands.
Overall, differentiation indices (GST and D) were moderate and revealed that gene flow is ongoing within each genetic group or taxon. They were lower (GST = 0.060 - 0.120) than those reported for the critically endangered M. coriacea (GST = 0.187) estimated from ISSR markers (Zhao et al. 2012). Similarly, genetic patterns in various Neotropical Magnolia species were recently investigated and population differentiation was considered moderate to high in most species (Veltjen et al. 2018). The Exact Text of population differentiation confirmed the genetic structure and relationships observed in the M. pacifica complex. Genetic differentiation in Magnolia might be favored by limitations of pollination mechanisms, seed dispersal and initial establishment of seedlings (Loveless & Hamrick 1984). It is known that the transfer of pollen by beetles in small and isolated populations of Magnolia is not very efficient (Huang & Guo 2002, Chen et al. 2016), and flowering times do not always correspond with the emergence of pollinators (Dieringer & Espinosa 1994). Habitat fragmentation and other anthropogenic factors may also be troublesome for Magnolia seed dispersal by birds (Oppel & Mack 2010, Chen et al. 2016); even if dispersal occurs, seeds may not germinate or seedlings may not survive in places without enough humidity all year round (José et al. 2011). The limited gene flow between M. pugana and the Barrancae group or M. pacifica s.s. and M. vallartensis can exacerbate genetic differentiation and loss of genetic variation. Finally, results from the Mantel test applied to all localities of the complex and to each main genetic cluster suggest that geographic distance may have contributed to differentiation between population lineages but not within them.
Results from this study provide relevant and significant information for conservation strategies of M. pugana, M. pacifica s.s. and M. vallartensis. These western Mexican species are in an endangered status because of their highly fragmented populations, surrounded by seasonal drier habitats (Fick & Hijmans 2017), lower genetic diversity and narrower extent of occurrence compared with other threatened Magnolia species (Newton et al. 2008). EOO and AOO estimates for all three taxa are in accordance with the IUCN criteria for endangered species: B1 (EOO < 5,000 km2) and B2 (AOO < 500 km2). In addition, populations of these taxa are severely fragmented and present a continuing decline in the area, quality of habitat, and number of mature individuals [IUCN criteria B1ab (iii,v) + B2ab (iii,v)] (IUCN 2019). Measures of genetic diversity are not explicitly considered in the IUCN Red List of threatened species, however, this criterion is taken into account in the MER of the Mexican Official Norm 059 (SEMARNAT 2010). The genetic diversity of the three taxa studied here is considered lower than the average (H = 0.22) for diverse plant species with dominant markers (Nybom 2004), and it is even lower than that of M. sharpii (I = 0.56), a species categorized as endangered because of its narrow EOO (2,228 km2), severely fragmented and degraded habitats, and the fact that it is known from only five locations (Newton et al. 2008, Rivers et al. 2016). The low genetic diversity of the three taxa of the M. pacifica complex is in accordance with criterion C-2 (intrinsic biological vulnerability: genetic variability factor) of the MER for endangered species (SEMARNAT 2010).
In order to maintain the genetic diversity within the M. pacifica species complex as high as possible (Avise & Hamrick 1996, Xiao et al. 2004), the two subpopulations of M. pugana, and the subpopulations of M. pacifica s.s. and M. vallartensis should be considered as separate units of conservation. The implementation of in situ and ex situ conservation programs should protect and preserve at least one locality of each subpopulation of M. pugana and the most divergent localities of M. pacifica s.s. (BA and the type locality SS) and M. vallartensis (PV and the type locality APM). The BA locality of M. pacifica s.s. presented the highest genetic variability, which coincides with being the most important forest (maple forest) in terms of woody species richness, vascular plant endemism and floristic composition for western Mexico (Vázquez-García et al. 2000, Vargas-Rodríguez et al. 2010). Concerted efforts must be made to conserve all of these taxa, and a greater focus is required to protect M. pugana characterized by low levels of genetic variation. Moreover, M. pacifica s.s. and M. vallartensis must also be protected, because of their current threatened status according to the IUCN criteria (IUCN Standards and Petitions Committee 2019). These three taxa must be fully assessed with the MER method for inclusion in the list of endangered species in the official Mexican Norm NOM-059. Ecological studies in relation to future climatic scenarios should also be carried out in order to assess their impact on endangered species of Magnolia in western Mexico. Management plans and education in conservation issues, linked to workshops to rural and indigenous communities are badly needed to train people in propagation techniques of magnolias in order to decrease their threatened status and raise awareness of the fate of these relict species in extinction risk.
Supplemental data
Supplemental material for this article can be accessed here: https://doi.org/10.17129/botsci.2551