Lay mistletoes on the Yucatán Peninsula: post-glacial expansion and genetic differentiation of Psittacanthus mayanus (Loranthaceae)

Lay mistletoes on the Yucatán Peninsula: post-glacial expansion and genetic differentiation of... Abstract We evaluate phylogeographical patterns of Psittacanthus mayanus (Loranthaceae), a widely distributed mistletoe species on the Yucatán Peninsula and Chiapas, using nuclear and plastid markers. An analysis of phylogeographic population structure and demography and Bayesian inference methods were conducted on P. mayanus from 16 localities across the range of the species. To assess the historical processes and changes through the Pleistocene climatic cycles that may have affected the distribution and demographic history of the species, the current potential distribution of the species was modelled using ecological niche modelling (ENM) and projected onto the mid Holocene (MH), Last Glacial Maximum (LGM) and Last Inter Glacial (LIG). Psittacanthus mayanus exhibited higher haplotype and nucleotide diversity in the southern part of its distribution than in the northern part. Two divergent lineages were revealed in the phylogenetic and phylogeographical analyses of ribotypes and haplotypes. Populations from the northern portion of the Yucatán Peninsula probably experienced a recent population growth, whereas populations from the southern portion of the Peninsula and Chiapas are marked by historical demographic stability and isolation. Approximate Bayesian computation (ABC) analyses strongly supported a scenario of post-glacial population growth. ENM results indicate that the distribution of P. mayanus expanded and was connected into the southern portion of the Yucatán Peninsula during LGM conditions and colonized the northern portion of the Peninsula and fragmented during MH conditions. According to the ABC and ENM results, the genetic differentiation of P. mayanus may be linked to the effects of the glacial/interglacial cycles and environmental factors, in which populations of P. mayanus expanded during the LGM into the Petén province currently covered with semideciduous tropical rain forest followed by northward expansion, southwards contractions and post-glacial colonization of the northernmost portion of the Yucatán province, a region currently covered with seasonally dry tropical deciduous forest. Our results highlight the influence of Pleistocene events in shaping the geographical distribution of genetic variation in Neotropical lowland forest. The phylogeographic and environmental patterns in P. mayanus provide an opportunity to investigate further the evolution of Mexican lowland forest biodiversity. INTRODUCTION The Yucatán Peninsula has witnessed numerous marine transgressions over the last 65 Myr, submerging under warm tropical waters what is now the Peninsula (Vázquez-Domínguez & Arita, 2010). Today, the Yucatán Peninsula occupies almost half the portion of the peninsula platform, including the Mexican states of Yucatán, Quintana Roo, Campeche and adjacent portions of Chiapas and Tabasco, and the northern portions of Guatemala and Belize (c. 181 200 km2; Lundell, 1934; Espadas-Manrique, Durán & Argáez, 2003). Submerged for millions of years, the Peninsula is basically a large limestone slab that emerged slowly from south to north. The ecological features and biogeography of the Yucatán Peninsula show the climatic effects of its 65-Myr history, with a rain gradient from humid in the south/south-east to dry in the north/north-west (Espadas-Manrique et al., 2003; Orellana, Islebe & Espadas, 2003; Torrescano-Valle & Islebe, 2015). The vegetation also follows the south-east to north-west rain gradient, from tropical rainforests in the Petén region to tropical scrubland in the extreme north-west part of the Yucatán peninsula, with extensive areas covered with deciduous or semideciduous tropical forests between these two extremes, and a decrease in species diversity from the humid communities in the south to the northern counterparts (Simpson’s peninsula effect; Simpson, 1964; Ramírez-Barahona et al., 2009; Vázquez-Domínguez & Arita, 2010). This distinctiveness has prompted most biogeographers to consider the Yucatán Peninsula a region with two provinces, the Petén province in the south and the dry Yucatán province in the north (e.g. Espadas-Manrique et al., 2003; Morrone, 2005; Ramírez-Barahona et al., 2009; Vázquez-Domínguez & Arita, 2010). Palaeogeographic data indicate that the Yucatán Peninsula was emergent until the Triassic–Jurassic, with several marine transgressions having occurred to the Pliocene–Miocene, and the present shape of the peninsula was completed near the end of the Pliocene and continued its formation into the Quaternary (López Ramos, 1975; Padilla y Sánchez, 2007). During most of the Pleistocene, the climate of the Yucatán Peninsula was dry and 6ºC cooler than current temperatures. Savanna and Juniperus-matorral covered most of the region until the early Holocene (Islebe et al., 1996; Orellana et al., 2003), in which the vegetation of the Peninsula changed to a more humid and tropical-type forest (Brosimum Sw.) established by 6800 years before present (BP) followed by a dry period between 6000 and 5000 BP (Leyden, 1984; Metcalfe et al., 2000). Wetter conditions were then re-established until c. 3500 BP, to be followed by drying and a more open canopy forest and more grasses until c. 1500 BP (Metcalfe et al., 2000). This dry period corresponds to the time of the ‘collapse’ of the lowland Maya civilization (Metcalfe et al., 2000; Haug et al., 2003; Vázquez-Domínguez & Arita, 2010). The biogeography of the region is not only influenced by the effects of the formation and emergence of the Yucatán Peninsula and the effects of Pleistocene climate changes but also by those related to the geology of and tectonic activity on the Maya block, between the biogeographic breaks of the Isthmus of Tehuantepec and the Motagua–Polochic–Jocotán fault system (Gutiérrez-García & Vázquez-Domínguez, 2012). In the Maya block, the most conspicuous and dynamic mountain ranges include the modern Chiapas volcanic arc (Altos de Chiapas) and the Chiapas massif (Sierra Madre de Chiapas) separated by the Central Depression of Chiapas and the Maya mountains of Belize (Gutiérrez-García & Vázquez-Domínguez, 2012). Although these geological and topographical features and the isolation of these mountain ranges certainly influenced the evolutionary history of Central American lowland tree species (Cavers, Navarro & Lowe, 2003; Poelchau & Hamrick, 2012, 2013a) and widespread Mesoamerican vertebrates (e.g. Guevara-Chumacero et al., 2010; González, Ornelas & Gutiérrez-Rodríguez, 2011; Smith et al., 2011; Gutiérrez-García & Vázquez-Domínguez, 2012; Rodríguez-Gómez & Ornelas, 2014; Jiménez & Ornelas, 2016; Williford et al., 2016), few studies have included phylogeographies of terrestrial species restricted to the Maya block (in the Yucatán Peninsula). Gutiérrez-García & Vázquez-Domínguez (2012) described patterns of historical divergence and population genetic differentiation of mitochondrial DNA in the rodent Ototylomys phyllotis distributed from the Isthmus of Tehuantepec (Mexico) to Costa Rica. They found that the Yucatán haplotypes formed the most derived lineage and hypothesized that O. phyllotis arrived to Mexico from the south (Central America). The colonization of the Yucatán Peninsula occurred c. 0.8 million years ago (Mya), from individuals that persisted in Pleistocene refugia on the Maya mountains of Belize that expanded into the Yucatán Peninsula. More recently, Oyama et al. (2016) assessed genetic differentiation of Guaiacum sanctum L. (Zygophyllaceae) trees on the Yucatán Peninsula, and found evidence of fragmentation and little genetic differentiation between two groups of populations: northern populations of Yucatán on porous karst terrain soils and southern populations from Campeche in seasonally flooded terrains. In this study, the demographic history, population structure and genetic differentiation of Psittacanthus mayanus Standl. & Steyerm. (Loranthaceae) populations are investigated using nuclear DNA (nrDNA) and plastid DNA sequences. We also used ecological niche modelling, isolation-with-migration and Bayesian inference methods to provide insight into the structuring of genetic variation of these populations. Specifically, we asked the following two questions. (1) What geographical features and climatic oscillations in the Pleistocene affected the genetic diversity and population genetic structure of P. mayanus? (2) Based on the observed genetic diversity, did populations expand across the flat Yucatán Peninsula and diverge according to the rain gradient? The distribution of P. mayanus populations in different habitats across the peninsula (Petén and Yucatán provinces) and Chiapas facilitates the testing of hypotheses on the relative importance of divergence and gene flow in the evolutionary process and colonization history of this species, in which both physical and ecological barriers might limit gene flow. Given that Quaternary climate fluctuations, particularly between now and the Last Glacial Maximum (LGM), may have driven distributional shifts in tree species of the Yucatán Peninsula, we predict dynamic shifts in the reconstructed distribution of P. mayanus along the Yucatán Peninsula as seasonally dry tropical deciduous forest contracted and expanded many times during the Pleistocene (Metcalfe et al., 2000). In this scenario, we assume that (1) the populations in Chiapas and those in the Petén and Yucatán provinces experienced cycling periods of isolation during LGM contraction of the seasonally dry tropical deciduous forest followed by secondary contact and gene flow during the warmer interglacial cycles; (2) populations of P. mayanus in the seasonally dry tropical deciduous forest of the Yucatán province contracted southwards into refugia during the LGM in the semideciduous tropical rain forest of the Petén province followed by northward expansion and recolonization of the seasonally dry tropical deciduous forest of the Yucatán province during warmer periods. The absence of shifts in suitable habitat between the LGM and the present day for several tree species of the seasonally dry tropical deciduous forest (Cavers et al., 2003; Poelchau & Hamrick, 2013a, b) would support an alternative hypothesis of long-term persistence and stability in climatically mediated distributions of P. mayanus with floristic affinities to seasonally dry tropical deciduous forest (Kuijt, 2009). This hypothesis assumes that the Petén and Yucatán populations of P. mayanus experienced genetic divergence in isolation. MATERIAL AND METHODS Study Organism Psittacanthus mayanus commonly grows in seasonally dry deciduous and semideciduous forests to seasonally dry deciduous forest with columnar cacti. The species frequently parasitizes leguminous trees (Fabaceae), mainly Lysiloma Benth., Lonchocarpus Kunth, Acacia Mill., Pithecellobium Mart., Inga Mill. and Haematoxylum Gronov., and less frequently Coccoloba L. (Polygonaceae), Conocarpus L. (Combretaceae) Crescentia L. (Bignoniaceae) and Bursera Jacq. ex L. (Burseraceae) trees (Kuijt, 2009). It is distributed in Mexico, Guatemala and Belize, mostly on the Yucatán Peninsula and Chiapas, extending somewhat to the north and south from almost sea level up to 500 m a.s.l. (Kuijt, 2009). Psittacanthus mayanus is easily recognizable from other Psittacanthus Mart. in the region by the presence of prophylls (two minute foliar organs, transversely placed in a leaf axil), also detected in P. elegans Kuijt and conspicuous in P. breedlovei Kuijt and P. sonorae (Watson) Kuijt, four or five cotyledons, non-foliaceous lowest triad bracts and straight, thinner flowers. In its southern range, however, P. mayanus is sometimes difficult to separate from other Psittacanthus species presumably occurring in the tropical dry forest, namely P. schiedeanus (Cham. & Schltdl.) G.Don, P. calyculatus G.Don, P. breedlovei Kuijt and P. minor Kuijt (Kuijt, 2009). However, P. breedlovei has extremely narrow, linear-lanceolate leaves, <1 cm wide (P. mayanus has leaves >2 cm wide, ovate, ovate-lanceolate) and P. calyculatus is a robust plant with denser inflorescences and buds to 6 cm long (P. mayanus is a slender plant with few flowers and buds to 3 cm long); P. minor is endemic to Nicaragua (Kuijt, 2009). The few-flowered inflorescences of P. mayanus produce red-orange flowers mainly visited by Campylopterus curvipennis pampa and Amazilia hummingbirds and its ripe, purplish-black fruits are consumed by several bird species, mainly Patagioenas cayennensis, Myiozetetes similis, Pitangus sulphuratus, Megarhynchus pitangua and Piranga roseogularis (Howell, 1972; A. Vargas-Varguéz, unpubl. data). Population Sampling Leaf tissue samples were collected from 87 individuals in 16 localities across most of the species range of P. mayanus (Table 1). Samples from Quintana Roo, Belize and Guatemala were not available for the present study. Target sampling localities were chosen based on localities of occurrence data acquired from Kuijt (2009). Only one plant of P. mayanus was sampled per individual host tree, with care taken to collect only samples from adult flowering individuals and individuals present in the same location were treated as a population, ranging from one to 14 individuals per population. Most populations collected have accompanying vouchers that were deposited at the XAL herbarium of the Instituto de Ecología, A.C. (INECOL; for voucher information, see Table 1). Table 1. Geographical information and sample sizes of the 16 Psittacanthus mayanus localities sampled in this study and voucher information. Pop. no.  Location  G2  G3  GP  n  Elevation (m)  Latitude (N)  Longitude (W)  Voucher (XAL)  1  Yucatán, Dzibilchaltún  YUC  YUC  dryYUC  3  13  21° 05′ 50′′  89° 34′ 55′′  E. Gándara 3131–3133  2  Yucatán, Hunucmá, Sisal  YUC  YUC  dryYUC  14  10  21° 02′ 58′′  89° 54′ 38′′  –  3  Yucatán, Celestún  YUC  YUC  dryYUC  9  8  21° 02′ 11′′  89° 48′ 07′′  E. Gándara 3137–3145  4  Yucatán, Mérida, Cuxtal  YUC  YUC  dryYUC  13  11  20° 54′ 37′′  89° 37′ 15′′  –  5  Yucatán, Maxcanú, Límite Estatal  YUC  YUC  dryYUC  8  11  20° 36′ 31′′  89° 57′ 54′′  Y. Licona-Vera 149–156  6  Campeche, Champotón, Seybaplaya  YUC  PET  dryYUC  10  18  19° 38′ 35′′  90° 39′ 56′′  E. Gándara 3121–3130  7  Campeche, Escárcega, Sabancuy  YUC  PET  humYUC  10  7  18° 54′ 19′′  91° 03′ 31′′  E. Gándara 3111–3120  8  Chiapas, Bochil, Km 205 Pto. Caté  CHI  CHI  humYUC  3  1158  16° 59′ 19′′  92° 52′ 21′′  Y. Licona-Vera 023-025  9  Chiapas, Bochil, Sn Cristobalito  CHI  CHI  humYUC  1  1134  16° 59′ 23′′  92° 52′ 59′′  Y. Licona-Vera 026  10  Chiapas, Sima Las Cotorras  CHI  CHI  dryCHI  6  814  16° 49′ 52′′  93° 27′ 30′′  Y. Licona-Vera 007-012  11  Chiapas, Ocozocuautla, El Aguacero  CHI  CHI  dryCHI  1  602  16° 45′ 50′′  93° 31′ 50′′  Y. Licona-Vera 001  12  Chiapas, Ocozocuautla  CHI  CHI  dryCHI  1  1078  16º 45′ 31′′  93º 19′ 34′′  E. Gándara 3095  13  Chiapas, Arriaga, Km 89–92  CHI  CHI  dryCHI  2  759  16° 44′ 03′′  93° 24′ 47′′  E. Ruiz 259  14  Chiapas, Km 9 Arriaga-Tuxtla  CHI  CHI  dryCHI  1  705  16º 29′ 23′′  93º 50′ 55′′  E. Gándara 3088  15  Chiapas, Arriaga, La Aurora  CHI  CHI  dryCHI  1  699  16° 28′ 18′′  93° 49′ 55′′  E. Gándara 3090  16  Chiapas, Arriaga, Km 28.3  CHI  CHI  dryCHI  6  667  16° 24′ 24′′  93° 48′ 22′′  E. Ruiz 342  Pop. no.  Location  G2  G3  GP  n  Elevation (m)  Latitude (N)  Longitude (W)  Voucher (XAL)  1  Yucatán, Dzibilchaltún  YUC  YUC  dryYUC  3  13  21° 05′ 50′′  89° 34′ 55′′  E. Gándara 3131–3133  2  Yucatán, Hunucmá, Sisal  YUC  YUC  dryYUC  14  10  21° 02′ 58′′  89° 54′ 38′′  –  3  Yucatán, Celestún  YUC  YUC  dryYUC  9  8  21° 02′ 11′′  89° 48′ 07′′  E. Gándara 3137–3145  4  Yucatán, Mérida, Cuxtal  YUC  YUC  dryYUC  13  11  20° 54′ 37′′  89° 37′ 15′′  –  5  Yucatán, Maxcanú, Límite Estatal  YUC  YUC  dryYUC  8  11  20° 36′ 31′′  89° 57′ 54′′  Y. Licona-Vera 149–156  6  Campeche, Champotón, Seybaplaya  YUC  PET  dryYUC  10  18  19° 38′ 35′′  90° 39′ 56′′  E. Gándara 3121–3130  7  Campeche, Escárcega, Sabancuy  YUC  PET  humYUC  10  7  18° 54′ 19′′  91° 03′ 31′′  E. Gándara 3111–3120  8  Chiapas, Bochil, Km 205 Pto. Caté  CHI  CHI  humYUC  3  1158  16° 59′ 19′′  92° 52′ 21′′  Y. Licona-Vera 023-025  9  Chiapas, Bochil, Sn Cristobalito  CHI  CHI  humYUC  1  1134  16° 59′ 23′′  92° 52′ 59′′  Y. Licona-Vera 026  10  Chiapas, Sima Las Cotorras  CHI  CHI  dryCHI  6  814  16° 49′ 52′′  93° 27′ 30′′  Y. Licona-Vera 007-012  11  Chiapas, Ocozocuautla, El Aguacero  CHI  CHI  dryCHI  1  602  16° 45′ 50′′  93° 31′ 50′′  Y. Licona-Vera 001  12  Chiapas, Ocozocuautla  CHI  CHI  dryCHI  1  1078  16º 45′ 31′′  93º 19′ 34′′  E. Gándara 3095  13  Chiapas, Arriaga, Km 89–92  CHI  CHI  dryCHI  2  759  16° 44′ 03′′  93° 24′ 47′′  E. Ruiz 259  14  Chiapas, Km 9 Arriaga-Tuxtla  CHI  CHI  dryCHI  1  705  16º 29′ 23′′  93º 50′ 55′′  E. Gándara 3088  15  Chiapas, Arriaga, La Aurora  CHI  CHI  dryCHI  1  699  16° 28′ 18′′  93° 49′ 55′′  E. Gándara 3090  16  Chiapas, Arriaga, Km 28.3  CHI  CHI  dryCHI  6  667  16° 24′ 24′′  93° 48′ 22′′  E. Ruiz 342  IDs reported below refer to accession numbers in the Instituto de Ecología, AC (XAL) herbarium. G2 = two groups by geography (YUC = Yucatán Peninsula, CHI = Chiapas), G3 = three groups by geography (YUC = Yucatan Province, PET = Petén Province, CHI = Chiapas), GP = three groups by precipitation gradient (dryYUC= dry Yucatán, humYUC = humid Yucatán, dryCHI = dry Central Depression of Chiapas). View Large Leaf tissue samples of P. mayanus were preserved in silica gel desiccant until DNA extractions were performed. We additionally included multiple DNA sequences of Psittacanthus schiedeanus (accession numbers ITS: KU922961–KU923004, KU923036–KU923037; trnL-F: KU923270–KU923278, KU923310–KU923311), P. calyculatus (accession numbers ITS: KX640930–KX640946, KX241619–KX241734, trnL-F: KX640947–KX640960, KX241735–KX241823) and 96 sequences of ten other Psittacanthus spp. (29 samples) and 18 other representatives of Loranthaceae (19 samples) downloaded from the GenBank (accession numbers and voucher information in Appendix 1) to be used as outgroups according to Wilson & Calvin (2006), Amico, Vidal-Russell & Nickrent (2007), Vidal-Russell & Nickrent (2007, 2008a, b), Díaz Infante et al. (2016), Ornelas et al. (2016), and Pérez-Crespo et al. (2017). DNA sequencing The internal transcribed spacer (ITS) region of the nuclear ribosomal DNA (nrDNA) and the trnL-F (cpDNA) intergenic spacer from the plastid genome were amplified by PCR and sequenced for 87 individuals. These markers have shown the appropriate level of variation for mistletoe studies at the species and population levels (e.g. Amico et al., 2012; Ornelas et al., 2016; Pérez-Crespo et al., 2017). As the plastid genome is maternally inherited in most angiosperms, plastid DNA sequences may provide a picture of evolutionary haplotype relationships of seed-mediated gene flow compared to the nrDNA sequences reflecting a combination of pollen and seed movement (Ornelas et al., 2016). Total genomic DNA was extracted from silica-dried material using a modified 2× cetyl trimethyl ammonium bromide (CTAB) protocol (Doyle & Doyle, 1987) or the DNeasy Plant Mini kit (Qiagen, Valencia, CA, USA) using the manufacturer protocol. Amplification of the nrDNA ITS region was conducted with the primers ITS-F2-Psitta (5′-TCGCAGTATGCTCCGTATTG-3′) and ITS-R2-Psitta (5′-TCGTAACAAGGTTTCCGTAGG-3′) designed for this project (Ornelas et al., 2016). For the trnL-F region we used the universal primers ‘e’ and ‘f’ (Taberlet et al., 1991). The 25-µL ITS PCR mix contained 5 µL 5× buffer (Promega, Madison, WI, USA), 2.0 μL MgCl2 (25 mM), 2 µL dNTPs mix (8 mM), 0.82 µL each primer (10 µM), 0.27 µL Taq polymerase (5U/µL) (Promega), 2 µL dimethyl sulphoxide (Sigma), 1.5–3.0 µL template DNA and dH2O to volume. The 25–µL trnL-F PCR mix contained 3.6 µL 5× buffer (Promega, Madison, WI, USA), 3 μL MgCl2 (25 mM), 1.8 µL dNTPs mix (0.57 mM), 0.40 µL each primer (0.16 µM), 0.20 µL Taq polymerase (5U/µL) (Promega), 1.5–5 µL template DNA, and dH2O added to volume. PCRs of ITS consisted of an initial denaturation at 94ºC for 5 min, followed by five cycles of 94ºC for 30 s, 50ºC for 30 s, 72ºC for 1 min, 35 cycles of 94ºC for 30 s, 48–52ºC for 30 s, 72ºC for 1 min and a final step of 72ºC for 5 min. Plastid (trnL-F) amplifications used the following profile: initial denaturation at 95ºC for 5 min, 40 cycles of 95ºC for 10 s, 55–57ºC for 1 min, 72ºC for 20 s and a final step of 72ºC for 4 min. PCR products were purified with the QIAquick kit (Qiagen) and sequenced in both directions to check the validity of the sequence data using the BigDye Terminator Cycle Sequencing kit (Applied Biosystems, Foster City, CA, USA). The products were analysed on a 310 automated DNA sequencer (Applied Biosystems) at the University of Washington High Throughput Genomics Unit, Seattle, Washington. Finally, assembled sequences of the two regions examined were edited and manually aligned with BIO-EDIT ver. 7.2.5 (Hall, 1999) and SE-AL ver. 2.0a111 (http://tree.bio.ed.ac.uk/software/seal). All unique sequences were submitted to GenBank (Accession numbers ITS: KY446904–KY446922, trnL-F: KY446923–KY446942). Divergence time estimation and genealogical relationships We estimated divergence time under a Bayesian approach using the combined ITS and plastid trnL-F sequence data set. The ingroup comprised all sequences of P. mayanus; sequences of other representatives of Loranthaceae including other Psittacanthus spp. downloaded from GenBank from Wilson & Calvin (2006), Amico, Vidal-Russell & Nickrent (2007), Vidal-Russell & Nickrent (2007, 2008a, b), Díaz Infante et al. (2016), Ornelas et al. (2016) and Pérez-Crespo et al. (2017) were used as outgroups. Phylogenetic reconstruction and divergence time estimation was performed with BEAST ver. 1.6.1 (Drummond & Rambaut, 2007) using the uncorrelated lognormal relaxed molecular clock and the closest nucleotide substitution model under the Bayesian information criterion (BIC), GTR+G+I for the combined data set, suggested by jMODELTEST ver. 0.1.1 (Posada, 2008). The tree prior model was set using a coalescent approach assuming constant population size. In this analysis, we constrained Nuytsia floribunda (Labill.) G.Don as sister to the aerial parasites based on Vidal-Russell & Nickrent (2008b) to calibrate the root node. The divergence time between Nuytsia and the aerial parasite clade was used as secondary calibration, approximating a median age of 28 Myr (normal distribution, mean 28, SD 3.7, range 34–21 Myr). The geometric mean of 5.798 × 10–9 substitutions per neutral site per year (s/s/y) was used to calibrate the tree based on the mean mutation rates of 4.130 × 10–9 s/s/y for ITS of herbaceous annual/perennial plants (Kay, Whittall & Hodges, 2006) and 8.240 × 10–9 s/s/y for trnL-F estimated for annual or perennial herbs (Richardson et al., 2001). For divergence time estimation, BEAST was run two times for 30 million generations, sampling every 3000 steps. We combined log and trees files from each independent run using LOGCOMBINER ver. 1.8.0 (Drummond & Rambaut, 2007), then viewed the combined log file in TRACER ver. 1.6 (http://tree.bio.ed.ac.uk/software/tracer/) to ensure that ESSs 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 10% trees were discarded as burn-in, and the remaining trees were annotated and summarized as a maximum clade credibility tree with mean divergence times and 95% highest posterior density (HPD) intervals of age estimates using TREEANNOTATOR ver. 1.8.0 (Drummond & Rambaut, 2007) and visualized in FIGTREE ver. 1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/). To infer genealogical relationships among haplotypes, statistical parsimony networks for ITS and trnL-F data sets were constructed using TCS ver. 1.2.1 (Clement, Posada & Crandall, 2000), with polymorphic sites coded as ambiguous (ITS) or gaps treated as single evolutionary events (trnL-F) and a connection limit set to 95%. Loops were resolved following the criteria given by Pfenninger & Posada (2002). Genetic diversity and population structure Molecular diversity indices including haplotype diversity (h) and nucleotide diversity (π), pairwise comparisons of FST values between genetic groups with 1600 permutations, and analyses of molecular variance (AMOVAs; Excoffier, Smouse & Quattro, 1992) were performed in ARLEQUIN ver. 3.01 (Excoffier, Laval & Schneider, 2005). Populations with fewer than three samples were combined with closest populations for this analysis (Table 1). To explore whether populations are structured along the Yucatán Peninsula and Chiapas, three AMOVAs were performed based on pairwise differences using ARLEQUIN with populations treated as (1) a single group to determine the amount of variation partitioned among and within populations, and grouped by geography into two (2) Yucatán Peninsula and Chiapas or three lineages (3) corresponding to the southern (PET = Petén province) or northern (YUC = Yucatán province) distribution of the Yucatán Peninsula and Chiapas (CHI) or (4) grouped by precipitation into three lineages corresponding to dry (dryYUC or dryCHI) and humid (humYUC) locations (Table 1). The AMOVAs were performed using the Kimura 2P model for ITS and the Jukes & Cantor model for trnL-F and 16 000 permutations to determine the significance of each AMOVA. The program PERMUT ver. 1.0 (Pons & Petit, 1996) was used to calculate for each marker within- population-diversity (hS), total diversity (hT), geographical total haplotype diversity (vT), geographical average haplotype diversity (vS), level of population differentiation at the species level (GST) and an estimate of population subdivisions for phylogenetically ordered alleles (NST). We further tested for phylogeographic structure at the species range between NST and GST parameters using 10 000 permutations and the U-statistic. An NST value significantly higher than the GST value means that two alleles sampled from the same region are phylogenetically more closely related than two alleles sampled from different regions, providing evidence for the presence of a significant phylogeographic signal in the data (Pons & Petit, 1996). Demographic history The demographic patterns of P. mayanus populations were assessed using different approaches. Fu’s Fs (Fu, 1997) and Tajima’s D (Tajima, 1989) tests were computed for each marker with 25 000 permutations to judge significance to test whether populations departed from a mutation-drift equilibrium model. A significantly negative value is often taken as an indicator of population expansion. In addition, mismatch distributions (Harpending, 1994) were simulated under the sudden expansion model of Schneider & Excoffier (1999) with 16 000 bootstrap replicates. The validity of the sudden expansion assumption was determined using the sum of squares differences (SSD) and Harpending’s raggedness index (Hri; Harpending, 1994), both of which are higher in stable, non-expanding populations (Rogers & Harpending, 1992). The SSD between the observed and the expected mismatch was calculated, and a significant SSD value means the observed value does not fit sudden expansion model. All tests were performed using ARLEQUIN and DnaSP (Librado & Rozas, 2009). The demographic changes in effective population size (Ne) through time of each P. mayanus group (Petén and Yucatán) were also inferred with Bayesian skyline plots (Drummond et al., 2005) using BEAST. We chose a HKY substitution model with empirical base frequencies, a strict clock model and a piecewise-linear coalescent Bayesian skyline tree prior with five starting groups. Two independent runs of 20 million generations each were run, with trees and parameters sampled every 1000-iterations and a burn-in of 10%. Results of each run were visualized using TRACER to ensure that stationarity and convergence had been reached, and that the ESS were higher than 200. We used the software IMa (Hey & Nielsen, 2007) on genetic groups in P. mayanus (Yucatán Peninsula = 63 samples, Chiapas = 24 samples) to estimate the effective population size of the ancestral (θa) and the two descendant populations (θ1 and θ2), effective migration rates in both directions (m1-to-2 and m2-to-1) and time since divergence (t) at which the ancestral population gave rise to the descendant populations. The isolation-with-migration (IM) model (Hey & Nielsen, 2004) is appropriate to estimate parameters for two descendant populations that have diverged recently from the ancestral population and that may be sharing haplotypes as a result of gene interchange. We began with multiple short runs of 10 000 steps (following 100 000 iterations as burn-in) to assess mixing and to fine-tune the parameter space. Once the six parameters were optimized (several replicates converged on the same answer), a long run of MCMC simulations was performed for one million generations and 30 million steps, under the HKY model of sequence evolution. Two additional independent runs were performed with different seed numbers to guarantee convergence of samples (Hey & Nielsen, 2004). We considered that the analyses had converged upon a stationary distribution if the independent runs had similar posterior distributions and the ESS for each parameter was at least 50. We report the mean parameter estimates of three runs and the 95% highest posterior densities (HPDs) intervals of each parameter. To convert raw parameter estimates into demographic quantities, we used the geometric mean of 4.58 × 10–5 of the median mutation rates of both loci reported for P. schiedeanus (Ornelas et al., 2016) and P. calyculatus (Pérez-Crespo et al., 2017). These mutation rates were converted to per-locus mutation rate (μ) by multiplying by the fragment length in base pairs for conversion, as required by IMa (Hey & Nielsen, 2007), to compute the divergence time (t) by using T = t/μ. We used an 11-year generation time (G) based on the observation that the age at maturity (seed production) begins c. 2 years after seed germination and an assumed annual survivorship of 0.9 based on seedling survivorship in the sister species P. schiedeanus (López de Buen & Ornelas, 2002; Ornelas et al., 2016) to convert the effective populations size estimates. The approximate average generation time (T) is calculated according to T = a + [s ⁄ (1–s)] (Lande, Engen & Sæther, 2003), where a is the time to maturity and s is the adult annual survival rate. Based on this, the estimate for T was 11 years. To convert the IMa time since divergence parameter to years, t, we divided the time parameter (B) by the mutation rate per year (U) converted to per locus rate by multiplying by the fragment length in base pairs, and calculated for the geometric mean rate described before. To calculate effective population size (N), we used N = θ/(4Gμ), where the generation time (G) is 11 years. To estimate the population migration rate, we used 2Nm = θ * m/2. Population expansion and approximate Bayesian computation (ABC) To assess glacial and post-glacial changes in population size within the Yucatán Peninsula that could explain population sizes and genetic diversity, we applied approximate Bayesian computation (ABC; Beaumont, Zhang & Balding, 2002) implemented in DIYABC (Cornuet et al., 2014). Two different models for population size changes were assumed, with population growth or population decline. We assumed either 6000 or 21 000 yr bp for the beginning of the population size change, which correspond to the MH or the LGM, respectively. For this analysis, we grouped all populations into three groups of samples (dryCHI, humYUC, dryYUC): scenario 1, population expansion from 21 000–10 000 yr bp after the split between the dryCHIS and the YUC groups; scenario 2, population expansion from 1000–8000 yr bp after the split between the dryCHIS and the YUC groups. As little is known about the effective population size of P. mayanus in the Yucatán Peninsula, we defined a Ne prior with a uniform distribution for the YUC population based on previous estimates of closely related Psittacanthus spp. (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Furthermore, we assumed 11 years as the generation time and mutation rates based on estimates for other Psittacanthus spp. (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Comparison of scenarios was implemented in the DIYABC, with two million coalescent-based simulated datasets under each evolutionary scenario. The posterior probability of scenarios was assessed using a logistic regression on the 1% of simulated datasets closest to the observed data (Fontaine et al., 2013). For the best-supported scenario, we performed a model checking procedure, as implemented in DIYABC, and assessed the confidence in scenario choice by simulating 500 pseudo-observed datasets (PODs) under each scenario to estimate Type I and Type II error rates (Robert et al., 2011). Finally, point estimates for demographic and temporal parameters were obtained by local linear regression on the 1% of simulations closest to the observed dataset for the best-supported scenario (Cornuet et al., 2008, 2014). Palaeodistribution modelling and environmental variation We used ecological niche modelling (ENM; Elith et al., 2011) to predict the distribution of P. mayanus at the MH (c. 6000 years ago), LGM (c. 20 000 years ago) and at the LIG (c. 120 000–140 000 years ago) and to examine whether populations of P. mayanus (1) persisted in the southern portion of the Yucatán Peninsula (in the wetter Petén province) in the semideciduous tropical rain forest during the cold glacial period and expanded northwards to the seasonally dry tropical deciduous forest after the LGM or (2) contracted into southern refugia during the LGM and expanded northwards only after this period. Mistletoe occurrence data were assembled mainly from herbarium records from the Global Biodiversity Information Facility (GBIF; https://www.gbif.org/species/4003770), and supplemented with records from Kuijt (2009) and our geo-referenced records from field collection. After careful verification of data for every location, excluding duplicate occurrence records or in close proximity to each other (c. 1 km2) to reduce the effects of spatial autocorrelation, we restricted the data set to 70 unique presence records for the analysis. Distributional records were input and analysed to infer an ENM with the maximum entropy algorithm in MAXENT ver. 3.2.2 (Phillips, Anderson & Schapire, 2006) and using ARCVIEW ver. 3.2 (ESRI, Redlands, CA, USA) to extract the data. Present-day and Last Interglacial (LIG) temperature and precipitation data (BIO 1–19 variables) were drawn as climate layers from the WorldClim database at a spatial resolution of c. 1 km2 in each raster (Hijmans et al., 2005; http://www.worldclim.org). We also obtained data for past climate layers for two MH (2.5 arc-minutes, pixel size c. 21.62 km2) and LGM (2.5 arc-minutes) past conditions based on CCSM4 and MIROC-ESM global models (Hijmans et al., 2005; Braconnot et al., 2007). We first constructed a model to the present using all climatic variables with 70% of the locality records as training data and the other 30% as testing data in addition to response curves, jackknife tests, logistic output format and random seed parameters. We ran 1000 iterations with no extrapolation in order to avoid artificial extrapolations from extreme values of the ecological variables; as such parameters are biased towards the environmental envelope of background points and occurrence data (Phillips et al., 2006). All other parameters in MAXENT were maintained at default settings. The logistic format was used to obtain the values of habitat suitability and the variable importance was determined comparing percentage contribution values and jackknife plots. Then, correlation coefficients between present variables were calculated using PAST ver. 2.12 (Hammer, 2011) with the objective to identify and remove the variables that were highly correlated (correlation values ≥ 0.8) and less explanatory, having as a rule of decision, the relative contributions of each of them in the MAXENT model. After removing the highly correlated variables, six variables were used in the final analysis for the mistletoe (BIO1 = annual mean temperature, BIO3 = isothermality, BIO4 = temperature seasonality, BIO6 = minimum temperature of coldest month, BIO7 = temperature annual range, BIO18 = precipitation of warmest quarter). The models were calibrated according the accessibility area of P. mayanus; this was estimated using a geographical clipping based on biogeographic provinces (Morrone, 2005, 2014) and elevational range limits (Kuijt, 2009), which represent the species historical barriers to dispersal (Kuijt, 2009). The resulting species distribution under current climate conditions was then projected onto the LIG and onto past climate reconstructions for the MH and LGM under two general atmospheric circulation models (CCSM4 and MIROC-ESM). Final models were constructed with ten cross-validation replicates without clamping and extrapolation and considering the average output grids as the final predictive models. The area under the curve (AUC) values of receiver operating characteristic (ROC) curves was used to evaluate the statistical performance of ENMs, a single measure of model performance independently of thresholds (Phillips et al., 2006), where 1 is the maximum prediction and 0.5 a random prediction. In addition, we also used the true skill statistic (TSS), which is a threshold-dependent measure of model performance, to evaluate the accuracy of predictive maps generated by presence-only data (Allouche, Tsoar & Kadmon, 2006; Liu et al., 2005), where TSS values ranging between 0.4–0.8 are considered useful (Landis & Koch, 1977; Fielding & Bell, 1997). For each replicate, the TSS was calculated using the ten-percentile training presence logistic threshold, and then TSS values were averaged among replicates. Measurements of temperature and precipitation (BIO 1–19 variables) were extracted for each of the sampling locations and for those with mistletoe occurrence data from WorldClim (Hijmans et al., 2005) as explained before. We performed principal components analysis (PCA) using SPSS Statistics for Mac ver. 17 (SPSS, Inc., Chicago, IL, USA) to reduce BIO variables to a smaller system of uncorrelated, independent variables and then used the PC loadings to examine habitat and ecological variation potentially related to population divergence. Differences in PC scores between genetic groups were tested with one-way ANOVAs in SPSS. RESULTS Divergence time estimation and genealogical relationships The aligned sequences of ITS and trnL-F were 640 and 381 base pairs (bp) long, respectively. Out of 1021 bp of the concatenated sequences (n = 197), 590 sites were variable and out of these 464 were parsimony informative. In P. mayanus samples, we detected polymorphic sites in several samples of the ITS region, with 17 of the 564 sites (3%). These were coded as ambiguous. The BEAST analyses placed the origin of the P. mayanus clade in the mid Pleistocene, 1.747 Mya (95% HPD 3.172–0.519 Mya; PP = 0.99; Fig. 1). Although a split estimated at 0.962 Mya (95% HPD 1.738–0.263 Mya) between samples of P. mayanus from Chiapas and the Petén province and samples from the Yucatán province is recovered with strong support (PP = 0.99), the relationships in each geographic region were poorly resolved, with samples from Chiapas intermingled among the Yucatán province clade (Fig. 1). Figure 1. View largeDownload slide (A) Psittacanthus mayanus. Photograph by Andrés Ortiz Rodríguez. (B) Chronogram based on a Bayesian approach using a coalescent prior under an uncorrelated lognormal relaxed clock model and assuming constant population size of P. mayanus DNA sequences in BEAST. Blue bars indicate 95% highest posterior density (HPD) intervals for nodes of particular interest. These nodes all have posterior probabilities >0.9. Figure 1. View largeDownload slide (A) Psittacanthus mayanus. Photograph by Andrés Ortiz Rodríguez. (B) Chronogram based on a Bayesian approach using a coalescent prior under an uncorrelated lognormal relaxed clock model and assuming constant population size of P. mayanus DNA sequences in BEAST. Blue bars indicate 95% highest posterior density (HPD) intervals for nodes of particular interest. These nodes all have posterior probabilities >0.9. The aligned ITS dataset for 82 samples of P. mayanus (564 bp) yielded 14 ribotypes (Table 2). In the ribotype network (Fig. 2), ribotype R1 was found in 53 of the 82 samples for ten of the 15 sampled populations and forms the core of the network. This widespread and most frequent ribotype (64.6% of the individuals and 66.6% of the sampled populations), inferred as ancestral based on its frequency and position in the network (Fig. 2), was separated by one or two mutations of most ribotypes in the network. The second most frequent ribotypes, R11 and R8 (8 and 6 samples, respectively) were mainly found in populations of Campeche and Chiapas and shared by one sample from southern Yucatán (Table 2). Some ribotypes had certain geographical structure and were exclusively distributed in localities of the Yucatán Peninsula (R2–R7) or Chiapas (R9–10, R12–R14) (Fig. 2). Table 2. Number of genetically analysed samples (n) for each molecular marker (ITS and trnL-F) and number of distinct ribotypes (R) and haplotypes (H) found in sampled Psittacanthus mayanus individuals, and the number of individuals per ribotype or haplotype in parentheses. Population number  Location  n  ITS ribotype  n  trnL-F haplotype  1  Yucatán, Dzibilchaltún  3  R1(2), R3(1)  3  H1(3)  2  Yucatán, Hunucmá, Sisal  14  R1(14)  13  H1(12), H3(1)  3  Yucatán, Celestún  9  R1(7), R2(2)  9  H1(9)  4  Yucatán, Mérida, Cuxtal  13  R1(11), R4(1), R6(1)  11  H1(10), H2(1)  5  Yucatán, Maxcanú, Límite Estatal  6  R1(5), R8(1)  7  H1(7)  6  Campeche, Champotón, Seybaplaya  9  R1(6), R3(1), R5(1), R7(1)  10  H1(10)  7  Campeche, Escárcega, Sabancuy  10  R1(3), R8(4), R11(3)  10  H1(10)  8  Chiapas, Bochil, Km 205 Pto. Caté  1  R11(1)  3  H1(3)  9  Chiapas, Bochil, Sn Cristobalito  1  R12(1)  3  H6(2), H7(1)  10  Chiapas, Sima de las Cotorras  6  R1(3), R10(1), R11(1), R14(1)  6  H1(6)  11  Chiapas, Ocozocuautla, El Aguacero  1  R1(1)  1  H4(1)  12  Chiapas, Ocozocuautla  1  R14(1)  1  H1(1)  13  Chiapas, Arriaga, Km 89–92  2  R9(2)  2  H1(1), H5(1)  14  Chiapas, Km 9 Arriaga-Tuxtla    -  1  H1(1)  15  Chiapas, Arriaga, La Aurora  1  R1(1)  1  H1(1)  16  Chiapas, Arriaga, Km 28.3  5  R8(1), R11(3), R13(1)  6  H1(6)  Population number  Location  n  ITS ribotype  n  trnL-F haplotype  1  Yucatán, Dzibilchaltún  3  R1(2), R3(1)  3  H1(3)  2  Yucatán, Hunucmá, Sisal  14  R1(14)  13  H1(12), H3(1)  3  Yucatán, Celestún  9  R1(7), R2(2)  9  H1(9)  4  Yucatán, Mérida, Cuxtal  13  R1(11), R4(1), R6(1)  11  H1(10), H2(1)  5  Yucatán, Maxcanú, Límite Estatal  6  R1(5), R8(1)  7  H1(7)  6  Campeche, Champotón, Seybaplaya  9  R1(6), R3(1), R5(1), R7(1)  10  H1(10)  7  Campeche, Escárcega, Sabancuy  10  R1(3), R8(4), R11(3)  10  H1(10)  8  Chiapas, Bochil, Km 205 Pto. Caté  1  R11(1)  3  H1(3)  9  Chiapas, Bochil, Sn Cristobalito  1  R12(1)  3  H6(2), H7(1)  10  Chiapas, Sima de las Cotorras  6  R1(3), R10(1), R11(1), R14(1)  6  H1(6)  11  Chiapas, Ocozocuautla, El Aguacero  1  R1(1)  1  H4(1)  12  Chiapas, Ocozocuautla  1  R14(1)  1  H1(1)  13  Chiapas, Arriaga, Km 89–92  2  R9(2)  2  H1(1), H5(1)  14  Chiapas, Km 9 Arriaga-Tuxtla    -  1  H1(1)  15  Chiapas, Arriaga, La Aurora  1  R1(1)  1  H1(1)  16  Chiapas, Arriaga, Km 28.3  5  R8(1), R11(3), R13(1)  6  H1(6)  Codes are from networks in Figure 3 View Large Figure 2. View largeDownload slide Geographical distribution and statistical parsimony networks of Psittacanthus mayanus overlaid on a map of the Yucatán Peninsula and adjacent states of southern Mexico. (A) Fourteen ITS ribotypes. (B) Seven trnL-F haplotypes. Ribotype or haplotype designations in the network correspond to those in Table 2, the size of the circles is proportional to the frequency of each ribotype or haplotype, and small black-filled circles are non-sampled ribotypes or haplotypes. Each ribotype or haplotype is coded with a different colour. The numbers in the ribotypes or haplotypes indicate the number of individuals that share the ribotype or haplotype. Pie charts represent ribotypes or haplotypes found in each sampling locality. The size of sections of the pie charts corresponds to the proportion of individuals with a given ribotype or haplotype. Figure 2. View largeDownload slide Geographical distribution and statistical parsimony networks of Psittacanthus mayanus overlaid on a map of the Yucatán Peninsula and adjacent states of southern Mexico. (A) Fourteen ITS ribotypes. (B) Seven trnL-F haplotypes. Ribotype or haplotype designations in the network correspond to those in Table 2, the size of the circles is proportional to the frequency of each ribotype or haplotype, and small black-filled circles are non-sampled ribotypes or haplotypes. Each ribotype or haplotype is coded with a different colour. The numbers in the ribotypes or haplotypes indicate the number of individuals that share the ribotype or haplotype. Pie charts represent ribotypes or haplotypes found in each sampling locality. The size of sections of the pie charts corresponds to the proportion of individuals with a given ribotype or haplotype. The aligned trnL-F dataset for 87 samples of P. mayanus (275 bp) yielded seven haplotypes (Fig. 2; Table 2). Haplotype H1, the most frequent (80 samples, 91.9% of the individuals), was retrieved from most sampled populations located both along the Yucatán Peninsula and in Chiapas (14 of the 16 sampled populations). Haplotypes H2 and H3 were singletons one step separated from H1 and located in northern Yucatán, and haplotypes H4–H7, only located in Chiapas, were also singletons or shared by two individuals (Fig. 2). Genetic diversity and population structure Differentiation among populations based on ITS variation indicated genetic subdivision (GST = 0.210, SE = 0.0604). Genetic diversity across all populations (hT = 0.667, SE = 0.1004; vT = 0.671, SE = 0.1046) was higher than within-population genetic diversity (hS = 0.527, SE = 0.0838; vS = 0.486, SE = 0.0884). PERMUT analysis showed that NST (0.276, SE = 0.0610) was not significantly higher (p > 0.05) than GST, indicating no phylogeographical structure. Similarly, phylogeographical structure was not observed based on trnL-F variation (GST = 0.461, SE = not calculated, NST = 0.609, SE = not calculated; NST = GST, p > 0.05), with genetic diversity across all populations (hT = 0.245, SE = 0.1457; vT = 0.248, SE = 0.1781) higher than within-populations (hS = 0.132, SE = 0.0612; vS = 0.097, SE = 0.0594). Pairwise comparisons of FST values were low but significant for both ITS and trnL-F when two groups of populations (Yucatán and Chiapas) or the most distant groups (YUC and CHI or dryYUC and dryCHI) in the three-group comparisons were compared (Table 3). Table 3. Pairwise comparisons of FST values of ITS (above the diagonal) and trnL-F (below the diagonal) among groups of Psittacanthus mayanus populations. Grouping  A. Two groups by geography  Yucatan  Chiapas    Yucatán  –  0.13088**    Chiapas  0.14846**  –    B. Three groups by geography  YUC  PET  CHI  YUC  _  0.07831*  0.21593***  PET  −0.02018  –  0.01033  CHI  0.11320*  0.11439  –  C. Three groups by precipitation  dryYUC  humYUC  dryCHI  dryYUC  –  0.21518**  0.23890***  humYUC  0.15614*  –  −0.04345  dryCHI  0.00791  −0.00197  –  Grouping  A. Two groups by geography  Yucatan  Chiapas    Yucatán  –  0.13088**    Chiapas  0.14846**  –    B. Three groups by geography  YUC  PET  CHI  YUC  _  0.07831*  0.21593***  PET  −0.02018  –  0.01033  CHI  0.11320*  0.11439  –  C. Three groups by precipitation  dryYUC  humYUC  dryCHI  dryYUC  –  0.21518**  0.23890***  humYUC  0.15614*  –  −0.04345  dryCHI  0.00791  −0.00197  –  Significant values at P < 0.001 are highlighted in bold. Group abbreviations are as follows: Yucatán = Yucatan Peninsula; Chiapas = Estado de Chiapas; YUC = Yucatán Province; PET = Petén Province; dryYUC = dry Yucatán Peninsula, humYUC = humid Yucatán Peninsula; dryCHI = dry Central Depression of Chiapas. * P < 0.01, ** P < 0.001, *** P < 0.0001. View Large The AMOVAs showed that 16.5% of genetic variation for ITS and 18.4% for trnL-F was explained by differences within populations and 83.5% and 81.6%, respectively, by differences between populations when all locations were treated as a single group (Table 4). The AMOVAs for ITS and trnL-F revealed population structure, with high and significant FCT values obtained when populations are grouped into two lineages corresponding to those along the Yucatán Peninsula and in Chiapas (Table 4). When sampling sites are grouped into tree lineages separated by geography or precipitation, a significant proportion of the ITS variation was also attributed to differences between groups (Table 4). Nucleotide diversity (π) was low for both datasets (Table 5). Table 4. AMOVA results for plausible Psittacanthus mayanus populations with no groups defined a priori (A), and grouping based on geography (two geographical groups: Yucatán Peninsula and Chiapas; three geographic groups: YUC = Yucatán province, PET = Petén province, CHI = Chiapas) (B and C) or based on the precipitation gradient (dryYUC, humYUC, dryCHI) (D). Source of variation  ITS  trnL-F  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  A. Single group   Among populations  10  10.417  0.08467  16.54    11  1.338  0.01054  18.38     Within populations  71  30.342  0.42735  83.46  FST = 0.16**  75  3.513  0.04684  81.62  FST = 0.18*   Total  81  40.759  0.51202      86  4.851  0.05738      B. Two geographic groups   Among groups  1  4.290  0.12918  21.85  FCT = 0.22**  1  0.393  0.00894  14.32  FCT = 0.14**   Among pop. within groups  9  6.127  0.03464  5.86  FSC = 0.07 ns  10  0.945  0.00667  10.68  FSC = 0.12 ns   Within populations  71  30.342  0.42735  72.29  FST = 0.28**  75  3.513  0.04684  75.00  FST = 0.25**  Total  81  40.759  0.59118      86  4.851  0.06245      C. Three geographic groups   Among groups  2  5.468  0.08492  15.74  FCT = 0.16*  2  0.408  0.00324  5.55  FCT = 0.05 ns   Among pop. within groups  8  4.949  0.02738  5.07  FSC = 0.06 ns  9  0.931  0.00830  14.21  FSC = 0.15 ns   Within populations  71  30.342  0.42735  79.19  FST = 0.21**  75  3.513  0.04684  80.23  FST = 0.20**   Total  81  40.759  0.53965      86  4.851  0.05838      D. Three precipitation groups   Among groups  2  6.967  0.14687  25.55  FCT = 0.25***  2  0.311  0.00158  2.73  FCT = 0.03 ns   Among pop. within groups  8  3.450  0.00055  0.10  FSC = 0.00 ns  9  1.027  0.00959  16.52  FSC = 0.17*   Within populations  71  30.342  0.42735  74.35  FST = 0.26**  75  3.513  0.04684  80.74  FST = 0.19**   Total  81  40.759  0.57477      86  4.851  0.05801      Source of variation  ITS  trnL-F  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  A. Single group   Among populations  10  10.417  0.08467  16.54    11  1.338  0.01054  18.38     Within populations  71  30.342  0.42735  83.46  FST = 0.16**  75  3.513  0.04684  81.62  FST = 0.18*   Total  81  40.759  0.51202      86  4.851  0.05738      B. Two geographic groups   Among groups  1  4.290  0.12918  21.85  FCT = 0.22**  1  0.393  0.00894  14.32  FCT = 0.14**   Among pop. within groups  9  6.127  0.03464  5.86  FSC = 0.07 ns  10  0.945  0.00667  10.68  FSC = 0.12 ns   Within populations  71  30.342  0.42735  72.29  FST = 0.28**  75  3.513  0.04684  75.00  FST = 0.25**  Total  81  40.759  0.59118      86  4.851  0.06245      C. Three geographic groups   Among groups  2  5.468  0.08492  15.74  FCT = 0.16*  2  0.408  0.00324  5.55  FCT = 0.05 ns   Among pop. within groups  8  4.949  0.02738  5.07  FSC = 0.06 ns  9  0.931  0.00830  14.21  FSC = 0.15 ns   Within populations  71  30.342  0.42735  79.19  FST = 0.21**  75  3.513  0.04684  80.23  FST = 0.20**   Total  81  40.759  0.53965      86  4.851  0.05838      D. Three precipitation groups   Among groups  2  6.967  0.14687  25.55  FCT = 0.25***  2  0.311  0.00158  2.73  FCT = 0.03 ns   Among pop. within groups  8  3.450  0.00055  0.10  FSC = 0.00 ns  9  1.027  0.00959  16.52  FSC = 0.17*   Within populations  71  30.342  0.42735  74.35  FST = 0.26**  75  3.513  0.04684  80.74  FST = 0.19**   Total  81  40.759  0.57477      86  4.851  0.05801      ns, not significant (P > 0.05), * P < 0.05, ** P < 0.001, *** P < 0.0001 View Large Table 5. Summary statistics of neutrality tests and demographic analysis of Psittacanthus mayanus samples by genetic groups to infer demographic range expansion. Parameter  Two groups  Three groups by geography  Three groups by precipitation  Yucatán  Chiapas  YUC  PET  CHI  dryYUC  humYUC  dryCHI  ITS  N  64  18  45  19  18  54  12  16  NH  8  6  6  5  6  8  3  5  h, gene diversity  0.3616  0.6797  0.2505  0.5789  0.6797  0.3068  0.5909  0.6667  π, nucleotide diversity  0.001871  0.001818  0.001164  0.003088  0.001818  0.001704  0.002055  0.001788  Tajima’s D  –1.467*  –0.420  –2.044**  –0.767  –0.420  –1.769*  0.472  –0.576  Fu’s Fs  –2.473  –2.449*  –2.497*  0.183  –2.449*  –3.086*  0.951  –1.469  SDD  0.0215  0.0051  0.0068  0.4324***  0.0051  0.0123  0.0880  0.0066  Hri  0.3059  0.0745  0.4096  0.2321  0.0744  0.3366  0.3289  0.0725  trnL-F  N  63  24  43  20  24  53  16  18  NH  3  5  3  1  5  3  3  3  h, gene diversity  0.0620  0.3913  0.0919  0.0000  0.3768  0.0747  0.3417  0.2157  π, nucleotide diversity  0.000227  0.002314  0.000338  0.000000  0.002227  0.000274  0.001692  0.000808  Tajima’s D  –1.437*  –0.317  –1.479*  n.a.  –0.317  –1.458*  0.000  –1.165  Fu’s Fs  –3.356**  –2.187*  –2.827**  n.a.  –2.187*  –3.113**  –0.571  –1.743*  SDD  0.0002  0.2109***  0.0001  n.a.  0.2109***  0.0000  0.1937  0.0024  Hri  0.7687  0.2669  0.6759  n.a.  0.2669  0.7299  0.2079  0.3719  Parameter  Two groups  Three groups by geography  Three groups by precipitation  Yucatán  Chiapas  YUC  PET  CHI  dryYUC  humYUC  dryCHI  ITS  N  64  18  45  19  18  54  12  16  NH  8  6  6  5  6  8  3  5  h, gene diversity  0.3616  0.6797  0.2505  0.5789  0.6797  0.3068  0.5909  0.6667  π, nucleotide diversity  0.001871  0.001818  0.001164  0.003088  0.001818  0.001704  0.002055  0.001788  Tajima’s D  –1.467*  –0.420  –2.044**  –0.767  –0.420  –1.769*  0.472  –0.576  Fu’s Fs  –2.473  –2.449*  –2.497*  0.183  –2.449*  –3.086*  0.951  –1.469  SDD  0.0215  0.0051  0.0068  0.4324***  0.0051  0.0123  0.0880  0.0066  Hri  0.3059  0.0745  0.4096  0.2321  0.0744  0.3366  0.3289  0.0725  trnL-F  N  63  24  43  20  24  53  16  18  NH  3  5  3  1  5  3  3  3  h, gene diversity  0.0620  0.3913  0.0919  0.0000  0.3768  0.0747  0.3417  0.2157  π, nucleotide diversity  0.000227  0.002314  0.000338  0.000000  0.002227  0.000274  0.001692  0.000808  Tajima’s D  –1.437*  –0.317  –1.479*  n.a.  –0.317  –1.458*  0.000  –1.165  Fu’s Fs  –3.356**  –2.187*  –2.827**  n.a.  –2.187*  –3.113**  –0.571  –1.743*  SDD  0.0002  0.2109***  0.0001  n.a.  0.2109***  0.0000  0.1937  0.0024  Hri  0.7687  0.2669  0.6759  n.a.  0.2669  0.7299  0.2079  0.3719  N = number of individuals, NH = number of ribotypes or haplotypes, SDD = differences in the sum of squares or mismatch distribution, Hri = Harpending’s raggedness index. DT and FS positive values are indicative of mutation-drift-equilibrium, which is typical of stable populations, and negative values that result from an excess of rare haplotypes indicate that populations have undergone recent expansions, often preceded by a bottleneck. Significantly negative values (at the 0.05 level) reveal in both tests historic demographic expansion events. Significant (P ≤ 0.05) SSD and Hri values indicate deviations from the sudden expansion model. n.a. = not available; *P < 0.05; **P < 0.01; ***P < 0.001. Group abbreviations are as follows: PET = Petén province, YUC = Yucatán province, CHI = Chiapas, dryYUC = dry Yucatán, humYUC = humid Yucatán, dryCHI = dry Central Depression of Chiapas View Large Demographic history When populations were analysed as two (Yucatán Peninsula and Chiapas) or three groups based on geography (YUC, PET, CHI) or precipitation (dryYUC, humYUC, dryCHI), negative and significant values for neutrality tests (Fu’s Fs, Tajima’s D) were observed (Table 5), indicating a departure from neutrality (probably an excess of rare alleles or new mutations from population expansion) for northern peninsular P. mayanus (Yucatán, YUC, dryYUC) populations in most cases. The model of sudden demographic expansion for most cases of Yucatán groups was not rejected by the generalized least square procedure (SSD) or by the raggedness index of the distribution (Hri), except for three cases, the Chiapas and CHI groups in the trnL-F dataset and PET in the ITS dataset (Table 5). The Bayesian skyline plots suggested that Chiapas and Yucatán Peninsula populations underwent a constant effective population size, with a marginal increase in population size occurring after the MH (Fig. 3). When groups of populations were pooled, the increase in population size occurring after the MH was more pronounced (Fig. 3). Figure 3. View largeDownload slide Bayesian skyline plots showing historical demographic trends of Psittacanthus mayanus for YUC (A), CHI (B), and for all the P. mayanus (C) populations using the concatenated data set. The y-axis is the product between effective population size and the generation time and the x-axis is time in thousands of years. Solid lines represent mean estimates and shaded areas represent 95% confidence intervals. YUC = Yucatán Peninsula, CHI = Chiapas. Figure 3. View largeDownload slide Bayesian skyline plots showing historical demographic trends of Psittacanthus mayanus for YUC (A), CHI (B), and for all the P. mayanus (C) populations using the concatenated data set. The y-axis is the product between effective population size and the generation time and the x-axis is time in thousands of years. Solid lines represent mean estimates and shaded areas represent 95% confidence intervals. YUC = Yucatán Peninsula, CHI = Chiapas. According to IMa results (summarized in Table 6), the peninsular (YUC) population size (Ne1) is estimated to be larger than population size of ancestral (Nea) and Chiapas (CHI, Ne2) descendant population. Divergence time between YUC and CHI populations was estimated at c. 66 kya. A peak at the left-hand side of the posterior distribution of t is observed, with a long right-hand tail over which the probability distribution remains flat. This means that this parameter value needs to be interpreted with caution because parameter values have essentially zero probability for most of the posterior distribution and the only values that are well supported are those less than 33 kya. When testing for migration during the time since population splitting, migration among populations occurred in one direction, moving from CHI to YUC (Table 6). Table 6. Results of isolation-with-migration model (IMa) for the split between peninsular (YUC) and Chiapas (CHI) populations of Psittacanthus mayanus.   Model parameter estimates  θ1  θ2  θa  t  m2-to-1  m1-to-2  YUC vs. CHI  Mean  1.7503  0.1994  0.9657  3.029  67.43  0.483  HPD95Lo  0.1553  0.0037  0.0479  0.199  17.31  0.024  HPD95Hi  3.4506  0.8537  1.9043  5.851  98.75  0.945    Demographic parameter estimates  Ne1  Ne2  Nea  t  Nm2-to-1  Nm1-to-2  YUC vs. CHI  Mean  868.487  98.939  479.165  66,137  59.012  0.048  HPD95Lo  77.074  1.819  23.751  4,344  1.344  0.000  HPD95Hi  1712.132  423.592  944.883  127,739  170.373  0.403    Model parameter estimates  θ1  θ2  θa  t  m2-to-1  m1-to-2  YUC vs. CHI  Mean  1.7503  0.1994  0.9657  3.029  67.43  0.483  HPD95Lo  0.1553  0.0037  0.0479  0.199  17.31  0.024  HPD95Hi  3.4506  0.8537  1.9043  5.851  98.75  0.945    Demographic parameter estimates  Ne1  Ne2  Nea  t  Nm2-to-1  Nm1-to-2  YUC vs. CHI  Mean  868.487  98.939  479.165  66,137  59.012  0.048  HPD95Lo  77.074  1.819  23.751  4,344  1.344  0.000  HPD95Hi  1712.132  423.592  944.883  127,739  170.373  0.403  Model parameters indicate estimates without use of molecular rate of evolution for six parameters (IMa output values). Values are averages of two runs of mean parameter estimates and the 95% highest posterior densities (HPD) intervals of each parameter: effective population size of the ancestral (θa) and descendant populations (θ1, θ2), effective migration rates in both directions (m1-to-2 and m2-to-1), and time since divergence (t) at which the ancestral population gave rise to the descendant populations. Demographic rates represent parameters scaled to rates of molecular evolution: effective population sizes (Ne, individuals), population migration rates (Nm, migrants per generation) and divergence time (t, years). Population size (Ne) based on the average generation time (T) of 11 years for a high (0.9) annual adult survival rate. View Large Population expansion and approximate Bayesian computation (ABC) DIYABC analysis indicated that population growth of the YUC population (scenario 1) is the best-supported scenario (Fig. 4), with a higher posterior probability value and 95% confidence intervals (PP, 95% CI; scenario 1: 0.882, 0.869–0.896) that did not overlap with those obtained for scenario 2 (0.117, 0.104–0.131). Analyses to estimate confidence in scenario choice indicated that type I and type II errors for the best-supported scenario were low (0.354 and 0.383, respectively). Under scenario 1, posterior mean parameter estimates indicated that the divergence between dryCHI and YUC groups (t3) occurred 43 900 years before present (CI: 61.8–40.1 kya) and the split between humYUC and dryYUC (t1) was dated to 1510 (4.15–1.02 kya) (Appendix 2). In accordance with the population expansion model, the DIYABC analysis supported a scenario of a YUC population size increase around the LGM (t2 =15 600 years before present, 20.8–10.3 kya; Appendix 2). Estimated mean mutation rates of nrDNA and cpDNA were estimated to be 2.03 × 10–7 and 3.25 × 10–7, respectively (Appendix 2). Figure 4. View largeDownload slide (A) Schematic representation of two demographic scenarios of Psittacanthus mayanus tested with DIYABC and the prior values used in DIYABC simulations. N represents effective population sizes (Ne) and t indicates divergence time. Timeline from present to past indicates Mid Holocene (MH), Last Glacial Maximum (LGM) and Last Interglacial (LIG). (B) The posterior probability of scenarios was assessed using a weighted logistic regression on the 1% of simulated datasets closest to the observed data and, for the best-supported scenario (scenario 1). (C) A principal components analysis was used to test the goodness-of-fit of each evaluated scenario against the simulated data in DIYABC. Figure 4. View largeDownload slide (A) Schematic representation of two demographic scenarios of Psittacanthus mayanus tested with DIYABC and the prior values used in DIYABC simulations. N represents effective population sizes (Ne) and t indicates divergence time. Timeline from present to past indicates Mid Holocene (MH), Last Glacial Maximum (LGM) and Last Interglacial (LIG). (B) The posterior probability of scenarios was assessed using a weighted logistic regression on the 1% of simulated datasets closest to the observed data and, for the best-supported scenario (scenario 1). (C) A principal components analysis was used to test the goodness-of-fit of each evaluated scenario against the simulated data in DIYABC. Palaeodistribution modelling and environmental variation The ENM yielded a good fit for the current geographical distribution of P. mayanus. The AUC value for the replicate runs was 0.993 ± 0.002 (Fig. 5). Determination of the threshold probability for predicted presence using TSS resulted in a mean proportion of correctly classified training observations of 0.510 ± 0.229. The palaeodistribution modelling revealed that suitable habitat for P. mayanus was out of its current distribution during the LIG and restricted to small areas in southern Guatemala with moderate-to-high suitability (Fig. 5). However, the predictions for the LGM applying both CCSM and MIROC simulations revealed that conditions of suitable habitat potentially expanded the distribution of P. mayanus, particularly in Guatemala and Chiapas, and an isolated portion of Quintana Roo with high probability according to the CCSM model and, with a high probability, a more continuous distribution in Guatemala, Chiapas, Belize and Quintana Roo according to the MIROC model (Fig. 5). Conditions of suitable habitat were contracted during the MH, with suitable habitat restricted to small areas between Belize and Quintana Roo, Central Depression of Chiapas (CCSM), and the northernmost part of the Yucatán Peninsula. Overall, the comparison between present and past climatic conditions predicted by the models suggests that climatic conditions for P. mayanus suitable habitat experienced range shifts from southern Guatemala into the Yucatán Peninsula, with expansion from LIG to LGM, contraction from LGM to MH and colonization of the northernmost portion of the Yucatán Peninsula and expansion from MH to present. Figure 5. View largeDownload slide Results from the MAXENT analyses showing binary transformation of ecological niche models for Psittacanthus mayanus at Last Inter Glacial (LIG, 140–120 kya), Last Glacial Maximum (LGM, CCSM and MIROC scenarios, 21 kya), Mid-Holocene (MH, CCSM and MIROC scenarios, 6 kya), and at present. Vegetation map on the top right indicates the biogeographic regions discussed here (Morrone, 2005, 2014). Yellow indicates binary transformation of occurrence using the threshold that maximized the true skill statistics (TSS). Yellow dots shown on the right map below indicate occurrence records. Figure 5. View largeDownload slide Results from the MAXENT analyses showing binary transformation of ecological niche models for Psittacanthus mayanus at Last Inter Glacial (LIG, 140–120 kya), Last Glacial Maximum (LGM, CCSM and MIROC scenarios, 21 kya), Mid-Holocene (MH, CCSM and MIROC scenarios, 6 kya), and at present. Vegetation map on the top right indicates the biogeographic regions discussed here (Morrone, 2005, 2014). Yellow indicates binary transformation of occurrence using the threshold that maximized the true skill statistics (TSS). Yellow dots shown on the right map below indicate occurrence records. Four principal components (PC) explained 91.4% of the total environmental variance in P. mayanus (Appendix 3; Fig. 6). The PC1 (45.7% of total variation) was positively associated with temperature (BIO1, BIO6, BIO9–BIO11) and precipitation variables (BIO14, BIO17, BIO19) and negatively with precipitation seasonality (BIO15), whereas PC2 (23.9%) was positively associated with mean temperature of wettest quarter (BIO8) and negatively associated with precipitation variables (BIO12–BIO13, BIO18) (Appendix 3). Univariate ANOVAs of PC scores, using grouping by precipitation as a fixed factor, showed strong significant environmental differences among groups in PC1 scores (Appendix 3, Fig. 6) and differences between dryYUC and the other groups (dryCHI and humYUC) in PC2 scores (Appendix 3, Fig. 6). Figure 6. View largeDownload slide (A) Map relating to annual precipitation (BIO12) and location of precipitation differences on the Yucatán Peninsula. The approximate divisions of the dry Yucatán (dryYUC), humid Yucatán (humYUC) and dry Central Depression of Chiapas (dryCHI) are indicated with different colours. White dots shown on the map indicate collection localities and coloured dots indicate occurrence records. (B) Principal components analysis (PCA) on the 19 BIOCLIM variables with the first principal component (PC1) largely a measure of temperature differences, and the second principal component (PC2) mainly determined by precipitation measures. Symbols correspond to populations of Psittacanthus mayanus on the dryYUC (blue), humYUC (green), and dryCHI populations (yellow). Means (± 1 SE) of PC1 (C) and PC2 (D) scores with different superscript letters are significantly different between groups (P < 0.05). Figure 6. View largeDownload slide (A) Map relating to annual precipitation (BIO12) and location of precipitation differences on the Yucatán Peninsula. The approximate divisions of the dry Yucatán (dryYUC), humid Yucatán (humYUC) and dry Central Depression of Chiapas (dryCHI) are indicated with different colours. White dots shown on the map indicate collection localities and coloured dots indicate occurrence records. (B) Principal components analysis (PCA) on the 19 BIOCLIM variables with the first principal component (PC1) largely a measure of temperature differences, and the second principal component (PC2) mainly determined by precipitation measures. Symbols correspond to populations of Psittacanthus mayanus on the dryYUC (blue), humYUC (green), and dryCHI populations (yellow). Means (± 1 SE) of PC1 (C) and PC2 (D) scores with different superscript letters are significantly different between groups (P < 0.05). DISCUSSION Our survey of molecular sequence data in populations of P. mayanus identified two genetic groups in the Yucatán Peninsula. Phylogeographic and distributional projections derived from ecological niche modelling are consistent with the hypothesis that the directionality of colonization was from southern Chiapas to the north on the Yucatán Peninsula and that the Altos de Chiapas mountain areas may have restricted northward gene flow. Together, our results suggest a recent geographical northern expansion into the seasonally tropical deciduous forest, consistent with a colonization model and a pattern of population expansions/contractions for species responding to the Pleistocene glacial cycles and reduced gene flow since by both physical and environmental barriers. Given that this region is poorly studied at the phylogeographic level, the phylogeography of P. mayanus in combination with a niche-modelling approach is particularly important in understanding the evolution of the biodiversity of Mexican lowland forests. Genetic diversity, structure and phylogeographic analysis Genetic diversity analysis for P. mayanus populations revealed a higher gene diversity (h) and nucleotide diversity (π) in Chiapas than in the northern part of the distribution on the Yucatán Peninsula (Table 5 and Fig. 2). The genetic diversity within populations and total gene diversity of P. mayanus showed that the species has low levels of genetic diversity relative to that observed in other mistletoe species (Zuber & Widmer, 2009; Amico et al., 2014; Nyagumbo, Ndagurwa & Mushonga, 2017) or that of closely related and more widespread congeners (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Because the fruits of P. mayanus are consumed and dispersed by several bird species, including the endemic (to the Yucatán Peninsula) Patagioenas cayennensis and Piranga roseogularis and the widespread Myiozetetes similis Pitangus sulphuratus and Megarhynchus pitangua, it is possible that the observed genetic structure and genetic differentiation have been eroded by large-scale seed recruitment and high historical rates of gene flow (Ornelas et al., 2016). The haplotype network and AMOVA results suggest genetic differentiation between two groups of populations for P. mayanus, from Chiapas and the Yucatán Peninsula. Similar results were also observed for Guaiacum sanctum L. (Zygophyllaceae) trees using nuclear microsatellites (Oyama et al., 2016). The strongest genetic differentiation between the groups of P. mayanus populations in Chiapas, located approximately between 16 and 17ºN, and the Yucatán Peninsula, located between 18 and 21ºN, corresponds to historically different biogeographic regions, and occurs between 600 and 1000 m a.s.l. to nearly above sea level, respectively (Table 1). The population groups of P. mayanus showed, however, haplotype sharing and short-term divergence (Figs 1 and 2). The widespread distribution and higher frequency of haplotype H1 (trnL-F) is not surprising considering that several widespread bird species mediate mistletoe seed dispersal. Migration among groups of populations occurred in one direction, moving from Chiapas to the Yucatán Peninsula, suggesting that the barriers between these groups of populations associated with different geographical regions and environmental conditions are permeable. A lower and similar genetic differentiation among groups of populations (haplotype network and AMOVA results) was found in the ITS when populations of P. mayanus were separated by geography and precipitation into three groups, respectively (Tables 1, 4–5, Fig. 2). It is possible that these habitat differences are related to the formation of genetic groups. PC1 and PC2 explained 70% of the total variance in current climatic conditions among groups of populations, respectively. PC1 was interpretable as an index of precipitation variability. For example, PC1 exhibited high positive loading for annual precipitation variables and low or negative loading for precipitation seasonality and temperature variables. Thus, high scores of PC1 are indicative of wet climatic conditions and low or negative PC1 scores of seasonal climatic conditions. Notably, the climatic PC1 scores differed between groups. Thus, in this dataset at least, the geographical position is relatively dependent of precipitation variability and seasonality. The fact that there are some environmental differences among areas does not necessarily point to adaptation as the driver of genetic differentiation, and thus these interpretations are correlative. Instead, the genetic structure may be a consequence of historical fragmentation in the region, generating responses to contractions and expansions of the semideciduous tropical dry forest along the Yucatán Peninsula during Pleistocene climatic cycles as suggested by the ENMs. Demographic history of P. mayanus Despite the observed genetic differentiation between groups of populations, the most frequent ribotypes and haplotypes of P. mayanus populations are widespread, suggesting effective gene flow between peninsular and Chiapas populations and signs of recent demographic expansion for populations from the northern part of the Yucatán Peninsula. Group of populations from Chiapas (CHI) showed no major changes in their effective population sizes, exhibiting a high genetic diversity; whereas group of populations of the northern portion of the Yucatán Peninsula (YUC) exhibited lower genetic diversity, no changes in effective population sizes and strong signatures of recent demographic expansion (Table 5, Fig. 3). Note that the historical population size changes or the spatial population expansion can be influenced by several sources of error (e.g. sample size, sampling scheme, levels of polymorphism, mutation rates) or by genetic hitchhiking of polymorphism linked to genes under positive selection (e.g. Grant, 2015). The levels of genetic differentiation found in P. mayanus reveal the presence of two genetic groups that corresponds with the distinction between samples from the Yucatán Peninsula and Chiapas, and that the pattern of differentiation is associated with ecological/elevational isolation. Our divergence time estimates indicate that the split between these groups of populations occurred during the Pleistocene, in agreement with divergence time estimates for lineages of other Psittacanthus spp. (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Therefore, Pleistocene climatic oscillations are invoked as the main factor driving mistletoe diversification in the region. The ENMs and gene flow estimates using IMa suggest that dispersal patterns among populations of P. mayanus occurred asymmetrically; the directionality of gene flow indicates that populations from Chiapas are the source of migrants into the Yucatán Peninsula. The facts that several samples from Chiapas clustered sparsely with samples of P. mayanus from the Yucatán Peninsula and no samples from the Yucatán Peninsula of P. mayanus clustered with samples of Chiapas is evidence of gene flow and suggests that the directionality of colonization was from southern Chiapas to the north of the Yucatán Peninsula and that the Altos de Chiapas mountain areas may have restricted northward gene flow. Approximate Bayesian computation (ABC) analyses strongly supported a scenario of post-glacial population growth and that the northern part of the Yucatán Peninsula was invaded by ancestral P. mayanus individuals from Chiapas. This scenario suggests that species-specific traits, availability of host tree species and/or variation of suitable habitat may be causing the phylogeographical pattern influenced by the Pleistocene glacial-interglacial cycles. A recent phylogeographic study of the bird-dispersed mistletoe Phoradendron californicum Nutt. (Santalaceae) in the Sonoran Desert showed that the overall plastid DNA variation was best explained by large-scale historical geographical events and changes in vegetation limiting seed dispersal at deeper evolutionary timescales but host species poorly explained plastid DNA variation (Lira-Noriega et al., 2015). More recent phylogeographic studies of Psittacanthus spp. using nrDNA and plastid DNA molecular markers also provide support for the predominant role of isolation and environmental factors in driving genetic differentiation between populations of P. schiedeanus from the cloud forests of eastern Mexico (Ornelas et al., 2016) and P. calyculatus from the temperate forests along the Trans-Mexican Volcanic Belt (Pérez-Crespo et al., 2017). Using nuclear microsatellite data, Ramírez-Barahona et al. (2017) showed that genetic differentiation among populations of the species complex composed of P. schiedeanus, P. calyculatus, P. breedlovei and P. angustifolius Kuijt was associated with environmental and host preferences, independent of geography. However, the observed genetic structure was best explained by environmental predictors rather than by host preferences, supporting the hypothesis that the occurrence of these parasites has been largely determined by its own climatic niche and to a lesser degree by host specificity. The palaeodistribution of P. mayanus on the Yucatán Peninsula The ENM predicted fragmentation in the southern regions of P. mayanus distribution and recent colonization of the northernmost part of the Yucatán Peninsula (Fig. 5). ENM results indicate that the range of P. mayanus expanded from the LIG to the LGM; it was wider and connected during the LGM and contracted during the MH and it is expanding under the warmer current post-glacial period (Fig. 5). The predicted LIG area is small. It is possible that the distribution of P. mayanus during this period was limited, but also that the species occurred in environments without a present analogue. The LIG palaeodistribution model shows a lack of suitable habitat in most of the Yucatán Peninsula 140 000–120 000 years ago. Based on the phylogenetic evidence and divergence time estimation suggesting that the species originated before the LIG, we do not interpret this ENM result as a hypothesis of complete absence of the species from the northern Central America. In the absence of any evidence to suggest a South American origin for the species ancestor we offer as an alternative hypothesis that the distribution of P. mayanus in the LIG was quite small or restricted to a different habitat in northern Central America. Other Psittacanthus spp. generally associated with habitats at higher elevation show a wider distribution during the LIG (Ornelas et al., 2016). Comparisons between lowland and highland mistletoe species will be an area for future palaeodistribution modelling in which development of additional LGM climate models may produce novel hypotheses. The two global circulation model simulations of the LGM climate (CCSM and MIROC) revealed dissimilar inferences regarding the potential palaeodistribution of P. mayanus. The CCSM model inferred two geographically separated regions with high suitability expanding northwards into partially the current distribution, one in the Pacific coast and lowlands of Chiapas and Guatemala, in the borderlands of southern Oaxaca and southern Guatemala, and the other one in the eastern side of the Yucatán Peninsula, stretching from the north to the south Caribbean slope of Quintana Roo. In contrast, the MIROC model inferred a single large area in the southern part of the Yucatán Peninsula including most of Chiapas, southern Tabasco, southern Quintana Roo, northern Guatemala and Belize. A common feature of both models, however, is the inference of suitable habitat in Chiapas and southern Quintana Roo (Fig. 5). The CCSM and MIROC models differ in temperature and precipitation patterns. The LGM climate as simulated by CCSM is colder and dryer than that simulated by MIROC. The CCSM model is consistent with the tropical refugia theory (Haffer, 1969), in which the cooler glacial climates characterized by a drastic reduction in precipitation and widespread aridity during glacial times caused tropical forests to become contracted in their distribution and fragmented into isolates of tropical forest patches separated by intervening non-forest vegetation (e.g. Haffer, 1969; Carnaval & Moritz, 2008). The knowledge on precipitation during the LGM, and hence on the distribution of Neotropical forests, is largely based on limited and conflicting palaeoecological data (Ramírez-Barahona & Eguiarte, 2013). It has been suggested that the summer precipitation regime was collapsed, resulting in arid conditions extending over much of Mexico, Central America and the Caribbean (Ramírez-Barahona & Eguiarte, 2013 and references therein). Accordingly, Anselmetti et al. (2006) suggested a reduction in wet season precipitation for upper Central America based on a reduction in the water volume of Lake Petén-Itzá, Guatemala. In contrast, Mueller et al. (2010) suggested that there was no effective decrease in humidity in upper Central America based on pollen records from Lake Petén-Itzá showing a rich clay deposit and a pollen diagram indicative of a montane forest community thriving at lower elevations under cool and moist conditions (Hodell et al., 2008; Bush et al., 2009). Both LGM models predict expanded suitable habitat compared with LIG prediction. Based on palynological evidence for the Petén region (Leyden, 1984), xeric vegetation was widely distributed during the LGM, which might have facilitated expansion of P. mayanus to coastal and central depression of Chiapas and Quintana Roo, and then it was supplanted by more mesic vegetation towards the Holocene transition and contraction of P. mayanus. Therefore, it is possible that the uncovered distribution on Quintana Roo, Guatemala and Belize had provided suitable conditions for P. mayanus during the LGM, MH and present. According to the work of Carnaval et al. (2009) on other tropical regions, this would correspond to an area of long-term population persistence, and thus a source of genetic variability during population expansion. Although we did not sample these areas for the present study, the two models identified similar areas of suitable habitat and the use of the two different climate models enabled us to identify the possibility of an area of long-term population persistence and to account for modelling uncertainty due to LGM climate based on limited and conflicting palaeoecological data (Ramírez-Barahona & Eguiarte, 2013). The distribution of potential suitable habitat appears to affect the distribution of P. mayanus, causing population differentiation. Several populations from both the northern and southern distribution of P. mayanus possess unique ribotypes and haplotypes, suggesting that the population of P. mayanus persisted in the southern region of the Yucatán Peninsula under LGM conditions where the ecological conditions were appropriate, and that the northern part of the Yucatán Peninsula was colonized more recently (MH) and contracted its distribution particularly in the Central Depression of Chiapas. The IMa results revealed that the divergence or splitting time (t) between descendant YUC and CHI populations overlaps with the tMRCA estimated using DIYABC. When testing for migration, migration among populations occurred predominantly from Chiapas to the Yucatán Peninsula. These estimates suggest a recent divergence scenario, in which haplotype sharing between regions results from retention of ancestral polymorphisms with not enough time for lineage sorting, ongoing gene flow across habitats or complete haplotype turnover. Along with these results, the haplotype network, FST pairwise comparisons and AMOVAs suggest moderate levels of genetic diversity among populations and the existence of three genetic groups when populations are grouped by the geography (YUC, PET, CHI) or by the precipitation gradient (dryYUC, humYUC, dryCHI), especially when considering the nrDNA (Tables 3 and 4). Altogether these results suggest that range shifts and changes in the habitat distribution that occurred during the Pleistocene climatic cycles might have fragmented the distribution of P. mayanus. According to the ENM, IMa and ABC results, the genetic differentiation of P. mayanus may be linked to the effects of the glacial/interglacial cycles and environmental factors. Under these scenarios, ancestral populations of P. mayanus colonized the semideciduous tropical rain forest of the Petén province during the LGM followed by northward expansion and colonization of the northern Yucatán province during the MH, currently covered by seasonally dry tropical deciduous forests, when the drier and cooler climate of the Peninsula and the savanna and Juniperus-matorral covering most of the region until the early Holocene (Islebe et al., 1996; Orellana et al., 2003) changed to a more humid and tropical-type forest (Metcalfe et al., 2000). Interpreting the resultant P. mayanus range dynamics through time is difficult without considering contemporaneous demography/occurrence of its principal host populations. Specifically, the putative recent expansion of this mistletoe into northern Yucatán might simply reflect broader changes in the host plants upon which it ultimately depends. If so, each of the genetic groups needs to be modelled by adding distribution data of each of the several host species to examine direct versus indirect determinants of distribution. Although such an approach would be the norm in phylogeographic studies of host-specialized parasitic plants, no data are available on the prevalence or specificity of these mistletoes for each of the recorded host species (Kuijt, 2009). Considering that parasitic plant distributions ultimately depend on host availability, future ENM should evaluate the effects of including the host distribution to predict distribution shifts of the mistletoe under past climates comparing models that use only climate variables with models that also take into account the biotic interactions with the hosts. CONCLUSIONS Our results revealed the overall phylogeographical patterns in P. mayanus, with two main haplogroups comprising groups of populations from the Yucatán Peninsula and Chiapas and low genetic differentiation between the Yucatán and Petén provinces. Despite the genetic differentiation of some P. mayanus populations, there are widespread ribotypes and haplotypes in the ITS and trnL-F markers, respectively. In particular, the low level of variation found in the ITS and trnL-F data and the presence of widespread ribotypes and haplotypes on the Yucatán Peninsula suggest effective nuclear and plastid gene flow, and negative values of Tajima’s D, Fu’s Fs and the results of mismatch distributions and BSPs are suggestive of recent demographic expansion for populations of P. mayanus in the region, most likely during the LGM. In contrast, the population differentiation of P. mayanus from Chiapas shown by ITS and plastid trnL-F datasets suggests that gene flow has been restricted. The ABC analysis strongly supported a scenario of population growth. Under this scenario, climatic fluctuations throughout the Pleistocene would have altered the distribution of suitable habitat for mistletoes along the Yucatán Peninsula leading to genetic differentiation of P. mayanus comprising moderate-to-low haplotype diversity and hinting widespread ribotypes and haplotypes lay on the Yucatán Peninsula involving geographical northern expansion into the seasonally tropical deciduous forest, consistent with a colonization model and a pattern of population expansions/contractions for species responding to the Pleistocene glacial cycles. ACKNOWLEDGEMENTS We thank C. Bárcenas, E. Gándara, E. Ruiz Sánchez, H. Gómez Domínguez and A. J. Vargas who helped in obtaining samples and/or generated sequences for this work. Permission to conduct our fieldwork was granted by the Mexican Government (Instituto Nacional de Ecología, Secretaría del Medio Ambiente y Recursos Naturales, SGPA/DGGFS/712/1299/12). We also thank C. Palma-Silva and two anonymous reviewers for their valuable criticisms and suggestions. This project was funded by a competitive grant (155686) from the Consejo Nacional de Ciencia y Tecnología (CONACyT; http://www.conacyt.mx) and research funds from the Departamento de Biología Evolutiva, Instituto de Ecología, A.C. (20030/10563) awarded to Juan Francisco Ornelas. A doctoral scholarship was granted from CONACyT to Y.L.V. (155686) and A.E.O.R. (262563). REFERENCES Allouche O, Tsoar A, Kadmon R. 2006. Assessing the accuracy of species distribution models: prevalence, kappa, and the true skill statistic (TSS). Journal of Applied Ecology  43: 1223– 1232. Google Scholar CrossRef Search ADS   Amico GC, Vidal-Russell R, Nickrent DL. 2007. Phylogenetic relationships and ecological speciation in the mistletoe Tristerix (Loranthaceae): the influence of pollinators, dispersers, and hosts. American Journal of Botany  94: 558– 567. Google Scholar CrossRef Search ADS PubMed  Amico GC, Vidal-Russell R, García MA, Nickrent DL. 2012. Evolutionary history of the South American mistletoe Tripodanthus (Loranthaceae) using nuclear and plastid markers. Systematic Botany  37: 218– 225. Google Scholar CrossRef Search ADS   Amico GC, Vidal-Russell R, Aizen MA, Nickrent D. 2014. Genetic diversity and population structure of the mistletoe Tristerix corymbosus (Loranthaceae). Plant Systematics and Evolution  300: 153– 162. Google Scholar CrossRef Search ADS   Anselmetti F, Ariztegui D, Hodell DA, Hillesheim MB, Brenner M, Gilli A, McKenzie JA, Mueller AD. 2006. Late quaternary climate-induced lake level variations in lake Petén Itzá, Guatemala, inferred from seismic stratigraphic analysis. Palaeogeography, Palaeoclimatology, Palaeoecology  230: 52– 69. Google Scholar CrossRef Search ADS   Beaumont MA, Zhang W, Balding DJ. 2002. Approximate Bayesian computation in population genetics. Genetics  162: 2025– 2035. Google Scholar PubMed  Braconnot P, Otto-Bliesner B, Harrison S, Joussaume S, Peterschmitt JY, Abe-Ouchi A, Crucifix M, Driesschaert E, Fichefet T, Hewitt CD, Kageyama M, Kitoh A, Loutre MF, Marti O, Merkel U, Ramstein G, Valdes P, Weber L, Yu Y, Zhao Y. 2007. Results of PMIP2 coupled simulations of the mid-holocene and last glacial maximum – part 2: feedbacks with emphasis on the location of the ITCZ and mid- and high latitudes heat budget. Climate of the Past  3: 279– 296. Google Scholar CrossRef Search ADS   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: Vimeuz F, Sylvestre F, Khodri M, eds. Past climate variability in South America and surrounding regions, from the Last Glacial Maximum to the Holocene , Vol. 14. Houten: Springer, 113– 128. Google Scholar CrossRef Search ADS   Carnaval AC, Moritz C. 2008. Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest. Journal of Biogeography  35: 1187– 1201. Google Scholar CrossRef Search ADS   Carnaval AC, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C. 2009. Stability predicts genetic diversity in the Brazilian Atlantic Forest hotspot. Science  23: 785– 789. Google Scholar CrossRef Search ADS   Cavers S, Navarro C, Lowe AJ. 2003. Chloroplast DNA phylogeography reveals colonization history of a Neotropical tree, Cedrela odorata L., in Mesoamerica. Molecular Ecology  12: 1451– 1460. Google Scholar CrossRef Search ADS PubMed  Clement M, Posada D, Crandall KA. 2000. TCS: a computer program to estimate gene genealogies. Molecular Ecology  9: 1657– 1659. Google Scholar CrossRef Search ADS PubMed  Cornuet JM, Santos F, Beaumont MA, Robert CP, Marin JM, Balding DJ, Guillemaud T, Estoup A. 2008. Inferring population history with DIYABC: a user-friendly approach to approximate Bayesian computation. Bioinformatics  24: 2713– 2719. Google Scholar CrossRef Search ADS PubMed  Cornuet JM, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, Leblois R, Marin JM, Estoup A. 2014. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics  30: 1187– 1189. Google Scholar CrossRef Search ADS PubMed  Díaz Infante S, Lara C, Arizmendi MC, Eguiarte LE, Ornelas JF. 2016. Reproductive ecology and isolation of Psittacanthus calyculatus and P. auriculatus mistletoes (Loranthaceae). PeerJ  4: e2491. Google Scholar CrossRef Search ADS PubMed  Doyle JJ, Doyle JL. 1987. A rapid DNA isolation procedure from small quantities of fresh leaf tissue. Phytochemical Bulletin, Botanical Society of America  19: 11– 15. Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology  7: 214. Google Scholar CrossRef Search ADS PubMed  Drummond AJ, Rambaut A, Shapiro B, Pybus OG. 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution  22: 1185– 1192. Google Scholar CrossRef Search ADS PubMed  Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ. 2011. A statistical explanation of MaxEnt for ecologists. Diversity and Distributions  17: 43– 57. Google Scholar CrossRef Search ADS   Espadas-Manrique C, Durán R, Argáez J. 2003. Phytogeographic analysis of taxa endemic to the Yucatán Peninsula using geographic information systems, the domain heuristic method and parsimony analysis of endemicity. Diversity and Distributions  9: 313– 330. Google Scholar CrossRef Search ADS   Excoffier L, Smouse PE, Quattro JM. 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics  131: 479– 491. Google Scholar PubMed  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. Fielding AH, Bell JF. 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation  24: 38– 49. Google Scholar CrossRef Search ADS   Fontaine MC, Austerlitz F, Giraud T, Labbé F, Papura D, Richard-Cervera S, Delmotte F. 2013. Genetic signature of a range expansion and leap-frog event after the recent invasion of Europe by the grapevine downy mildew pathogen Plasmopara viticola. Molecular Ecology  22: 2771– 2786. Google Scholar CrossRef Search ADS PubMed  Fu YX. 1997. Statistical neutrality of mutations against population growth, hitchhiking and background selection. Genetics  147: 915– 925. Google Scholar PubMed  González C, Ornelas JF, Gutiérrez-Rodríguez C. 2011. Selection and geographic isolation influence hummingbird speciation: genetic, acoustic and morphological divergence in the wedge-tailed sabrewing (Campylopterus curvipennis). BMC Evolutionary Biology  11: 38. Google Scholar CrossRef Search ADS PubMed  Grant WS. 2015. Problems and cautions with sequence mismatch analysis and bayesian skyline plots to infer historical demography. The Journal of Heredity  106: 333– 346. Google Scholar CrossRef Search ADS PubMed  Guevara-Chumacero LM, López-Wilchis R, Pedroche FF, Juste J, Ibáñez C, Barriga-Sosa IDLA. 2010. Molecular phylogeography of Pteronotus davyi (Chiroptera: Mormoopidae) in Mexico. Journal of Mammalogy  91: 220– 232. Google Scholar CrossRef Search ADS   Gutiérrez-García TA, Vázquez-Domínguez E. 2012. Biogeographically dynamic genetic structure bridging two continents in the monotypic Central American rodent Ototylomys phyllotis. Biological Journal of the Linnean Society  107: 593– 610. Google Scholar CrossRef Search ADS   Haffer J. 1969. Speciation in Amazonian forest birds. Science  165: 131– 137. Google Scholar CrossRef Search ADS PubMed  Hall TA. 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium Series  41: 95– 98. Hammer O. 2011. PAST (Paleontological Statistics) . Natural History Museum, University of Oslo. Available at: http://folk.uio.no/ohammer/past/index.html. Harpending RC. 1994. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Human Biology  66: 591– 600. Google Scholar PubMed  Haug GH, Günther D, Peterson LC, Sigman DM, Hughen KA, Aeschlimann B. 2003. Climate and the collapse of Maya civilization. Science  299: 1731– 1735. Google Scholar CrossRef Search ADS PubMed  Hey J, Nielsen R. 2004. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics  167: 747– 760. Google Scholar CrossRef Search ADS PubMed  Hey J, Nielsen R. 2007. Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proceedings of the National Academy of Sciences USA  104: 2785– 2790. Google Scholar CrossRef Search ADS   Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high resolution interpolated climate surfaces from global land areas. International Journal of Climatology  25: 1965– 1978. Google Scholar CrossRef Search ADS   Hodell DA, Anselmetti FS, Ariztegui D, Brenner M, Curtis JH, Gilli A, Grzesik DA, Guilderson TJ, Müller AD, Bush MB, Correa-Metrio A, Escobar J, Kutterolf S. 2008. An 85-ka record of climate change in lowland Central America. Quaternary Science Reviews  27: 1152– 1165. Google Scholar CrossRef Search ADS   Howell TR. 1972. Birds of the lowland pine savannah of northeastern Nicaragua. The Condor  74: 316– 340. Google Scholar CrossRef Search ADS   Islebe GA, Hooghiemstra H, Brenner M, Curtis JH, Hodell D. 1996. A Holocene vegetation history from lowland Guatemala. The Holocene  6: 265– 271. Google Scholar CrossRef Search ADS   Jiménez RA, Ornelas JF. 2016. Historical and current introgression in a Mesoamerican hummingbird species complex: a biogeographic perspective. PeerJ  4: e1556. Google Scholar CrossRef Search ADS PubMed  Kay KM, Whittall JB, Hodges SA. 2006. A survey of nuclear ribosomal internal transcribed spacer substitution rates across angiosperms: an approximate molecular clock with life history effects. BMC Evolutionary Biology  6: 36. Google Scholar CrossRef Search ADS PubMed  Kuijt J. 2009. Monograph of Psittacanthus (Loranthaceae). Systematic Botany Monographs  86: 1– 361. Lande R, Engen S, Sæther BE. 2003. Stochastic population dynamics in ecology and conservation . Oxford: Oxford University Press. Google Scholar CrossRef Search ADS   Landis JR, Koch GG. 1977. The measurement of observer agreement for categorical data. Biometrics  33: 159– 174. Google Scholar CrossRef Search ADS PubMed  Leyden BW. 1984. Guatemalan forest synthesis after Pleistocene aridity. Proceedings of the National Academy of Sciences, USA  81: 4856– 4859. Google Scholar CrossRef Search ADS   Li W, Wang Z, Ma Z, Tang H. 1997. A regression model for the spatial distribution of red-crown crane in Yancheng Biosphere Reserve, China. Ecological Modelling  103: 115– 121. Google Scholar CrossRef Search ADS   Librado P, Rozas J. 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics (Oxford, England)  25: 1451– 1452. Google Scholar CrossRef Search ADS PubMed  Lira-Noriega A, Toro-Núñez O, Oaks JR, Mort ME. 2015. The roles of history and ecology in chloroplast phylogeographic patterns of the bird-dispersed plant parasite Phoradendron californicum (Viscaceae) in the Sonoran Desert. American Journal of Botany  102: 149– 164. Google Scholar CrossRef Search ADS PubMed  Liu C, Berry PM, Dawson TP, Pearson RG. 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography  28: 385– 393. Google Scholar CrossRef Search ADS   López de Buen L, Ornelas JF. 2002. Host compatibility of the cloud forest mistletoe Psittacanthus schiedeanus (Loranthaceae) in central Veracruz, Mexico. American Journal of Botany  89: 95– 102. Google Scholar CrossRef Search ADS PubMed  López Ramos E. 1975. Geological summary of the Yucatán Peninsula. In: Nairin AEM, Stehli FG, eds. The ocean basis and margins. Volume 3. The Gulf of Mexico and the Caribbean . New York: Plenum Press, 257– 282. Lundell CL. 1934. Preliminary sketch of the phytogeography of the Yucatán Peninsula. Contributions to American Archaeology  12: 257– 321. 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. Google Scholar CrossRef Search ADS   Morrone JJ. 2005. Hacia una síntesis biogeográfica de México. Revista Mexicana de Biodiversidad  76: 207– 252. Morrone JJ. 2014. Biogeographical regionalisation of the Neotropical region. Zootaxa  3782: 1– 110. Google Scholar CrossRef Search ADS PubMed  Mueller AD, Anselmetii FS, Ariztequi D, Brenner M, Hodell DA, Curtis JH, Escobar J, Gilli A, Grzesik DA, Guilderson TP, Kutterolf S, Plötze M. 2010. Late Quaternary palaeoenvironment of northern Guatemala: evidence from deep drill cores and seismic statigraphy of Lake Petén Itzá. Sedimentology  57: 1220– 1245. Nyagumbo E, Ndagurwa HGT, Mushonga K. 2017. High genetic variation and low gene flow among populations of Viscum verrucosum in semi-arid savannah, southwest Zimbabwe. Ecological Genetics and Genomics 3–5 : 41– 46. 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. Ornelas JF, Gándara E, Vásquez-Aguilar AA, Ramírez-Barahona S, Ortiz-Rodriguez AE, González C, Mejía S aules MT, Ruiz-Sanchez E. 2016. A mistletoe tale: postglacial invasion of Psittacanthus schiedeanus (Loranthaceae) to Mesoamerican cloud forests revealed by molecular data and species distribution modeling. BMC Evolutionary Biology  16: 78. Google Scholar CrossRef Search ADS PubMed  Oyama K, Martínez-Ramos M, Peñaloza-Ramírez JM, Rocha-Ramírez V, Armenta-Medina EG, Hernández-Soto P. 2016. Population genetic structure of an extremely logged tree species Guaiacum sanctum L. in the Yucatán Peninsula, Mexico. Botanical Sciences  94: 345– 356. Google Scholar CrossRef Search ADS   Padilla y Sánchez RJ. 2007. Evolución geológica del sureste mexicano desde el Mesozoico al presente en el context regional del Golfo de México. Boletín de la Sociedad Geológica Mexicana  59: 19– 42. Google Scholar CrossRef Search ADS   Pérez-Crespo MJ, Ornelas JF, González-Rodríguez A, Ruiz-Sanchez E, Vásquez-Aguilar AA, Ramírez-Barahona S. 2017. Phylogeography and population differentiation in the Psittacanthus calyculatus (Loranthaceae) mistletoe: a complex scenario of climate-volcanism interaction along the Trans-Mexican Volcanic Belt. Journal of Biogeography , 44: 2501– 2514. Google Scholar CrossRef Search ADS   Pfenninger M, Posada D. 2002. Phylogeographic history of the land snail Candidula unifasciata (Helicellinae, Stylommatophora): fragmentation, corridor migration, and secondary contact. Evolution; International Journal of Organic Evolution  56: 1776– 1788. Google Scholar CrossRef Search ADS PubMed  Phillips SJ, Anderson RP, Schapire RE. 2006. Maximum entropy modelling of species geographic distributions. Ecological Modelling  190: 231– 259. Google Scholar CrossRef Search ADS   Poelchau MF, Hamrick JL. 2012. Differential effects of landscape-level environmental features on genetic structure in three codistributed tree species in Central America. Molecular Ecology  21: 4970– 4982. Google Scholar CrossRef Search ADS PubMed  Poelchau MF, Hamrick JL. 2013a. Comparative phylogeography of three common Neotropical tree species. Journal of Biogeography  40: 618– 631. Google Scholar CrossRef Search ADS   Poelchau MF, Hamrick JL. 2013b. Palaeodistribution modelling does not support disjunct Pleistocene refugia in several Central American plant taxa. Journal of Biogeography  40: 662– 675. Google Scholar CrossRef Search ADS   Pons O, Petit RJ. 1996. Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics  144: 1237– 1245. Google Scholar PubMed  Posada D. 2008. jModelTest: phylogenetic model averaging. Molecular Biology and Evolution  25: 1253– 1256. Google Scholar CrossRef Search ADS PubMed  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 and Evolution  3: 725– 738. Google Scholar CrossRef Search ADS PubMed  Ramírez-Barahona S, González C, González-Rodríguez A, Ornelas JF. 2017. The influence of climatic niche preferences on the population genetic structure of a mistletoe species complex. New Phytologist  214: 1751– 1761. Google Scholar CrossRef Search ADS PubMed  Ramírez-Barahona S, Torres-Miranda A, Palacios-Ríos M, Luna-Vega I. 2009. Historical biogeography of the Yucatán Peninsula, Mexico: a perspective from ferns (Monilophyta) and lycopods (Lycophyta). Biological Journal of the Linnean Society  98: 775– 786. Google Scholar CrossRef Search ADS   Richardson JE, Pennington RT, Pennington TD, Hollingsworth PM. 2001. Rapid diversification of a species-rich genus of Neotropical rain forest trees. Science (New York, N.Y.)  293: 2242– 2245. Google Scholar CrossRef Search ADS PubMed  Robert CP, Cornuet JM, Marin JM, Pillai NS. 2011. Lack of confidence in approximate Bayesian computation model choice. Proceedings of the National Academy of Sciences USA  108: 15112– 15117. Google Scholar CrossRef Search ADS   Rodríguez-Gómez F, Ornelas JF. 2014. Genetic divergence of the Mesoamerican azure-crowned hummingbird (Amazilia cyanocephala) across the Motagua–Polochic–Jocotán fault system. Journal of Zoological Systematics and Evolutionary Research  52: 142– 153. Google Scholar CrossRef Search ADS   Rogers AR, Harpending H. 1992. Population growth makes waves in the distribution of pairwise genetic differences. Molecular Biology and Evolution  9: 552– 569. Google Scholar PubMed  Simpson GG. 1964. Species density of North American recent mammals. Systematic Zoology  12: 57– 73. Google Scholar CrossRef Search ADS   Smith BT, Escalante P, Hernández B años BE, Navarro-Sigüenza AG, Rohwer S, Klicka J. 2011. The role of historical and contemporary processes on phylogeographic structure and genetic diversity in the Northern Cardinal, Cardinalis cardinalis. BMC Evolutionary Biology  11: 136. Google Scholar CrossRef Search ADS PubMed  Schneider S, Excoffier L. 1999. Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA. Genetics  152: 1079– 1089. Google Scholar PubMed  Taberlet P, Gielly L, Pautou G, Bouvet J. 1991. Universal primers for amplification of three non-coding regions of chloroplast DNA. Plant Molecular Biology  17: 1105– 1109. Google Scholar CrossRef Search ADS PubMed  Tajima F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics  123: 585– 595. Google Scholar PubMed  Vázquez-Domínguez E, Arita HT. 2010. The Yucatán peninsula: biogeographical history 65 million years in the making. Ecography  33: 212– 219. Torrescano-Valle N, Islebe GA. 2015. Holocene paleoecology, climate history and human influence in the southwestern Yucatán Peninsula. Review of Paleobotany and Palynology  217: 1– 8. Google Scholar CrossRef Search ADS   Vidal-Russell R, Nickrent DL. 2007. A molecular phylogeny of the feathery mistletoe Misodendrum. Systematic Botany  32: 560– 568. Google Scholar CrossRef Search ADS   Vidal-Russell R, Nickrent DL. 2008a. Evolutionary relationships in the showy mistletoe family (Loranthaceae). American Journal of Botany  95: 1015– 1029. Google Scholar CrossRef Search ADS   Vidal-Russell R, Nickrent DL. 2008b. The first mistletoes: origins of aerial parasitism in Santalales. Molecular Phylogenetics and Evolution  47: 523– 537. Google Scholar CrossRef Search ADS   Williford D, Deyoung RW, Honeycutt RL, Brennan LA, Hernández F. 2016. Phylogeography of the bobwhite (Colinus) quails. Wildlife Monographs  193: 1– 49. Google Scholar CrossRef Search ADS   Wilson CA, Calvin CL. 2006. An origin of aerial branch parasitism in the mistletoe family, Loranthaceae. American Journal of Botany  93: 787– 796. Google Scholar CrossRef Search ADS PubMed  Zuber D, Widmer A. 2009. Phylogeography and host race differentiation in the European mistletoe (Viscum album L.). Molecular Ecology  18: 1946– 1962. Google Scholar CrossRef Search ADS PubMed  APPENDICES Appendix 1. GenBank accession numbers of DNA sequences of other Psittacanthus species and representatives of Loranthaceae from Wilson & Calvin (2006), Amico, Vidal-Russell & Nickrent (2007), Vidal-Russell & Nickrent (2007, 2008a, b), Díaz Infante et al. (2016), Ornelas et al. (2016) and Pérez-Crespo et al. (2017) used in this study. Species  Voucher  Country  ITS  trnL-F  Alepis flavida Tiegh.  Calvin and Wilson NZ98- 04 (NA)  New Zealand  DQ333847  DQ340598  Atkinsonia ligustrina F.Muell.  Calvin and Wilson AU00-01 (NA)  Australia  DQ333865  DQ788714  Cladocolea cupulata Kuijt  Calvin and Wilson MX03- 08 (NA)  Mexico  DQ333861  DQ340612  Cladocolea mcvaughii Kuijt  Calvin and Wilson MX03- 09 (NA)  Mexico  DQ333860  DQ340611  Desmaria mutabilis (Poepp. & Endl.) Tiegh. ex B.D.Jacks.  Calvin and Wilson CL03-07 (NA)  Chile  DQ333852  DQ340603  Gaiadendron punctatum (Ruiz & Pav.) G.Don.  Calvin and Wilson CR01- 08 (NA)  Costa Rica  DQ333866  DQ340617  Ligaria cuneifolia (Ruiz & Pav.) Tiegh.  Calvin and Wilson CL03-01 (NA)  Chile  DQ333853  DQ340604  Notanthera heterophylla (Ruiz & Pav.) G.Don.  Calvin and Wilson Cl03-03 (NA)  Chile  DQ333855  DQ340606  Nuytsia floribunda R.Br.  Nickrent 2747 (NA)  Australia  DQ788705  -  Nuytsia floribunda R.Br.  Nickrent 3080 (NA)  Australia  -  DQ788716  Nuytsia floribunda R.Br.  Calvin and Wilson AU01-22 (NA)  Australia  DQ333867  DQ340618  Oryctanthus occidentalis (L.) Eichler  Calvin and Wilson CR01- 11 (NA)  Costa Rica  DQ333862  DQ340613  Peraxilla tetrapetala (L.f.) Tiegh.  Calvin and Wilson NZ98- 03 (NA)  New Zealand  DQ333846  DQ340597  Phragmanthera regularis (Steud. ex Sprague) M.G.Gilbert  Calvin and Wilson Y88-01 (NA)  NA  DQ333830  DQ340579  Phthirusa pyrifolia (Kunth) Eichler  Calvin and Wilson CR01- 03 (NA)  Costa Rica  DQ333857  -  Phthirusa pyrifolia (Kunth) Eichler  Nickrent 2762 (NA)  NA  -  EU544504  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670A (USP)  Brazil  KU923005  KU923279  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670B (USP)  Brazil  KU923006  KU923280  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 004 (XAL)  Mexico  KU923007  KU923281  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 005 (XAL)  Mexico  KU923008  KU923282  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez (NA)  Mexico  KU923009  KU923283  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3785 (USP)  Brazil  KU923010  KU923284  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3786 (USP)  Brazil  KU923011  KU923285  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3671 (USP)  Brazil  KU923019  KU923293  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3672 (USP)  Brazil  KU923020  KU923294  Psittacanthus macrantherus Eichler  E. Ruiz-Sanchez 348 (XAL)  Mexico  KU923021  KU923295  Psittacanthus macrantherus Eichler  F. Rodriguez NA (XAL)  Mexico  KU923022  KU923296  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923029  KU923303  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923030  KU923304  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923031  KU923305  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 013 (XAL)  Mexico  KU923032  KU923306  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 014 (XAL)  Mexico  KU923033  KU923307  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 015 (XAL  Mexico  KU923034  KU923308  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 016 (XAL)  Mexico  KU923035  KU923309  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923036  KU923310  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923037  KU923311  Psittacanthus rhynchanthus (Benth.) Kuijt  P. Carrillo Reyes 5372 (XAL)  Guatemala  KU923038  KU923312  Psittacanthus rhynchanthus (Benth.) Kuijt  E. Ruiz-Sanchez 417 (XAL)  Mexico  KU923040  KU923314  Psittacanthus rhynchanthus (Benth.) Kuijt  T. Mejia-Saules 2048 (XAL)  Mexico  KU923041  KU923315  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3588 (USP)  Brazil  KU923042  KU923316  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3589 (USP)  Brazil  KU923043  KU923317  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3596 (USP)  Brazil  KU923044  KU923318  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923045  KU923319  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923046  KU923320  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923047  KU923321  Struthanthus orbicularis (Kunth) Blume  Calvin and Wilson CR01- 02 (NA)  Costa Rica  DQ333856  DQ340607  Tripodanthus acutifolius (Ruiz & Pav.) Tiegh.  Calvin and Wilson AR03-04 (NA)  Argentina  DQ333864  DQ340615  Tristerix aphyllus (Miers ex DC.) Barlow & Wiens  G. Amico 166 (BRCU)  Chile  DQ442966  DQ442919  Tristerix corymbosus (L.) Kuijt  Calvin and Wilson CL03-02 (NA)  Chile  DQ333854  DQ340605  Tupeia antarctica (G.Forst.) Cham. & Schltdl.  Calvin and Wilson NZ98- 02 (NA)  New Zealand  DQ333850  DQ340601  Species  Voucher  Country  ITS  trnL-F  Alepis flavida Tiegh.  Calvin and Wilson NZ98- 04 (NA)  New Zealand  DQ333847  DQ340598  Atkinsonia ligustrina F.Muell.  Calvin and Wilson AU00-01 (NA)  Australia  DQ333865  DQ788714  Cladocolea cupulata Kuijt  Calvin and Wilson MX03- 08 (NA)  Mexico  DQ333861  DQ340612  Cladocolea mcvaughii Kuijt  Calvin and Wilson MX03- 09 (NA)  Mexico  DQ333860  DQ340611  Desmaria mutabilis (Poepp. & Endl.) Tiegh. ex B.D.Jacks.  Calvin and Wilson CL03-07 (NA)  Chile  DQ333852  DQ340603  Gaiadendron punctatum (Ruiz & Pav.) G.Don.  Calvin and Wilson CR01- 08 (NA)  Costa Rica  DQ333866  DQ340617  Ligaria cuneifolia (Ruiz & Pav.) Tiegh.  Calvin and Wilson CL03-01 (NA)  Chile  DQ333853  DQ340604  Notanthera heterophylla (Ruiz & Pav.) G.Don.  Calvin and Wilson Cl03-03 (NA)  Chile  DQ333855  DQ340606  Nuytsia floribunda R.Br.  Nickrent 2747 (NA)  Australia  DQ788705  -  Nuytsia floribunda R.Br.  Nickrent 3080 (NA)  Australia  -  DQ788716  Nuytsia floribunda R.Br.  Calvin and Wilson AU01-22 (NA)  Australia  DQ333867  DQ340618  Oryctanthus occidentalis (L.) Eichler  Calvin and Wilson CR01- 11 (NA)  Costa Rica  DQ333862  DQ340613  Peraxilla tetrapetala (L.f.) Tiegh.  Calvin and Wilson NZ98- 03 (NA)  New Zealand  DQ333846  DQ340597  Phragmanthera regularis (Steud. ex Sprague) M.G.Gilbert  Calvin and Wilson Y88-01 (NA)  NA  DQ333830  DQ340579  Phthirusa pyrifolia (Kunth) Eichler  Calvin and Wilson CR01- 03 (NA)  Costa Rica  DQ333857  -  Phthirusa pyrifolia (Kunth) Eichler  Nickrent 2762 (NA)  NA  -  EU544504  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670A (USP)  Brazil  KU923005  KU923279  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670B (USP)  Brazil  KU923006  KU923280  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 004 (XAL)  Mexico  KU923007  KU923281  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 005 (XAL)  Mexico  KU923008  KU923282  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez (NA)  Mexico  KU923009  KU923283  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3785 (USP)  Brazil  KU923010  KU923284  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3786 (USP)  Brazil  KU923011  KU923285  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3671 (USP)  Brazil  KU923019  KU923293  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3672 (USP)  Brazil  KU923020  KU923294  Psittacanthus macrantherus Eichler  E. Ruiz-Sanchez 348 (XAL)  Mexico  KU923021  KU923295  Psittacanthus macrantherus Eichler  F. Rodriguez NA (XAL)  Mexico  KU923022  KU923296  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923029  KU923303  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923030  KU923304  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923031  KU923305  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 013 (XAL)  Mexico  KU923032  KU923306  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 014 (XAL)  Mexico  KU923033  KU923307  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 015 (XAL  Mexico  KU923034  KU923308  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 016 (XAL)  Mexico  KU923035  KU923309  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923036  KU923310  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923037  KU923311  Psittacanthus rhynchanthus (Benth.) Kuijt  P. Carrillo Reyes 5372 (XAL)  Guatemala  KU923038  KU923312  Psittacanthus rhynchanthus (Benth.) Kuijt  E. Ruiz-Sanchez 417 (XAL)  Mexico  KU923040  KU923314  Psittacanthus rhynchanthus (Benth.) Kuijt  T. Mejia-Saules 2048 (XAL)  Mexico  KU923041  KU923315  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3588 (USP)  Brazil  KU923042  KU923316  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3589 (USP)  Brazil  KU923043  KU923317  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3596 (USP)  Brazil  KU923044  KU923318  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923045  KU923319  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923046  KU923320  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923047  KU923321  Struthanthus orbicularis (Kunth) Blume  Calvin and Wilson CR01- 02 (NA)  Costa Rica  DQ333856  DQ340607  Tripodanthus acutifolius (Ruiz & Pav.) Tiegh.  Calvin and Wilson AR03-04 (NA)  Argentina  DQ333864  DQ340615  Tristerix aphyllus (Miers ex DC.) Barlow & Wiens  G. Amico 166 (BRCU)  Chile  DQ442966  DQ442919  Tristerix corymbosus (L.) Kuijt  Calvin and Wilson CL03-02 (NA)  Chile  DQ333854  DQ340605  Tupeia antarctica (G.Forst.) Cham. & Schltdl.  Calvin and Wilson NZ98- 02 (NA)  New Zealand  DQ333850  DQ340601  NA = Not available. View Large Appendix 2. Posterior parameter estimates for the best-supported scenario (scenario 1) considering the three Psittacanthus mayanus groups (dryCHI, humYUC, dryYUC). Parameter  mean  median  mode  q025  q050  q250  q750  q950  q975  N1  6.02e+02  6.31e+02  6.87e+02  1.92e+02  2.64e+02  4.92e+02  7.41e+02  8.37e+02  8.56e+02  N2b  3.99e+03  4.13e+03  4.84e+03  2.24e+03  2.57e+03  3.53e+03  4.57e+03  4.91e+03  4.95e+03  N3b  4.36e+03  4.49e+03  4.79e+03  3.00e+03  3.28e+03  4.10e+03  4.76e+03  4.95e+03  4.97e+03  N2a2  4.33e+03  4.52e+03  4.97e+03  2.61e+03  3.00e+03  4.05e+03  4.80e+03  4.96e+03  4.98e+03  N2b2  2.65e+03  2.65e+03  3.04e+03  1.10e+03  1.19e+03  1.89e+03  3.39e+03  4.16e+03  4.35e+03  Na  7.72e+02  8.09e+02  8.94e+02  4.50e+02  5.22e+02  7.20e+02  8.62e+02  8.93e+02  8.97e+02  t1  1.51e+03  1.24e+03  1.08e+03  1.02e+03  1.03e+03  1.11e+03  1.53e+03  3.02e+03  4.15e+03  t2  1.56e+04  1.57e+04  2.07e+04  1.03e+04  1.06e+04  1.29e+04  1.84e+04  2.05e+04  2.08e+04  t3  4.39e+04  4.19e+04  4.04e+04  4.01e+04  4.02e+04  4.07e+04  4.45e+04  5.51e+04  6.18e+04  μnuDNA  2.03e-07  1.92e-07  1.80e-07  1.08e-07  1.19e-07  1.57e-07  2.38e-07  3.24e-07  3.56e-07  μcpDNA  3.25e-07  3.07e-07  2.65e-07  1.06e-07  1.38e-07  2.34e-07  3.96e-07  5.71e-07  6.39e-07  Parameter  mean  median  mode  q025  q050  q250  q750  q950  q975  N1  6.02e+02  6.31e+02  6.87e+02  1.92e+02  2.64e+02  4.92e+02  7.41e+02  8.37e+02  8.56e+02  N2b  3.99e+03  4.13e+03  4.84e+03  2.24e+03  2.57e+03  3.53e+03  4.57e+03  4.91e+03  4.95e+03  N3b  4.36e+03  4.49e+03  4.79e+03  3.00e+03  3.28e+03  4.10e+03  4.76e+03  4.95e+03  4.97e+03  N2a2  4.33e+03  4.52e+03  4.97e+03  2.61e+03  3.00e+03  4.05e+03  4.80e+03  4.96e+03  4.98e+03  N2b2  2.65e+03  2.65e+03  3.04e+03  1.10e+03  1.19e+03  1.89e+03  3.39e+03  4.16e+03  4.35e+03  Na  7.72e+02  8.09e+02  8.94e+02  4.50e+02  5.22e+02  7.20e+02  8.62e+02  8.93e+02  8.97e+02  t1  1.51e+03  1.24e+03  1.08e+03  1.02e+03  1.03e+03  1.11e+03  1.53e+03  3.02e+03  4.15e+03  t2  1.56e+04  1.57e+04  2.07e+04  1.03e+04  1.06e+04  1.29e+04  1.84e+04  2.05e+04  2.08e+04  t3  4.39e+04  4.19e+04  4.04e+04  4.01e+04  4.02e+04  4.07e+04  4.45e+04  5.51e+04  6.18e+04  μnuDNA  2.03e-07  1.92e-07  1.80e-07  1.08e-07  1.19e-07  1.57e-07  2.38e-07  3.24e-07  3.56e-07  μcpDNA  3.25e-07  3.07e-07  2.65e-07  1.06e-07  1.38e-07  2.34e-07  3.96e-07  5.71e-07  6.39e-07  Estimates are based on 1% of simulated datasets closest to the observed values. Simulations and approximate Bayesian computation (ABC) analyses were performed considering both nrDNA and plastid DNA sequences. Parameters are N = effective population size for dryCHI (pop1, N1), humYUC (Pop2, N2), dryYUC (Pop3, N3), and the ancestral populations (Na), t = time since divergence (years), and μ = mutation rate for nuclear and plastid DNA, respectively. Group abbreviations are as follows: dryYUC = dry Yucatán, humYUC = humid Yucatán, dryCHI = dry Central Depression of Chiapas. View Large Appendix 3. Factor loadings from the principal components analysis of Psittacanthus mayanus in Mexico on temperature and precipitation variables from WorldClim. Variables  Factor loading    PC1  PC2  PC3  PC4  Annual mean temperature (BIO1)  0.809  0.558  0.051  0.084  Mean diurnal range (BIO2)  –0.621  0.234  0.675  0.320  Isothermality (BIO3)  –0.645  0.057  0.086  0.651  Temperature seasonality (BIO4)  0.427  0.476  0.467  –0.300  Max. temperature of warmest month (BIO5)  0.310  0.653  0.674  0.102  Min. temperature of coldest month (BIO6)  0.869  0.218  –0.391  –0.096  Temperature annual range (BIO7)  –0.506  0.271  0.791  0.150  Mean temperature of wettest quarter (BIO8)  0.668  0.689  –0.053  0.070  Mean temperature of driest quarter (BIO9)  0.826  0.434  0.040  –0.141  Mean temperature of warmest quarter (BIO10)  0.799  0.574  0.140  0.012  Mean temperature of coldest quarter (BIO11)  0.786  0.476  –0.099  0.187  Annual precipitation (BIO12)  0.640  –0.645  0.388  –0.044  Precipitation of wettest month (BIO13)  0.508  –0.611  0.524  –0.216  Precipitation of driest month (BIO14)  0.777  –0.414  0.155  0.358  Precipitation seasonality (BIO15)  –0.765  0.060  0.292  –0.393  Precipitation of wettest quarter (BIO16)  0.396  –0.589  0.606  –0.237  Precipitation of driest quarter (BIO17)  0.825  –0.397  0.067  0.289  Precipitation of warmest quarter (BIO18)  0.459  –0.649  –0.065  0.314  Precipitation of coldest quarter (BIO19)  0.799  –0.538  0.071  –0.011  Proportion of variance explained  45.66%  23.91%  15.08%  6.82%  ANOVA  F2,103 = 77.8 P < 0.0001  F2,103 = 39.9 P < 0.0001  F2,103 = 4.2 P = 0.017  F2,103 = 5.3 P = 0.006  Variables  Factor loading    PC1  PC2  PC3  PC4  Annual mean temperature (BIO1)  0.809  0.558  0.051  0.084  Mean diurnal range (BIO2)  –0.621  0.234  0.675  0.320  Isothermality (BIO3)  –0.645  0.057  0.086  0.651  Temperature seasonality (BIO4)  0.427  0.476  0.467  –0.300  Max. temperature of warmest month (BIO5)  0.310  0.653  0.674  0.102  Min. temperature of coldest month (BIO6)  0.869  0.218  –0.391  –0.096  Temperature annual range (BIO7)  –0.506  0.271  0.791  0.150  Mean temperature of wettest quarter (BIO8)  0.668  0.689  –0.053  0.070  Mean temperature of driest quarter (BIO9)  0.826  0.434  0.040  –0.141  Mean temperature of warmest quarter (BIO10)  0.799  0.574  0.140  0.012  Mean temperature of coldest quarter (BIO11)  0.786  0.476  –0.099  0.187  Annual precipitation (BIO12)  0.640  –0.645  0.388  –0.044  Precipitation of wettest month (BIO13)  0.508  –0.611  0.524  –0.216  Precipitation of driest month (BIO14)  0.777  –0.414  0.155  0.358  Precipitation seasonality (BIO15)  –0.765  0.060  0.292  –0.393  Precipitation of wettest quarter (BIO16)  0.396  –0.589  0.606  –0.237  Precipitation of driest quarter (BIO17)  0.825  –0.397  0.067  0.289  Precipitation of warmest quarter (BIO18)  0.459  –0.649  –0.065  0.314  Precipitation of coldest quarter (BIO19)  0.799  –0.538  0.071  –0.011  Proportion of variance explained  45.66%  23.91%  15.08%  6.82%  ANOVA  F2,103 = 77.8 P < 0.0001  F2,103 = 39.9 P < 0.0001  F2,103 = 4.2 P = 0.017  F2,103 = 5.3 P = 0.006  View Large © 2018 The Linnean Society of London, Botanical Journal of the Linnean Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Botanical Journal of the Linnean Society Oxford University Press

Lay mistletoes on the Yucatán Peninsula: post-glacial expansion and genetic differentiation of Psittacanthus mayanus (Loranthaceae)

Loading next page...
 
/lp/ou_press/lay-mistletoes-on-the-yucat-n-peninsula-post-glacial-expansion-and-XLPr478N5V
Publisher
The Linnean Society of London
Copyright
© 2018 The Linnean Society of London, Botanical Journal of the Linnean Society
ISSN
0024-4074
eISSN
1095-8339
D.O.I.
10.1093/botlinnean/box098
Publisher site
See Article on Publisher Site

Abstract

Abstract We evaluate phylogeographical patterns of Psittacanthus mayanus (Loranthaceae), a widely distributed mistletoe species on the Yucatán Peninsula and Chiapas, using nuclear and plastid markers. An analysis of phylogeographic population structure and demography and Bayesian inference methods were conducted on P. mayanus from 16 localities across the range of the species. To assess the historical processes and changes through the Pleistocene climatic cycles that may have affected the distribution and demographic history of the species, the current potential distribution of the species was modelled using ecological niche modelling (ENM) and projected onto the mid Holocene (MH), Last Glacial Maximum (LGM) and Last Inter Glacial (LIG). Psittacanthus mayanus exhibited higher haplotype and nucleotide diversity in the southern part of its distribution than in the northern part. Two divergent lineages were revealed in the phylogenetic and phylogeographical analyses of ribotypes and haplotypes. Populations from the northern portion of the Yucatán Peninsula probably experienced a recent population growth, whereas populations from the southern portion of the Peninsula and Chiapas are marked by historical demographic stability and isolation. Approximate Bayesian computation (ABC) analyses strongly supported a scenario of post-glacial population growth. ENM results indicate that the distribution of P. mayanus expanded and was connected into the southern portion of the Yucatán Peninsula during LGM conditions and colonized the northern portion of the Peninsula and fragmented during MH conditions. According to the ABC and ENM results, the genetic differentiation of P. mayanus may be linked to the effects of the glacial/interglacial cycles and environmental factors, in which populations of P. mayanus expanded during the LGM into the Petén province currently covered with semideciduous tropical rain forest followed by northward expansion, southwards contractions and post-glacial colonization of the northernmost portion of the Yucatán province, a region currently covered with seasonally dry tropical deciduous forest. Our results highlight the influence of Pleistocene events in shaping the geographical distribution of genetic variation in Neotropical lowland forest. The phylogeographic and environmental patterns in P. mayanus provide an opportunity to investigate further the evolution of Mexican lowland forest biodiversity. INTRODUCTION The Yucatán Peninsula has witnessed numerous marine transgressions over the last 65 Myr, submerging under warm tropical waters what is now the Peninsula (Vázquez-Domínguez & Arita, 2010). Today, the Yucatán Peninsula occupies almost half the portion of the peninsula platform, including the Mexican states of Yucatán, Quintana Roo, Campeche and adjacent portions of Chiapas and Tabasco, and the northern portions of Guatemala and Belize (c. 181 200 km2; Lundell, 1934; Espadas-Manrique, Durán & Argáez, 2003). Submerged for millions of years, the Peninsula is basically a large limestone slab that emerged slowly from south to north. The ecological features and biogeography of the Yucatán Peninsula show the climatic effects of its 65-Myr history, with a rain gradient from humid in the south/south-east to dry in the north/north-west (Espadas-Manrique et al., 2003; Orellana, Islebe & Espadas, 2003; Torrescano-Valle & Islebe, 2015). The vegetation also follows the south-east to north-west rain gradient, from tropical rainforests in the Petén region to tropical scrubland in the extreme north-west part of the Yucatán peninsula, with extensive areas covered with deciduous or semideciduous tropical forests between these two extremes, and a decrease in species diversity from the humid communities in the south to the northern counterparts (Simpson’s peninsula effect; Simpson, 1964; Ramírez-Barahona et al., 2009; Vázquez-Domínguez & Arita, 2010). This distinctiveness has prompted most biogeographers to consider the Yucatán Peninsula a region with two provinces, the Petén province in the south and the dry Yucatán province in the north (e.g. Espadas-Manrique et al., 2003; Morrone, 2005; Ramírez-Barahona et al., 2009; Vázquez-Domínguez & Arita, 2010). Palaeogeographic data indicate that the Yucatán Peninsula was emergent until the Triassic–Jurassic, with several marine transgressions having occurred to the Pliocene–Miocene, and the present shape of the peninsula was completed near the end of the Pliocene and continued its formation into the Quaternary (López Ramos, 1975; Padilla y Sánchez, 2007). During most of the Pleistocene, the climate of the Yucatán Peninsula was dry and 6ºC cooler than current temperatures. Savanna and Juniperus-matorral covered most of the region until the early Holocene (Islebe et al., 1996; Orellana et al., 2003), in which the vegetation of the Peninsula changed to a more humid and tropical-type forest (Brosimum Sw.) established by 6800 years before present (BP) followed by a dry period between 6000 and 5000 BP (Leyden, 1984; Metcalfe et al., 2000). Wetter conditions were then re-established until c. 3500 BP, to be followed by drying and a more open canopy forest and more grasses until c. 1500 BP (Metcalfe et al., 2000). This dry period corresponds to the time of the ‘collapse’ of the lowland Maya civilization (Metcalfe et al., 2000; Haug et al., 2003; Vázquez-Domínguez & Arita, 2010). The biogeography of the region is not only influenced by the effects of the formation and emergence of the Yucatán Peninsula and the effects of Pleistocene climate changes but also by those related to the geology of and tectonic activity on the Maya block, between the biogeographic breaks of the Isthmus of Tehuantepec and the Motagua–Polochic–Jocotán fault system (Gutiérrez-García & Vázquez-Domínguez, 2012). In the Maya block, the most conspicuous and dynamic mountain ranges include the modern Chiapas volcanic arc (Altos de Chiapas) and the Chiapas massif (Sierra Madre de Chiapas) separated by the Central Depression of Chiapas and the Maya mountains of Belize (Gutiérrez-García & Vázquez-Domínguez, 2012). Although these geological and topographical features and the isolation of these mountain ranges certainly influenced the evolutionary history of Central American lowland tree species (Cavers, Navarro & Lowe, 2003; Poelchau & Hamrick, 2012, 2013a) and widespread Mesoamerican vertebrates (e.g. Guevara-Chumacero et al., 2010; González, Ornelas & Gutiérrez-Rodríguez, 2011; Smith et al., 2011; Gutiérrez-García & Vázquez-Domínguez, 2012; Rodríguez-Gómez & Ornelas, 2014; Jiménez & Ornelas, 2016; Williford et al., 2016), few studies have included phylogeographies of terrestrial species restricted to the Maya block (in the Yucatán Peninsula). Gutiérrez-García & Vázquez-Domínguez (2012) described patterns of historical divergence and population genetic differentiation of mitochondrial DNA in the rodent Ototylomys phyllotis distributed from the Isthmus of Tehuantepec (Mexico) to Costa Rica. They found that the Yucatán haplotypes formed the most derived lineage and hypothesized that O. phyllotis arrived to Mexico from the south (Central America). The colonization of the Yucatán Peninsula occurred c. 0.8 million years ago (Mya), from individuals that persisted in Pleistocene refugia on the Maya mountains of Belize that expanded into the Yucatán Peninsula. More recently, Oyama et al. (2016) assessed genetic differentiation of Guaiacum sanctum L. (Zygophyllaceae) trees on the Yucatán Peninsula, and found evidence of fragmentation and little genetic differentiation between two groups of populations: northern populations of Yucatán on porous karst terrain soils and southern populations from Campeche in seasonally flooded terrains. In this study, the demographic history, population structure and genetic differentiation of Psittacanthus mayanus Standl. & Steyerm. (Loranthaceae) populations are investigated using nuclear DNA (nrDNA) and plastid DNA sequences. We also used ecological niche modelling, isolation-with-migration and Bayesian inference methods to provide insight into the structuring of genetic variation of these populations. Specifically, we asked the following two questions. (1) What geographical features and climatic oscillations in the Pleistocene affected the genetic diversity and population genetic structure of P. mayanus? (2) Based on the observed genetic diversity, did populations expand across the flat Yucatán Peninsula and diverge according to the rain gradient? The distribution of P. mayanus populations in different habitats across the peninsula (Petén and Yucatán provinces) and Chiapas facilitates the testing of hypotheses on the relative importance of divergence and gene flow in the evolutionary process and colonization history of this species, in which both physical and ecological barriers might limit gene flow. Given that Quaternary climate fluctuations, particularly between now and the Last Glacial Maximum (LGM), may have driven distributional shifts in tree species of the Yucatán Peninsula, we predict dynamic shifts in the reconstructed distribution of P. mayanus along the Yucatán Peninsula as seasonally dry tropical deciduous forest contracted and expanded many times during the Pleistocene (Metcalfe et al., 2000). In this scenario, we assume that (1) the populations in Chiapas and those in the Petén and Yucatán provinces experienced cycling periods of isolation during LGM contraction of the seasonally dry tropical deciduous forest followed by secondary contact and gene flow during the warmer interglacial cycles; (2) populations of P. mayanus in the seasonally dry tropical deciduous forest of the Yucatán province contracted southwards into refugia during the LGM in the semideciduous tropical rain forest of the Petén province followed by northward expansion and recolonization of the seasonally dry tropical deciduous forest of the Yucatán province during warmer periods. The absence of shifts in suitable habitat between the LGM and the present day for several tree species of the seasonally dry tropical deciduous forest (Cavers et al., 2003; Poelchau & Hamrick, 2013a, b) would support an alternative hypothesis of long-term persistence and stability in climatically mediated distributions of P. mayanus with floristic affinities to seasonally dry tropical deciduous forest (Kuijt, 2009). This hypothesis assumes that the Petén and Yucatán populations of P. mayanus experienced genetic divergence in isolation. MATERIAL AND METHODS Study Organism Psittacanthus mayanus commonly grows in seasonally dry deciduous and semideciduous forests to seasonally dry deciduous forest with columnar cacti. The species frequently parasitizes leguminous trees (Fabaceae), mainly Lysiloma Benth., Lonchocarpus Kunth, Acacia Mill., Pithecellobium Mart., Inga Mill. and Haematoxylum Gronov., and less frequently Coccoloba L. (Polygonaceae), Conocarpus L. (Combretaceae) Crescentia L. (Bignoniaceae) and Bursera Jacq. ex L. (Burseraceae) trees (Kuijt, 2009). It is distributed in Mexico, Guatemala and Belize, mostly on the Yucatán Peninsula and Chiapas, extending somewhat to the north and south from almost sea level up to 500 m a.s.l. (Kuijt, 2009). Psittacanthus mayanus is easily recognizable from other Psittacanthus Mart. in the region by the presence of prophylls (two minute foliar organs, transversely placed in a leaf axil), also detected in P. elegans Kuijt and conspicuous in P. breedlovei Kuijt and P. sonorae (Watson) Kuijt, four or five cotyledons, non-foliaceous lowest triad bracts and straight, thinner flowers. In its southern range, however, P. mayanus is sometimes difficult to separate from other Psittacanthus species presumably occurring in the tropical dry forest, namely P. schiedeanus (Cham. & Schltdl.) G.Don, P. calyculatus G.Don, P. breedlovei Kuijt and P. minor Kuijt (Kuijt, 2009). However, P. breedlovei has extremely narrow, linear-lanceolate leaves, <1 cm wide (P. mayanus has leaves >2 cm wide, ovate, ovate-lanceolate) and P. calyculatus is a robust plant with denser inflorescences and buds to 6 cm long (P. mayanus is a slender plant with few flowers and buds to 3 cm long); P. minor is endemic to Nicaragua (Kuijt, 2009). The few-flowered inflorescences of P. mayanus produce red-orange flowers mainly visited by Campylopterus curvipennis pampa and Amazilia hummingbirds and its ripe, purplish-black fruits are consumed by several bird species, mainly Patagioenas cayennensis, Myiozetetes similis, Pitangus sulphuratus, Megarhynchus pitangua and Piranga roseogularis (Howell, 1972; A. Vargas-Varguéz, unpubl. data). Population Sampling Leaf tissue samples were collected from 87 individuals in 16 localities across most of the species range of P. mayanus (Table 1). Samples from Quintana Roo, Belize and Guatemala were not available for the present study. Target sampling localities were chosen based on localities of occurrence data acquired from Kuijt (2009). Only one plant of P. mayanus was sampled per individual host tree, with care taken to collect only samples from adult flowering individuals and individuals present in the same location were treated as a population, ranging from one to 14 individuals per population. Most populations collected have accompanying vouchers that were deposited at the XAL herbarium of the Instituto de Ecología, A.C. (INECOL; for voucher information, see Table 1). Table 1. Geographical information and sample sizes of the 16 Psittacanthus mayanus localities sampled in this study and voucher information. Pop. no.  Location  G2  G3  GP  n  Elevation (m)  Latitude (N)  Longitude (W)  Voucher (XAL)  1  Yucatán, Dzibilchaltún  YUC  YUC  dryYUC  3  13  21° 05′ 50′′  89° 34′ 55′′  E. Gándara 3131–3133  2  Yucatán, Hunucmá, Sisal  YUC  YUC  dryYUC  14  10  21° 02′ 58′′  89° 54′ 38′′  –  3  Yucatán, Celestún  YUC  YUC  dryYUC  9  8  21° 02′ 11′′  89° 48′ 07′′  E. Gándara 3137–3145  4  Yucatán, Mérida, Cuxtal  YUC  YUC  dryYUC  13  11  20° 54′ 37′′  89° 37′ 15′′  –  5  Yucatán, Maxcanú, Límite Estatal  YUC  YUC  dryYUC  8  11  20° 36′ 31′′  89° 57′ 54′′  Y. Licona-Vera 149–156  6  Campeche, Champotón, Seybaplaya  YUC  PET  dryYUC  10  18  19° 38′ 35′′  90° 39′ 56′′  E. Gándara 3121–3130  7  Campeche, Escárcega, Sabancuy  YUC  PET  humYUC  10  7  18° 54′ 19′′  91° 03′ 31′′  E. Gándara 3111–3120  8  Chiapas, Bochil, Km 205 Pto. Caté  CHI  CHI  humYUC  3  1158  16° 59′ 19′′  92° 52′ 21′′  Y. Licona-Vera 023-025  9  Chiapas, Bochil, Sn Cristobalito  CHI  CHI  humYUC  1  1134  16° 59′ 23′′  92° 52′ 59′′  Y. Licona-Vera 026  10  Chiapas, Sima Las Cotorras  CHI  CHI  dryCHI  6  814  16° 49′ 52′′  93° 27′ 30′′  Y. Licona-Vera 007-012  11  Chiapas, Ocozocuautla, El Aguacero  CHI  CHI  dryCHI  1  602  16° 45′ 50′′  93° 31′ 50′′  Y. Licona-Vera 001  12  Chiapas, Ocozocuautla  CHI  CHI  dryCHI  1  1078  16º 45′ 31′′  93º 19′ 34′′  E. Gándara 3095  13  Chiapas, Arriaga, Km 89–92  CHI  CHI  dryCHI  2  759  16° 44′ 03′′  93° 24′ 47′′  E. Ruiz 259  14  Chiapas, Km 9 Arriaga-Tuxtla  CHI  CHI  dryCHI  1  705  16º 29′ 23′′  93º 50′ 55′′  E. Gándara 3088  15  Chiapas, Arriaga, La Aurora  CHI  CHI  dryCHI  1  699  16° 28′ 18′′  93° 49′ 55′′  E. Gándara 3090  16  Chiapas, Arriaga, Km 28.3  CHI  CHI  dryCHI  6  667  16° 24′ 24′′  93° 48′ 22′′  E. Ruiz 342  Pop. no.  Location  G2  G3  GP  n  Elevation (m)  Latitude (N)  Longitude (W)  Voucher (XAL)  1  Yucatán, Dzibilchaltún  YUC  YUC  dryYUC  3  13  21° 05′ 50′′  89° 34′ 55′′  E. Gándara 3131–3133  2  Yucatán, Hunucmá, Sisal  YUC  YUC  dryYUC  14  10  21° 02′ 58′′  89° 54′ 38′′  –  3  Yucatán, Celestún  YUC  YUC  dryYUC  9  8  21° 02′ 11′′  89° 48′ 07′′  E. Gándara 3137–3145  4  Yucatán, Mérida, Cuxtal  YUC  YUC  dryYUC  13  11  20° 54′ 37′′  89° 37′ 15′′  –  5  Yucatán, Maxcanú, Límite Estatal  YUC  YUC  dryYUC  8  11  20° 36′ 31′′  89° 57′ 54′′  Y. Licona-Vera 149–156  6  Campeche, Champotón, Seybaplaya  YUC  PET  dryYUC  10  18  19° 38′ 35′′  90° 39′ 56′′  E. Gándara 3121–3130  7  Campeche, Escárcega, Sabancuy  YUC  PET  humYUC  10  7  18° 54′ 19′′  91° 03′ 31′′  E. Gándara 3111–3120  8  Chiapas, Bochil, Km 205 Pto. Caté  CHI  CHI  humYUC  3  1158  16° 59′ 19′′  92° 52′ 21′′  Y. Licona-Vera 023-025  9  Chiapas, Bochil, Sn Cristobalito  CHI  CHI  humYUC  1  1134  16° 59′ 23′′  92° 52′ 59′′  Y. Licona-Vera 026  10  Chiapas, Sima Las Cotorras  CHI  CHI  dryCHI  6  814  16° 49′ 52′′  93° 27′ 30′′  Y. Licona-Vera 007-012  11  Chiapas, Ocozocuautla, El Aguacero  CHI  CHI  dryCHI  1  602  16° 45′ 50′′  93° 31′ 50′′  Y. Licona-Vera 001  12  Chiapas, Ocozocuautla  CHI  CHI  dryCHI  1  1078  16º 45′ 31′′  93º 19′ 34′′  E. Gándara 3095  13  Chiapas, Arriaga, Km 89–92  CHI  CHI  dryCHI  2  759  16° 44′ 03′′  93° 24′ 47′′  E. Ruiz 259  14  Chiapas, Km 9 Arriaga-Tuxtla  CHI  CHI  dryCHI  1  705  16º 29′ 23′′  93º 50′ 55′′  E. Gándara 3088  15  Chiapas, Arriaga, La Aurora  CHI  CHI  dryCHI  1  699  16° 28′ 18′′  93° 49′ 55′′  E. Gándara 3090  16  Chiapas, Arriaga, Km 28.3  CHI  CHI  dryCHI  6  667  16° 24′ 24′′  93° 48′ 22′′  E. Ruiz 342  IDs reported below refer to accession numbers in the Instituto de Ecología, AC (XAL) herbarium. G2 = two groups by geography (YUC = Yucatán Peninsula, CHI = Chiapas), G3 = three groups by geography (YUC = Yucatan Province, PET = Petén Province, CHI = Chiapas), GP = three groups by precipitation gradient (dryYUC= dry Yucatán, humYUC = humid Yucatán, dryCHI = dry Central Depression of Chiapas). View Large Leaf tissue samples of P. mayanus were preserved in silica gel desiccant until DNA extractions were performed. We additionally included multiple DNA sequences of Psittacanthus schiedeanus (accession numbers ITS: KU922961–KU923004, KU923036–KU923037; trnL-F: KU923270–KU923278, KU923310–KU923311), P. calyculatus (accession numbers ITS: KX640930–KX640946, KX241619–KX241734, trnL-F: KX640947–KX640960, KX241735–KX241823) and 96 sequences of ten other Psittacanthus spp. (29 samples) and 18 other representatives of Loranthaceae (19 samples) downloaded from the GenBank (accession numbers and voucher information in Appendix 1) to be used as outgroups according to Wilson & Calvin (2006), Amico, Vidal-Russell & Nickrent (2007), Vidal-Russell & Nickrent (2007, 2008a, b), Díaz Infante et al. (2016), Ornelas et al. (2016), and Pérez-Crespo et al. (2017). DNA sequencing The internal transcribed spacer (ITS) region of the nuclear ribosomal DNA (nrDNA) and the trnL-F (cpDNA) intergenic spacer from the plastid genome were amplified by PCR and sequenced for 87 individuals. These markers have shown the appropriate level of variation for mistletoe studies at the species and population levels (e.g. Amico et al., 2012; Ornelas et al., 2016; Pérez-Crespo et al., 2017). As the plastid genome is maternally inherited in most angiosperms, plastid DNA sequences may provide a picture of evolutionary haplotype relationships of seed-mediated gene flow compared to the nrDNA sequences reflecting a combination of pollen and seed movement (Ornelas et al., 2016). Total genomic DNA was extracted from silica-dried material using a modified 2× cetyl trimethyl ammonium bromide (CTAB) protocol (Doyle & Doyle, 1987) or the DNeasy Plant Mini kit (Qiagen, Valencia, CA, USA) using the manufacturer protocol. Amplification of the nrDNA ITS region was conducted with the primers ITS-F2-Psitta (5′-TCGCAGTATGCTCCGTATTG-3′) and ITS-R2-Psitta (5′-TCGTAACAAGGTTTCCGTAGG-3′) designed for this project (Ornelas et al., 2016). For the trnL-F region we used the universal primers ‘e’ and ‘f’ (Taberlet et al., 1991). The 25-µL ITS PCR mix contained 5 µL 5× buffer (Promega, Madison, WI, USA), 2.0 μL MgCl2 (25 mM), 2 µL dNTPs mix (8 mM), 0.82 µL each primer (10 µM), 0.27 µL Taq polymerase (5U/µL) (Promega), 2 µL dimethyl sulphoxide (Sigma), 1.5–3.0 µL template DNA and dH2O to volume. The 25–µL trnL-F PCR mix contained 3.6 µL 5× buffer (Promega, Madison, WI, USA), 3 μL MgCl2 (25 mM), 1.8 µL dNTPs mix (0.57 mM), 0.40 µL each primer (0.16 µM), 0.20 µL Taq polymerase (5U/µL) (Promega), 1.5–5 µL template DNA, and dH2O added to volume. PCRs of ITS consisted of an initial denaturation at 94ºC for 5 min, followed by five cycles of 94ºC for 30 s, 50ºC for 30 s, 72ºC for 1 min, 35 cycles of 94ºC for 30 s, 48–52ºC for 30 s, 72ºC for 1 min and a final step of 72ºC for 5 min. Plastid (trnL-F) amplifications used the following profile: initial denaturation at 95ºC for 5 min, 40 cycles of 95ºC for 10 s, 55–57ºC for 1 min, 72ºC for 20 s and a final step of 72ºC for 4 min. PCR products were purified with the QIAquick kit (Qiagen) and sequenced in both directions to check the validity of the sequence data using the BigDye Terminator Cycle Sequencing kit (Applied Biosystems, Foster City, CA, USA). The products were analysed on a 310 automated DNA sequencer (Applied Biosystems) at the University of Washington High Throughput Genomics Unit, Seattle, Washington. Finally, assembled sequences of the two regions examined were edited and manually aligned with BIO-EDIT ver. 7.2.5 (Hall, 1999) and SE-AL ver. 2.0a111 (http://tree.bio.ed.ac.uk/software/seal). All unique sequences were submitted to GenBank (Accession numbers ITS: KY446904–KY446922, trnL-F: KY446923–KY446942). Divergence time estimation and genealogical relationships We estimated divergence time under a Bayesian approach using the combined ITS and plastid trnL-F sequence data set. The ingroup comprised all sequences of P. mayanus; sequences of other representatives of Loranthaceae including other Psittacanthus spp. downloaded from GenBank from Wilson & Calvin (2006), Amico, Vidal-Russell & Nickrent (2007), Vidal-Russell & Nickrent (2007, 2008a, b), Díaz Infante et al. (2016), Ornelas et al. (2016) and Pérez-Crespo et al. (2017) were used as outgroups. Phylogenetic reconstruction and divergence time estimation was performed with BEAST ver. 1.6.1 (Drummond & Rambaut, 2007) using the uncorrelated lognormal relaxed molecular clock and the closest nucleotide substitution model under the Bayesian information criterion (BIC), GTR+G+I for the combined data set, suggested by jMODELTEST ver. 0.1.1 (Posada, 2008). The tree prior model was set using a coalescent approach assuming constant population size. In this analysis, we constrained Nuytsia floribunda (Labill.) G.Don as sister to the aerial parasites based on Vidal-Russell & Nickrent (2008b) to calibrate the root node. The divergence time between Nuytsia and the aerial parasite clade was used as secondary calibration, approximating a median age of 28 Myr (normal distribution, mean 28, SD 3.7, range 34–21 Myr). The geometric mean of 5.798 × 10–9 substitutions per neutral site per year (s/s/y) was used to calibrate the tree based on the mean mutation rates of 4.130 × 10–9 s/s/y for ITS of herbaceous annual/perennial plants (Kay, Whittall & Hodges, 2006) and 8.240 × 10–9 s/s/y for trnL-F estimated for annual or perennial herbs (Richardson et al., 2001). For divergence time estimation, BEAST was run two times for 30 million generations, sampling every 3000 steps. We combined log and trees files from each independent run using LOGCOMBINER ver. 1.8.0 (Drummond & Rambaut, 2007), then viewed the combined log file in TRACER ver. 1.6 (http://tree.bio.ed.ac.uk/software/tracer/) to ensure that ESSs 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 10% trees were discarded as burn-in, and the remaining trees were annotated and summarized as a maximum clade credibility tree with mean divergence times and 95% highest posterior density (HPD) intervals of age estimates using TREEANNOTATOR ver. 1.8.0 (Drummond & Rambaut, 2007) and visualized in FIGTREE ver. 1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/). To infer genealogical relationships among haplotypes, statistical parsimony networks for ITS and trnL-F data sets were constructed using TCS ver. 1.2.1 (Clement, Posada & Crandall, 2000), with polymorphic sites coded as ambiguous (ITS) or gaps treated as single evolutionary events (trnL-F) and a connection limit set to 95%. Loops were resolved following the criteria given by Pfenninger & Posada (2002). Genetic diversity and population structure Molecular diversity indices including haplotype diversity (h) and nucleotide diversity (π), pairwise comparisons of FST values between genetic groups with 1600 permutations, and analyses of molecular variance (AMOVAs; Excoffier, Smouse & Quattro, 1992) were performed in ARLEQUIN ver. 3.01 (Excoffier, Laval & Schneider, 2005). Populations with fewer than three samples were combined with closest populations for this analysis (Table 1). To explore whether populations are structured along the Yucatán Peninsula and Chiapas, three AMOVAs were performed based on pairwise differences using ARLEQUIN with populations treated as (1) a single group to determine the amount of variation partitioned among and within populations, and grouped by geography into two (2) Yucatán Peninsula and Chiapas or three lineages (3) corresponding to the southern (PET = Petén province) or northern (YUC = Yucatán province) distribution of the Yucatán Peninsula and Chiapas (CHI) or (4) grouped by precipitation into three lineages corresponding to dry (dryYUC or dryCHI) and humid (humYUC) locations (Table 1). The AMOVAs were performed using the Kimura 2P model for ITS and the Jukes & Cantor model for trnL-F and 16 000 permutations to determine the significance of each AMOVA. The program PERMUT ver. 1.0 (Pons & Petit, 1996) was used to calculate for each marker within- population-diversity (hS), total diversity (hT), geographical total haplotype diversity (vT), geographical average haplotype diversity (vS), level of population differentiation at the species level (GST) and an estimate of population subdivisions for phylogenetically ordered alleles (NST). We further tested for phylogeographic structure at the species range between NST and GST parameters using 10 000 permutations and the U-statistic. An NST value significantly higher than the GST value means that two alleles sampled from the same region are phylogenetically more closely related than two alleles sampled from different regions, providing evidence for the presence of a significant phylogeographic signal in the data (Pons & Petit, 1996). Demographic history The demographic patterns of P. mayanus populations were assessed using different approaches. Fu’s Fs (Fu, 1997) and Tajima’s D (Tajima, 1989) tests were computed for each marker with 25 000 permutations to judge significance to test whether populations departed from a mutation-drift equilibrium model. A significantly negative value is often taken as an indicator of population expansion. In addition, mismatch distributions (Harpending, 1994) were simulated under the sudden expansion model of Schneider & Excoffier (1999) with 16 000 bootstrap replicates. The validity of the sudden expansion assumption was determined using the sum of squares differences (SSD) and Harpending’s raggedness index (Hri; Harpending, 1994), both of which are higher in stable, non-expanding populations (Rogers & Harpending, 1992). The SSD between the observed and the expected mismatch was calculated, and a significant SSD value means the observed value does not fit sudden expansion model. All tests were performed using ARLEQUIN and DnaSP (Librado & Rozas, 2009). The demographic changes in effective population size (Ne) through time of each P. mayanus group (Petén and Yucatán) were also inferred with Bayesian skyline plots (Drummond et al., 2005) using BEAST. We chose a HKY substitution model with empirical base frequencies, a strict clock model and a piecewise-linear coalescent Bayesian skyline tree prior with five starting groups. Two independent runs of 20 million generations each were run, with trees and parameters sampled every 1000-iterations and a burn-in of 10%. Results of each run were visualized using TRACER to ensure that stationarity and convergence had been reached, and that the ESS were higher than 200. We used the software IMa (Hey & Nielsen, 2007) on genetic groups in P. mayanus (Yucatán Peninsula = 63 samples, Chiapas = 24 samples) to estimate the effective population size of the ancestral (θa) and the two descendant populations (θ1 and θ2), effective migration rates in both directions (m1-to-2 and m2-to-1) and time since divergence (t) at which the ancestral population gave rise to the descendant populations. The isolation-with-migration (IM) model (Hey & Nielsen, 2004) is appropriate to estimate parameters for two descendant populations that have diverged recently from the ancestral population and that may be sharing haplotypes as a result of gene interchange. We began with multiple short runs of 10 000 steps (following 100 000 iterations as burn-in) to assess mixing and to fine-tune the parameter space. Once the six parameters were optimized (several replicates converged on the same answer), a long run of MCMC simulations was performed for one million generations and 30 million steps, under the HKY model of sequence evolution. Two additional independent runs were performed with different seed numbers to guarantee convergence of samples (Hey & Nielsen, 2004). We considered that the analyses had converged upon a stationary distribution if the independent runs had similar posterior distributions and the ESS for each parameter was at least 50. We report the mean parameter estimates of three runs and the 95% highest posterior densities (HPDs) intervals of each parameter. To convert raw parameter estimates into demographic quantities, we used the geometric mean of 4.58 × 10–5 of the median mutation rates of both loci reported for P. schiedeanus (Ornelas et al., 2016) and P. calyculatus (Pérez-Crespo et al., 2017). These mutation rates were converted to per-locus mutation rate (μ) by multiplying by the fragment length in base pairs for conversion, as required by IMa (Hey & Nielsen, 2007), to compute the divergence time (t) by using T = t/μ. We used an 11-year generation time (G) based on the observation that the age at maturity (seed production) begins c. 2 years after seed germination and an assumed annual survivorship of 0.9 based on seedling survivorship in the sister species P. schiedeanus (López de Buen & Ornelas, 2002; Ornelas et al., 2016) to convert the effective populations size estimates. The approximate average generation time (T) is calculated according to T = a + [s ⁄ (1–s)] (Lande, Engen & Sæther, 2003), where a is the time to maturity and s is the adult annual survival rate. Based on this, the estimate for T was 11 years. To convert the IMa time since divergence parameter to years, t, we divided the time parameter (B) by the mutation rate per year (U) converted to per locus rate by multiplying by the fragment length in base pairs, and calculated for the geometric mean rate described before. To calculate effective population size (N), we used N = θ/(4Gμ), where the generation time (G) is 11 years. To estimate the population migration rate, we used 2Nm = θ * m/2. Population expansion and approximate Bayesian computation (ABC) To assess glacial and post-glacial changes in population size within the Yucatán Peninsula that could explain population sizes and genetic diversity, we applied approximate Bayesian computation (ABC; Beaumont, Zhang & Balding, 2002) implemented in DIYABC (Cornuet et al., 2014). Two different models for population size changes were assumed, with population growth or population decline. We assumed either 6000 or 21 000 yr bp for the beginning of the population size change, which correspond to the MH or the LGM, respectively. For this analysis, we grouped all populations into three groups of samples (dryCHI, humYUC, dryYUC): scenario 1, population expansion from 21 000–10 000 yr bp after the split between the dryCHIS and the YUC groups; scenario 2, population expansion from 1000–8000 yr bp after the split between the dryCHIS and the YUC groups. As little is known about the effective population size of P. mayanus in the Yucatán Peninsula, we defined a Ne prior with a uniform distribution for the YUC population based on previous estimates of closely related Psittacanthus spp. (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Furthermore, we assumed 11 years as the generation time and mutation rates based on estimates for other Psittacanthus spp. (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Comparison of scenarios was implemented in the DIYABC, with two million coalescent-based simulated datasets under each evolutionary scenario. The posterior probability of scenarios was assessed using a logistic regression on the 1% of simulated datasets closest to the observed data (Fontaine et al., 2013). For the best-supported scenario, we performed a model checking procedure, as implemented in DIYABC, and assessed the confidence in scenario choice by simulating 500 pseudo-observed datasets (PODs) under each scenario to estimate Type I and Type II error rates (Robert et al., 2011). Finally, point estimates for demographic and temporal parameters were obtained by local linear regression on the 1% of simulations closest to the observed dataset for the best-supported scenario (Cornuet et al., 2008, 2014). Palaeodistribution modelling and environmental variation We used ecological niche modelling (ENM; Elith et al., 2011) to predict the distribution of P. mayanus at the MH (c. 6000 years ago), LGM (c. 20 000 years ago) and at the LIG (c. 120 000–140 000 years ago) and to examine whether populations of P. mayanus (1) persisted in the southern portion of the Yucatán Peninsula (in the wetter Petén province) in the semideciduous tropical rain forest during the cold glacial period and expanded northwards to the seasonally dry tropical deciduous forest after the LGM or (2) contracted into southern refugia during the LGM and expanded northwards only after this period. Mistletoe occurrence data were assembled mainly from herbarium records from the Global Biodiversity Information Facility (GBIF; https://www.gbif.org/species/4003770), and supplemented with records from Kuijt (2009) and our geo-referenced records from field collection. After careful verification of data for every location, excluding duplicate occurrence records or in close proximity to each other (c. 1 km2) to reduce the effects of spatial autocorrelation, we restricted the data set to 70 unique presence records for the analysis. Distributional records were input and analysed to infer an ENM with the maximum entropy algorithm in MAXENT ver. 3.2.2 (Phillips, Anderson & Schapire, 2006) and using ARCVIEW ver. 3.2 (ESRI, Redlands, CA, USA) to extract the data. Present-day and Last Interglacial (LIG) temperature and precipitation data (BIO 1–19 variables) were drawn as climate layers from the WorldClim database at a spatial resolution of c. 1 km2 in each raster (Hijmans et al., 2005; http://www.worldclim.org). We also obtained data for past climate layers for two MH (2.5 arc-minutes, pixel size c. 21.62 km2) and LGM (2.5 arc-minutes) past conditions based on CCSM4 and MIROC-ESM global models (Hijmans et al., 2005; Braconnot et al., 2007). We first constructed a model to the present using all climatic variables with 70% of the locality records as training data and the other 30% as testing data in addition to response curves, jackknife tests, logistic output format and random seed parameters. We ran 1000 iterations with no extrapolation in order to avoid artificial extrapolations from extreme values of the ecological variables; as such parameters are biased towards the environmental envelope of background points and occurrence data (Phillips et al., 2006). All other parameters in MAXENT were maintained at default settings. The logistic format was used to obtain the values of habitat suitability and the variable importance was determined comparing percentage contribution values and jackknife plots. Then, correlation coefficients between present variables were calculated using PAST ver. 2.12 (Hammer, 2011) with the objective to identify and remove the variables that were highly correlated (correlation values ≥ 0.8) and less explanatory, having as a rule of decision, the relative contributions of each of them in the MAXENT model. After removing the highly correlated variables, six variables were used in the final analysis for the mistletoe (BIO1 = annual mean temperature, BIO3 = isothermality, BIO4 = temperature seasonality, BIO6 = minimum temperature of coldest month, BIO7 = temperature annual range, BIO18 = precipitation of warmest quarter). The models were calibrated according the accessibility area of P. mayanus; this was estimated using a geographical clipping based on biogeographic provinces (Morrone, 2005, 2014) and elevational range limits (Kuijt, 2009), which represent the species historical barriers to dispersal (Kuijt, 2009). The resulting species distribution under current climate conditions was then projected onto the LIG and onto past climate reconstructions for the MH and LGM under two general atmospheric circulation models (CCSM4 and MIROC-ESM). Final models were constructed with ten cross-validation replicates without clamping and extrapolation and considering the average output grids as the final predictive models. The area under the curve (AUC) values of receiver operating characteristic (ROC) curves was used to evaluate the statistical performance of ENMs, a single measure of model performance independently of thresholds (Phillips et al., 2006), where 1 is the maximum prediction and 0.5 a random prediction. In addition, we also used the true skill statistic (TSS), which is a threshold-dependent measure of model performance, to evaluate the accuracy of predictive maps generated by presence-only data (Allouche, Tsoar & Kadmon, 2006; Liu et al., 2005), where TSS values ranging between 0.4–0.8 are considered useful (Landis & Koch, 1977; Fielding & Bell, 1997). For each replicate, the TSS was calculated using the ten-percentile training presence logistic threshold, and then TSS values were averaged among replicates. Measurements of temperature and precipitation (BIO 1–19 variables) were extracted for each of the sampling locations and for those with mistletoe occurrence data from WorldClim (Hijmans et al., 2005) as explained before. We performed principal components analysis (PCA) using SPSS Statistics for Mac ver. 17 (SPSS, Inc., Chicago, IL, USA) to reduce BIO variables to a smaller system of uncorrelated, independent variables and then used the PC loadings to examine habitat and ecological variation potentially related to population divergence. Differences in PC scores between genetic groups were tested with one-way ANOVAs in SPSS. RESULTS Divergence time estimation and genealogical relationships The aligned sequences of ITS and trnL-F were 640 and 381 base pairs (bp) long, respectively. Out of 1021 bp of the concatenated sequences (n = 197), 590 sites were variable and out of these 464 were parsimony informative. In P. mayanus samples, we detected polymorphic sites in several samples of the ITS region, with 17 of the 564 sites (3%). These were coded as ambiguous. The BEAST analyses placed the origin of the P. mayanus clade in the mid Pleistocene, 1.747 Mya (95% HPD 3.172–0.519 Mya; PP = 0.99; Fig. 1). Although a split estimated at 0.962 Mya (95% HPD 1.738–0.263 Mya) between samples of P. mayanus from Chiapas and the Petén province and samples from the Yucatán province is recovered with strong support (PP = 0.99), the relationships in each geographic region were poorly resolved, with samples from Chiapas intermingled among the Yucatán province clade (Fig. 1). Figure 1. View largeDownload slide (A) Psittacanthus mayanus. Photograph by Andrés Ortiz Rodríguez. (B) Chronogram based on a Bayesian approach using a coalescent prior under an uncorrelated lognormal relaxed clock model and assuming constant population size of P. mayanus DNA sequences in BEAST. Blue bars indicate 95% highest posterior density (HPD) intervals for nodes of particular interest. These nodes all have posterior probabilities >0.9. Figure 1. View largeDownload slide (A) Psittacanthus mayanus. Photograph by Andrés Ortiz Rodríguez. (B) Chronogram based on a Bayesian approach using a coalescent prior under an uncorrelated lognormal relaxed clock model and assuming constant population size of P. mayanus DNA sequences in BEAST. Blue bars indicate 95% highest posterior density (HPD) intervals for nodes of particular interest. These nodes all have posterior probabilities >0.9. The aligned ITS dataset for 82 samples of P. mayanus (564 bp) yielded 14 ribotypes (Table 2). In the ribotype network (Fig. 2), ribotype R1 was found in 53 of the 82 samples for ten of the 15 sampled populations and forms the core of the network. This widespread and most frequent ribotype (64.6% of the individuals and 66.6% of the sampled populations), inferred as ancestral based on its frequency and position in the network (Fig. 2), was separated by one or two mutations of most ribotypes in the network. The second most frequent ribotypes, R11 and R8 (8 and 6 samples, respectively) were mainly found in populations of Campeche and Chiapas and shared by one sample from southern Yucatán (Table 2). Some ribotypes had certain geographical structure and were exclusively distributed in localities of the Yucatán Peninsula (R2–R7) or Chiapas (R9–10, R12–R14) (Fig. 2). Table 2. Number of genetically analysed samples (n) for each molecular marker (ITS and trnL-F) and number of distinct ribotypes (R) and haplotypes (H) found in sampled Psittacanthus mayanus individuals, and the number of individuals per ribotype or haplotype in parentheses. Population number  Location  n  ITS ribotype  n  trnL-F haplotype  1  Yucatán, Dzibilchaltún  3  R1(2), R3(1)  3  H1(3)  2  Yucatán, Hunucmá, Sisal  14  R1(14)  13  H1(12), H3(1)  3  Yucatán, Celestún  9  R1(7), R2(2)  9  H1(9)  4  Yucatán, Mérida, Cuxtal  13  R1(11), R4(1), R6(1)  11  H1(10), H2(1)  5  Yucatán, Maxcanú, Límite Estatal  6  R1(5), R8(1)  7  H1(7)  6  Campeche, Champotón, Seybaplaya  9  R1(6), R3(1), R5(1), R7(1)  10  H1(10)  7  Campeche, Escárcega, Sabancuy  10  R1(3), R8(4), R11(3)  10  H1(10)  8  Chiapas, Bochil, Km 205 Pto. Caté  1  R11(1)  3  H1(3)  9  Chiapas, Bochil, Sn Cristobalito  1  R12(1)  3  H6(2), H7(1)  10  Chiapas, Sima de las Cotorras  6  R1(3), R10(1), R11(1), R14(1)  6  H1(6)  11  Chiapas, Ocozocuautla, El Aguacero  1  R1(1)  1  H4(1)  12  Chiapas, Ocozocuautla  1  R14(1)  1  H1(1)  13  Chiapas, Arriaga, Km 89–92  2  R9(2)  2  H1(1), H5(1)  14  Chiapas, Km 9 Arriaga-Tuxtla    -  1  H1(1)  15  Chiapas, Arriaga, La Aurora  1  R1(1)  1  H1(1)  16  Chiapas, Arriaga, Km 28.3  5  R8(1), R11(3), R13(1)  6  H1(6)  Population number  Location  n  ITS ribotype  n  trnL-F haplotype  1  Yucatán, Dzibilchaltún  3  R1(2), R3(1)  3  H1(3)  2  Yucatán, Hunucmá, Sisal  14  R1(14)  13  H1(12), H3(1)  3  Yucatán, Celestún  9  R1(7), R2(2)  9  H1(9)  4  Yucatán, Mérida, Cuxtal  13  R1(11), R4(1), R6(1)  11  H1(10), H2(1)  5  Yucatán, Maxcanú, Límite Estatal  6  R1(5), R8(1)  7  H1(7)  6  Campeche, Champotón, Seybaplaya  9  R1(6), R3(1), R5(1), R7(1)  10  H1(10)  7  Campeche, Escárcega, Sabancuy  10  R1(3), R8(4), R11(3)  10  H1(10)  8  Chiapas, Bochil, Km 205 Pto. Caté  1  R11(1)  3  H1(3)  9  Chiapas, Bochil, Sn Cristobalito  1  R12(1)  3  H6(2), H7(1)  10  Chiapas, Sima de las Cotorras  6  R1(3), R10(1), R11(1), R14(1)  6  H1(6)  11  Chiapas, Ocozocuautla, El Aguacero  1  R1(1)  1  H4(1)  12  Chiapas, Ocozocuautla  1  R14(1)  1  H1(1)  13  Chiapas, Arriaga, Km 89–92  2  R9(2)  2  H1(1), H5(1)  14  Chiapas, Km 9 Arriaga-Tuxtla    -  1  H1(1)  15  Chiapas, Arriaga, La Aurora  1  R1(1)  1  H1(1)  16  Chiapas, Arriaga, Km 28.3  5  R8(1), R11(3), R13(1)  6  H1(6)  Codes are from networks in Figure 3 View Large Figure 2. View largeDownload slide Geographical distribution and statistical parsimony networks of Psittacanthus mayanus overlaid on a map of the Yucatán Peninsula and adjacent states of southern Mexico. (A) Fourteen ITS ribotypes. (B) Seven trnL-F haplotypes. Ribotype or haplotype designations in the network correspond to those in Table 2, the size of the circles is proportional to the frequency of each ribotype or haplotype, and small black-filled circles are non-sampled ribotypes or haplotypes. Each ribotype or haplotype is coded with a different colour. The numbers in the ribotypes or haplotypes indicate the number of individuals that share the ribotype or haplotype. Pie charts represent ribotypes or haplotypes found in each sampling locality. The size of sections of the pie charts corresponds to the proportion of individuals with a given ribotype or haplotype. Figure 2. View largeDownload slide Geographical distribution and statistical parsimony networks of Psittacanthus mayanus overlaid on a map of the Yucatán Peninsula and adjacent states of southern Mexico. (A) Fourteen ITS ribotypes. (B) Seven trnL-F haplotypes. Ribotype or haplotype designations in the network correspond to those in Table 2, the size of the circles is proportional to the frequency of each ribotype or haplotype, and small black-filled circles are non-sampled ribotypes or haplotypes. Each ribotype or haplotype is coded with a different colour. The numbers in the ribotypes or haplotypes indicate the number of individuals that share the ribotype or haplotype. Pie charts represent ribotypes or haplotypes found in each sampling locality. The size of sections of the pie charts corresponds to the proportion of individuals with a given ribotype or haplotype. The aligned trnL-F dataset for 87 samples of P. mayanus (275 bp) yielded seven haplotypes (Fig. 2; Table 2). Haplotype H1, the most frequent (80 samples, 91.9% of the individuals), was retrieved from most sampled populations located both along the Yucatán Peninsula and in Chiapas (14 of the 16 sampled populations). Haplotypes H2 and H3 were singletons one step separated from H1 and located in northern Yucatán, and haplotypes H4–H7, only located in Chiapas, were also singletons or shared by two individuals (Fig. 2). Genetic diversity and population structure Differentiation among populations based on ITS variation indicated genetic subdivision (GST = 0.210, SE = 0.0604). Genetic diversity across all populations (hT = 0.667, SE = 0.1004; vT = 0.671, SE = 0.1046) was higher than within-population genetic diversity (hS = 0.527, SE = 0.0838; vS = 0.486, SE = 0.0884). PERMUT analysis showed that NST (0.276, SE = 0.0610) was not significantly higher (p > 0.05) than GST, indicating no phylogeographical structure. Similarly, phylogeographical structure was not observed based on trnL-F variation (GST = 0.461, SE = not calculated, NST = 0.609, SE = not calculated; NST = GST, p > 0.05), with genetic diversity across all populations (hT = 0.245, SE = 0.1457; vT = 0.248, SE = 0.1781) higher than within-populations (hS = 0.132, SE = 0.0612; vS = 0.097, SE = 0.0594). Pairwise comparisons of FST values were low but significant for both ITS and trnL-F when two groups of populations (Yucatán and Chiapas) or the most distant groups (YUC and CHI or dryYUC and dryCHI) in the three-group comparisons were compared (Table 3). Table 3. Pairwise comparisons of FST values of ITS (above the diagonal) and trnL-F (below the diagonal) among groups of Psittacanthus mayanus populations. Grouping  A. Two groups by geography  Yucatan  Chiapas    Yucatán  –  0.13088**    Chiapas  0.14846**  –    B. Three groups by geography  YUC  PET  CHI  YUC  _  0.07831*  0.21593***  PET  −0.02018  –  0.01033  CHI  0.11320*  0.11439  –  C. Three groups by precipitation  dryYUC  humYUC  dryCHI  dryYUC  –  0.21518**  0.23890***  humYUC  0.15614*  –  −0.04345  dryCHI  0.00791  −0.00197  –  Grouping  A. Two groups by geography  Yucatan  Chiapas    Yucatán  –  0.13088**    Chiapas  0.14846**  –    B. Three groups by geography  YUC  PET  CHI  YUC  _  0.07831*  0.21593***  PET  −0.02018  –  0.01033  CHI  0.11320*  0.11439  –  C. Three groups by precipitation  dryYUC  humYUC  dryCHI  dryYUC  –  0.21518**  0.23890***  humYUC  0.15614*  –  −0.04345  dryCHI  0.00791  −0.00197  –  Significant values at P < 0.001 are highlighted in bold. Group abbreviations are as follows: Yucatán = Yucatan Peninsula; Chiapas = Estado de Chiapas; YUC = Yucatán Province; PET = Petén Province; dryYUC = dry Yucatán Peninsula, humYUC = humid Yucatán Peninsula; dryCHI = dry Central Depression of Chiapas. * P < 0.01, ** P < 0.001, *** P < 0.0001. View Large The AMOVAs showed that 16.5% of genetic variation for ITS and 18.4% for trnL-F was explained by differences within populations and 83.5% and 81.6%, respectively, by differences between populations when all locations were treated as a single group (Table 4). The AMOVAs for ITS and trnL-F revealed population structure, with high and significant FCT values obtained when populations are grouped into two lineages corresponding to those along the Yucatán Peninsula and in Chiapas (Table 4). When sampling sites are grouped into tree lineages separated by geography or precipitation, a significant proportion of the ITS variation was also attributed to differences between groups (Table 4). Nucleotide diversity (π) was low for both datasets (Table 5). Table 4. AMOVA results for plausible Psittacanthus mayanus populations with no groups defined a priori (A), and grouping based on geography (two geographical groups: Yucatán Peninsula and Chiapas; three geographic groups: YUC = Yucatán province, PET = Petén province, CHI = Chiapas) (B and C) or based on the precipitation gradient (dryYUC, humYUC, dryCHI) (D). Source of variation  ITS  trnL-F  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  A. Single group   Among populations  10  10.417  0.08467  16.54    11  1.338  0.01054  18.38     Within populations  71  30.342  0.42735  83.46  FST = 0.16**  75  3.513  0.04684  81.62  FST = 0.18*   Total  81  40.759  0.51202      86  4.851  0.05738      B. Two geographic groups   Among groups  1  4.290  0.12918  21.85  FCT = 0.22**  1  0.393  0.00894  14.32  FCT = 0.14**   Among pop. within groups  9  6.127  0.03464  5.86  FSC = 0.07 ns  10  0.945  0.00667  10.68  FSC = 0.12 ns   Within populations  71  30.342  0.42735  72.29  FST = 0.28**  75  3.513  0.04684  75.00  FST = 0.25**  Total  81  40.759  0.59118      86  4.851  0.06245      C. Three geographic groups   Among groups  2  5.468  0.08492  15.74  FCT = 0.16*  2  0.408  0.00324  5.55  FCT = 0.05 ns   Among pop. within groups  8  4.949  0.02738  5.07  FSC = 0.06 ns  9  0.931  0.00830  14.21  FSC = 0.15 ns   Within populations  71  30.342  0.42735  79.19  FST = 0.21**  75  3.513  0.04684  80.23  FST = 0.20**   Total  81  40.759  0.53965      86  4.851  0.05838      D. Three precipitation groups   Among groups  2  6.967  0.14687  25.55  FCT = 0.25***  2  0.311  0.00158  2.73  FCT = 0.03 ns   Among pop. within groups  8  3.450  0.00055  0.10  FSC = 0.00 ns  9  1.027  0.00959  16.52  FSC = 0.17*   Within populations  71  30.342  0.42735  74.35  FST = 0.26**  75  3.513  0.04684  80.74  FST = 0.19**   Total  81  40.759  0.57477      86  4.851  0.05801      Source of variation  ITS  trnL-F  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  d.f.  Sum of squares  Estimated variance  % of variation  Fixation indices  A. Single group   Among populations  10  10.417  0.08467  16.54    11  1.338  0.01054  18.38     Within populations  71  30.342  0.42735  83.46  FST = 0.16**  75  3.513  0.04684  81.62  FST = 0.18*   Total  81  40.759  0.51202      86  4.851  0.05738      B. Two geographic groups   Among groups  1  4.290  0.12918  21.85  FCT = 0.22**  1  0.393  0.00894  14.32  FCT = 0.14**   Among pop. within groups  9  6.127  0.03464  5.86  FSC = 0.07 ns  10  0.945  0.00667  10.68  FSC = 0.12 ns   Within populations  71  30.342  0.42735  72.29  FST = 0.28**  75  3.513  0.04684  75.00  FST = 0.25**  Total  81  40.759  0.59118      86  4.851  0.06245      C. Three geographic groups   Among groups  2  5.468  0.08492  15.74  FCT = 0.16*  2  0.408  0.00324  5.55  FCT = 0.05 ns   Among pop. within groups  8  4.949  0.02738  5.07  FSC = 0.06 ns  9  0.931  0.00830  14.21  FSC = 0.15 ns   Within populations  71  30.342  0.42735  79.19  FST = 0.21**  75  3.513  0.04684  80.23  FST = 0.20**   Total  81  40.759  0.53965      86  4.851  0.05838      D. Three precipitation groups   Among groups  2  6.967  0.14687  25.55  FCT = 0.25***  2  0.311  0.00158  2.73  FCT = 0.03 ns   Among pop. within groups  8  3.450  0.00055  0.10  FSC = 0.00 ns  9  1.027  0.00959  16.52  FSC = 0.17*   Within populations  71  30.342  0.42735  74.35  FST = 0.26**  75  3.513  0.04684  80.74  FST = 0.19**   Total  81  40.759  0.57477      86  4.851  0.05801      ns, not significant (P > 0.05), * P < 0.05, ** P < 0.001, *** P < 0.0001 View Large Table 5. Summary statistics of neutrality tests and demographic analysis of Psittacanthus mayanus samples by genetic groups to infer demographic range expansion. Parameter  Two groups  Three groups by geography  Three groups by precipitation  Yucatán  Chiapas  YUC  PET  CHI  dryYUC  humYUC  dryCHI  ITS  N  64  18  45  19  18  54  12  16  NH  8  6  6  5  6  8  3  5  h, gene diversity  0.3616  0.6797  0.2505  0.5789  0.6797  0.3068  0.5909  0.6667  π, nucleotide diversity  0.001871  0.001818  0.001164  0.003088  0.001818  0.001704  0.002055  0.001788  Tajima’s D  –1.467*  –0.420  –2.044**  –0.767  –0.420  –1.769*  0.472  –0.576  Fu’s Fs  –2.473  –2.449*  –2.497*  0.183  –2.449*  –3.086*  0.951  –1.469  SDD  0.0215  0.0051  0.0068  0.4324***  0.0051  0.0123  0.0880  0.0066  Hri  0.3059  0.0745  0.4096  0.2321  0.0744  0.3366  0.3289  0.0725  trnL-F  N  63  24  43  20  24  53  16  18  NH  3  5  3  1  5  3  3  3  h, gene diversity  0.0620  0.3913  0.0919  0.0000  0.3768  0.0747  0.3417  0.2157  π, nucleotide diversity  0.000227  0.002314  0.000338  0.000000  0.002227  0.000274  0.001692  0.000808  Tajima’s D  –1.437*  –0.317  –1.479*  n.a.  –0.317  –1.458*  0.000  –1.165  Fu’s Fs  –3.356**  –2.187*  –2.827**  n.a.  –2.187*  –3.113**  –0.571  –1.743*  SDD  0.0002  0.2109***  0.0001  n.a.  0.2109***  0.0000  0.1937  0.0024  Hri  0.7687  0.2669  0.6759  n.a.  0.2669  0.7299  0.2079  0.3719  Parameter  Two groups  Three groups by geography  Three groups by precipitation  Yucatán  Chiapas  YUC  PET  CHI  dryYUC  humYUC  dryCHI  ITS  N  64  18  45  19  18  54  12  16  NH  8  6  6  5  6  8  3  5  h, gene diversity  0.3616  0.6797  0.2505  0.5789  0.6797  0.3068  0.5909  0.6667  π, nucleotide diversity  0.001871  0.001818  0.001164  0.003088  0.001818  0.001704  0.002055  0.001788  Tajima’s D  –1.467*  –0.420  –2.044**  –0.767  –0.420  –1.769*  0.472  –0.576  Fu’s Fs  –2.473  –2.449*  –2.497*  0.183  –2.449*  –3.086*  0.951  –1.469  SDD  0.0215  0.0051  0.0068  0.4324***  0.0051  0.0123  0.0880  0.0066  Hri  0.3059  0.0745  0.4096  0.2321  0.0744  0.3366  0.3289  0.0725  trnL-F  N  63  24  43  20  24  53  16  18  NH  3  5  3  1  5  3  3  3  h, gene diversity  0.0620  0.3913  0.0919  0.0000  0.3768  0.0747  0.3417  0.2157  π, nucleotide diversity  0.000227  0.002314  0.000338  0.000000  0.002227  0.000274  0.001692  0.000808  Tajima’s D  –1.437*  –0.317  –1.479*  n.a.  –0.317  –1.458*  0.000  –1.165  Fu’s Fs  –3.356**  –2.187*  –2.827**  n.a.  –2.187*  –3.113**  –0.571  –1.743*  SDD  0.0002  0.2109***  0.0001  n.a.  0.2109***  0.0000  0.1937  0.0024  Hri  0.7687  0.2669  0.6759  n.a.  0.2669  0.7299  0.2079  0.3719  N = number of individuals, NH = number of ribotypes or haplotypes, SDD = differences in the sum of squares or mismatch distribution, Hri = Harpending’s raggedness index. DT and FS positive values are indicative of mutation-drift-equilibrium, which is typical of stable populations, and negative values that result from an excess of rare haplotypes indicate that populations have undergone recent expansions, often preceded by a bottleneck. Significantly negative values (at the 0.05 level) reveal in both tests historic demographic expansion events. Significant (P ≤ 0.05) SSD and Hri values indicate deviations from the sudden expansion model. n.a. = not available; *P < 0.05; **P < 0.01; ***P < 0.001. Group abbreviations are as follows: PET = Petén province, YUC = Yucatán province, CHI = Chiapas, dryYUC = dry Yucatán, humYUC = humid Yucatán, dryCHI = dry Central Depression of Chiapas View Large Demographic history When populations were analysed as two (Yucatán Peninsula and Chiapas) or three groups based on geography (YUC, PET, CHI) or precipitation (dryYUC, humYUC, dryCHI), negative and significant values for neutrality tests (Fu’s Fs, Tajima’s D) were observed (Table 5), indicating a departure from neutrality (probably an excess of rare alleles or new mutations from population expansion) for northern peninsular P. mayanus (Yucatán, YUC, dryYUC) populations in most cases. The model of sudden demographic expansion for most cases of Yucatán groups was not rejected by the generalized least square procedure (SSD) or by the raggedness index of the distribution (Hri), except for three cases, the Chiapas and CHI groups in the trnL-F dataset and PET in the ITS dataset (Table 5). The Bayesian skyline plots suggested that Chiapas and Yucatán Peninsula populations underwent a constant effective population size, with a marginal increase in population size occurring after the MH (Fig. 3). When groups of populations were pooled, the increase in population size occurring after the MH was more pronounced (Fig. 3). Figure 3. View largeDownload slide Bayesian skyline plots showing historical demographic trends of Psittacanthus mayanus for YUC (A), CHI (B), and for all the P. mayanus (C) populations using the concatenated data set. The y-axis is the product between effective population size and the generation time and the x-axis is time in thousands of years. Solid lines represent mean estimates and shaded areas represent 95% confidence intervals. YUC = Yucatán Peninsula, CHI = Chiapas. Figure 3. View largeDownload slide Bayesian skyline plots showing historical demographic trends of Psittacanthus mayanus for YUC (A), CHI (B), and for all the P. mayanus (C) populations using the concatenated data set. The y-axis is the product between effective population size and the generation time and the x-axis is time in thousands of years. Solid lines represent mean estimates and shaded areas represent 95% confidence intervals. YUC = Yucatán Peninsula, CHI = Chiapas. According to IMa results (summarized in Table 6), the peninsular (YUC) population size (Ne1) is estimated to be larger than population size of ancestral (Nea) and Chiapas (CHI, Ne2) descendant population. Divergence time between YUC and CHI populations was estimated at c. 66 kya. A peak at the left-hand side of the posterior distribution of t is observed, with a long right-hand tail over which the probability distribution remains flat. This means that this parameter value needs to be interpreted with caution because parameter values have essentially zero probability for most of the posterior distribution and the only values that are well supported are those less than 33 kya. When testing for migration during the time since population splitting, migration among populations occurred in one direction, moving from CHI to YUC (Table 6). Table 6. Results of isolation-with-migration model (IMa) for the split between peninsular (YUC) and Chiapas (CHI) populations of Psittacanthus mayanus.   Model parameter estimates  θ1  θ2  θa  t  m2-to-1  m1-to-2  YUC vs. CHI  Mean  1.7503  0.1994  0.9657  3.029  67.43  0.483  HPD95Lo  0.1553  0.0037  0.0479  0.199  17.31  0.024  HPD95Hi  3.4506  0.8537  1.9043  5.851  98.75  0.945    Demographic parameter estimates  Ne1  Ne2  Nea  t  Nm2-to-1  Nm1-to-2  YUC vs. CHI  Mean  868.487  98.939  479.165  66,137  59.012  0.048  HPD95Lo  77.074  1.819  23.751  4,344  1.344  0.000  HPD95Hi  1712.132  423.592  944.883  127,739  170.373  0.403    Model parameter estimates  θ1  θ2  θa  t  m2-to-1  m1-to-2  YUC vs. CHI  Mean  1.7503  0.1994  0.9657  3.029  67.43  0.483  HPD95Lo  0.1553  0.0037  0.0479  0.199  17.31  0.024  HPD95Hi  3.4506  0.8537  1.9043  5.851  98.75  0.945    Demographic parameter estimates  Ne1  Ne2  Nea  t  Nm2-to-1  Nm1-to-2  YUC vs. CHI  Mean  868.487  98.939  479.165  66,137  59.012  0.048  HPD95Lo  77.074  1.819  23.751  4,344  1.344  0.000  HPD95Hi  1712.132  423.592  944.883  127,739  170.373  0.403  Model parameters indicate estimates without use of molecular rate of evolution for six parameters (IMa output values). Values are averages of two runs of mean parameter estimates and the 95% highest posterior densities (HPD) intervals of each parameter: effective population size of the ancestral (θa) and descendant populations (θ1, θ2), effective migration rates in both directions (m1-to-2 and m2-to-1), and time since divergence (t) at which the ancestral population gave rise to the descendant populations. Demographic rates represent parameters scaled to rates of molecular evolution: effective population sizes (Ne, individuals), population migration rates (Nm, migrants per generation) and divergence time (t, years). Population size (Ne) based on the average generation time (T) of 11 years for a high (0.9) annual adult survival rate. View Large Population expansion and approximate Bayesian computation (ABC) DIYABC analysis indicated that population growth of the YUC population (scenario 1) is the best-supported scenario (Fig. 4), with a higher posterior probability value and 95% confidence intervals (PP, 95% CI; scenario 1: 0.882, 0.869–0.896) that did not overlap with those obtained for scenario 2 (0.117, 0.104–0.131). Analyses to estimate confidence in scenario choice indicated that type I and type II errors for the best-supported scenario were low (0.354 and 0.383, respectively). Under scenario 1, posterior mean parameter estimates indicated that the divergence between dryCHI and YUC groups (t3) occurred 43 900 years before present (CI: 61.8–40.1 kya) and the split between humYUC and dryYUC (t1) was dated to 1510 (4.15–1.02 kya) (Appendix 2). In accordance with the population expansion model, the DIYABC analysis supported a scenario of a YUC population size increase around the LGM (t2 =15 600 years before present, 20.8–10.3 kya; Appendix 2). Estimated mean mutation rates of nrDNA and cpDNA were estimated to be 2.03 × 10–7 and 3.25 × 10–7, respectively (Appendix 2). Figure 4. View largeDownload slide (A) Schematic representation of two demographic scenarios of Psittacanthus mayanus tested with DIYABC and the prior values used in DIYABC simulations. N represents effective population sizes (Ne) and t indicates divergence time. Timeline from present to past indicates Mid Holocene (MH), Last Glacial Maximum (LGM) and Last Interglacial (LIG). (B) The posterior probability of scenarios was assessed using a weighted logistic regression on the 1% of simulated datasets closest to the observed data and, for the best-supported scenario (scenario 1). (C) A principal components analysis was used to test the goodness-of-fit of each evaluated scenario against the simulated data in DIYABC. Figure 4. View largeDownload slide (A) Schematic representation of two demographic scenarios of Psittacanthus mayanus tested with DIYABC and the prior values used in DIYABC simulations. N represents effective population sizes (Ne) and t indicates divergence time. Timeline from present to past indicates Mid Holocene (MH), Last Glacial Maximum (LGM) and Last Interglacial (LIG). (B) The posterior probability of scenarios was assessed using a weighted logistic regression on the 1% of simulated datasets closest to the observed data and, for the best-supported scenario (scenario 1). (C) A principal components analysis was used to test the goodness-of-fit of each evaluated scenario against the simulated data in DIYABC. Palaeodistribution modelling and environmental variation The ENM yielded a good fit for the current geographical distribution of P. mayanus. The AUC value for the replicate runs was 0.993 ± 0.002 (Fig. 5). Determination of the threshold probability for predicted presence using TSS resulted in a mean proportion of correctly classified training observations of 0.510 ± 0.229. The palaeodistribution modelling revealed that suitable habitat for P. mayanus was out of its current distribution during the LIG and restricted to small areas in southern Guatemala with moderate-to-high suitability (Fig. 5). However, the predictions for the LGM applying both CCSM and MIROC simulations revealed that conditions of suitable habitat potentially expanded the distribution of P. mayanus, particularly in Guatemala and Chiapas, and an isolated portion of Quintana Roo with high probability according to the CCSM model and, with a high probability, a more continuous distribution in Guatemala, Chiapas, Belize and Quintana Roo according to the MIROC model (Fig. 5). Conditions of suitable habitat were contracted during the MH, with suitable habitat restricted to small areas between Belize and Quintana Roo, Central Depression of Chiapas (CCSM), and the northernmost part of the Yucatán Peninsula. Overall, the comparison between present and past climatic conditions predicted by the models suggests that climatic conditions for P. mayanus suitable habitat experienced range shifts from southern Guatemala into the Yucatán Peninsula, with expansion from LIG to LGM, contraction from LGM to MH and colonization of the northernmost portion of the Yucatán Peninsula and expansion from MH to present. Figure 5. View largeDownload slide Results from the MAXENT analyses showing binary transformation of ecological niche models for Psittacanthus mayanus at Last Inter Glacial (LIG, 140–120 kya), Last Glacial Maximum (LGM, CCSM and MIROC scenarios, 21 kya), Mid-Holocene (MH, CCSM and MIROC scenarios, 6 kya), and at present. Vegetation map on the top right indicates the biogeographic regions discussed here (Morrone, 2005, 2014). Yellow indicates binary transformation of occurrence using the threshold that maximized the true skill statistics (TSS). Yellow dots shown on the right map below indicate occurrence records. Figure 5. View largeDownload slide Results from the MAXENT analyses showing binary transformation of ecological niche models for Psittacanthus mayanus at Last Inter Glacial (LIG, 140–120 kya), Last Glacial Maximum (LGM, CCSM and MIROC scenarios, 21 kya), Mid-Holocene (MH, CCSM and MIROC scenarios, 6 kya), and at present. Vegetation map on the top right indicates the biogeographic regions discussed here (Morrone, 2005, 2014). Yellow indicates binary transformation of occurrence using the threshold that maximized the true skill statistics (TSS). Yellow dots shown on the right map below indicate occurrence records. Four principal components (PC) explained 91.4% of the total environmental variance in P. mayanus (Appendix 3; Fig. 6). The PC1 (45.7% of total variation) was positively associated with temperature (BIO1, BIO6, BIO9–BIO11) and precipitation variables (BIO14, BIO17, BIO19) and negatively with precipitation seasonality (BIO15), whereas PC2 (23.9%) was positively associated with mean temperature of wettest quarter (BIO8) and negatively associated with precipitation variables (BIO12–BIO13, BIO18) (Appendix 3). Univariate ANOVAs of PC scores, using grouping by precipitation as a fixed factor, showed strong significant environmental differences among groups in PC1 scores (Appendix 3, Fig. 6) and differences between dryYUC and the other groups (dryCHI and humYUC) in PC2 scores (Appendix 3, Fig. 6). Figure 6. View largeDownload slide (A) Map relating to annual precipitation (BIO12) and location of precipitation differences on the Yucatán Peninsula. The approximate divisions of the dry Yucatán (dryYUC), humid Yucatán (humYUC) and dry Central Depression of Chiapas (dryCHI) are indicated with different colours. White dots shown on the map indicate collection localities and coloured dots indicate occurrence records. (B) Principal components analysis (PCA) on the 19 BIOCLIM variables with the first principal component (PC1) largely a measure of temperature differences, and the second principal component (PC2) mainly determined by precipitation measures. Symbols correspond to populations of Psittacanthus mayanus on the dryYUC (blue), humYUC (green), and dryCHI populations (yellow). Means (± 1 SE) of PC1 (C) and PC2 (D) scores with different superscript letters are significantly different between groups (P < 0.05). Figure 6. View largeDownload slide (A) Map relating to annual precipitation (BIO12) and location of precipitation differences on the Yucatán Peninsula. The approximate divisions of the dry Yucatán (dryYUC), humid Yucatán (humYUC) and dry Central Depression of Chiapas (dryCHI) are indicated with different colours. White dots shown on the map indicate collection localities and coloured dots indicate occurrence records. (B) Principal components analysis (PCA) on the 19 BIOCLIM variables with the first principal component (PC1) largely a measure of temperature differences, and the second principal component (PC2) mainly determined by precipitation measures. Symbols correspond to populations of Psittacanthus mayanus on the dryYUC (blue), humYUC (green), and dryCHI populations (yellow). Means (± 1 SE) of PC1 (C) and PC2 (D) scores with different superscript letters are significantly different between groups (P < 0.05). DISCUSSION Our survey of molecular sequence data in populations of P. mayanus identified two genetic groups in the Yucatán Peninsula. Phylogeographic and distributional projections derived from ecological niche modelling are consistent with the hypothesis that the directionality of colonization was from southern Chiapas to the north on the Yucatán Peninsula and that the Altos de Chiapas mountain areas may have restricted northward gene flow. Together, our results suggest a recent geographical northern expansion into the seasonally tropical deciduous forest, consistent with a colonization model and a pattern of population expansions/contractions for species responding to the Pleistocene glacial cycles and reduced gene flow since by both physical and environmental barriers. Given that this region is poorly studied at the phylogeographic level, the phylogeography of P. mayanus in combination with a niche-modelling approach is particularly important in understanding the evolution of the biodiversity of Mexican lowland forests. Genetic diversity, structure and phylogeographic analysis Genetic diversity analysis for P. mayanus populations revealed a higher gene diversity (h) and nucleotide diversity (π) in Chiapas than in the northern part of the distribution on the Yucatán Peninsula (Table 5 and Fig. 2). The genetic diversity within populations and total gene diversity of P. mayanus showed that the species has low levels of genetic diversity relative to that observed in other mistletoe species (Zuber & Widmer, 2009; Amico et al., 2014; Nyagumbo, Ndagurwa & Mushonga, 2017) or that of closely related and more widespread congeners (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Because the fruits of P. mayanus are consumed and dispersed by several bird species, including the endemic (to the Yucatán Peninsula) Patagioenas cayennensis and Piranga roseogularis and the widespread Myiozetetes similis Pitangus sulphuratus and Megarhynchus pitangua, it is possible that the observed genetic structure and genetic differentiation have been eroded by large-scale seed recruitment and high historical rates of gene flow (Ornelas et al., 2016). The haplotype network and AMOVA results suggest genetic differentiation between two groups of populations for P. mayanus, from Chiapas and the Yucatán Peninsula. Similar results were also observed for Guaiacum sanctum L. (Zygophyllaceae) trees using nuclear microsatellites (Oyama et al., 2016). The strongest genetic differentiation between the groups of P. mayanus populations in Chiapas, located approximately between 16 and 17ºN, and the Yucatán Peninsula, located between 18 and 21ºN, corresponds to historically different biogeographic regions, and occurs between 600 and 1000 m a.s.l. to nearly above sea level, respectively (Table 1). The population groups of P. mayanus showed, however, haplotype sharing and short-term divergence (Figs 1 and 2). The widespread distribution and higher frequency of haplotype H1 (trnL-F) is not surprising considering that several widespread bird species mediate mistletoe seed dispersal. Migration among groups of populations occurred in one direction, moving from Chiapas to the Yucatán Peninsula, suggesting that the barriers between these groups of populations associated with different geographical regions and environmental conditions are permeable. A lower and similar genetic differentiation among groups of populations (haplotype network and AMOVA results) was found in the ITS when populations of P. mayanus were separated by geography and precipitation into three groups, respectively (Tables 1, 4–5, Fig. 2). It is possible that these habitat differences are related to the formation of genetic groups. PC1 and PC2 explained 70% of the total variance in current climatic conditions among groups of populations, respectively. PC1 was interpretable as an index of precipitation variability. For example, PC1 exhibited high positive loading for annual precipitation variables and low or negative loading for precipitation seasonality and temperature variables. Thus, high scores of PC1 are indicative of wet climatic conditions and low or negative PC1 scores of seasonal climatic conditions. Notably, the climatic PC1 scores differed between groups. Thus, in this dataset at least, the geographical position is relatively dependent of precipitation variability and seasonality. The fact that there are some environmental differences among areas does not necessarily point to adaptation as the driver of genetic differentiation, and thus these interpretations are correlative. Instead, the genetic structure may be a consequence of historical fragmentation in the region, generating responses to contractions and expansions of the semideciduous tropical dry forest along the Yucatán Peninsula during Pleistocene climatic cycles as suggested by the ENMs. Demographic history of P. mayanus Despite the observed genetic differentiation between groups of populations, the most frequent ribotypes and haplotypes of P. mayanus populations are widespread, suggesting effective gene flow between peninsular and Chiapas populations and signs of recent demographic expansion for populations from the northern part of the Yucatán Peninsula. Group of populations from Chiapas (CHI) showed no major changes in their effective population sizes, exhibiting a high genetic diversity; whereas group of populations of the northern portion of the Yucatán Peninsula (YUC) exhibited lower genetic diversity, no changes in effective population sizes and strong signatures of recent demographic expansion (Table 5, Fig. 3). Note that the historical population size changes or the spatial population expansion can be influenced by several sources of error (e.g. sample size, sampling scheme, levels of polymorphism, mutation rates) or by genetic hitchhiking of polymorphism linked to genes under positive selection (e.g. Grant, 2015). The levels of genetic differentiation found in P. mayanus reveal the presence of two genetic groups that corresponds with the distinction between samples from the Yucatán Peninsula and Chiapas, and that the pattern of differentiation is associated with ecological/elevational isolation. Our divergence time estimates indicate that the split between these groups of populations occurred during the Pleistocene, in agreement with divergence time estimates for lineages of other Psittacanthus spp. (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Therefore, Pleistocene climatic oscillations are invoked as the main factor driving mistletoe diversification in the region. The ENMs and gene flow estimates using IMa suggest that dispersal patterns among populations of P. mayanus occurred asymmetrically; the directionality of gene flow indicates that populations from Chiapas are the source of migrants into the Yucatán Peninsula. The facts that several samples from Chiapas clustered sparsely with samples of P. mayanus from the Yucatán Peninsula and no samples from the Yucatán Peninsula of P. mayanus clustered with samples of Chiapas is evidence of gene flow and suggests that the directionality of colonization was from southern Chiapas to the north of the Yucatán Peninsula and that the Altos de Chiapas mountain areas may have restricted northward gene flow. Approximate Bayesian computation (ABC) analyses strongly supported a scenario of post-glacial population growth and that the northern part of the Yucatán Peninsula was invaded by ancestral P. mayanus individuals from Chiapas. This scenario suggests that species-specific traits, availability of host tree species and/or variation of suitable habitat may be causing the phylogeographical pattern influenced by the Pleistocene glacial-interglacial cycles. A recent phylogeographic study of the bird-dispersed mistletoe Phoradendron californicum Nutt. (Santalaceae) in the Sonoran Desert showed that the overall plastid DNA variation was best explained by large-scale historical geographical events and changes in vegetation limiting seed dispersal at deeper evolutionary timescales but host species poorly explained plastid DNA variation (Lira-Noriega et al., 2015). More recent phylogeographic studies of Psittacanthus spp. using nrDNA and plastid DNA molecular markers also provide support for the predominant role of isolation and environmental factors in driving genetic differentiation between populations of P. schiedeanus from the cloud forests of eastern Mexico (Ornelas et al., 2016) and P. calyculatus from the temperate forests along the Trans-Mexican Volcanic Belt (Pérez-Crespo et al., 2017). Using nuclear microsatellite data, Ramírez-Barahona et al. (2017) showed that genetic differentiation among populations of the species complex composed of P. schiedeanus, P. calyculatus, P. breedlovei and P. angustifolius Kuijt was associated with environmental and host preferences, independent of geography. However, the observed genetic structure was best explained by environmental predictors rather than by host preferences, supporting the hypothesis that the occurrence of these parasites has been largely determined by its own climatic niche and to a lesser degree by host specificity. The palaeodistribution of P. mayanus on the Yucatán Peninsula The ENM predicted fragmentation in the southern regions of P. mayanus distribution and recent colonization of the northernmost part of the Yucatán Peninsula (Fig. 5). ENM results indicate that the range of P. mayanus expanded from the LIG to the LGM; it was wider and connected during the LGM and contracted during the MH and it is expanding under the warmer current post-glacial period (Fig. 5). The predicted LIG area is small. It is possible that the distribution of P. mayanus during this period was limited, but also that the species occurred in environments without a present analogue. The LIG palaeodistribution model shows a lack of suitable habitat in most of the Yucatán Peninsula 140 000–120 000 years ago. Based on the phylogenetic evidence and divergence time estimation suggesting that the species originated before the LIG, we do not interpret this ENM result as a hypothesis of complete absence of the species from the northern Central America. In the absence of any evidence to suggest a South American origin for the species ancestor we offer as an alternative hypothesis that the distribution of P. mayanus in the LIG was quite small or restricted to a different habitat in northern Central America. Other Psittacanthus spp. generally associated with habitats at higher elevation show a wider distribution during the LIG (Ornelas et al., 2016). Comparisons between lowland and highland mistletoe species will be an area for future palaeodistribution modelling in which development of additional LGM climate models may produce novel hypotheses. The two global circulation model simulations of the LGM climate (CCSM and MIROC) revealed dissimilar inferences regarding the potential palaeodistribution of P. mayanus. The CCSM model inferred two geographically separated regions with high suitability expanding northwards into partially the current distribution, one in the Pacific coast and lowlands of Chiapas and Guatemala, in the borderlands of southern Oaxaca and southern Guatemala, and the other one in the eastern side of the Yucatán Peninsula, stretching from the north to the south Caribbean slope of Quintana Roo. In contrast, the MIROC model inferred a single large area in the southern part of the Yucatán Peninsula including most of Chiapas, southern Tabasco, southern Quintana Roo, northern Guatemala and Belize. A common feature of both models, however, is the inference of suitable habitat in Chiapas and southern Quintana Roo (Fig. 5). The CCSM and MIROC models differ in temperature and precipitation patterns. The LGM climate as simulated by CCSM is colder and dryer than that simulated by MIROC. The CCSM model is consistent with the tropical refugia theory (Haffer, 1969), in which the cooler glacial climates characterized by a drastic reduction in precipitation and widespread aridity during glacial times caused tropical forests to become contracted in their distribution and fragmented into isolates of tropical forest patches separated by intervening non-forest vegetation (e.g. Haffer, 1969; Carnaval & Moritz, 2008). The knowledge on precipitation during the LGM, and hence on the distribution of Neotropical forests, is largely based on limited and conflicting palaeoecological data (Ramírez-Barahona & Eguiarte, 2013). It has been suggested that the summer precipitation regime was collapsed, resulting in arid conditions extending over much of Mexico, Central America and the Caribbean (Ramírez-Barahona & Eguiarte, 2013 and references therein). Accordingly, Anselmetti et al. (2006) suggested a reduction in wet season precipitation for upper Central America based on a reduction in the water volume of Lake Petén-Itzá, Guatemala. In contrast, Mueller et al. (2010) suggested that there was no effective decrease in humidity in upper Central America based on pollen records from Lake Petén-Itzá showing a rich clay deposit and a pollen diagram indicative of a montane forest community thriving at lower elevations under cool and moist conditions (Hodell et al., 2008; Bush et al., 2009). Both LGM models predict expanded suitable habitat compared with LIG prediction. Based on palynological evidence for the Petén region (Leyden, 1984), xeric vegetation was widely distributed during the LGM, which might have facilitated expansion of P. mayanus to coastal and central depression of Chiapas and Quintana Roo, and then it was supplanted by more mesic vegetation towards the Holocene transition and contraction of P. mayanus. Therefore, it is possible that the uncovered distribution on Quintana Roo, Guatemala and Belize had provided suitable conditions for P. mayanus during the LGM, MH and present. According to the work of Carnaval et al. (2009) on other tropical regions, this would correspond to an area of long-term population persistence, and thus a source of genetic variability during population expansion. Although we did not sample these areas for the present study, the two models identified similar areas of suitable habitat and the use of the two different climate models enabled us to identify the possibility of an area of long-term population persistence and to account for modelling uncertainty due to LGM climate based on limited and conflicting palaeoecological data (Ramírez-Barahona & Eguiarte, 2013). The distribution of potential suitable habitat appears to affect the distribution of P. mayanus, causing population differentiation. Several populations from both the northern and southern distribution of P. mayanus possess unique ribotypes and haplotypes, suggesting that the population of P. mayanus persisted in the southern region of the Yucatán Peninsula under LGM conditions where the ecological conditions were appropriate, and that the northern part of the Yucatán Peninsula was colonized more recently (MH) and contracted its distribution particularly in the Central Depression of Chiapas. The IMa results revealed that the divergence or splitting time (t) between descendant YUC and CHI populations overlaps with the tMRCA estimated using DIYABC. When testing for migration, migration among populations occurred predominantly from Chiapas to the Yucatán Peninsula. These estimates suggest a recent divergence scenario, in which haplotype sharing between regions results from retention of ancestral polymorphisms with not enough time for lineage sorting, ongoing gene flow across habitats or complete haplotype turnover. Along with these results, the haplotype network, FST pairwise comparisons and AMOVAs suggest moderate levels of genetic diversity among populations and the existence of three genetic groups when populations are grouped by the geography (YUC, PET, CHI) or by the precipitation gradient (dryYUC, humYUC, dryCHI), especially when considering the nrDNA (Tables 3 and 4). Altogether these results suggest that range shifts and changes in the habitat distribution that occurred during the Pleistocene climatic cycles might have fragmented the distribution of P. mayanus. According to the ENM, IMa and ABC results, the genetic differentiation of P. mayanus may be linked to the effects of the glacial/interglacial cycles and environmental factors. Under these scenarios, ancestral populations of P. mayanus colonized the semideciduous tropical rain forest of the Petén province during the LGM followed by northward expansion and colonization of the northern Yucatán province during the MH, currently covered by seasonally dry tropical deciduous forests, when the drier and cooler climate of the Peninsula and the savanna and Juniperus-matorral covering most of the region until the early Holocene (Islebe et al., 1996; Orellana et al., 2003) changed to a more humid and tropical-type forest (Metcalfe et al., 2000). Interpreting the resultant P. mayanus range dynamics through time is difficult without considering contemporaneous demography/occurrence of its principal host populations. Specifically, the putative recent expansion of this mistletoe into northern Yucatán might simply reflect broader changes in the host plants upon which it ultimately depends. If so, each of the genetic groups needs to be modelled by adding distribution data of each of the several host species to examine direct versus indirect determinants of distribution. Although such an approach would be the norm in phylogeographic studies of host-specialized parasitic plants, no data are available on the prevalence or specificity of these mistletoes for each of the recorded host species (Kuijt, 2009). Considering that parasitic plant distributions ultimately depend on host availability, future ENM should evaluate the effects of including the host distribution to predict distribution shifts of the mistletoe under past climates comparing models that use only climate variables with models that also take into account the biotic interactions with the hosts. CONCLUSIONS Our results revealed the overall phylogeographical patterns in P. mayanus, with two main haplogroups comprising groups of populations from the Yucatán Peninsula and Chiapas and low genetic differentiation between the Yucatán and Petén provinces. Despite the genetic differentiation of some P. mayanus populations, there are widespread ribotypes and haplotypes in the ITS and trnL-F markers, respectively. In particular, the low level of variation found in the ITS and trnL-F data and the presence of widespread ribotypes and haplotypes on the Yucatán Peninsula suggest effective nuclear and plastid gene flow, and negative values of Tajima’s D, Fu’s Fs and the results of mismatch distributions and BSPs are suggestive of recent demographic expansion for populations of P. mayanus in the region, most likely during the LGM. In contrast, the population differentiation of P. mayanus from Chiapas shown by ITS and plastid trnL-F datasets suggests that gene flow has been restricted. The ABC analysis strongly supported a scenario of population growth. Under this scenario, climatic fluctuations throughout the Pleistocene would have altered the distribution of suitable habitat for mistletoes along the Yucatán Peninsula leading to genetic differentiation of P. mayanus comprising moderate-to-low haplotype diversity and hinting widespread ribotypes and haplotypes lay on the Yucatán Peninsula involving geographical northern expansion into the seasonally tropical deciduous forest, consistent with a colonization model and a pattern of population expansions/contractions for species responding to the Pleistocene glacial cycles. ACKNOWLEDGEMENTS We thank C. Bárcenas, E. Gándara, E. Ruiz Sánchez, H. Gómez Domínguez and A. J. Vargas who helped in obtaining samples and/or generated sequences for this work. Permission to conduct our fieldwork was granted by the Mexican Government (Instituto Nacional de Ecología, Secretaría del Medio Ambiente y Recursos Naturales, SGPA/DGGFS/712/1299/12). We also thank C. Palma-Silva and two anonymous reviewers for their valuable criticisms and suggestions. This project was funded by a competitive grant (155686) from the Consejo Nacional de Ciencia y Tecnología (CONACyT; http://www.conacyt.mx) and research funds from the Departamento de Biología Evolutiva, Instituto de Ecología, A.C. (20030/10563) awarded to Juan Francisco Ornelas. A doctoral scholarship was granted from CONACyT to Y.L.V. (155686) and A.E.O.R. (262563). REFERENCES Allouche O, Tsoar A, Kadmon R. 2006. Assessing the accuracy of species distribution models: prevalence, kappa, and the true skill statistic (TSS). Journal of Applied Ecology  43: 1223– 1232. Google Scholar CrossRef Search ADS   Amico GC, Vidal-Russell R, Nickrent DL. 2007. Phylogenetic relationships and ecological speciation in the mistletoe Tristerix (Loranthaceae): the influence of pollinators, dispersers, and hosts. American Journal of Botany  94: 558– 567. Google Scholar CrossRef Search ADS PubMed  Amico GC, Vidal-Russell R, García MA, Nickrent DL. 2012. Evolutionary history of the South American mistletoe Tripodanthus (Loranthaceae) using nuclear and plastid markers. Systematic Botany  37: 218– 225. Google Scholar CrossRef Search ADS   Amico GC, Vidal-Russell R, Aizen MA, Nickrent D. 2014. Genetic diversity and population structure of the mistletoe Tristerix corymbosus (Loranthaceae). Plant Systematics and Evolution  300: 153– 162. Google Scholar CrossRef Search ADS   Anselmetti F, Ariztegui D, Hodell DA, Hillesheim MB, Brenner M, Gilli A, McKenzie JA, Mueller AD. 2006. Late quaternary climate-induced lake level variations in lake Petén Itzá, Guatemala, inferred from seismic stratigraphic analysis. Palaeogeography, Palaeoclimatology, Palaeoecology  230: 52– 69. Google Scholar CrossRef Search ADS   Beaumont MA, Zhang W, Balding DJ. 2002. Approximate Bayesian computation in population genetics. Genetics  162: 2025– 2035. Google Scholar PubMed  Braconnot P, Otto-Bliesner B, Harrison S, Joussaume S, Peterschmitt JY, Abe-Ouchi A, Crucifix M, Driesschaert E, Fichefet T, Hewitt CD, Kageyama M, Kitoh A, Loutre MF, Marti O, Merkel U, Ramstein G, Valdes P, Weber L, Yu Y, Zhao Y. 2007. Results of PMIP2 coupled simulations of the mid-holocene and last glacial maximum – part 2: feedbacks with emphasis on the location of the ITCZ and mid- and high latitudes heat budget. Climate of the Past  3: 279– 296. Google Scholar CrossRef Search ADS   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: Vimeuz F, Sylvestre F, Khodri M, eds. Past climate variability in South America and surrounding regions, from the Last Glacial Maximum to the Holocene , Vol. 14. Houten: Springer, 113– 128. Google Scholar CrossRef Search ADS   Carnaval AC, Moritz C. 2008. Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest. Journal of Biogeography  35: 1187– 1201. Google Scholar CrossRef Search ADS   Carnaval AC, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C. 2009. Stability predicts genetic diversity in the Brazilian Atlantic Forest hotspot. Science  23: 785– 789. Google Scholar CrossRef Search ADS   Cavers S, Navarro C, Lowe AJ. 2003. Chloroplast DNA phylogeography reveals colonization history of a Neotropical tree, Cedrela odorata L., in Mesoamerica. Molecular Ecology  12: 1451– 1460. Google Scholar CrossRef Search ADS PubMed  Clement M, Posada D, Crandall KA. 2000. TCS: a computer program to estimate gene genealogies. Molecular Ecology  9: 1657– 1659. Google Scholar CrossRef Search ADS PubMed  Cornuet JM, Santos F, Beaumont MA, Robert CP, Marin JM, Balding DJ, Guillemaud T, Estoup A. 2008. Inferring population history with DIYABC: a user-friendly approach to approximate Bayesian computation. Bioinformatics  24: 2713– 2719. Google Scholar CrossRef Search ADS PubMed  Cornuet JM, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, Leblois R, Marin JM, Estoup A. 2014. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics  30: 1187– 1189. Google Scholar CrossRef Search ADS PubMed  Díaz Infante S, Lara C, Arizmendi MC, Eguiarte LE, Ornelas JF. 2016. Reproductive ecology and isolation of Psittacanthus calyculatus and P. auriculatus mistletoes (Loranthaceae). PeerJ  4: e2491. Google Scholar CrossRef Search ADS PubMed  Doyle JJ, Doyle JL. 1987. A rapid DNA isolation procedure from small quantities of fresh leaf tissue. Phytochemical Bulletin, Botanical Society of America  19: 11– 15. Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology  7: 214. Google Scholar CrossRef Search ADS PubMed  Drummond AJ, Rambaut A, Shapiro B, Pybus OG. 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution  22: 1185– 1192. Google Scholar CrossRef Search ADS PubMed  Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ. 2011. A statistical explanation of MaxEnt for ecologists. Diversity and Distributions  17: 43– 57. Google Scholar CrossRef Search ADS   Espadas-Manrique C, Durán R, Argáez J. 2003. Phytogeographic analysis of taxa endemic to the Yucatán Peninsula using geographic information systems, the domain heuristic method and parsimony analysis of endemicity. Diversity and Distributions  9: 313– 330. Google Scholar CrossRef Search ADS   Excoffier L, Smouse PE, Quattro JM. 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics  131: 479– 491. Google Scholar PubMed  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. Fielding AH, Bell JF. 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation  24: 38– 49. Google Scholar CrossRef Search ADS   Fontaine MC, Austerlitz F, Giraud T, Labbé F, Papura D, Richard-Cervera S, Delmotte F. 2013. Genetic signature of a range expansion and leap-frog event after the recent invasion of Europe by the grapevine downy mildew pathogen Plasmopara viticola. Molecular Ecology  22: 2771– 2786. Google Scholar CrossRef Search ADS PubMed  Fu YX. 1997. Statistical neutrality of mutations against population growth, hitchhiking and background selection. Genetics  147: 915– 925. Google Scholar PubMed  González C, Ornelas JF, Gutiérrez-Rodríguez C. 2011. Selection and geographic isolation influence hummingbird speciation: genetic, acoustic and morphological divergence in the wedge-tailed sabrewing (Campylopterus curvipennis). BMC Evolutionary Biology  11: 38. Google Scholar CrossRef Search ADS PubMed  Grant WS. 2015. Problems and cautions with sequence mismatch analysis and bayesian skyline plots to infer historical demography. The Journal of Heredity  106: 333– 346. Google Scholar CrossRef Search ADS PubMed  Guevara-Chumacero LM, López-Wilchis R, Pedroche FF, Juste J, Ibáñez C, Barriga-Sosa IDLA. 2010. Molecular phylogeography of Pteronotus davyi (Chiroptera: Mormoopidae) in Mexico. Journal of Mammalogy  91: 220– 232. Google Scholar CrossRef Search ADS   Gutiérrez-García TA, Vázquez-Domínguez E. 2012. Biogeographically dynamic genetic structure bridging two continents in the monotypic Central American rodent Ototylomys phyllotis. Biological Journal of the Linnean Society  107: 593– 610. Google Scholar CrossRef Search ADS   Haffer J. 1969. Speciation in Amazonian forest birds. Science  165: 131– 137. Google Scholar CrossRef Search ADS PubMed  Hall TA. 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium Series  41: 95– 98. Hammer O. 2011. PAST (Paleontological Statistics) . Natural History Museum, University of Oslo. Available at: http://folk.uio.no/ohammer/past/index.html. Harpending RC. 1994. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Human Biology  66: 591– 600. Google Scholar PubMed  Haug GH, Günther D, Peterson LC, Sigman DM, Hughen KA, Aeschlimann B. 2003. Climate and the collapse of Maya civilization. Science  299: 1731– 1735. Google Scholar CrossRef Search ADS PubMed  Hey J, Nielsen R. 2004. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics  167: 747– 760. Google Scholar CrossRef Search ADS PubMed  Hey J, Nielsen R. 2007. Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proceedings of the National Academy of Sciences USA  104: 2785– 2790. Google Scholar CrossRef Search ADS   Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high resolution interpolated climate surfaces from global land areas. International Journal of Climatology  25: 1965– 1978. Google Scholar CrossRef Search ADS   Hodell DA, Anselmetti FS, Ariztegui D, Brenner M, Curtis JH, Gilli A, Grzesik DA, Guilderson TJ, Müller AD, Bush MB, Correa-Metrio A, Escobar J, Kutterolf S. 2008. An 85-ka record of climate change in lowland Central America. Quaternary Science Reviews  27: 1152– 1165. Google Scholar CrossRef Search ADS   Howell TR. 1972. Birds of the lowland pine savannah of northeastern Nicaragua. The Condor  74: 316– 340. Google Scholar CrossRef Search ADS   Islebe GA, Hooghiemstra H, Brenner M, Curtis JH, Hodell D. 1996. A Holocene vegetation history from lowland Guatemala. The Holocene  6: 265– 271. Google Scholar CrossRef Search ADS   Jiménez RA, Ornelas JF. 2016. Historical and current introgression in a Mesoamerican hummingbird species complex: a biogeographic perspective. PeerJ  4: e1556. Google Scholar CrossRef Search ADS PubMed  Kay KM, Whittall JB, Hodges SA. 2006. A survey of nuclear ribosomal internal transcribed spacer substitution rates across angiosperms: an approximate molecular clock with life history effects. BMC Evolutionary Biology  6: 36. Google Scholar CrossRef Search ADS PubMed  Kuijt J. 2009. Monograph of Psittacanthus (Loranthaceae). Systematic Botany Monographs  86: 1– 361. Lande R, Engen S, Sæther BE. 2003. Stochastic population dynamics in ecology and conservation . Oxford: Oxford University Press. Google Scholar CrossRef Search ADS   Landis JR, Koch GG. 1977. The measurement of observer agreement for categorical data. Biometrics  33: 159– 174. Google Scholar CrossRef Search ADS PubMed  Leyden BW. 1984. Guatemalan forest synthesis after Pleistocene aridity. Proceedings of the National Academy of Sciences, USA  81: 4856– 4859. Google Scholar CrossRef Search ADS   Li W, Wang Z, Ma Z, Tang H. 1997. A regression model for the spatial distribution of red-crown crane in Yancheng Biosphere Reserve, China. Ecological Modelling  103: 115– 121. Google Scholar CrossRef Search ADS   Librado P, Rozas J. 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics (Oxford, England)  25: 1451– 1452. Google Scholar CrossRef Search ADS PubMed  Lira-Noriega A, Toro-Núñez O, Oaks JR, Mort ME. 2015. The roles of history and ecology in chloroplast phylogeographic patterns of the bird-dispersed plant parasite Phoradendron californicum (Viscaceae) in the Sonoran Desert. American Journal of Botany  102: 149– 164. Google Scholar CrossRef Search ADS PubMed  Liu C, Berry PM, Dawson TP, Pearson RG. 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography  28: 385– 393. Google Scholar CrossRef Search ADS   López de Buen L, Ornelas JF. 2002. Host compatibility of the cloud forest mistletoe Psittacanthus schiedeanus (Loranthaceae) in central Veracruz, Mexico. American Journal of Botany  89: 95– 102. Google Scholar CrossRef Search ADS PubMed  López Ramos E. 1975. Geological summary of the Yucatán Peninsula. In: Nairin AEM, Stehli FG, eds. The ocean basis and margins. Volume 3. The Gulf of Mexico and the Caribbean . New York: Plenum Press, 257– 282. Lundell CL. 1934. Preliminary sketch of the phytogeography of the Yucatán Peninsula. Contributions to American Archaeology  12: 257– 321. 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. Google Scholar CrossRef Search ADS   Morrone JJ. 2005. Hacia una síntesis biogeográfica de México. Revista Mexicana de Biodiversidad  76: 207– 252. Morrone JJ. 2014. Biogeographical regionalisation of the Neotropical region. Zootaxa  3782: 1– 110. Google Scholar CrossRef Search ADS PubMed  Mueller AD, Anselmetii FS, Ariztequi D, Brenner M, Hodell DA, Curtis JH, Escobar J, Gilli A, Grzesik DA, Guilderson TP, Kutterolf S, Plötze M. 2010. Late Quaternary palaeoenvironment of northern Guatemala: evidence from deep drill cores and seismic statigraphy of Lake Petén Itzá. Sedimentology  57: 1220– 1245. Nyagumbo E, Ndagurwa HGT, Mushonga K. 2017. High genetic variation and low gene flow among populations of Viscum verrucosum in semi-arid savannah, southwest Zimbabwe. Ecological Genetics and Genomics 3–5 : 41– 46. 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. Ornelas JF, Gándara E, Vásquez-Aguilar AA, Ramírez-Barahona S, Ortiz-Rodriguez AE, González C, Mejía S aules MT, Ruiz-Sanchez E. 2016. A mistletoe tale: postglacial invasion of Psittacanthus schiedeanus (Loranthaceae) to Mesoamerican cloud forests revealed by molecular data and species distribution modeling. BMC Evolutionary Biology  16: 78. Google Scholar CrossRef Search ADS PubMed  Oyama K, Martínez-Ramos M, Peñaloza-Ramírez JM, Rocha-Ramírez V, Armenta-Medina EG, Hernández-Soto P. 2016. Population genetic structure of an extremely logged tree species Guaiacum sanctum L. in the Yucatán Peninsula, Mexico. Botanical Sciences  94: 345– 356. Google Scholar CrossRef Search ADS   Padilla y Sánchez RJ. 2007. Evolución geológica del sureste mexicano desde el Mesozoico al presente en el context regional del Golfo de México. Boletín de la Sociedad Geológica Mexicana  59: 19– 42. Google Scholar CrossRef Search ADS   Pérez-Crespo MJ, Ornelas JF, González-Rodríguez A, Ruiz-Sanchez E, Vásquez-Aguilar AA, Ramírez-Barahona S. 2017. Phylogeography and population differentiation in the Psittacanthus calyculatus (Loranthaceae) mistletoe: a complex scenario of climate-volcanism interaction along the Trans-Mexican Volcanic Belt. Journal of Biogeography , 44: 2501– 2514. Google Scholar CrossRef Search ADS   Pfenninger M, Posada D. 2002. Phylogeographic history of the land snail Candidula unifasciata (Helicellinae, Stylommatophora): fragmentation, corridor migration, and secondary contact. Evolution; International Journal of Organic Evolution  56: 1776– 1788. Google Scholar CrossRef Search ADS PubMed  Phillips SJ, Anderson RP, Schapire RE. 2006. Maximum entropy modelling of species geographic distributions. Ecological Modelling  190: 231– 259. Google Scholar CrossRef Search ADS   Poelchau MF, Hamrick JL. 2012. Differential effects of landscape-level environmental features on genetic structure in three codistributed tree species in Central America. Molecular Ecology  21: 4970– 4982. Google Scholar CrossRef Search ADS PubMed  Poelchau MF, Hamrick JL. 2013a. Comparative phylogeography of three common Neotropical tree species. Journal of Biogeography  40: 618– 631. Google Scholar CrossRef Search ADS   Poelchau MF, Hamrick JL. 2013b. Palaeodistribution modelling does not support disjunct Pleistocene refugia in several Central American plant taxa. Journal of Biogeography  40: 662– 675. Google Scholar CrossRef Search ADS   Pons O, Petit RJ. 1996. Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics  144: 1237– 1245. Google Scholar PubMed  Posada D. 2008. jModelTest: phylogenetic model averaging. Molecular Biology and Evolution  25: 1253– 1256. Google Scholar CrossRef Search ADS PubMed  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 and Evolution  3: 725– 738. Google Scholar CrossRef Search ADS PubMed  Ramírez-Barahona S, González C, González-Rodríguez A, Ornelas JF. 2017. The influence of climatic niche preferences on the population genetic structure of a mistletoe species complex. New Phytologist  214: 1751– 1761. Google Scholar CrossRef Search ADS PubMed  Ramírez-Barahona S, Torres-Miranda A, Palacios-Ríos M, Luna-Vega I. 2009. Historical biogeography of the Yucatán Peninsula, Mexico: a perspective from ferns (Monilophyta) and lycopods (Lycophyta). Biological Journal of the Linnean Society  98: 775– 786. Google Scholar CrossRef Search ADS   Richardson JE, Pennington RT, Pennington TD, Hollingsworth PM. 2001. Rapid diversification of a species-rich genus of Neotropical rain forest trees. Science (New York, N.Y.)  293: 2242– 2245. Google Scholar CrossRef Search ADS PubMed  Robert CP, Cornuet JM, Marin JM, Pillai NS. 2011. Lack of confidence in approximate Bayesian computation model choice. Proceedings of the National Academy of Sciences USA  108: 15112– 15117. Google Scholar CrossRef Search ADS   Rodríguez-Gómez F, Ornelas JF. 2014. Genetic divergence of the Mesoamerican azure-crowned hummingbird (Amazilia cyanocephala) across the Motagua–Polochic–Jocotán fault system. Journal of Zoological Systematics and Evolutionary Research  52: 142– 153. Google Scholar CrossRef Search ADS   Rogers AR, Harpending H. 1992. Population growth makes waves in the distribution of pairwise genetic differences. Molecular Biology and Evolution  9: 552– 569. Google Scholar PubMed  Simpson GG. 1964. Species density of North American recent mammals. Systematic Zoology  12: 57– 73. Google Scholar CrossRef Search ADS   Smith BT, Escalante P, Hernández B años BE, Navarro-Sigüenza AG, Rohwer S, Klicka J. 2011. The role of historical and contemporary processes on phylogeographic structure and genetic diversity in the Northern Cardinal, Cardinalis cardinalis. BMC Evolutionary Biology  11: 136. Google Scholar CrossRef Search ADS PubMed  Schneider S, Excoffier L. 1999. Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA. Genetics  152: 1079– 1089. Google Scholar PubMed  Taberlet P, Gielly L, Pautou G, Bouvet J. 1991. Universal primers for amplification of three non-coding regions of chloroplast DNA. Plant Molecular Biology  17: 1105– 1109. Google Scholar CrossRef Search ADS PubMed  Tajima F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics  123: 585– 595. Google Scholar PubMed  Vázquez-Domínguez E, Arita HT. 2010. The Yucatán peninsula: biogeographical history 65 million years in the making. Ecography  33: 212– 219. Torrescano-Valle N, Islebe GA. 2015. Holocene paleoecology, climate history and human influence in the southwestern Yucatán Peninsula. Review of Paleobotany and Palynology  217: 1– 8. Google Scholar CrossRef Search ADS   Vidal-Russell R, Nickrent DL. 2007. A molecular phylogeny of the feathery mistletoe Misodendrum. Systematic Botany  32: 560– 568. Google Scholar CrossRef Search ADS   Vidal-Russell R, Nickrent DL. 2008a. Evolutionary relationships in the showy mistletoe family (Loranthaceae). American Journal of Botany  95: 1015– 1029. Google Scholar CrossRef Search ADS   Vidal-Russell R, Nickrent DL. 2008b. The first mistletoes: origins of aerial parasitism in Santalales. Molecular Phylogenetics and Evolution  47: 523– 537. Google Scholar CrossRef Search ADS   Williford D, Deyoung RW, Honeycutt RL, Brennan LA, Hernández F. 2016. Phylogeography of the bobwhite (Colinus) quails. Wildlife Monographs  193: 1– 49. Google Scholar CrossRef Search ADS   Wilson CA, Calvin CL. 2006. An origin of aerial branch parasitism in the mistletoe family, Loranthaceae. American Journal of Botany  93: 787– 796. Google Scholar CrossRef Search ADS PubMed  Zuber D, Widmer A. 2009. Phylogeography and host race differentiation in the European mistletoe (Viscum album L.). Molecular Ecology  18: 1946– 1962. Google Scholar CrossRef Search ADS PubMed  APPENDICES Appendix 1. GenBank accession numbers of DNA sequences of other Psittacanthus species and representatives of Loranthaceae from Wilson & Calvin (2006), Amico, Vidal-Russell & Nickrent (2007), Vidal-Russell & Nickrent (2007, 2008a, b), Díaz Infante et al. (2016), Ornelas et al. (2016) and Pérez-Crespo et al. (2017) used in this study. Species  Voucher  Country  ITS  trnL-F  Alepis flavida Tiegh.  Calvin and Wilson NZ98- 04 (NA)  New Zealand  DQ333847  DQ340598  Atkinsonia ligustrina F.Muell.  Calvin and Wilson AU00-01 (NA)  Australia  DQ333865  DQ788714  Cladocolea cupulata Kuijt  Calvin and Wilson MX03- 08 (NA)  Mexico  DQ333861  DQ340612  Cladocolea mcvaughii Kuijt  Calvin and Wilson MX03- 09 (NA)  Mexico  DQ333860  DQ340611  Desmaria mutabilis (Poepp. & Endl.) Tiegh. ex B.D.Jacks.  Calvin and Wilson CL03-07 (NA)  Chile  DQ333852  DQ340603  Gaiadendron punctatum (Ruiz & Pav.) G.Don.  Calvin and Wilson CR01- 08 (NA)  Costa Rica  DQ333866  DQ340617  Ligaria cuneifolia (Ruiz & Pav.) Tiegh.  Calvin and Wilson CL03-01 (NA)  Chile  DQ333853  DQ340604  Notanthera heterophylla (Ruiz & Pav.) G.Don.  Calvin and Wilson Cl03-03 (NA)  Chile  DQ333855  DQ340606  Nuytsia floribunda R.Br.  Nickrent 2747 (NA)  Australia  DQ788705  -  Nuytsia floribunda R.Br.  Nickrent 3080 (NA)  Australia  -  DQ788716  Nuytsia floribunda R.Br.  Calvin and Wilson AU01-22 (NA)  Australia  DQ333867  DQ340618  Oryctanthus occidentalis (L.) Eichler  Calvin and Wilson CR01- 11 (NA)  Costa Rica  DQ333862  DQ340613  Peraxilla tetrapetala (L.f.) Tiegh.  Calvin and Wilson NZ98- 03 (NA)  New Zealand  DQ333846  DQ340597  Phragmanthera regularis (Steud. ex Sprague) M.G.Gilbert  Calvin and Wilson Y88-01 (NA)  NA  DQ333830  DQ340579  Phthirusa pyrifolia (Kunth) Eichler  Calvin and Wilson CR01- 03 (NA)  Costa Rica  DQ333857  -  Phthirusa pyrifolia (Kunth) Eichler  Nickrent 2762 (NA)  NA  -  EU544504  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670A (USP)  Brazil  KU923005  KU923279  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670B (USP)  Brazil  KU923006  KU923280  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 004 (XAL)  Mexico  KU923007  KU923281  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 005 (XAL)  Mexico  KU923008  KU923282  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez (NA)  Mexico  KU923009  KU923283  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3785 (USP)  Brazil  KU923010  KU923284  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3786 (USP)  Brazil  KU923011  KU923285  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3671 (USP)  Brazil  KU923019  KU923293  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3672 (USP)  Brazil  KU923020  KU923294  Psittacanthus macrantherus Eichler  E. Ruiz-Sanchez 348 (XAL)  Mexico  KU923021  KU923295  Psittacanthus macrantherus Eichler  F. Rodriguez NA (XAL)  Mexico  KU923022  KU923296  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923029  KU923303  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923030  KU923304  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923031  KU923305  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 013 (XAL)  Mexico  KU923032  KU923306  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 014 (XAL)  Mexico  KU923033  KU923307  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 015 (XAL  Mexico  KU923034  KU923308  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 016 (XAL)  Mexico  KU923035  KU923309  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923036  KU923310  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923037  KU923311  Psittacanthus rhynchanthus (Benth.) Kuijt  P. Carrillo Reyes 5372 (XAL)  Guatemala  KU923038  KU923312  Psittacanthus rhynchanthus (Benth.) Kuijt  E. Ruiz-Sanchez 417 (XAL)  Mexico  KU923040  KU923314  Psittacanthus rhynchanthus (Benth.) Kuijt  T. Mejia-Saules 2048 (XAL)  Mexico  KU923041  KU923315  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3588 (USP)  Brazil  KU923042  KU923316  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3589 (USP)  Brazil  KU923043  KU923317  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3596 (USP)  Brazil  KU923044  KU923318  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923045  KU923319  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923046  KU923320  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923047  KU923321  Struthanthus orbicularis (Kunth) Blume  Calvin and Wilson CR01- 02 (NA)  Costa Rica  DQ333856  DQ340607  Tripodanthus acutifolius (Ruiz & Pav.) Tiegh.  Calvin and Wilson AR03-04 (NA)  Argentina  DQ333864  DQ340615  Tristerix aphyllus (Miers ex DC.) Barlow & Wiens  G. Amico 166 (BRCU)  Chile  DQ442966  DQ442919  Tristerix corymbosus (L.) Kuijt  Calvin and Wilson CL03-02 (NA)  Chile  DQ333854  DQ340605  Tupeia antarctica (G.Forst.) Cham. & Schltdl.  Calvin and Wilson NZ98- 02 (NA)  New Zealand  DQ333850  DQ340601  Species  Voucher  Country  ITS  trnL-F  Alepis flavida Tiegh.  Calvin and Wilson NZ98- 04 (NA)  New Zealand  DQ333847  DQ340598  Atkinsonia ligustrina F.Muell.  Calvin and Wilson AU00-01 (NA)  Australia  DQ333865  DQ788714  Cladocolea cupulata Kuijt  Calvin and Wilson MX03- 08 (NA)  Mexico  DQ333861  DQ340612  Cladocolea mcvaughii Kuijt  Calvin and Wilson MX03- 09 (NA)  Mexico  DQ333860  DQ340611  Desmaria mutabilis (Poepp. & Endl.) Tiegh. ex B.D.Jacks.  Calvin and Wilson CL03-07 (NA)  Chile  DQ333852  DQ340603  Gaiadendron punctatum (Ruiz & Pav.) G.Don.  Calvin and Wilson CR01- 08 (NA)  Costa Rica  DQ333866  DQ340617  Ligaria cuneifolia (Ruiz & Pav.) Tiegh.  Calvin and Wilson CL03-01 (NA)  Chile  DQ333853  DQ340604  Notanthera heterophylla (Ruiz & Pav.) G.Don.  Calvin and Wilson Cl03-03 (NA)  Chile  DQ333855  DQ340606  Nuytsia floribunda R.Br.  Nickrent 2747 (NA)  Australia  DQ788705  -  Nuytsia floribunda R.Br.  Nickrent 3080 (NA)  Australia  -  DQ788716  Nuytsia floribunda R.Br.  Calvin and Wilson AU01-22 (NA)  Australia  DQ333867  DQ340618  Oryctanthus occidentalis (L.) Eichler  Calvin and Wilson CR01- 11 (NA)  Costa Rica  DQ333862  DQ340613  Peraxilla tetrapetala (L.f.) Tiegh.  Calvin and Wilson NZ98- 03 (NA)  New Zealand  DQ333846  DQ340597  Phragmanthera regularis (Steud. ex Sprague) M.G.Gilbert  Calvin and Wilson Y88-01 (NA)  NA  DQ333830  DQ340579  Phthirusa pyrifolia (Kunth) Eichler  Calvin and Wilson CR01- 03 (NA)  Costa Rica  DQ333857  -  Phthirusa pyrifolia (Kunth) Eichler  Nickrent 2762 (NA)  NA  -  EU544504  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670A (USP)  Brazil  KU923005  KU923279  Psittacanthus acinarius (Mart.) Mart.  G. Ceccantini 3670B (USP)  Brazil  KU923006  KU923280  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 004 (XAL)  Mexico  KU923007  KU923281  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez 005 (XAL)  Mexico  KU923008  KU923282  Psittacanthus auriculatus (Oliv.) Eichler  M. J. Perez (NA)  Mexico  KU923009  KU923283  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3785 (USP)  Brazil  KU923010  KU923284  Psittacanthus biternatus (Hoffmanns.) G.Don  G. Ceccantini 3786 (USP)  Brazil  KU923011  KU923285  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3671 (USP)  Brazil  KU923019  KU923293  Psittacanthus cordatus (Hoffmanns.) G.Don  G. Ceccantini 3672 (USP)  Brazil  KU923020  KU923294  Psittacanthus macrantherus Eichler  E. Ruiz-Sanchez 348 (XAL)  Mexico  KU923021  KU923295  Psittacanthus macrantherus Eichler  F. Rodriguez NA (XAL)  Mexico  KU923022  KU923296  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923029  KU923303  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923030  KU923304  Psittacanthus palmeri (S.Watson) Barlow & Wiens  C. Soberanes (NA)  Mexico  KU923031  KU923305  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 013 (XAL)  Mexico  KU923032  KU923306  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 014 (XAL)  Mexico  KU923033  KU923307  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 015 (XAL  Mexico  KU923034  KU923308  Psittacanthus ramiflorus (DC.) G.Don  Y. Licona Vera 016 (XAL)  Mexico  KU923035  KU923309  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923036  KU923310  Psittacanthus ramiflorus (DC.) G.Don  NA (NA)  Panama  KU923037  KU923311  Psittacanthus rhynchanthus (Benth.) Kuijt  P. Carrillo Reyes 5372 (XAL)  Guatemala  KU923038  KU923312  Psittacanthus rhynchanthus (Benth.) Kuijt  E. Ruiz-Sanchez 417 (XAL)  Mexico  KU923040  KU923314  Psittacanthus rhynchanthus (Benth.) Kuijt  T. Mejia-Saules 2048 (XAL)  Mexico  KU923041  KU923315  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3588 (USP)  Brazil  KU923042  KU923316  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3589 (USP)  Brazil  KU923043  KU923317  Psittacanthus robustus (Mart.) Mart.  G. Ceccantini 3596 (USP)  Brazil  KU923044  KU923318  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923045  KU923319  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923046  KU923320  Psittacanthus sonorae (S.Watson) Kuijt  T. S. Rodríguez (NA)  Mexico  KU923047  KU923321  Struthanthus orbicularis (Kunth) Blume  Calvin and Wilson CR01- 02 (NA)  Costa Rica  DQ333856  DQ340607  Tripodanthus acutifolius (Ruiz & Pav.) Tiegh.  Calvin and Wilson AR03-04 (NA)  Argentina  DQ333864  DQ340615  Tristerix aphyllus (Miers ex DC.) Barlow & Wiens  G. Amico 166 (BRCU)  Chile  DQ442966  DQ442919  Tristerix corymbosus (L.) Kuijt  Calvin and Wilson CL03-02 (NA)  Chile  DQ333854  DQ340605  Tupeia antarctica (G.Forst.) Cham. & Schltdl.  Calvin and Wilson NZ98- 02 (NA)  New Zealand  DQ333850  DQ340601  NA = Not available. View Large Appendix 2. Posterior parameter estimates for the best-supported scenario (scenario 1) considering the three Psittacanthus mayanus groups (dryCHI, humYUC, dryYUC). Parameter  mean  median  mode  q025  q050  q250  q750  q950  q975  N1  6.02e+02  6.31e+02  6.87e+02  1.92e+02  2.64e+02  4.92e+02  7.41e+02  8.37e+02  8.56e+02  N2b  3.99e+03  4.13e+03  4.84e+03  2.24e+03  2.57e+03  3.53e+03  4.57e+03  4.91e+03  4.95e+03  N3b  4.36e+03  4.49e+03  4.79e+03  3.00e+03  3.28e+03  4.10e+03  4.76e+03  4.95e+03  4.97e+03  N2a2  4.33e+03  4.52e+03  4.97e+03  2.61e+03  3.00e+03  4.05e+03  4.80e+03  4.96e+03  4.98e+03  N2b2  2.65e+03  2.65e+03  3.04e+03  1.10e+03  1.19e+03  1.89e+03  3.39e+03  4.16e+03  4.35e+03  Na  7.72e+02  8.09e+02  8.94e+02  4.50e+02  5.22e+02  7.20e+02  8.62e+02  8.93e+02  8.97e+02  t1  1.51e+03  1.24e+03  1.08e+03  1.02e+03  1.03e+03  1.11e+03  1.53e+03  3.02e+03  4.15e+03  t2  1.56e+04  1.57e+04  2.07e+04  1.03e+04  1.06e+04  1.29e+04  1.84e+04  2.05e+04  2.08e+04  t3  4.39e+04  4.19e+04  4.04e+04  4.01e+04  4.02e+04  4.07e+04  4.45e+04  5.51e+04  6.18e+04  μnuDNA  2.03e-07  1.92e-07  1.80e-07  1.08e-07  1.19e-07  1.57e-07  2.38e-07  3.24e-07  3.56e-07  μcpDNA  3.25e-07  3.07e-07  2.65e-07  1.06e-07  1.38e-07  2.34e-07  3.96e-07  5.71e-07  6.39e-07  Parameter  mean  median  mode  q025  q050  q250  q750  q950  q975  N1  6.02e+02  6.31e+02  6.87e+02  1.92e+02  2.64e+02  4.92e+02  7.41e+02  8.37e+02  8.56e+02  N2b  3.99e+03  4.13e+03  4.84e+03  2.24e+03  2.57e+03  3.53e+03  4.57e+03  4.91e+03  4.95e+03  N3b  4.36e+03  4.49e+03  4.79e+03  3.00e+03  3.28e+03  4.10e+03  4.76e+03  4.95e+03  4.97e+03  N2a2  4.33e+03  4.52e+03  4.97e+03  2.61e+03  3.00e+03  4.05e+03  4.80e+03  4.96e+03  4.98e+03  N2b2  2.65e+03  2.65e+03  3.04e+03  1.10e+03  1.19e+03  1.89e+03  3.39e+03  4.16e+03  4.35e+03  Na  7.72e+02  8.09e+02  8.94e+02  4.50e+02  5.22e+02  7.20e+02  8.62e+02  8.93e+02  8.97e+02  t1  1.51e+03  1.24e+03  1.08e+03  1.02e+03  1.03e+03  1.11e+03  1.53e+03  3.02e+03  4.15e+03  t2  1.56e+04  1.57e+04  2.07e+04  1.03e+04  1.06e+04  1.29e+04  1.84e+04  2.05e+04  2.08e+04  t3  4.39e+04  4.19e+04  4.04e+04  4.01e+04  4.02e+04  4.07e+04  4.45e+04  5.51e+04  6.18e+04  μnuDNA  2.03e-07  1.92e-07  1.80e-07  1.08e-07  1.19e-07  1.57e-07  2.38e-07  3.24e-07  3.56e-07  μcpDNA  3.25e-07  3.07e-07  2.65e-07  1.06e-07  1.38e-07  2.34e-07  3.96e-07  5.71e-07  6.39e-07  Estimates are based on 1% of simulated datasets closest to the observed values. Simulations and approximate Bayesian computation (ABC) analyses were performed considering both nrDNA and plastid DNA sequences. Parameters are N = effective population size for dryCHI (pop1, N1), humYUC (Pop2, N2), dryYUC (Pop3, N3), and the ancestral populations (Na), t = time since divergence (years), and μ = mutation rate for nuclear and plastid DNA, respectively. Group abbreviations are as follows: dryYUC = dry Yucatán, humYUC = humid Yucatán, dryCHI = dry Central Depression of Chiapas. View Large Appendix 3. Factor loadings from the principal components analysis of Psittacanthus mayanus in Mexico on temperature and precipitation variables from WorldClim. Variables  Factor loading    PC1  PC2  PC3  PC4  Annual mean temperature (BIO1)  0.809  0.558  0.051  0.084  Mean diurnal range (BIO2)  –0.621  0.234  0.675  0.320  Isothermality (BIO3)  –0.645  0.057  0.086  0.651  Temperature seasonality (BIO4)  0.427  0.476  0.467  –0.300  Max. temperature of warmest month (BIO5)  0.310  0.653  0.674  0.102  Min. temperature of coldest month (BIO6)  0.869  0.218  –0.391  –0.096  Temperature annual range (BIO7)  –0.506  0.271  0.791  0.150  Mean temperature of wettest quarter (BIO8)  0.668  0.689  –0.053  0.070  Mean temperature of driest quarter (BIO9)  0.826  0.434  0.040  –0.141  Mean temperature of warmest quarter (BIO10)  0.799  0.574  0.140  0.012  Mean temperature of coldest quarter (BIO11)  0.786  0.476  –0.099  0.187  Annual precipitation (BIO12)  0.640  –0.645  0.388  –0.044  Precipitation of wettest month (BIO13)  0.508  –0.611  0.524  –0.216  Precipitation of driest month (BIO14)  0.777  –0.414  0.155  0.358  Precipitation seasonality (BIO15)  –0.765  0.060  0.292  –0.393  Precipitation of wettest quarter (BIO16)  0.396  –0.589  0.606  –0.237  Precipitation of driest quarter (BIO17)  0.825  –0.397  0.067  0.289  Precipitation of warmest quarter (BIO18)  0.459  –0.649  –0.065  0.314  Precipitation of coldest quarter (BIO19)  0.799  –0.538  0.071  –0.011  Proportion of variance explained  45.66%  23.91%  15.08%  6.82%  ANOVA  F2,103 = 77.8 P < 0.0001  F2,103 = 39.9 P < 0.0001  F2,103 = 4.2 P = 0.017  F2,103 = 5.3 P = 0.006  Variables  Factor loading    PC1  PC2  PC3  PC4  Annual mean temperature (BIO1)  0.809  0.558  0.051  0.084  Mean diurnal range (BIO2)  –0.621  0.234  0.675  0.320  Isothermality (BIO3)  –0.645  0.057  0.086  0.651  Temperature seasonality (BIO4)  0.427  0.476  0.467  –0.300  Max. temperature of warmest month (BIO5)  0.310  0.653  0.674  0.102  Min. temperature of coldest month (BIO6)  0.869  0.218  –0.391  –0.096  Temperature annual range (BIO7)  –0.506  0.271  0.791  0.150  Mean temperature of wettest quarter (BIO8)  0.668  0.689  –0.053  0.070  Mean temperature of driest quarter (BIO9)  0.826  0.434  0.040  –0.141  Mean temperature of warmest quarter (BIO10)  0.799  0.574  0.140  0.012  Mean temperature of coldest quarter (BIO11)  0.786  0.476  –0.099  0.187  Annual precipitation (BIO12)  0.640  –0.645  0.388  –0.044  Precipitation of wettest month (BIO13)  0.508  –0.611  0.524  –0.216  Precipitation of driest month (BIO14)  0.777  –0.414  0.155  0.358  Precipitation seasonality (BIO15)  –0.765  0.060  0.292  –0.393  Precipitation of wettest quarter (BIO16)  0.396  –0.589  0.606  –0.237  Precipitation of driest quarter (BIO17)  0.825  –0.397  0.067  0.289  Precipitation of warmest quarter (BIO18)  0.459  –0.649  –0.065  0.314  Precipitation of coldest quarter (BIO19)  0.799  –0.538  0.071  –0.011  Proportion of variance explained  45.66%  23.91%  15.08%  6.82%  ANOVA  F2,103 = 77.8 P < 0.0001  F2,103 = 39.9 P < 0.0001  F2,103 = 4.2 P = 0.017  F2,103 = 5.3 P = 0.006  View Large © 2018 The Linnean Society of London, Botanical Journal of the Linnean Society

Journal

Botanical Journal of the Linnean SocietyOxford University Press

Published: Mar 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off