SciELO - Scientific Electronic Library Online

 
vol.97 número4Fenología, sincronía floral y éxito reproductivo de Neolloydia conoidea (Cactaceae)Sinorhizobium formadores de nódulos y hongos micorrízicos arbusculares (HMA) mejoran el crecimiento de Acacia farnesiana (Fabaceae): una alternativa para la reforestación del Cerro de la Estrella, México índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

Links relacionados

  • No hay artículos similaresSimilares en SciELO

Compartir


Botanical Sciences

versión On-line ISSN 2007-4476versión impresa ISSN 2007-4298

Bot. sci vol.97 no.4 México oct./dic. 2019  Epub 04-Feb-2020

https://doi.org/10.17129/botsci.2195 

Ecología

Impact of Late Pleistocene-Holocene climatic fluctuations on the phylogeographic structure and historical demography of Zamia prasina (Cycadales: Zamiaceae)

Impacto de las fluctuaciones climáticas en el Pleistoceno tardío-Holoceno sobre la estructura filogeográfica y demografía histórica de Zamia prasina (Cycadales: Zamiaceae)

Grecia Montalvo-Fernández1 
http://orcid.org/0000-0002-0627-3347

Lorenzo Felipe Sánchez-Teyer1 

Germán Carnevali1  2 

Andrew P. Vovides3 

Ricardo Gaytán-Legaria4 

Matilde Margarita Ortiz-García1 

Jaime Alejandro Muñoz-López1 

Jaime Martínez-Castillo1  * 
http://orcid.org/0000-0002-6484-6406

1 Centro de Investigación Científica de Yucatán, A. C, Mérida, Yucatán, México.

2 Orchid Herbarium of Oakes Ames, Harvard University Herbaria, Cambridge, Massachusetts, United States of America.

3 Instituto de Ecología, A.C., Xalapa, Veracruz, México.

4 Instituto de Investigaciones en Ecosistemas y Sustentabilidad, Universidad Nacional Autónoma de México, Morelia, Michoacán, México.


Abstract:

Background:

Glacial periods during the Pleistocene have been hypothesized to have greatly influenced geographical patterns of genetic structure and demography of many tropical species. The Glacial Refugium Hypothesis proposes that, during cold, dry glacial periods, populations of moisture-affinities tropical species were restricted to sheltered, humid areas and that, during warmer and more humid interglacial periods, these populations expanded. Some mountain regions in the tropics acted as refugia during the cold, dry periods of the Pleistocene for several temperate forest taxa, which recolonized the humid areas farther north during the interglacial periods.

Questions:

(1) Did Late Pleistocene-Holocene climate changes affect the historical demophraphy of Zamia prasina? (2) Does the historical distribution of Zamia prasina agree with the Glacial Refugium Hypothesis?

Study species:

Zamia prasina W.Bull. (Zamiaceae), the only cycad native to the Yucatan Peninsula Biotic Province (YPBP).

Methods:

Five individuals were collected in 23 populations and characterized using two DNA regions: plastid atpF-atpH, and nuclear ITS2. Genetic diversity, phylogeographic structure, historical demography, and potential distributions were assessed.

Results:

Our results showed moderately high genetic diversity and low, but significant, phylogeographic structure. Two genetic groups were identified, one in the eastern part of the Peninsula, the other in the western. The changes in historical demography suggest that Z. prasina experienced a population expansion following the warm conditions of the Holocene.

Conclusions:

The population dynamics of Zamia prasina are in accordance with the Glacial Refugium Hypothesis.

Key words: Cycads; demography expansion; diversity and genetic differentiation; glacial refugium hypothesis; Yucatan Peninsula Biotic Province (YPBP)

Resumen:

Antecedentes:

Los períodos glaciales ocurridos en el Pleistoceno tuvieron gran influencia en la estructura genética y demografía de muchas especies tropicales. La Hipótesis del Refugio Glacial propone que, durante los períodos glaciales secos y fríos, las poblaciones de especies tropicales se restringieron a áreas húmedas y que, durante los períodos interglaciales más cálidos y húmedos, estas poblaciones se expandieron. Algunas regiones montañosas de los trópicos actuaron como refugios durante los períodos fríos y secos del Pleistoceno para varios taxones de bosques templados, que recolonizaron las áreas húmedas más al norte durante los períodos interglaciales.

Preguntas:

1) ¿Los cambios climáticos del Pleistoceno tardío-Holoceno afectaron la demografía histórica de Zamia prasina? 2) ¿La distribución histórica de Zamia prasina responde a la Hipótesis de Refugio Glacial?

Especies de estudio:

Zamia prasina W. Bull. (Zamiaceae), la única cícada nativa de la Provincia Biótica peninsula de Yucatán (PBPY).

Métodos:

Se colectaron cinco individuos en 23 poblaciones y se caracterizaron utilizando dos regiones de ADN: cloroplasto (atpF-atpH) y nuclear (ITS2). Se evaluó la diversidad genética, estructura filogeográfica, demografía histórica y distribución potencial.

Resultados:

Nuestros resultados mostraron una diversidad genética moderadamente alta y una estructura filogeográfica baja, pero significativa. Se identificaron dos grupos de poblaciones diferenciados geneticamente, uno ubicado en el sureste y el otro al noroeste de la Península. Los cambios en la demografía histórica sugieren que Z. prasina experimentó una expansión poblacional en las condiciones cálidas del Holoceno.

Conclusiones:

La dinámica poblacional de Zamia prasina se ajustan a la Hipótesis de Refugio Glacial.

Palabras clave: Cicadas; expansión demográfica; diversidad y diferenciación genética; hipótesis de refugio glacial; Provincia Biótica Península de Yucatán (PBPY)

Throughout the Pleistocene (≈ 2.5 Ma), extended blankets of glacial ice covered the highest latitudinal zones of the planet, giving rise to glacial periods, which alternated with interglacial periods, when these areas were partially free of ice (Svensson et al. 2005). The glaciation period had a greater impact on the northern hemisphere biota than on the southern hemisphere given that the large land masses very close to the Arctic could channel the glaciers toward the south, whereas an ocean separated the Antarctic from the southern continents (Svensson et al. 2005). The glacial periods of the Pleistocene also greatly influenced the geographic distribution, the genetic structure and demography of tropical species (Stewart et al. 2010, Ramírez-Barahona & Eguiarte 2013). Some the mountain regions located in the tropics acted as areas of refugium during the cold, dry periods of the Pleistocene for several temperate forest taxa, which recolonized the humid areas farther north during the interglacial periods (Ramírez-Barahona & Eguiarte 2013). These Neotropical refugia might have remained continuously humid, while the savannahs and dry tropical forests expanded, leaving populations featuring high genetic structure in different areas of refugium (Haffer 1969, Graham 1973, Toledo 1982, Myers1982, Pennington et al. 2000). The Holocene was also a period of climatic fluctuations. The Holocene Climate Optimum (9-5 Ka) was a warm period; the global climate was most likely between 0.5-3 °C warmer than it is today. Then the temperatures decreased progressively with cyclic heating/cooling periods until now (Walker et al. 2012).

Tropical forests and their biota have a complex evolutionary history (Ornelas et al. 2013); for example, Dioon edule (Zamiaceae) (González-Astorga et al. 2003) and Zamia paucijuga (Zamiaceae) (Nolasco-Soto et al. 2015) experienced demography changes in their populations influenced by the Pleistocene glaciations. Historic demography changes in plant species could be consistent with either the Glacial Refugium Hypothesis (GRH), which propose that populations contracted to one or more southerly refugia during the cold-dry glacial and expanded out from them in warm-humidity interglacials. This hypothesis is widely accepted for temperate species such as subtropical columnar cacti in the mid latitudes of the Northern Hemisphere (Soltis et al. 2006). In contrast, the Interglacial Refugium Hypothesis (IRH) suggests that in the intertropical open dry vegetation of South America (Caatinga, Cerrado and Chaco biomes), populations of some species contracted to warm and humid refugia during the interglacials and expanded outward under the cold/dry climate of the Ultimate Glacial Maximum (LGM) (≈ 80 Ka) (Cornejo-Romero et al. 2017). Several authors have proposed a refugium theory as the underlying model for glacial and postglacial population dynamics of tropical species during the LGM (Farrera et al. 1999, Ramírez-Barahona & Eguiarte 2013, Cornejo-Romero et al. 2017).

The degree of aridity in the tropics generated by the reduced precipitation during the LGM and its influence on the distribution of the Neotropical forests have been controversial topics (Ramírez-Barahona & Eguiarte 2013). In one of the pioneer publications on the impact of the glaciations of the Pleistocene in the Neotropics, Van der Hammen (1961) postulated two opposing hypotheses: (1) During the glaciations, the lower temperatures generated a reduction in precipitation and, consequently, an increase in aridity in the Neotropics. (2) The cold phases of the glaciations were accompanied by an increase in precipitation, and the warm phases were much drier; therefore, the glaciations did not lead to a reduction in the precipitation and did not significantly affect the continuous and stable distribution of the tropical rainforests (Farrera et al. 1999). However, the information available to date has been based largely on limited and conflicting paleoecological data (Ramírez-Barahona & Eguiarte 2013). In general, the existing records indicate greater tropical aridity, with low lake levels in regions such as tropical Africa (Caballero et al. 1999). These records have been interpreted as an indication that, during the LMG, climates were drier (Kutzbach et al. 1993). In the Trans-Mexican Volcanic Belt, several records of the internal basins suggest greater aridity (Lozano-García et al. 2005, Caballero et al. 2010), while west-central Mexico have been proposed to have had relatively humid conditions (Bush et al. 2009, Bradbury 2000).

During the Last Interglacial periods (120 Ka), the climate of the YPBP became warmer and more humid (Metcalfe et al. 2000). During most of the LGM, the climate of the YPBP was dry, and the temperature was approximately 6 ºC colder than at present. Savannahs and scrublands covered most of the region up to the Early Holocene (10 Ka) (Orellana et al. 2003). In an analysis of Holocene fossil pollen samples from the YPBP, to reconstruct the vegetation and to develop a precipitation record for the last 7,900 years, a gradual increase in precipitation and expansion of the tropical forest during the Middle Holocene was shown (Vela-Pelaez et al. 2018).

Zamia prasina W. Bull. (Zamiaceae), distributed in the tropical rainforests of the central and southeast areas of the YPBP, is the only cycad native to this region. According to the International Union for the Conservation of Nature (IUCN), it is severely threatened (Vovides & Nicolalde-Morejón 2010). The populations of Z. prasina almost certainly were affected by the climate changes of the Pleistocene, which have been proposed as the driving force of the diversification in cycads (Nagalingum et al. 2011). The cycads are a species group with a long evolutionary history originating in the late Permian (Rull 2012). However, molecular data indicate that current species had a more recent origin (Late Miocene, ≈10-5 Ma) (Nagalingum et al. 2011). The almost simultaneous initiation of the diversification of six of the living genera of the cycads (in Australia, Africa, South-east Asia and tropical South America) suggests a single promotor may have been responsible for their diversification and that event may have been global climate change (Nagalingum et al. 2011). The climate fluctuations during the Pleistocene significantly influenced the current distribution of the cycads (González & Vovides 2002) and that, for the cycads of Mexico, the refugium areas of the LGM were fundamental in the definition of their phylogeographic dynamics (Contreras-Medina & Luna-Vega 2002). Consequently, we sought to determine whether the Late Pleistocene-Holocene (80-10 Ka) climate changes affected the historical demography of Zamia prasina and whether the historical distribution of Zamia prasina fits with the scenario proposed by the glacial refugium hypothesis. Thus, if this species underwent a contraction during the glaciations (LGM), then we would expect a reduction in the environmentally suitable area and in the effective size of the population during the LGM. Conversely, the effective size of the population and suitable area would increase during an interglacial period. To answer these questions, we used a phylogeographic approach, Bayesian analysis, and niche model. This study represents one of the first phylogeographic studies in the YPBP. We addressed the following objectives: (1) to determine the genetic diversity and phylogeographic structure of Z. prasina, (2) to evaluate the influence of the climate changes during the Late Pleistocene-Holocene on the historical demography of Z. prasina and (3) to determine whether the climate changes of the Late Pleistocene-Holocene modified the potential distribution range of Z. prasina.

Materials and Methods

Study area, populations and samples. This study was conducted in the Yucatan Peninsula Biotic Province (YPBP), which includes the Mexican states of Campeche, Quintana Roo, Yucatan and part of the states of Tabasco and Chiapas, northern districts of Belize (Orange Walk, Corozal, and Belize) and the department of Peten in Guatemala (Carnevali et al. 2010). The selection of populations sampled was based on information obtained from Herbarium CICY and from previous studies, Velasco (2015). Herbarium vouchers for each population were deposited at Herbarium CICY. We collected samples from 21 natural populations of Z. prasina in the Mexican part of the YPBP and five samples from two populations in Belize (two in Cayo District and three in Belize District) that were donated by Michael Calonje (Montgomery Botanical Center) (Figure 1, Table 1). It was not possible to obtain samples of Z. prasina from Guatemala; therefore, any reference to the YPBP here does not include Guatemala.

Figure 1 Location of the 23 populations of Zamia prasina sampled in the Yucatan Peninsula Biotic Province (YPBP). Biotic Province (YPBP). 1: Valladolid, 2: Peto, 3: Becanchén, 4: Rancho Duarte, 5: Kaxil Kiuic, 6: Champotón, 7: Calakmul, 8: Escárcega, 9: Hormiguero, 10: Virgensita, 11: Tenosique, 12: Palenque, 13: Puerto Morelos, 14: Cobá, 15: Carrillo Puerto, 16: Xhazil, 17: José María Morelos, 18: Pedro A Santos, 19: Bacalar, 20: Pantoha, 21: Nachicocom, 22: Belize-Belize 23: Belize-Cayo. 

Table 1 Population, Number of individuals for population (N), latitude, longitude and respective haplotypes for populations of Zamia prasina studied. 

Population N N. Latitude W. Longitude Haplotype atpF-atpH Haplotype ITS2
1 Valladolid 5 20.633342 -88.343722 A1 H1, H2
2 Peto 4 20.125567 -88.980083 A1, A2 H1
3 Becanchén 5 20.065992 -89.104694 A1, A2 H1, H3, H4, H5
4 Rancho Duarte 4 20.028056 -89.169531 A1, A2 H1, H3
5 Kaxil Kiuic 3 20.086944 -89.551528 A3 H1, H6
6 Champotón 5 19.04675 -90.405556 A1 H1, H3
7 Calakmul 5 18.409306 -89.899444 A1 H1, H10
8 Escárcega 5 18.617028 -90.850889 A1 H1, H3, H7
9 Hormiguero 3 18.359667 -89.491722 A1, A3, A7 H8
10 Virgensita 5 18.266389 -91.554 A3, A4, A5, A6 H6
11 Tenosique 5 17.591472 -91.554 A1 H1, H3, H9
12 Palenque 5 17.539278 -91.959889 A1 H1, H3, H7
13 Puerto Morelos 5 20.862083 -87.035139 A3 H1, H3
14 Cobá 5 20.499972 -87.711806 A1, A3 H1, H3
15 Carrillo Puerto 4 19.501583 -88.034056 A1, A3 H3, H12, H13, H14
16 Xhazil 5 19.400417 -88.092278 A1, A3 H3, H8
17 José María Morelos 5 19.140333 -88.594694 A1, A3 H3, H8, H15, H16, H17
18 Pedro A.Santos 5 18.950778 -88.178833 A3, A8 H3, H8
19 Bacalar 5 18.913278 -88.2325 A3 -----
20 Panto-Ha 5 18.578861 -88.454167 A1, A3 H8, H11, H12
21 Nachicocom 4 18.481778 -88.787833 A1, A3 H8, H12
22 Belize-Belize 2 16.998917 -89.071533 A1, A3 H12, H18
23 Belice-Cayo 3 17.302667 -88.488167 A1, A9 H12

DNA extraction, amplification and sequencing. DNA was extracted from individuals from each population, using the CTAB protocol (Doyle et al. 1987). We amplified the chloroplast intergenic region atpF-atpH and nuclear region ITS2 selected by Nicolalde-Morejón et al. (2011) using the methods of Lahaye et al. (2008) and Baldwin (1992), respectively. These markers were previously used in genetic studies of the genus Zamia. The amplified products were visualized in 1 % agarose, and 110 samples for chloroplast and nuclear regions were sent to in Macrogen 2018 for unidirectional sequencing using the forward primer and Sanger sequencing. The sequences were edited in Sequencher v. 5 (Gene Codes Corp., Ann Arbor, MI, USA) and aligned in the software PhyDE® (Müller et al. 2010) with the multiple alignment tool Muscle, then they were examined manually. For our ITS2 data set, we confirmed that multiple copies were not present for the following reasons: (1) all sequences had a GC content > 65 %, (2) all BLAST results coincided with the number of base pairs for ITS2 sequences for Zamia in GenBank, (3) conserved regions were homologous to ITS2 sequences from Zamia in GenBank, and (4) a few nucleotides differed between our sequences and those in GenBank.

Genetic diversity and phylogeographic structure. The number of haplotypes (k) and nucleotide and haplotype diversities (π and h, respectively) were determined using the DnaSP v.5.10 (Librado & Rozas 2009); indels were not considered. The genetic diversity was determined for each DNA region, population, population groups (based on genetic differences) and specie (all populations). The genetic structure was determined by means of a non-hierarchical molecular analysis of variance (AMOVA) with 1,000 permutations, using paired genetic differences with the program Arlequin 3.5 (Excoffier et al. 2005). To infer the phylogeographic structure, the relationship of N ST/G ST was calculated. Both are estimators of genetic structure; N ST considers the nucleotide differences among the haplotypes, whereas G ST is based on the number of haplotypes and their frequency. The statistical analysis was performed in the PERMUT 2.0 program (Pons & Petit 1996) with 1,000 permutations where, if 5 % of the permuted values are lower than the N ST observed, then N ST > G ST, thereby indicating phylogeographic structure. The relationships among the haplotypes were determined through haplotype networks using the TCS network algorithm in the PopArt 1.7 program (Clement et al. 2002). This network shows the relationship among the haplotypes and the number of mutational steps separating them. In addition, the geographic distribution of the haplotypes was represented on a map to graphically record the aggregation of the haplotypes.

To confirm the existence of population groups genetically differentiated, an analysis was conducted to detect the geographical location of genetic discontinuities among the populations using Monmonier’s maximum difference algorithm implemented in the program BARRIER (Manni et al. 2004). From the pairwise genetic distances obtained in Arlequin 3.5 (Excoffier et al. 2005), BARRIER identifies the edges where there are greater genetic distances. For obtaining confidence levels for the barriers, 100 replicas of the genetic distance matrix were calculated using Program R ver 4.3.2.

To explore the influence of climatic variables on the genetic structure, a principal component analysis (PCA) was carried out using data for 19 environmental variables, derived from temperature and precipitation, obtained from WorldClim Global Climate Data V. 1.4 (http://www.worldclim.org/version1) (Hijmans et al. 2005) with a resolution of 1 km2. Any enviromental variables with a Spearman’s rank correlation higher than 0.8 between them were not used. A final set of 11 environmental varibles was extracted for the 23 populations of our study (Appendix 1): Mean Diurnal Range (Mean of monthly (max temp - min temp)) (BIO2), Isothermality (BIO3), Temperature Seasonality (BIO4), Max Temperature of Warmest Month (BIO5), Temperature Annual Range (BIO7), Mean Temperature of Coldest Quarter (BIO11), Annual Precipitation (BIO12), Precipitation Seasonality (BIO15), Precipitation of Driest Quarter (BIO17), Precipitation of Warmest Quarter (BIO18) and Precipitation of Coldest Querter (BIO19). Past software version 3 was used (Hammer et al. 2001).

Estimate of divergence times. The time of divergence for intraspecific diversification of Zamia prasina and a possible relationship with pre-Pleistocene and Pleistocene events was estimated using Bayesian inference (BI) implemented by the program BEAST 2 (Bouckaert et al. 2014). The model of sequence evolution HKY+G was employed for both regions according to the results of the AIC model selection from jMODELTEST (Posada 2008). An uncorrelated relaxed clock Log Normal model (UCLD) and a coalescent model assuming constant size were used to model the tree prior. The tree was calibrated using the 95 % highest posterior density (HPD) age reported by Calonje et al. (2019). Prior distributions for all calibrated nodes were conservatively set to uniform using the minimum and maximum age bounds outlined below.

The age intervals for divergence between Microcycas calocoma and Zamia genus was used for the tree root node (33 to 84.5 Ma) (I), the crown node of Zamia (9 to 22.1 Ma) (II). Divergence between Zamia prasina and Zamia variegata Warsz. (0 to 0.82 Ma) (III) and the divergence between the clades formed by Z. prasina plus Z. variegata and the clade Z. spartea Ac.DC. plus Z. furfuracea Aiton plus Z. loddiguesii Miq. (0.51 to 1.53 Ma) (III). Twelve species (i.e., Microcycas calocoma, Zamia variegata, Z. furfuracea, Z. paucijuga, Z. loddigesii, Z. spartea, Z. lacandona, Z. pseudoparasitica, Z. manicata, Z. inermis, Z. soconuscensis, Z. fischeri, were chosen as outgroups.

For divergence time estimations, Markov chain Monte Carlo (MCMC) were run for three independent 50 million generations, sampling every 5,000 generations. BI analyses were run using the CIPRES Science Gateway (Miller et al. 2010). We combined the log and trees files from each independent run using LOGCOMBINER 1.8.0 (Drummond & Rambaut 2007), then viewed the combined log file in TRACER 1.6 to ensure that effective sample sizes for all priors and the posterior distribution were > 200, making sure that parameter values were fluctuating at stable levels. Based on these results, the first 5,000 trees were discarded as burn-in, and the remaining samples were summarized as a maximum clade credibility tree with mean divergence times and 95 % highest posterior density (HPD) intervals of age estimates in TREEANNOTATOR. Finally, these results were summarized in a single tree visualized in FIGTREE ver 1.3.1 (Rambaut 2009).

Historical demography. The demography processes were analyzed at two levels: (1) species and (2) population groups. Three types of analyses were carried out. (1) Tajima’s D (Tajima 1989) and Fu’s F S (Fu 1997): neutrality tests were used to detect departures from a constant population size under the neutral model. Population growth was indicated by significant negative values (p < 0.05) using Arlequin with 10,000 permutations. (2) Mismatch distribution of pairwise nucleotide differences (Rogers & Harpending 1992) was calculated and compared with expected values for an expanding population using the Ramos-Onsins & Rozas (2010)R 2 statistic. This statistic considers the sample size, the number of singleton mutations in a sequence, the average number of nucleotide differences between two sequences, and the total number of segregating sites. Lower values of R 2 (< 0.05) are expected under a recent population growth event (Ramos-Onsins & Rozas 2010). A unimodal type graph shifted to the left indicates many comparisons where the differences between pairs of sequences are small, suggesting a recent expansion in the populations this was carried out in the program program DnaSP 5.10 (Librado & Rozas 2009) with 10,000 permutations. (3) Bayesian skyline plot (Drummond & Rambaut 2007) to infer the changes of the effective population size over time, allowing the use of mutation models and independent replacement rates for each region of DNA. The model HKY+G was used for both DNA regions, with a strict molecular clock model and coalescent model. This analysis was performed at three levels: 1) each DNA region, 2) species and, 3) populations groups. The number of substitutions per site per year (s/s/y) for atpF-atpH was 0.00056 and 0.00181 for ITS2 and were used to date the crown radiation of Zamia paucijuga (Nolasco-Soto et al. 2015), also published by Nagalingum et al. (2011) for the crown-age of Zamia spp. Thirty million permutations (MCMC) were carried out, and trees were collected every 3,000 generations, using the program BEAST ver 1.8.0. Outputs were visualized with TRACER ver 1.6 to assess stationarity of the MCMC (effective sample sizes > 200).

Potential distribution. The ecological niche of Zamia prasina was modelated based on 80 occurrences obtained from the Global Biodiversity Information Facility (GBIF, http://data.gbif.org/species/browse/taxon/) and 103 occurrences compiled from fieldwork. Duplicates ocurrences and ocurrences with a distance < 20 km among them were removed using the package spThin (Lammens et al. 2015) to reduce overfitting promoted by spatial aggregation of the occurrence. A total of 43 records were used. The same 11 environmental variables that were used for the PCA were used for the Ecological Niche Modelling (ENM). The 11 bioclimatic variables were masked to extend of ecoregions (WWF 2006) were Z. prasina occurrence exist as hypothesis of the accessible area (M).

The niche of Z. prasina was modeled using the maximum entrophy algorithm in MAXENT v. 3.3.3 (Phillips et al. 2006) with 100 replicates using bootstrapping as a resampling method with no extrapolation and no clamping. For model evaluation 20 % of the total records were used. Logistic output was selected to obtain maps with suitability values. The ENM was projected in the YPBP and Central America, a center of diveristy for the genera Zamia. The minimum training presence was used as threshold to obtain a binary map of presence-absence. Finally, the ENM obtained was transferred into three past climatic scenarios: the Last Interglacial (LIG ≈ 120 ka), the Last Glacial Maximum (LGM ≈ 21 ka) and the Middle Holocene (MH ≈ 6 ka). Two general models of global circulation (GCM, http://www.worldclim.org/paleo-climate1) were employed to transfer the ENM into the LGM and MH scenarios: the Community Climate System Model (CCSM, Collins 2004) and the Model of Interdisciplinary Research on Climate (MIROC, Hasumi & Emori 2004). Both models simulate the climatic conditions in the LGM and MH, with a stronger reduction in the temperature assumed in the CCSM model in comparison with the MIROC model (Otto-Bliesner et al. 2007).

Results

We analyzed 102 sequences of the atpF-atpH region, with a length of 462 bp and 90 sequences of the ITS2 region, with a length of 301 bp. For atpF-atpH, 23 populations were analyzed. For ITS2 one population was excluded because its sequences coincided with an endophytic fungus.

Genetic diversity. The genealogical relationships indicated that, for atpF-atpH, the most frequent haplotypes were A1 and A3, which are thus candidates as ancestral haplotypes given their internal position in the network and their higher number of connections with the other less-frequent haplotypes. Most haplotypes were differentiated in a single mutational step, except for A6 and A7 (Figure 2-A). For ITS2, haplotypes H1, H3 and H8 were the most frequent and thus possible ancestral haplotypes, given their internal position in the network. Most haplotypes were differentiated in only one mutational step. The exceptions were H13, H14, H16 and H17 (Figure 3-A). For atpF-atpH, nine haplotypes were obtained (Table 1). The overall haplotype diversity (h) was moderate, with a mean value of h = 0.564; the nucleotide diversity (π) was low, with a mean value of π = 0.00186. The populations of Valladolid, Kaxil Kiuic, Champotón, Escarcega, Calakmul, Palenque, Tenosique, Puerto Morelos y Bacalar had only one haplotype; so, h and π were equal to zero (Appendix 2). Four populations had exclusive haplotypes: Virgencita (A4, A5 and A6), Hormiguero (A7), Belize-Cayo (A9) and Pedro A. Santos (A8) (Figure 2-B). For ITS2, 18 haplotypes were obtained (Table 1). Haplotype diversity (h = 0.827) and nucleotide (π = 0.00852) values were higher than those obtained for atpF-atpH. The populations of Peto, Virgencita, Hormiguero and Belize-Cayo presented only one haplotype; therefore, h and π were equal to zero (Appendix 3). Five populations had exclusive haplotypes: Becanchén (H4 and H5), Calakmul (H10), Carrillo Puerto (H13 and H14) and José María Morelos (H15 and H16) (Figure 3-B). List of GenBank accession for haplotypes are in Appendix 4.

Figure 2 (A) Haplotype network obtained with chloroplast sequence atpF-atpH. Circle size is proportional to the frequency of each haplotype; lines between the haplotypes represent the mutational steps. (B) Geographic distribution of haplotypes in the 23 populations of Zamia prasina sampled in the YPBP. Pie charts represent the haplotypes found for each population; section size is proportional to the number of individuals with that haplotype. Gray dotted lines represent the location of the most probable barriers obtained with BARRIER. Numbers in the pie are population codes (see Figure 1 for locations). 

Figure 3 (A) Haplotype network obtained with nuclear sequence ITS2. Circle size is proportional to the frequency of each haplotype; lines between the haplotypes represent the mutational steps. (B) Geographic distribution of haplotypes in 22 populations of Zamia prasina in YPBP. Pie charts represent the haplotypes found for each population; section size is proportional to the number of individuals with that haplotype. Gray dotted lines represent the location of the most probable barriers obtained with BARRIER. Numbers in the pie are populations codes (see Figure 1 for locations). 

Phylogeographic structure. The non-hierarchical AMOVA showed similar results for atpF-atpH and ITS2 regions: 33.85 and 37.23 %, respectively, of the total variation found among populations. For chloroplast and nuclear regions, the values of F ST were statistically significant (F ST = 0.338, p < 0.001; F ST = 0.372, p < 0.001; respectively), indicating genetic differentiation among the populations studied (Table 2). For both regions, the greatest diversity is found within the populations, 66.15 % for atpF-atpH and 62.77 % for ITS2. For atpF-atpH, the N ST value observed was statistically higher in comparison with the value of GST (0.433 and 0.370, respectively; 96.1 % of permuted values of N ST were lower than observed N ST). For ITS2, the observed value of N ST was significantly higher than the value of G ST (0.362 and 0.257, respectively; 86.8 % of permuted values of N ST were less than the observed N ST). Both results indicate signatures of phylogeographic structure.

Table 2 Non-hierarchical Molecular variance analysis (AMOVA) of 23 populations of Zamia prasina using the chloroplast region atpF-atpH and the nuclear region ITS2. 

atfF-atpH ITS2
Source of variation Sum of squares Variance components Percent variation Fixation index Sum of squares Variance components Percent variation Fixation index
Among populations 23.989 0.17088 33.85 F ST = 0.338* 64.778 0.53492 37.23 F ST = 0.372**
Within populations 26.383 0.33397 66.15 61.333 0.90196 62.77
Total 50.373 0.50484 126.111 1.43688

The chloroplast and nuclear regions of DNA analyzed, the geographical distribution of the haplotypes suggested the existence of two groups, one in the eastern part of YPBP and the other in the western part (Figures 2-B and 3-B). For atpF-atpH, haplotype A1 was the most represented in the western region of the YPBP. On the other hand, haplotype A3 was more abundant in the populations of the eastern section of the YPBP. For ITS2, the haplotype H1 was the most represented in the populations of western YPBP, and haplotypes H3, H8 were much more frequent in eastern YPBP. BARRIER analysis showed, for both regions (atpF-atpH and ITS2), that the most probable genetic discontinuities among populations correspond to the east and west regions of the PBYP: the group in the east includes Puerto Morelos, Carrillo Puerto, Bacalar, José María Morelos, Cobá, Panto-Ha, X-Hazil, Pedro A. Santos, Nachi Cocom, Calakmul, Hormiguero and Virgencita; the group in the west includes Valladolid, Peto, Becanchén, Rancho Duarte, Kaxil Kiuic, Escarcega, Champotón, Palenque and Tenosique (Figures 2-B and 3-B). This analysis supported the geographic distribution of haplotypes mentioned before.

Considering the evidence obtained regarding the population groups in the eastern and western YPBP, we then determined their genetic diversity (Table 3). For atpF-atpH, the mean haplotype diversity (h) and nucleotide diversity (π) were low for both groups: h = 0.373 and π = 0.00080 for the group in the east; h = 0.406 and π = 0.00210 for the group in the west. For ITS2 in both groups, diversity values were moderate than for atpF-atpH (h = 0.566 and π = 0.00701 for the group in the east; h = 0.600 and π = 0.00281 for the group in the west.

Table 3 Average genetic diversity in the eastern and western groups using the chloroplast region atpF-atpH and nuclear region ITS2. 

Statistic atpF-atpH
Eastern group Western group
h 0.373 ± 0.0075 0.406 ± 0.1290 (p = 0.002))
π 0.0008 ± 0.0001 0.0021 ± 0.0010 (p < 0.0001)
ITS2
h 0.566 ± 0.100 0.600 ± 0.106 (p = 0.822)
π 0.007 ± 0.0039 0.002 ± 0.0005 (p < 0.0001)

The principal component analysis (PCA) shows two groups of populations, where the populations in the east of the YPBP are located and another group where most of the western populations of the YPBP are found (Appendix 5). Most of the variation was obtained in two components: PCA1 contributed 55.03 % of the variation and PCA2 contributed 21.22 % explaining in total 76.25 % of the variation. For PCA1, six variables were associated, with two temperature variables (BIO2 and BIO7) contributing the most; however, four of the variables in this component corresponded to variables related to precipitation (BIO15, BIO17, BIO18 and BIO19). Two variables were associated with PCA2, one related to temperature (BIO3) and another to precipitation (BIO12), which contributed the most to the component.

Divergence time. The comparison of all the coalescence analyses revealed high convergence between the inferred parameters, and the effective sample sizes were always greater than 200. The trees obtained for atpF-atpH and ITS2 (Figures 4 and 5, respectively) were congruent in the principal nodes and with the phylogeny of Zamia (Calonje et al. 2019): Microcycas A. DC. as sister group of Zamia; the clade formed by Z. furfuracea, Z paucijuga Wieland, Z. loddigesii and Z. spartea as a sister group of the clade formed by Z. variegata and Z. prasina; and Z. variegata as sister species of Z. prasina. The genealogy for atpF-atpH showed that the haplotypes of Z. prasina grouped into two clades: (1) haplotypes A1, A2, A7 and A9 and (2) A3, A4, A5, A6 and A8. The greatest divergence age between these two haplotypes groups was 1.37 (95 % HPD, 2.3-0.46 Ma), and the ancestral haplotype is A7 (Figure 4, Appendix 6), which was found exclusively in the Hormiguero population (Campeche) (Table 1).

Figure 4 (A) Plant of Zamia prasina in its natural habit. (B) Bayesian chronogram for chloroplast sequence atpF-atpH. The divergence time of Zamia prasina haplotypes in the YPBP and other cycad species is shown. The 95 % confidence intervals are shown with purple bars; the numbers at the nodes indicate the estimated age, the numbers below the branches indicate the posterior probability (Appendix 6). (C) Extract from the complete chronogram corresponding to the Zamia prasina haplotypes (marked with the red line), the colors of the squares correspond to the haplotypes shown in Figure 2-A

Figure 5 (A) Bayesian chronogram for nuclear sequence ITS2. The divergence time of Zamia prasina haplotypes in the YPBP and other cycad species is shown. The 95 % confidence intervals are shown with purple bars; the numbers at the nodes indicate the estimated age, the numbers below the branches indicate the posterior probability (Appendix 7). (B) Female cone of Zamia prasina. (C) Extract from the complete chronogram corresponding to the Zamia prasina haplotypes (marked with the red line), the colors of the squares correspond to the haplotypes shown in Figure 3-A

The topology obtained with ITS2 had several clades with short branches, indicating low differentiation between the haplotypes and that, based on molecular clock theory, the times of divergence must be recent. The haplotypes of Z. prasina grouped into six clades: (1) haplotypes H13, H14 and H16; (2) H17, 2, 6, 4, 1 and H10; (3) H3, H5, H7; (4) H18; (5) H8 and H9; and (6) H15, H11, H12. The greatest divergence age among group of haplotypes was 1.36 (95 % HPD, 2.23-0.48 Ma). The ancestral haplotype was H12 with a divergence time of 69,000 years (Figure 5, Appendix 7). This haplotype was present in individuals from Panto-Ha, Nachi Cocom, Carrillo Puerto and was the only haplotype present in the population Belize-Cayo (Table 1).

Historical demography. Tajima’s D and Fu’s F S were negative values at the species and population levels; however, these values were not significant (p > 0.05). This result was similar for atpF-atpH and for ITS2 (Table 4). The mismatch distribution analysis showed a distinctive unimodal pattern, indicating recent demography expansion (R 2= 0.0477 and R 2= 0.0455 for atpF-atpH and ITS2 respectively). The skyline plot for the species level indicated that the effective size of the population decreased, and much more so for ITS2, then increased recently, in the last 10,000 years approximately, which can be interpreted as population growth, which we expect would be followed by an expansion (Figures 6-A and 6-D). At the population group level, the skyline plot showed one pattern of demography expansion only for the population groups in the east for ITS2 (Figure 6-F).

Table 4 Average values of the neutrality indexes Tajima’s D and Fu’s F S , for the population groups in the eastern and western parts of the peninsula and for all populations. 

Statistic atpF-atpH
Eastern Western Total
Tajima’s D 0.01856 (p = 0.68700) -0.17476 (p = 0.83408) -0.05336 (p = 0.80242)
Fu’s F S 0.31829 (NA) 0.09316 (NA) 0.24214 (NA)
ITS2
Tajima’s D 0.0397 (p = 0.71136) -0.006129 (p = 0.68373) -0.03263 (p =0.69850)
Fu’s F S -0.05500 (NA) 0.05405 (NA) 0.11408 (NA)

Positive values for D and F S are indicative of mutation-drift equilibrium, typical of stable populations; negative and significant values (p < 0.05) resulting from an excess of rare haplotypes indicate that populations have experienced recent expansions, often preceded by a bottleneck.

Figure 6 Demography history of Zamia prasina in YPBP based on Bayesian skyline plot. Substitution rates reported by Nolasco-Soto et al. (2015) were used. The vertical axis corresponds to the effective size of the population (Ne), the black line indicates the trend for the median Ne over time; purple lines are the 95% confidence intervals. The horizontal axis represents time in thousands of years. Graphs on left atpF-atpH (A: All populations, B: populations of the east group of YPBP, C: populations of the west group). Graphs on right ITS2 (D: All populations, E: populations of the east group of the YPBP, F: populations of the west group). 

Potential distribution. The ENM performed well (AUC= 0.907 ± 0.018) and the predicted distribution (Figure 7) matched the actual known range of Z. prasina (Vovides & Nicolalde-Morejón 2010) associated with tropical evergreen and tropical deciduos forests in the YPBP and north of Nicaragua where other Zamia species are distributed. The ENM transferred into a past climatic scenarios shows drastic range shift since the LIG from present. During the last interglacial period (LIG ≈ 140-120 ka), only the northern area of Guatemala and Chiapas were suitable for Z. prasina. The suitable area during the Last Glacial Maximum (LGM ≈ 21-18 ka) for the two GCM were highly contracted. According to the CCSM model, during the LGM the species was possibly distributed on the eastern coast of the YPBP and in some areas of Central America. The MIROC model, on the other hand, indicated that conditions were suitable in a small area in the eastern part of the YPBP and western part of the Yucatan Peninsula. During the Middle Holocene (MH) for both models, the suitable area for Z. prasina increased, and expanded more to the northwest, so that it inhabited a large part of the YPBP and the southern coast of the Gulf of Mexico.

Figure 7 Ecological niche modelling for Zamia prasina. Current potential distribution. (LIG) Predicted distribution during the interglacial period (≈140-120 ka). (LMG-CCSM) Predicted distribution during the Last Glacial Maximum, CCSM model (≈ 21-18 ka). (LMG-MIROC) Predicted distribution during the Last Glacial Maximum, MIROC model. (HM-CCSM) Predicted distribution during the Middle Holocene, CCSM model (≈ 6 ka). (HM-MIROC) Predicted distribution during the Middle Holocene, MIROC model. Green: high suitability, light brown: low suitability. 

Discussion

Phylogeography. The moderate and high haplotype diversity for atpF-atpH and ITS2, respectively, and their low nucleotide diversity indicate the presence of many haplotypes at high frequency but with few nucleotide differences among them; thus, the diversification of these haplotypes is recent. Similar results have been reported by Nolasco-Soto et al. (2015) for Z. paucijuga, with h = 0.669 and π = 0.0013 for the chloroplast region psbK-psbI and h = 0.843 and π = 0.0063 for ITS2. For Dioon sonorense,Gutiérrez-Ortega et al. (2014) reported h = 0.629 and π = 0.0004 for the chloroplast region trnL-trnF. For the cycads in Mexico, there are reports of high genetic diversity in some species, with few populations and restricted habitats, as in the case of Zamia loddigesii Miq. (González-Astorga et al. 2006), Dioon caputoi De Luca, Sabato & Vázq. Torres, D. merolae De Luca, Sabato & Vázq. Torres (Cabrera-Toledo et al. 2010) and D. sonorense (De Luca, Sabato & Vázq.Torres) Chemnick, T. J. Greg. & Salas-Mor (González-Astorga et al. 2008).

Disagreements in the estimated genetic diversity values for atpF-atpH and ITS2 in Zamia prasina may be due to the differences between these two DNA types. cpDNA is inheritanced maternally in cycads (Cafasso et al. 2001) so there is no genetic recombination, whereas it is reasonable to assume that nDNA is inherited bipaternally. In plants, cpDNA markers display average evolutionary rates in the range of 10-9 substitutions/site/year, relatively slower than nDNA rates (Wolf et al. 1987); our results are consistent with this difference.

The moderately-high genetic diversity found in Z. prasina may be due to three factors. (1) Although its geographical distribution is restricted and its habitat perturbed by anthropogenic activities, this species is abundant, and there are still large extension of vegetation where the species is common (Vovides & Nicolalde-Morejón 2010). The Global Biodiversity Information Facility (GBIF, http://data.gbif.org/species/browse/taxon/) reports more than 100 records of Z. prasina. The majority of them are in Mexico and Belize, but there are still areas that have been little explored to the south of Peten in Guatemala (Districts of Izabal and Alta Verapaz), where the species can be found (M. Veliz, curator of the BIGU herbarium pers. comm). (2) The species is dioecious with a cross-mating system that is mediated by pollinators, which implies that genetic material is exchanged between individuals, thus generating greater genetic diversity (Vovides & Nicolalde-Morejón 2010). (3) Most likely, a recent population expansion increased the haplotype diversity with low values of nucleotide diversity; this pattern was reported for other zamias (i.e., Z. paucijuga) (Nolasco-Soto et al. 2015). As a population expands, it may come in contact with populations that had been isolated during the glaciation and subjected to different selection pressures. This contact between divergent populations can increase the diversity due to genetic flow between them.

The number of samples analysed also can influences the amount of genetic diversity detected. In cycads, a range of sample sizes have been used for phylogeographic studies, for example, one individual was used per population for Dioon sonorense (Zamiaceae), and moderate values of nucleotide diversity (h = 0.629) were obtained using the chloroplast region trnL-trnF (Gutiérrez-Ortega et al. 2014). For Zamia paucijuga (Zamiaceae) with 8-10 individuals per population, high haplotid diversity was found using ITS2 (h = 0.843) and moderate diversity for the chloroplast region psbK-psbI (h = 0.669) (Nolasco-Soto et al. 2015). Here we used a sample size between these two examples among these reported (five individuals per population) and still managed to detect a moderate to-high genetic diversity (atpF-atpH, h = 0.564; ITS2, h = 0.827) in agreement with the previous studies. Since we are working with non-coding regions that have low mutation rates, small sample sizes might not detected, the actual diversity.

According to our chronogram (Figures 4 and 5), the diversification within Z. prasina happened during the Pleistocene. During this period, neotropical montane forests experienced extremely complex glacial-interglacial dynamics and the effect of climatic fluctuations on the genetic structure and population history of species distributed in these habitats led to different outcomes, such as rapid radiation or local extinction (Ramírez-Barahona & Eguiarte 2013). (4) The moderate genetic diversity detected for Z. prasina might also be the result of a still unfinished process of habitat fragmentation, in which subsequent isolation has not yet affected the distribution of its molecular variants.

The greatest variation is within the populations and not between them may be because not enough time has elapsed for a significant divergence between them.The comparison of genetic diversity between the eastern and western population groups on peninsula indicates that there are no differences between these two. We found populations that stand out due to their high genetic diversity, for atpF-atpH in the Virgencita and Hormiguero populations and for ITS2 in the Jose Maria Morelos, Carrillo Puerto and Bekanche populations.

The east-west phylogeographic structure in the populations of Z. prasina (suggested by the geographic distribution of the haplotypes and BARRIER analysis) can be explained by climatic factors. The PCA supports the east - west divergence and indicates that the variables related to precipitation contributed greatly to PCA1 and PCA2; thus, the genetic differentiation between the group of populations from the eastern and western parts of the peninsula may be due to the precipitation gradient in the YPBP. There is a humid region in the southeast and a dry region in the north-northeast (Carnevali et al. 2010). The eastern part of the YPBP is humid with 1,200-1,500 mm of annual rainfall, and the basin of Laguna de Terminos in southeastern YPBP receives up to 1,400 mm, compared with 500-1,000 mm in the northwest part of the YPBP (Orellana & Espada 2009). This gradient of precipitation determines the type of vegetation; in the northeast, low-elevation deciduous forest predominates, in the middle is low-elevation deciduous forest, in the east is predominantly medium-elevation semideciduous forest, and in the southeast, high-elevation evergreen forest predominates (Carnevali et al. 2012). Thus, the geographic distribution of Zamia prasina in the YPBP could be determined by precipitation. Although no published studies have focused on the influence of precipitation, the drier northeastern strip constitutes a barrier to the distribution of this species, as indicated by the absence of historical records in this area (GBIF, http://data.gbif.org/species/browse/taxon/). The populations group in the west is found in the medium subdeciduous forest, bordering the dry northeast portion of the YPBP, while, the eastern population group is in a medium evergreen forest in a more humid area and thus more suitable for Z. prasina. Differences in geography, forest and climate may thus also contribute to the differentiation (Vovides & Nicolalde-Morejón 2010). To test this hypothesis, a more comprehensive study with a landscape genomics approach should be carried out.

Regarding the influence of the precipitation gradient of the YPBP on Zamia prasina, Vovides & Olivares (1996) analyzed the caryotypical variation in 11 individuals of Z. prasina in the Mexican part of the YPBP. They found that the plants collected in the northeastern part of this region, where habitats are drier, have higher chromosome numbers (2n = 24-27) than in plants collected in the southeastern region (2n = 17), where the habitats are more humid. Vovides & Olivares (1996) concluded that the stressful drought conditions present in the northeast generate chromosomal changes. The chromosomal changes in Zamia are attributed to chromosomal fusion or fission. The fusion occurs between two acrocentric chromosomes (chromosomes whose centromere is closer to one of the ends), giving rise to a metacentric chromosome (with the centromere in the center of the chromosome) and thus reducing the number of chromosomes. During fission, a metacentric chromosome arises by fission or rupture of a metacentric chromosome into two acrocentric chromosomes, in this case increasing the number of chromosomes. Chromosomal fissions might confer adaptive advantages to survive in dry environments. Similar results were reported for other Zamia species of the Caribbean, which live in more stable conditions than do the continental Zamia species, which are exposed to a wide range of ecological conditions, from humid to dry and semi-xeric (Balduzzi et al. 1982).

In a study of the influence of environmental conditions on phenotypic diversity of Z. prasina,Limón et al. (2016) found that in populations with low precipitation and high temperature, the plants have fewer but wider leafs; in contrast, in areas with greater precipitation, the plants have more but thinner leafs. These authors proposed that Z. prasina inhabits areas under the dense canopy of tropical rainforests and semi-deciduous forests; therefore, the leaflets may be an adaptation to facilitate evapotranspiration and light capture in the shade. According to Stevenson et al. (1996), individuals of Z. prasina in intermediate environments have intermediate morphological forms. Thus, evidence of phylogeographic structure, in karyotype and morphology indicate that the precipitation gradient in the YPBP greatly influences the populations of Z. prasina. A study on the influence of precipitation on the genetic diversity of this species is necessary.

The present study shows that the estimated time of divergence among the most ancestral haplotypes of Zamia prasina was approximately 1.37 Ma based on atpF-atpH and 1.36 Ma based on ITS2; although most of the haplotypes have a more recent origin, indicating a recent population expansion that coincides with the Pleistocene (1.8-0.10 Ma), a period of events relevant to the evolution of the cycads. During the climatic changes of the Pleistocene, the tropical forests were subjected to a complex dynamic, an aspect that favored the presentation of rapid radiations and local extinctions (Ramírez-Barahona & Eguiarte 2013), and in the case of Z. prasina, there was an increase in genetic diversification.

Based on molecular phylogeny, Nagalingum et al. (2011) pointed out that, diversification increased in the genus Zamia during the Pleistocene and that it was one of the most speciose rich genera of cycads at that time. Pleistocene speciation prevailed in the neotropics, mainly for Zamia and Cerotozamia with the glacial-interglacial cycles being the predominat environmental force. The climate during that period also became more seasonal. Nolasco-Soto et al. (2015), employing nuclear and chloroplast markers determined the origin of Z. paucijuga between 6.4-1.6 Ma. Nolasco-Soto et al. (2015) suggested that the historical factors that affected Z. paucijuga at the population level during the Late Pleistocene similar affected other Zamia species, as our results showed for Z. prasina.

The values of the Tajima’s D and the Fu’s F S were negative but not significant, indicating that there was no population expansion. Because Z. prasina is a perennial with a long generational cycle the mutation rate is very low (Loveless & Hamrick 1984), thus, the significance of these two statistics is difficult to assess. However, the present work generated several lines of evidence that suggest that populations of Z. prasina underwent a recent demography expansion. (1) In the populations sampled, haplotype diversity is high, and nucleotide diversity is low. The low genetic divergence between the haplotypes can be explained because not enough time has elapsed for a significant divergence between the haplotypes (Lavery et al. 1996). (2) The star-shaped network of haplotypes found for Z. prasina in atpF-atpH and ITS2 regions also indicate a population expansion and one or a few haplotypes were more frequent, from which many differing haplotypes were derived in one or a few mutational steps; moreover, these haplotypes are found in low frequency. (3) The skyline plot showed a recent population expantion at species level, the effective size of the populations underwent a recent increase in the late Pleistocene-Holocene. (4) The diversification of haplotypes increased in the Late Pleistocene-Holocene.

Current and palaeodistribution modelling. As mentioned earlier, the palaeoecological studies of fossil pollen in the YPBP that reconstructed past of vegetation distribution and climate (Carrillo-Bastos et al. 2010, Sánchez-Sánchez & Islebe 2002) showed that the distribution in the YPBP was governed mainly by precipitation levels. A clear example is given by the certain ones semi-evergreen forest species, which grow in areas with a precipitation between 1,000 and 1,500 mm/year (Sánchez-Sánchez & Islebe 2002). The mosaic of current vegetation formed during the early Holocene, which implies that the modern pattern of isohyets was also established then (Leyden 2002). In the present, the climate is more seasonal, with rainier summers and less-intense winters, which encourages the regeneration and survival of seedlings. Z. prasina is actually more widely distributed in the YPBP, in ecosystems with seasonal climates (evergreen, deciduous and subdeciduous tropical forest).

In a geospatial analysis study of pollen records from the YPBP, Carrillo-Bastos et al. (2012) found that changes in vegetation during the Holocene in the peninsula were caused by changes in climate. In addition to climate influences during the Classic period (2,500-1,200 BC), the activities carried out by the Mayan people placed pressure on the forest coverage. In spite of this pressure, the vegetation coverage was not severely diminished; thus, the use of forest resources did not result in total deforestation. Eventually, vegetation recovered in relation to an increase of precipitation that coincided with the Medieval warm period during the Holocene. The highest pollen percentages taxa of semi-evergreen forests (85-92 %) were predicted to be in the southern and central parts of the peninsula and the lowest in the northern part. This distribution is similar to the modern distribution according to the isohyets, with precipitation increasing toward the eastern coast (Carrillo-Bastos et al. (2012) and consistent with our results that, during the Middle Holocene (MH), the environmentally suitable area for Z. prasina increased and expanded from the southeast to the center, covering a large part of the YPBP.

Although the LGM was a very dry period, some pleaces have been proposed to be relatively wet, e.g., center-west of Mexico (Bradbury 2000) and south of the YPBP (Hodell et al. 2008, Bush et al. 2009). Drill cores obtained from Lake Petén Itzá, Petén, Guatemala, contain a ~85-kyr record of terrestrial climate from lowland Central America that was used to reconstruct hydrologic changes in the northern Neotropics during the last glaciation.

Sediments are composed of clay, reflecting a relatively wet climate, and pollen from the same period indicates vegetation consisted of a temperate pine-oak forest. This finding contradicts previous inferences that the climate was arid during the Last Glacial Maximum (LGM) at least in the south of the YPBP (Hodell et al. 2008). This result coincides with ours that, during the LGM, the environment was suitable for the establishment of Z. prasina was concentrated in the eastern YPBP where the species may have found refuge during the cold, dry periods of the glaciation.

We propose that the absence of environmental suitability in the YPBP for the establishment of Z. prasina in the Last Interglacial (≈ 120 ka) could be due to some Zamia species needing a seasonal climate, while others have a tree habit and are adapted to very humid, woodland habitats (Nicolalde-Morejón et al. 2011). It is possible that the last interglacial was warmer, humid and without seasonality, in comparison with the Holocene (Walker et al. 2012, Cornejo-Romero et al. 2017), and thus was not propitious for the growth of Z. prasina.

The predictions by theMIROC and the CCSM models for the paleodistribution of Z. prasina in the LGM are different. The CCSM model predicts environmental suitability areas on the eastern coast of the YPBP, by the MIROC model does not predict any favorable areas on the peninsula (Figure 7). This difference may be due to these models basing their simulations on different initial experimental conditions using in different algorithms (Harrison et al. 2016). CCSM model estimates a lower temperature and greater precipitation than the MIROC model does (Otto-Bliesner et al. 2007) although the differences vary depending on the area of ​​study (Taylor et al. 2012, Harrison et al. 2016).

The genetic data (phylogeographic structure and historical demography) suggest that the populations of Z. prasina expanded during the warmest periods of the Mid-Holocene, a pattern consistent with the glacial refugium hypothesis and validated by its paleodistribution. Post-Pleistocene expansion has been proposed for other species of the genus Zamia, such as Z. paucijuga (Nolasco-Soto et al. 2015); after the glaciations, its populations that originated in central Mexico became established to the south in the Pacific watershed.

In our study, the presence of ancestral haplotypes on the southeastern peninsula (based on ITS2 evidence) indicates possible center of origin, although based on atpF-atpH the center could be in Campeche. This difference may be due to the differing type of inheritance between the two markers. However, these results are only preliminary because they are based on information from two short genome regions. Future studies should address this incongruity by using sampling more populations and genomic regions.

In general, our results indicate that climatic fluctuations during the Pleistocene-Holocene influenced the evolutionary history of Z. prasina. The ecological niche modelling suggests that this species prevail during glacial events in restricted populations on the eastern coast of the peninsula, and expanded post-glacially into the northwest to occupy almost the entire territory of the YPBP. This result are congruent with the phylogeographic pattern of other Zamia species (Nolasco-Soto et al. 2015) and helps us understand the historical phylogeographic patterns of the Zamiaceae in the Neotropics. An analysis of the comparative phylogeography at the genus level for Zamia with greater genomic sampling will allow us to make inferences about the temporal congruence of species and about the influence of climate on the composition of the genus.

Acknowledgments

This paper is part of the research for the first author’s Ph. D. thesis at the Centro de Investigación Científica de Yucatán, A. C., postgraduate studies in Biological Sciences Option Natural Resources. The first author thanks the Consejo Nacional de Ciencia y Tecnología-Mexico for the scholarship awarded for her postgraduate studies. Authors thank to four anonymous reviewers for their critical review of a previous version of this manuscript. To Dr. Diego Angulo for his help with the edition of the figures. Dr. Michael Calonje for Belize samples, we thanks too MSc. Iván Tamayo and Dr. Michael Calonje for theirs help in the diversification timing analysis and Paulino Simá for his assistance in the field.

Literature cited

Balduzzi A, De Luca P, Sabato S. 1982. A phytogeographical approach to the New World cycads. Delpinoa 23-24: 185-202. [ Links ]

Baldwin BG. 1992. Phylogenetic Utility of the internal Transcribed Spacet of Nuclear Ribosomal DNA in Plants: An Example from the Compositae. Molecular Phylogenetic and Evolution 1: 3-16. DOI: https://doi.org/10.1016/1055-7903(92)90030-k [ Links ]

Bradbury JP. 2000. Limnologic history of Lago de Pátzcuaro, Michoacán, Mexico for the past 48000 years: impacts of climate and man. Palaeogeography, Palaeoclimatology, Palaeoecology 163: 65-95. DOI: https://doi.org/10.1016/s0031-0182(00)00146-2 [ Links ]

Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, Suchard MA, Rambaut A, Drummond AJ. 2014. BEAST 2: A Software Platform for Bayesian Evolutionary Analysis. PLOS Computational Biology 10: e1003537. DOI: https://doi.org/10.1371/journal.pcbi.1003537 [ Links ]

Bush MB, Correa Metrio AY, Hodell DA, Brenner M, Anselmetti FS, Ariztegui D, Mueller AD, Curtis JH, Grzesik DA, Burton C, Gilli A. 2009. Re-evaluation of climate change in lowland Central America during the Last Glacial Maximum using new sediment cores from lake Petén Itzá, Guatemala. In Vimeux F, Sylvestre F, Khodri M, eds. Past Climate Cariability in South America and Surrounding Regions: Developments in Paleoenvironmental Research 14: 113-128. DOI: https://doi.org/10.1007/978-90-481-2672-9_5 [ Links ]

Caballero M, Lozano GS, Ortega B, Urrutia J, Macías JL. 1999. Environmental characteristics of lake Tecocomulco, northern basin of Mexico, for the last 50,000 years. Journal of Paleolimnology 22: 399-411. DOI: https://doi.org/10.1023/A:1008012813412 [ Links ]

Caballero M, Lozano-García S, Vázquez-Selem L, Ortega B. 2010. Evidencias de cambio climático y ambiental en registros glaciales y en cuencas lacustres del centro de México durante el último máximo glacial. Boletín de la Sociedad Geológica Mexicana 62: 359-377. DOI: https://doi.org/10.18268/bsgm2010v62n3a4 [ Links ]

Cabrera-Toledo D, González-Astorga J, Nicolalde-Morejón F, Vergara-Silva F, Vovides AP. 2010. Allozyme diversity levels on two congeneric Dioon spp (Zamiaceae, Cycadales) with contrasting rarities. Plant Systematics and Evolution 290: 115-125. DOI: https://doi.org/10.1007/s00606-010-0354-6 [ Links ]

Cafasso D, Cozzolino S, Caputo P, De Luca P. 2001. Maternal inheritance of plastids in Encephalartos Lehm. (Zamiaceae, Cycadales). Genome 44: 239-241. [ Links ]

Calonje M, Meerow AW, Patrick-Griffith M, Salas-Leiva D, Vovides A P, Coiro M, Francisco-Ortega J. 2019. A time-calibrated species tree phylogeny of the New World cycad genus Zamia L. (Zamiaceae, Cycadales). International Journal of Plant Sciences 180: 286-314. DOI: https://doi.org/10.1086/702642 [ Links ]

Carnevali FCG, Tapia-Muñoz JL, Duno de Stefano R, Ramírez-Morillo I. 2010. Flora Ilustrada de la Península de Yucatán: Listado Florístico. Yucatán, México: Centro de Investigación Científica de Yucatán, A.C. ISBN: 9686077823070 [ Links ]

Carnevali FCG, Tapia-Muñoz JL, Duno de Stefano R, Ramírez-Morillo I, Can-Itzá L, Hernández-Aguilar S, Castillo A. 2012. La Flora de la Península de Yucatán Mexicana: 250 años de conocimiento florístico. Biodiversitas, 101: 6-10. [ Links ]

Carrillo-Bastos A, Islebe GA, Torrescano-Valle N. 2012. Geospatial analysis of pollen records from the Yucatan peninsula, Mexico. Vegetation History Archaeobotany 21: 429-437. DOI: https://doi.org/10.1007/s00334-012-0355-1 [ Links ]

Carrillo-Bastos A, Islebe GA, Torrescano-Valle N, Gonzalez NE. 2010. Holocene vegetation and climate history of central Quintana Roo, Yucatan Penınsula, Mexico. Review of Palaeobotanic and Palynology 160: 189-196. DOI: https://doi.org/10.1016/j.revpalbo.2010.02.013 [ Links ]

Clement M, Snell Q, Walker P, Posada D, Crandall K. 2002. TCS: Estimating gene genealogies. Parallel and Distributed Processing Symposium, International Proceedings 2: 184. DOI: https://doi.org/10.1109/ipdps.2002.1016585 [ Links ]

Collins WD, Bitz CM, Blackmon ML, Bonan GB, Bretherton CS, Carton JA , Chang P, Scott CD, James JH, Henderson TB, Kiehl JT, Large WG, McKenna DS, Santer BD and Smith RD. 2004. The community climate system model: CCSM3. Journal of Climate 19: 2122-2143. DOI: https://doi.org/10.1175/jcli3761.1 [ Links ]

Contreras-Medina R, Luna-Vega I. 2002. On the distribution of gymnosperm genera, their areas of endemism and cladistic biogeography. Australian Systematic Botany 15: 193-203. DOI: https://doi.org/10.1071/sb00004 [ Links ]

Cornejo-Romero A, Vargas-Mendoza CF, Aguilar-Martínez GF, Medina-Sánchez J, Rendón-Aguilar B, Valverde PL, Zavala-Hurtado JA, Serrato A, Rivas-Arancibia S, Pérez- Hernández MA, López-Ortega G, Jiménez-Sierra C. 2017. Alternative glacial-interglacial refugia demography hypotheses tested on Cephalocereus columna-trajani (Cactaceae) in the intertropical Mexican drylands. PLOS ONE 12: e0175905. DOI: https://doi.org/10.1371/journal.pone.0175905 [ Links ]

Doyle JJ, Doyle JL, Doyle JA, Doyle FJ. 1987. A rapid DNA isolation procedure from small quantities of fresh leaf tissues. Phytochemical Bulletin 19: 11-15. DOI: [ Links ]

Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC evolutionary biology 7: 214. DOI: https://doi.org/10.1186/1471-2148-7-214 [ Links ]

Excoffier L, Laval G, Schneider S. 2005. Arlequin ver. 3.0: an integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online 1: 47-50. DOI: https://doi.org/10.1177/117693430500100003 [ Links ]

Farrera I, Harrison SP, Prentice IC, Ramstein G, Guiot J, Bartlein PJ, Bonnefille R, Bush M, Cramer W, von Grafenstein U, Holmgren K, Hooghiemstra H, Hope G, Jolly D, Lauritzen SE, Ono Y, Pinot S, Stute M, Yu G. 1999. Tropical climates at the Last Glacial Maximum: a new synthesis of terrestrial palaeoclimate data. I. Vegetation, lake-levels and geochemistry. Climate Dynamics 15: 823-856. DOI: https://doi.org/10.1007/s003820050317 [ Links ]

Fu XY. 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915-925. [ Links ]

González-Astorga J, Vovides AP, Cabrera-Toledo D, Nicolalde-Morejón F. 2008. Diversity and genetic structure of the endangered Cycad Dioon sonorense (Zamiaceae) from Sonora, Mexico: evolutionary and conservation implications. Biochemical Systematics and Ecology 36: 891-899. DOI: https://doi.org/10.1016/j.bse.2008.11.006 [ Links ]

González-Astorga J, Vovides AP, Octavio-Aguilar P, Aguirre-Fey D, Nicolalde-Morejón F, Iglesias C. 2006. Genetic diversity and structure of the Cycad Zamia loddigesii Miq. (Zamiaceae): implications for evolution and conservation. Botanical Journal of the Linnean Society 152: 533-544. DOI: https://doi.org/10.1111/j.1095-8339.2006.00579.x [ Links ]

González-Astorga J, Vovides AP, Ferrer MM, Iglesias C. 2003. Population genetics of Dioon edule Lindl. (Zamiaceae, Cycadales): biogeographical and evolutionary implications. Biological Journal of the Linnean Society 80: 457-467. DOI: https://doi.org/10.1046/j.1095-8312.2003.00257.x [ Links ]

González D, Vovides AP. 2002. Low intralineage divergence in Ceratozamia (Zamiaceae) detected with nuclear ribosomal DNA ITS and chloroplast DNA trnL-F non-coding region. Systematic Botany 27: 654-661. DOI: https://doi.org/10.1043/0363-6445-27.4.654 [ Links ]

Graham A. 1973. History of the arborescent temperate element in the northern Latin American biota. In: Graham A, ed. Vegetation and vegetational history of northern Latin America. New York: Elsevier Scientific Publishing Co., pp. 301-314. ISBN-10: 0444410562; ISBN-13: 978-0444410566 [ Links ]

Gutiérrez-Ortega JS, Kajita T, Molina-Freaner FE. 2014. Conservation genetic of endangered Cycad, Dioon sonorense (Zamiaceae): Implication for variation of chloroplast DNA. Botanical Sciences 92: 441-451. DOI: https://doi.org/10.17129/botsci.112 [ Links ]

Haffer J. 1969. Speciation in Amazonian Forest Birds. Science 165: 131-136. https://doi.org/10.1126/science.165.3889.131 [ Links ]

Hammer Ø, Harper DAT, Ryan PD. 2001. Past: Paleontological Statistics Software Package for Education and Data Analysis. Palaeontologia Electronica 4: 1-9. [ Links ]

Hasumi H, Emori S. 2004. K-1 coupled GCM (MIROC) description. Center for Climate System Japan, Tokyo: Research, University of Tokyo. [ Links ]

Harrison SP, Bartlein PJ, Prentice IC. 2016. What have we learnt from palaeoclimate simulations? Journal of Quaternary Science, 31: 363-385. https://doi.org/10.1002/jqs.2842 [ Links ]

Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25: 1965-1978. DOI: https://doi.org/10.1002/joc.1276 [ Links ]

Hodell DA, Anselmetti FS, Ariztegui D, Brenner M, Curtis JH, Gilli A, Grzesik DA, Guilderson TJ, Müller AD, Bush MB, Correa-Metrio YA, Escobar J, Kutterolf S. 2008. An 85-ka Record of Climate Change in Lowland Central America. Quaternary Science Reviews 27: 1152-1165. DOI: https://doi.org/10.1016/j.quascirev.2008.02.008 [ Links ]

Kutzbach JE, Guetter PJ, Behling PJ, Selin R. 1993. Simulated climatic changes: results of the COHMAP climate-model experiments. In: Wright HE, Kutzbeach JE, Webb T, Ruddiman WF, StreetPerrott FA, Bartlein PJ, eds. Global Climates Since the Last Glacial Maximum. Minneapolis, Minnesota, EUA: University of Minnesota Press, pp. 24-93. ISBN-10: 0816621454: ISBN-13: 978-0816621453 [ Links ]

Lahaye R, Savolainen V, Duthoit S, Maurin O, Van der Bank M. 2008. A test of psbK-psbI and atpF-atpH as potential plant DNA barcodes using the flora of the Kruger National Park as a model system (South Africa). Nature Precedings. DOI: https://doi.org/10.1038/npre.2008.1896.1 [ Links ]

Lammens MEA, Boria RA, Radosavljevic A, Vilela B, Anderson RP. 2015. spThin: An R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography 38: 541-545. DOI: https://doi.org/10.1111/ecog.01132 [ Links ]

Lavery S, Moritz C, Fielder D. 1996. Genetic patterns suggest exponential population growth in a declining species. Molecular Biology and Evolution 13: 1106-113. DOI: https://doi.org/10.1093/oxfordjournals.molbev.a025672 [ Links ]

Leyden BW. 2002. Pollen evidence for climatic variability and cultural disturbance in the Maya lowlands. Ancient Mesoamerica 13: 85-101. DOI: https://doi.org/10.1017/s0956536102131099 [ Links ]

Librado P, Rozas J. 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451-1452. DOI: https://doi.org/10.1093/bioinformatics/btp187 [ Links ]

Limón F, González-Astorga J, Nicolalde-Morejón F, Roger-Guevara. 2016. Phenotypic variation of Zamia loddigesii Miq. and Z. prasina W. Bull. (Zamiaceae, Cycadales): the effect of environmental heterogeneity. Plant Systematics and Evolution 302: 1395-1404. DOI: https://doi.org/10.1007/s00606-016-1338-y [ Links ]

LoveLess MD, Hamrick JL. 1984. Ecological Determinants of Genetic Structure in Plant Populations. Annual Review of Ecology, Evolution and Systematics 15: 65-95. DOI: https://doi.org/10.1146/annurev.es.15.110184.000433 [ Links ]

Lozano-García MS, Vázquez-Selem L. 2005. A High Elevation Holocene Pollen Record from Iztaccihuatl volcano, Central México. The Holocene 15: 329-338. DOI: https://doi.org/10.1191/0959683605hl814rp [ Links ]

Macrogen. 2018. Macrogen Inc. < http://dna.macrogen.com > (accesssed Enero, 25, 2018). [ Links ]

Manni F, Guerard E, Heyer E. 2004. Geographic patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by using Monmoniers algorithm. Human Biology 76: 173-190. DOI: https://doi.org/10.1353/hub.2004.0034 [ Links ]

Metcalfe SE, O’Hara SL, Caballero M, Davies SJ. 2000. Records of late Pleistocene-Holocene climatic change in Mexico a review. Quaternary Science Reviews 19: 699-721. DOI: https://doi.org/10.1016/s0277-3791(99)00022-0 [ Links ]

Miller MA, Pfeiffer W, Schwartz T. 2010. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Proceedings of the Gateway Computing Environments Workshop 1-8. DOI: https://doi.org/10.1109/gce.2010.5676129 [ Links ]

Müller J, Müller K, Neinhuis C, Quandt D. 2010. PhyDE-Phylogenetic data editor. Version 0.9971. Program distributed by the author. < http://www.phyde.de > (accessed 15, September 2017). [ Links ]

Myers N. 1982. Biological Diversification in the Tropics. Environmental Conservation 10: 277-278. DOI: https://doi.org/10.1017/S037689290001290X [ Links ]

Nagalingum NS, Marshall CR, Quental TB, Rai HS, Little DP, Mathews S. 2011. Recent synchronous radiation of a living fossil. Science 334: 796-799. DOI: https://doi.org/10.1126/science.1209926 [ Links ]

Nicolalde-Morejón F, Vergara-Silva F, González-Astorga J, Stevenson DW, Vovides AP, Sosa V. 2011. A character-based approach in the Mexican cycads supports diverse multigene combinations for DNA barcoding. Cladistics 27: 150-164. DOI: https://doi.org/10.1111/j.1096-0031.2010.00321.x [ Links ]

Nolasco-Soto J, González-Astorga J, Nicolalde-Morejón F, Vergara-Silva F, Espinosa de los Monteros A, Medina-Villarreal A. 2015. Phylogeography and demography history of Zamia paucijuga Wieland (Zamiaceae), a cycad species from the Mexican Pacific slope. Plant Systematics and Evolution 301: 623-637. DOI: https://doi.org/10.1007/s00606-014-1101-1 [ Links ]

Orellana R, Islebe G, Espadas C. 2003. Presente, pasado y futuro de los climas de la Península de Yucatán. In: Colunga García Marín P, Larqué-Saavedra A, eds. Naturaleza y sociedad en el área maya, pasado, presente y futuro. Mexico City: Academia Mexicana de Ciencias, CICY 37-52. [ Links ]

Orellana R, Espada C. 2009. < http://www.ccpy.gob.mx/pdf/Regional/escenarios-cambio-climatico/precipitacion_total.pdf >. (accessed April 24, 2018). [ Links ]

Ornelas JF, Sosa V, Soltis DE, Daza JM, González C, Soltis PS, Gutiérrez-Rodríguez C, Espinosa de los Monteros A, Castoe TA, Bell C, Ruiz-Sánchez E. 2013. Comparative phylogeographic analyses illustrate the complex evolutionary history of threatened cloud forests of northern Mesoamerica. PLOS ONE 8: e56283. DOI: https://doi.org/10.1371/journal.pone.0056283 [ Links ]

Otto-Bliesner BL, Hewitt CD, Marchitto TM, Brady E, Abe-Ouchi A, Crucifix M, Murakami S, Weber SL. 2007. Last Glacial Maximum ocean thermohaline circulation: PMIP2 model intercomparisons and data constraints. Geophysical Research Letters 34: L12706. DOI: https://doi.org/10.1029/2007gl029475 [ Links ]

Pennington RT, Prado DE, Pendry CA. 2000. Neotropical seasonally dry forests and Quaternary vegetation changes. Journal of Biogeography 27: 261-273. DOI: https://doi.org/10.1046/j.1365-2699.2000.00397.x [ Links ]

Phillips SJ, Anderson RP, Schapire RE. 2006. Maximum entropy modeling of species geographic distributions. Ecological Modelling 190: 231-259. DOI: https://doi.org/10.1016/j.ecolmodel.2005.03.026 . [ Links ]

Pons O, Petit RJ. 1996. Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics 144: 1237-1245. [ Links ]

Posada D. 2008. jModelTest: Phylogenetic Model Averaging. Molecular Biology and Evolution 25: 1253-1256. DOI: https://doi.org/10.1093/molbev/msn083 [ Links ]

Rambaut A. 2009. < http://tree.bio.ed.ac.uk/software/figtree/ > (accessed October 12, 2017). [ Links ]

Ramírez-Barahona S, Eguiarte LE. 2013. The role of glacial cycles in promoting genetic diversity in the Neotropics: the case of cloud forests during the Last Glacial Maximum. Ecology Evolution 3: 725-738. DOI: https://doi.org/10.1002/ece3.483 [ Links ]

Ramos-Onsins. Rozas J. 2010. Statistical properties of new neutrality tests against population growth. Molecular Biology and Evolution 12: 2092-2100. DOI: https://doi.org/10.1093/oxfordjournals.molbev.a004034 [ Links ]

Rogers AR, Harpending H. 1992. Population growth makes waves in the distribution of pairwise differences. Molecular Biology and Evolution 9: 552-569. DOI: https://doi.org/10.1093/oxfordjournals.molbev.a040727 [ Links ]

Rull V. 2012. Cycad diversification and tropical biodiversity. Collectanea Botanica 31: 103-106. DOI: https://doi.org/10.3989/collectbot.2012.v31.008 [ Links ]

Sánchez-Sánchez O, Islebe GA. 2002. Tropical forest communities in southeastern Mexico. Plant Ecology 158: 183-200. https://doi.org/10.1023/A:1015509832734 [ Links ]

Stevenson DW, Moretti A, Gaudio l. 1996. A new species of Zamia (Zamiaceae) from Belize and the Yucatan Peninsula of Mexico. Delpinoa 37-38: 3-8. [ Links ]

Soltis DE, Morris AB, McLachlan JS, Manos PS, Soltis PS. 2006. Comparative phylogeography of unglaciated eastern North America. Molecular Ecology 15: 4261-4293. DOI: https://doi.org/10.1111/j.1365-294X.2006.03061.x [ Links ]

Stewart JR, Lister AM, Barnes I, Dalén L. 2010. Refugia revisited: individualistic responses of species in space and time. Proceedings of the Royal Society B: Biological Sciences 277: 661-671. DOI: https://doi.org/10.1098/rspb.2009.1272 [ Links ]

Svensson A, Nielsen SW, Kipfstuhl S, Johnsen SJ, Steffensen JP, Bigler M, Ruth U, Röthlisberger R. 2005. Visual stratigraphy of the North Greenland Ice Core Project (NorthGRIP) ice core during the last glacial period. Journal of Geophysical Research 110: (D02108). DOI: https://doi.org/10.1029/2004jd005134 [ Links ]

Tajima F. 1989. The effect of change in population size on DNA polymorphism. Genetics 123: 597-601. [ Links ]

Taylor KE, Stouffer RJ, Meehl GA. 2012. An overview of CMIP5 and the experiment design. Bulletin of the American Meteorological Society 93: 485-498. DOI: https://doi.org/10.1175/bams-d-11-00094.1 [ Links ]

TRACER 1.6. < http://tree. bio.ed.ac.uk/software/tracer/ > (accesssed Diciembre, 12, 2017). [ Links ]

Toledo VM. 1982. Pleistocene changes of vegetation in tropical Mexico. In: Prance GT, ed. Biological Diversification in the Tropics: Proceedings of the Fifth International Symposium of the Association for Tropical Biology, New York: Columbia University Press. 93-111. ISBN: 9780231048767 [ Links ]

van der Hammen Th. 1961. The Quaternary climatic changes of northern South America. Annals of the New York Academy of Sciences 95: 676-683. DOI: https://doi.org/10.1111/j.1749-6632.1961.tb50066.x [ Links ]

Vela-Pelaez AA, Torrescano-Valle N, Islebe GA, Mas JF, Weissenberger H. 2018. Holocene precipitation changes in the Maya forest, Yucatán peninsula, Mexico. Palaeogeography, Palaeoclimatology, Palaeoecology 505: 42-52. DOI: https://doi.org/10.1016/j.palaeo.2018.05.024 [ Links ]

Velasco-Martínez Y. 2015. Variación Genética de Zamia prasina W. BULL 1881 (Zamiaceae). BSc. Thesis. Universidad de Ciencias y Artes de Chiapas. [ Links ]

Vovides AP, Nicolalde-Morejón F. 2010. Ficha técnica de Zamia polymorpha. In: Vovides AP. (comp.) <) http://www.conabio.gob.mx/conocimiento/ise/fichasnom/Zamia%20polymorpha.pdf > (accesssed Marzo, 10, 2018). [ Links ]

Vovides AP, Olivares M. 1996. Karyotype polymorphism in the cycad Zamia loddigesii (Zamiaceae) of the Yucatán Peninsula, Mexico. Botanical Journal of the Linnean Society 120: 77-83. DOI: https://doi.org/10.1111/j.1095-8339.1996.tb00481.x [ Links ]

Walker MJ C, Berkelhammer M, Bjorck S, Cwynar LC, Fisher DA, Long AJ, Lowe JJ, Newnham RM, Rasmussen SO, Weiss H. 2012. Formal subdivision of the Holocene Series/Epoch: a Discussion Paper by a Working Group of Intimate (Integration of ice-core, marine and terrestrial records) and the Subcommission on Quaternary Stratigraphy (International Commission on Stratigraphy). Journal of Quaternary Science 27: 649-659. DOI: https://doi.org/10.1002/jqs.2565 [ Links ]

Wolf KH, Li WH, Sharp PM. 1987. Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear DNAs. Proceedings of the National Academy of Sciences of the United States of America 84: 9054-9058. DOI: https://doi.org/10.1073/pnas.84.24.9054 [ Links ]

1Associated editor: Monserrat Vázquez Sánchez

Appendix 1

Climatic variables per population 

Poblaciones Localidades bio02 bio03 bio04 bio05 bio07 bio11 bio12 bio15 bio17 bio18 bio19
1 Valladolid 130 68 1899 347 190 227 1220 65 103 452 103
2 Peto 135 68 1914 356 196 231 1109 62 95 296 111
3 Becanchén 134 67 1988 356 198 228 1071 62 92 283 103
4 Rancho Duarte 133 67 2011 353 197 226 1070 61 91 286 101
5 Kaxil Kiuic 135 66 2196 360 204 227 1029 65 81 277 89
6 Champotón 132 65 2222 357 203 226 1223 75 66 280 99
7 Calakmul 117 61 2198 344 189 221 1171 63 89 294 127
8 Escárcega 127 65 2111 357 194 231 1410 73 81 320 114
9 Hormiguero 112 61 2179 334 182 215 1171 57 98 322 142
10 Viergencita 107 62 1986 350 170 237 1624 68 104 352 169
11 Tenosique 110 63 2009 353 174 236 2169 54 206 473 328
12 Palenque 109 60 1981 357 181 238 2443 54 246 649 468
13 Puerto Morelos 100 67 1734 327 149 231 1152 50 129 354 185
14 Cobá 117 69 1755 336 169 229 1151 53 132 404 139
15 Carrillo Puerto 114 69 1787 336 165 232 1290 53 137 481 167
16 Xhazil 111 68 1733 334 161 233 1255 52 135 434 167
17 José María Morelos 110 66 1786 337 166 232 1087 51 118 351 146
18 Pedro A Santos 93 68 1552 321 136 235 1467 52 144 457 226
19 Bacalar 92 68 1528 320 135 236 1445 53 138 530 215
20 Panto-Ha 95 66 1593 324 142 234 1243 54 112 458 167
21 Nachicocom 104 65 1763 329 160 227 1338 59 113 512 166
22 Belize-Belize 85 59 1842 306 142 214 1510 47 155 318 241
23 Belize-Cayo 80 63 1553 310 126 230 2207 51 197 773 364

Appendix 2

Genetic diversity of the populations of Zamia prasina in the Yucatan Peninsula Biotic Province (YPBP). Biotic Province (YPBP) using the chloroplast DNA region atpF-atpH. Number of individuals (N), number of haplotypes (k), haplotype diversity (h) and nucleotide diversity (π). 

Population/state N k h π
Valladolid/Yucatán 5 1 0.000 0.00000
Peto/Yucatán 4 2 0.500 0.00110
Becanchén/Yucatán 5 2 0.400 0.00088
Rancho duarte/Yucatán 4 2 0.500 0.00110
Kaxil Kiuic/Yucatán 3 1 0.000 0.00000
Champotón/Campeche 5 1 0.000 0.00000
Escarcega/Campeche 5 1 0.000 0.00000
Virgencita/Campeche 5 4 0.900 0.00574
Calakmul/Campeche 5 1 0.000 0.00000
Hormiguero/Campeche 3 3 1.000 0.01177
Palenque/Chiapa 5 1 0.000 0.00000
Tenosique/Tabasco 5 1 0.000 0.00000
Xhazil/Quintana Roo 5 2 0.600 0.00132
Pedro A. Santos/Quintana Roo 5 2 0.400 0.00088
Panto Ha/Quintana Roo 5 2 0.600 0.00132
Nachi cocom/Quintana Roo 4 2 0.500 0.00110
Puerto Morelos/Quintana Roo 5 1 0.000 0.00000
Cobá/Quintana Roo 5 2 0.600 0.00132
Carrillo Puerto/Quintana Roo 4 2 0.500 0.00110
Bacalar/Quintana Roo 5 1 0.000 0.00000
José María Morelos/Quintana Roo 5 2 0.400 0.00088
Belize/Belize 2 2 1.000 0.00221
Belize/Cayo 3 2 0.666 0.00147
Total 102 9 0.564 0.00186

Appendix 3

Genetic diversity of populations of Zamia prasina in the Yucatan Peninsula Biotic Province (YPBP). Biotic Province (YPBP) using the nuclear DNA region ITS2. Number of individuals (N), number of haplotypes (k), haplotype diversity (h) and nucleotide diversity (π). 

Population/state N k h π
Valladolid/Yucatán 5 2 0.400 0.00135
Peto/Yucatán 5 1 0.000 0.00000
Becanchén/Yucatán 4 4 1.000 0.00617
Rancho duarte/Yucatán 5 2 0.600 0.00202
Kaxil Kiuic/Yucatán 3 2 0.666 0.00224
Champotón/Campeche 5 2 0.400 0.00135
Escárcega/Campeche 5 3 0.700 0.00337
Virgencita/Campeche 3 1 0.000 0.00000
Calakmul/Campeche 4 2 0.500 0.00168
Hormiguero/Campeche 3 1 0.000 0.00000
Palenque/Chiapa 5 4 0.900 0.00404
Tenosique/Tabasco 4 3 0.833 0.00505
Xhazil/Quintana Roo 5 2 0.600 0.00202
Pedro A. Santos/Quintana Roo 5 3 0.700 0.00269
Panto Ha/Quintana Roo 5 3 0.700 0.00337
Nachi co com/Quintana Roo 4 2 0.500 0.00168
Puerto Morelos/Quintana Roo 3 2 0.666 0.00224
Cobá/Quintana Roo 3 2 0.666 0.00224
Carrillo Puerto/Quintana Roo 4 4 1.000 0.04433
José M Morelos/Quintana Roo 5 5 1.000 0.01751
Belice/Belice 2 2 1.000 0.00337
Belice/Cayo 3 1 0.000 0.00000
Total 90 18 0.827 0.00852

Appendix 4

List of GenBank accession for haplotypes. 

Genome Code GenBank accession
Chloroplast (atpF-atpH) ZP-A1 MK513628
ZP-A2 MK513629
ZP-A3 MK513630
ZP-A4 MK513631
ZP-A5 MK513632
ZP-A6 MK513633
ZP-A7 MK513634
ZP-A8 MK513635
ZP-A9 MK513636
Nuclear (ITS2) ZP-H1 MK513637
ZP-H2 MK513638
ZP-H3 MK513639
ZP-H4 MK513640
ZP-H5 MK513641
ZP-H6 MK513642
ZP-H7 MK513643
ZP-H8 MK513644
ZP-H9 MK513645
ZP-H10 MK513646
ZP-H11 MK513647
ZP-H12 MK513648
ZP-H13 MK513649
ZP-H14 MK513650
ZP-H15 MK513651
ZP-H16 MK513652
ZP-H17 MK513653
ZP-H18 MK513654

Appendix 5

Principal component analysis of 11 climatic variables for 23 populations of Zamia prasina in the Biotic Province of the Yucatán Peninsula (BPYP) showing the presence of two groups of populations. Red triangles: eastern populations, green triangles: western group. 

Appendix 6

Bayesian chronogram for chloroplast sequence atpF-atpH. The divergence time of Zamia prasina haplotypes in the YPBP and other cycad species is shown. The 95 % confidence intervals are shown with purple bars; the numbers at the nodes indicate the estimated age, the numbers below the branches indicate the posterior probability. 

Appendix 7

Bayesian chronogram for nuclear sequence ITS2. The divergence time of Zamia prasina haplotypes in the YPBP and other cycad species is shown. The 95 % confidence intervals are shown with purple bars; the numbers at the nodes indicate the estimated age, the numbers below the branches indicate the posterior probability. 

Received: February 13, 2019; Accepted: July 24, 2019

* Corresponding author: Jaime Martínez-Castillo, e-mail: jmartinez@cicy.mx.

Author Contributions: GMF (https://orcid.org/0000-0002-0627-3347). Collected most of the botanical material, conceived, designed, and performed the experiments, analyzed the data, and drafted the first version of the paper as well as reviewed drafts of the manuscript. LFST Drafted the first version of the paper. Funded the sequences analyses. GCFC Conceived and designed parts of the experiments and analyses. Drafted the first version of the paper as well as reviewed all drafts of the manuscript. APV Drafted the first version of the paper as well as reviewed all drafts of the manuscript. RGL Carried out the analysis of potencial distribution with the tool of Ecological Niche Modelling. MMOG Participated in laboratory work, and reviewed drafts of the first version. JAML Participated in the botanical collections, and reviewed drafts of the first version. JMC (https://orcid.org/0000-0002-6484-6406). Conceived and designed parts of the experiments and reviewed all drafts of the manuscript. Funded the sequences analysis and the collecting trips.

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License