A mistletoe tale: postglacial invasion of Psittacanthus schiedeanus (Loranthaceae) to Mesoamerican cloud forests revealed by molecular data and species distribution modeling

Abstract

Background

Ecological adaptation to host taxa is thought to result in mistletoe speciation via race formation. However, historical and ecological factors could also contribute to explain genetic structuring particularly when mistletoe host races are distributed allopatrically. Using sequence data from nuclear (ITS) and chloroplast (trnL-F) DNA, we investigate the genetic differentiation of 31 Psittacanthus schiedeanus (Loranthaceae) populations across the Mesoamerican species range. We conducted phylogenetic, population and spatial genetic analyses on 274 individuals of P. schiedeanus to gain insight of the evolutionary history of these populations. Species distribution modeling, isolation with migration and Bayesian inference methods were used to infer the evolutionary transition of mistletoe invasion, in which evolutionary scenarios were compared through posterior probabilities.

Results

Our analyses revealed shallow levels of population structure with three genetic groups present across the sample area. Nine haplotypes were identified after sequencing the trnL-F intergenic spacer. These haplotypes showed phylogeographic structure, with three groups with restricted gene flow corresponding to the distribution of individuals/populations separated by habitat (cloud forest localities from San Luis Potosí to northwestern Oaxaca and Chiapas, localities with xeric vegetation in central Oaxaca, and localities with tropical deciduous forests in Chiapas), with post-glacial population expansions and potentially corresponding to post-glacial invasion types. Similarly, 44 ITS ribotypes suggest phylogeographic structure, despite the fact that most frequent ribotypes are widespread indicating effective nuclear gene flow via pollen. Gene flow estimates, a significant genetic signal of demographic expansion, and range shifts under past climatic conditions predicted by species distribution modeling suggest post-glacial invasion of P. schiedeanus mistletoes to cloud forests. However, Approximate Bayesian Computation (ABC) analyses strongly supported a scenario of simultaneous divergence among the three groups isolated recently.

Conclusions

Our results provide support for the predominant role of isolation and environmental factors in driving genetic differentiation of Mesoamerican parrot-flower mistletoes. The ABC results are consistent with a scenario of post-glacial mistletoe invasion, independent of host identity, and that habitat types recently isolated P. schiedeanus populations, accumulating slight phenotypic differences among genetic groups due to recent migration across habitats. Under this scenario, climatic fluctuations throughout the Pleistocene would have altered the distribution of suitable habitat for mistletoes throughout Mesoamerica leading to variation in population continuity and isolation. Our findings add to an understanding of the role of recent isolation and colonization in shaping cloud forest communities in the region.

Background

Most aerial Loranthaceae mistletoes use bird vectors for seed dispersal and, eventually, branch attachment and penetration of woody host plants [1]. Given their close affinities with their dispersal vectors and the connections with and dependence on their hosts [2, 3], the mistletoe geographical ranges are directly related to the morphological and behavioral adaptations of mistletoe fruit specialists and the availability of suitable host trees. However, most avian vectors of mistletoe parasites do not seem to be sufficiently specialized to isolate mistletoe populations on different hosts ([4] and references therein) and, therefore, the genetic structuring of mistletoe populations is likely more influenced by the geographic structuring of host populations [5–7].

Research on mistletoe evolution has highlighted the importance of host specialization during speciation (reviewed in [5]), which can occur in three ways: (1) mistletoe co-speciation with their hosts (dispersal hypothesis), (2) mistletoe speciation by changing host specificity (host-switching hypothesis), and (3) mistletoe speciation without changes in the host (ecological hypothesis). Each of these modes of mistletoe speciation could occur following the allopatric, peripheral isolates or sympatric models of speciation [8]. Psittacanthus Mart. (ca. 120 species; [3, 9]) mistletoes show a variety of patterns of host specificity from specialist to generalist [10]. For example, three species with allopatric distribution in Mesoamerica, P. angustifolius (Chiapas to Nicaragua), P. macrantherus (Trans-Mexican Volcanic Belt and Sierra Madre Occidental) and P. pinicola (Atlantic coast from Belize to Nicaragua), parasitize exclusively species of Pinus and Abies [11, 12], while P. sonorae (Sonora and Baja California Sur, Mexico), P. nudus (Yeguare Valley, Honduras) and P. palmeri (along the Trans-Mexican Volcanic Belt) parasitize almost exclusively species of Bursera, and P. auriculatus (Tehuacan-Cuicatlán Valley, Mexico) and P. breedlovei (Chiapas Central Depression, Mexico) parasitize Acacia species and other legumes [3, 13]. In contrast, other Psittacanthus species widely distributed in Mesoamerica often in sympatry parasitize distantly related host species (e.g., P. schiedeanus, P. mayanus, P. calyculatus, P. ramiflorus; [3, 14, 15]). That is, there is a tendency among Psittacanthus mistletoes that infect more than one host species to infect distantly related hosts and their ranges overlap, whereas Psittacanthus species that infect closely related host species have allopatric distributions.

The tendency for closely related mistletoes to infect closely related host species supports the model of mistletoe co-speciation with their hosts (dispersal hypothesis). In contrast, host specialization might play an important role in speciation to occur among species of mistletoes using several host species, and mistletoe speciation in sympatry would occur by changing host specificity (host-switching hypothesis) leading to specialization on different hosts and the formation of mistletoe races. If a mistletoe species utilizes two or more host species, the genetic pool of the mistletoe species may become differentiated into “host races,” which may eventually diverge into new species (i.e., follow a model of speciation by host race formation).

The existence of host-specific relationships through cross-dispersal experiments has been suggested as evidence for the presence of host-specific races in a number of mistletoe species (e.g., [15–21]), or by using genetic markers, including isozymes [22–26], AFLPs [27] and chloroplast DNA sequences [6, 7, 28, 29]. Cross-dispersal experiments have shown that seedling development is greatest when mistletoe seeds are placed on their source host species [15–17, 20, 21, 30, 31]. Using genetic data, however, the distinct host races identified were geographically isolated by some distance (except [25, 26]) or the phylogeographic groups were best explained by isolation linked to Pliocene geological events and Pleistocene climate changes [6, 7], which makes it difficult to distinguish between the potential effects of isolation and those linked to changes in host-tree availability on population divergence.

A phylogeographic approach is appropriate for evolutionary studies of mistletoes, because the assessment of geographic structure and genetic variation within and among populations might reflect underlying spatial and/or temporal discontinuities in available host tree species. Despite the appeal of the speciation model of host race formation by host switching, species formation in mistletoes that are transmitted between hosts by avian vectors is not understood. Speciation via race formation is thought to be the result of ecological adaptation to host taxa, and presumably the evolutionary pathway that eventually lead to mistletoe diversification [5]. However, other factors could also contribute to explain genetic structuring of mistletoe populations particularly when mistletoe “host races” are distributed allopatrically. It may occur when gene flow between conspecific populations is diminished because of the synergistic effects of both paleogeographical events and Pleistocene glacial-interglacial cycles that, coupled with the complex interactions between mistletoes and their tree hosts, pollinators and seed dispersers, could have shaped at different time scales the distribution and genetic structure of species in a region (e.g., [32]).

Here the evolutionary history of Psittacanthus schiedeanus (Schltdl. & Cham.) G. Don (Loranthaceae) is examined through species distribution modeling and phylogeographic and population genetic analyses of nuclear and chloroplast DNA. We used a broad geographical sampling of populations of P. schiedeanus DNA sequences from the internal transcribed spacer (ITS) region of the nuclear ribosomal DNA (nrDNA) in combination with sequences of the rapidly evolving trnL-F intergenic spacer region from the chloroplast genome (cpDNA). This multilocus approach allowed us to explicitly test the set of historical and climatic scenarios in the evolution of P. schiedeanus and identify patterns of population structure that may reflect the evolutionary history of the species. As the chloroplast genome is maternally inherited in most angiosperms, cpDNA sequences may provide a picture of evolutionary haplotype relationships of seed-mediated gene flow compared to the ITS sequences reflecting a combination of both pollen and seed movement [33, 34]. Specifically, we address the following questions: (1) Does the phylogeography of P. schiedeanus based on nuclear and chloroplast markers support a vicariance model in which major geographic disjunctions (in particular the Isthmus of Tehuantepec break) or historical fragmentation of cloud forest distribution correspond with genetic diversification? (2) Are the phylogeographic patterns concordant between the two types of markers? (3) Does niche modeling also support the vicariance model, and does the proposed patterns of suitable habitat distribution though Pleistocene glacial-interglacial cycles match those suggested by phylogeography? And (4), Does the phylogeography of P. schiedeanus support the model of ecological speciation due to host-specificity, and does this reflect the observed host population expansion/contraction events historically related to the availability of suitable host trees?

Methods

Study system

Psittacanthus schiedeanus mistletoes have yellow-to-orange, self-compatible hermaphroditic flowers pollinated mainly by hummingbirds [35], and ripe purplish-black, lipid-rich, fleshy (one-seed) fruits dispersed by a variety of birds [14, 20, 36, 37]. These mistletoes are characteristic of the canopy in the cloud forests edges in northern Mesoamerica, here delineated from northeastern Mexico to the Guatemalan highlands [38]. The northern boundary of the region comprises cloud forests along the Sierra Madre Oriental, from southern Tamaulipas to northern Oaxaca. The southern boundary includes cloud forests in the Sierra Madre de Chiapas, central highlands of Chiapas, and in the Sierra de Las Minas and Sierra de Cuchumatanes in Guatemala (Additional file 1). Between these boundaries, cloud forests encompass an extremely heterogeneous mixture of North American temperate tree species that are present in the region since the Tertiary and tropical species with South American origins during the early Miocene [32, 39–42].

Psittacanthus schiedeanus often parasitizes more than 20 host tree species, introduced to or native to cloud forests [43–45]. In central Veracruz, the most severe infections occur on Liquidambar styraciflua [14, 20, 45], a temperate tree species present in Mesoamerica as early as the Late Miocene [42, 46]. A phylogeographic study in L. styraciflua revealed breaks between populations in the deciduous forests of the southeastern USA and the Mesoamerican cloud forests, populations distributed north and south of the Trans-Mexican Volcanic Belt (TMVB), and private haplotypes for populations in the cloud forest of the Sierra de Los Tuxtlas [46]. The genetic divergence between L. styraciflua populations distributed north and south of the TMVB occurred 4.27–1.42 million years ago (MYA) during the formation of the Pliocene to Quaternary volcanic arc, one of the most recent episodes of the geologic evolution of the TMVB (at ca. 3.6 MYA), and subsequent climatic changes in the region. Ruiz-Sanchez and Ornelas [46] also found signals of demographic expansion during the Last Glacial Maximum (LGM, ca. 20,000 years ago) in the Mesoamerican populations of L. styraciflua that presumably connected and expanded to lower coastal areas. Given that the geographic structuring of host populations must influence the genetic structuring of mistletoe populations, the current distribution range and phylogeographic structure of P. schiedeanus populations should reflect the observed host population expansion/contraction events historically related to the availability of suitable host trees. However, this hemiparasitic mistletoe may not track their L. styraciflua hosts during range expansions and contractions linked to Pleistocene climate changes because it lives on several host tree species (non-obligate host-parasite relationship) and, as a consequence, its distribution range was not influenced by geographical host range changes during glacial/interglacial periods.

Samples and DNA sequencing

Leaf tissue samples were collected from 274 individuals on a variety of host tree species in 31 populations throughout the species range in Mexico (Additional files 1 and 2). Target sampling localities were chosen based on localities of occurrence data acquired from Kuijt [3]. Psittacanthus schiedeanus can be relatively easy to identify in the summer time because of its conspicuous and long flowers. However, it is an extremely variable species in flower and leaf morphology sometimes difficult to separate from sympatric P. calyculatus [3], even though its typical form west of the Isthmus of Tehuantepec has much longer and more slender flowers [35, 47] and P. angustifolius with narrow leaves similar to occasional individuals of P. schiedeanus but prefers Pinus chiapensis as host at higher elevations in Chiapas and Pinus oocarpa in Honduras [3, 48], and from P. mayanus (Yucatan Peninsula), P. breedlovei (Chiapas) and P. minor (northern Nicaragua) allopatrically distributed in tropical deciduous forests mainly on Acacia host trees at lower elevations [3]. Therefore, our geographic sampling scheme includes populations considered phenotypically P. schiedeanus along the cloud forest distribution and samples phenotypically similar to P. schiedeanus of adjacent populations regardless of taxonomy, focusing on those distributed in central Oaxaca (P. calyculatus) and Chiapas (P. breedlovei). Most populations collected have an accompanying voucher that is deposited at the XAL herbarium of the Instituto de Ecología, AC (INECOL) (Additional file 2). Only one P. schiedeanus plant was sampled per individual host tree.

Leaf tissue samples were preserved in silica gel desiccant until DNA extractions were performed. We additionally obtained leaf samples for DNA sequencing from other Psittacanthus species (Additional file 3) and downloaded DNA sequences from the GenBank of representatives of Loranthaceae tribes to be used as outgroups: Nuytsia floribunda (DQ333867, DQ788716), Atkinsonia ligustrina (DQ333865, DQ788714), Gaiadendron punctatum (DQ333866, DQ340617), Peraxilla tetrapetalla (DQ333846, DQ340597), Alepis flavida (DQ333847, DQ340598), Desmaria mutabilis (DQ333852, DQ340603), Phthirusa pyrifolia (DQ333857, EU544504), Oryctanthus occidentalis (DQ333862, DQ340613), Tupeia antarctica (DQ333850, DQ340601), Phragmanthera regularis (DQ333830, DQ340579), Tristerix corymbosus (DQ333854, DQ340605), Tristerix aphyllus (DQ442966, DQ442919), Tripodanthus acutifolius (DQ333864, DQ340615), Notanthera heterophylla (DQ333855, DQ340606), Ligaria cuneifolia (DQ333853, DQ340604), Cladocolea cupulata (DQ333861, DQ340612), Cladocolea mcvaughii (DQ333860, DQ340611) and Struthanthus orbicularis (DQ333856, DQ340607) from Wilson and Calvin [49], Amico et al. [50] and Vidal-Russell and Nickrent [51].

Total genomic DNA was extracted from silica-dried material using a modified 2 × cetyl trimethyl ammonium bromide (CTAB) protocol [52] 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 ITS5HP [53] and ITS4 [54], whereas for trnL-F intergenic spacer region we used the universal primers e and f [55]. For targeting successful of ITS region sequencing, we used the primers ITS-F2-Psitta (5′-TCGCAGTATGCTCCGTATTG-3′) and ITS-R2-Psitta (5′-TCGTAACAAGGTTTCCGTAGG-3′) designed for this study. The 25–μL ITS PCR mix contained 5 μL of 5 × buffer (Promega, Madison, WI, USA), 2.0 μL MgCl₂ (25 mM), 2 μL dNTPs mix (8 mM), 0.82 μL of each primer (10 μM), 0.27 μL Taq polymerase (5U/μL) (Promega), 2 μL of DMSO (Sigma), 1.5–3 μL of template DNA, and finally dH₂O added to bring to volume. The 25–μL trnL-F PCR mix contained 3.5 μL of 5 × buffer (Promega, Madison, WI, USA), 3.5 μL MgCl₂ (25 mM), 2 μL dNTPs mix (8 mM), 0.45 μL of each primer (10 μM), 0.20 μL Taq polymerase (5U/μL) (Promega), 1.5–3 μL of template DNA, and finally dH₂O added to bring to volume. PCRs of ITS consisted of an initial denaturation at 94 °C for 4 min, followed by 5 cycles of 94 °C for 1 min, 50 °C for 30 s, 72 °C for 1 min, followed by 30 cycles of 94 °C for 30 s, 45–52 °C for 30 s, 72 °C for 1 min, and a final step of 72 °C for 7 min. Chloroplast (trnL-F) amplifications used the following profile: initial denaturation at 95 °C for 5 min, followed by 40 cycles of 95 °C for 10 s, 55 °C for 1 min, 72 °C for 20 s, and a final step of 72 °C for 7 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, California, USA). The products were analyzed on a 310 automated DNA sequencer (Applied Biosystems) at the INECOL’s sequencing facility, or at University of Washington High Throughput Genomics Unit, Seattle, Washington. Finally, sequences were assembled using Sequencher v4.9 (Gene Codes, Ann Arbor, MI, USA) and then were manually aligned with SE-AL v2.0a111 (http://tree.bio.ed.ac.uk/software/seal). All unique sequences of P. schiedeanus and those new of other Psittacanthus species have been submitted to GenBank (Accession nos. ITS: KU922961–KU923004, KU923036–KU923037; trnL-F: KU923270–KU923278, KU923310–KU923311; Additional file 3).

Paleodistribution modeling and environmental variation

We used species distribution modeling (SDM; [56]) to predict where populations of P. schiedeanus resided at the LGM (ca. 20,000 years ago) and Last Inter Glacial (LIG; ca. 130,000 years ago), and whether range expansion and population connectivity are observed according to the predicted changes in the distribution of cloud forests [57, 58]. Coordinates of occurrence data were assembled for P. schiedeanus obtained from the Global Biodiversity Information Facility (GBIF; http://www.gbif.org/species/4003073), supplemented with records from field collection (see dataset on Dryad: doi:10.5061/dryad.t49h6). After careful verification of every data location and removing duplicate occurrence records, we restricted the data sets to unique records for the analyses, leaving 58 unique presence records. Distributional records were input into and analyzed to infer a SDM with the maximum entropy algorithm in MAXENT v3.2.2 [59] and using the ArcView v3.2 (ESRI, Redlands, CA, USA) to extract the GIS data. Present-day temperature and precipitation data (BIO 1–19 variables) were drawn as climate layers from the WorldClim database at a spatial resolution of ca. at 1 km² in each raster ([60]; http://www.worldclim.org). We first constructed a model to the present using all climatic variables with 80 % of the locality records as training data and the other 20 % as testing data in addition to response curves, jackknife tests, logistic output format and random seed parameters. Variable importance was determined comparing percentage contribution values and jackknife plots. Then, correlation coefficients between present variables were calculated using JMP v5.1 (SAS Institute Inc. Cary, NC, USA) 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, five were used in the final analyses: Temperature Seasonality (BIO4), Min Temperature of Coldest Month (BIO6), Temperature Annual Range (BIO7), Precipitation of Wettest Month (BIO13) and Precipitation Seasonality (BIO15). Final models were constructed with 10 cross-validation replicates without extrapolation and considering the average output grids as the final predictive models. The area under the receiver operating characteristic curve (AUC) was used to evaluate the prediction performance of the models, where 1 is the maximum prediction and 0.5 suggest a random prediction. Resulting species distributions under current climate conditions were projected onto past climate scenarios, at the LGM (at 2.5 arc-minutes) and LIG (at 30 arc-seconds). Past climate layers were also drawn from the WorldClim webpage for two LGM scenarios developed by the Palaeoclimate Modelling Intercomparison Project Phase II [61]: the Community Climate System Model (CCSM; [62]) and the Model for Interdisciplinary Research on Climate (MIROC; [63]), and the LIG [64]. The CCSM and MIROC climate models simulate different climate conditions, with cooler sea-surface temperature conditions assumed in CCSM than in MIROC, resulting in higher annual precipitation in CCSM than in MIROC [58, 65, 66].

Measurements of temperature and precipitation (BIO 1–19 variables) were extracted for each of the sampling locations from WorldClim [60] as explained above. We performed principal components analysis (PCA) using SPSS v17 (SPSS, Armonk, NY, USA) to reduce intercorrelated BIO variables to a smaller system of uncorrelated, independent variables, and used the resulting PC loadings to examine habitat and ecological variation potentially related to population divergence.

Phylogenetic analysis, divergence time estimation and haplotype relationships

Phylogenetic relationships among sequences based on Bayesian Inference (BI) were reconstructed using MRBAYES v3.12 [67, 68]. BI analyses were run using the CIPRES Science Gateway [69] for the ITS, trnL-F and the concatenated data sets (see ITS and trnL-F datasets on Dryad: doi:10.5061/dryad.t49h6). jMODELTEST v0.1.1 [70] was run to choose the model of molecular evolution that best fitted our sequence data under the Bayesian information criterion (BIC), GTR+G, for the ITS, trnL-F and for the concatenated data set. Two parallel Markov chain Monte Carlo (MCMC) analyses were executed simultaneously, and each was run for 10 million generations, sampling every 1000 generations. Bayesian posterior probability values were calculated from the sampled trees remaining after 2500 samples were discarded as burn-in [67] to only include trees after stationarity was reached. The remaining trees were used to generate a 50 % majority-rule consensus tree, showing nodes with a posterior probability (PP) of 0.5 or more. We consider nodes significantly supported if posterior probabilities were ≥0.95 [67]. While Struthanthus, Cladocolea or Aetanthus may be the closest relatives of Psittacanthus [49, 51], relationships among phenotypically similar Psittacanthus species remain uncertain [3]. Thus, we included several congeners to root P. schiedeanus sequences. The likelihood scores under the ITS, trnL-F and the concatenated data sets were compared with ln Bayes factors (BF) tests to determine which tree topology significantly improved explanation of the data.

To relate genetic differentiation found among P. schiedeanus ribotypes and haplotypes to pre-Pleistocene and Pleistocene events, we estimated divergence time under a Bayesian approach as implemented in BEAST v1.6.1 [71] using the combined ITS and plastid trnL-F data set. The ingroup comprised all nrDNA or plastid cpDNA sequences of P. schiedeanus and sequences downloaded from GenBank of species representing Loranthaceae tribes from Amico et al. [50] and Vidal-Russell and Nickrent [51, 72] used as multiple outgroups. Divergence time estimation was performed in BEAST using the uncorrelated lognormal relaxed molecular clock and the closest nucleotide substitution model under the Bayesian information criterion (BIC), GTR+G for the concatenated data set, suggested by jMODELTEST. The tree prior model was set using a coalescent approach assuming constant population size. In this analysis, we constrained Nuytsiae as the sister group of members of other subtribes and Psittacanthus to be monophyletic based on Vidal-Russell and Nickrent [72] and the BI tree topology.

To calibrate the root node, the divergence time between Nuytsia floribunda Australian root parasite and aerial loranth parasites clade [72] was used as secondary calibration, approximating a median age of 28 MYA (normal distribution, mean 28, SD 3.7, range 34–21 MYA). The geometric mean of 5.798 × 10⁻⁹ substitutions per neutral site per year (s/s/y) was used to calibrate the tree based on the mean mutation rates of 4.13 × 10⁻⁹ s/s/y for ITS of herbaceous annual/perennial plants [73] and 8.24 × 10⁻⁹ s/s/y for the trnL-F estimated for annual or perennial herbs [74]. For divergence time estimation, BEAST was run two times for 10 million generations, sampling every 1000 steps. We combined log and trees files from each independent run using LOGCOMBINER v1.8.0 [71], then viewed the combined log file in TRACER v1.6 (http://tree.bio.ed.ac.uk/software/tracer/) to ensure that effective sample size (ESS) 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 v1.8.0 [71] and visualized in FIGTREE v1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/).

To infer genealogical relationships among ribotypes (ITS) and haplotypes (trnL-F), statistical parsimony networks for both data sets were constructed as implemented in TCS v1.2.1 [75], with gaps treated as missing data and a connection limit set to 95 %. Loops were resolved following the criteria given by Pfenninger and Posada [76].

Population genetic diversity indices, geographic structure and relationships among populations

Population diversity for unordered (h _S, h _T) and ordered haplotypes (v _S, v _T) and differentiation (G _ST, N _ST) parameters were estimated using PERMUT v1.0 [77]. Significant differences between N _ST and G _ST parameters were tested with 10,000 permutations. If N _ST is significantly larger than G _ST, this could indicate that the haplotypes found in a given population are phylogenetically closely related [77]. We also calculated haplotype diversity (h), nucleotide diversity (π) and pairwise comparisons of F _ST values between groups (Additional file 2) in ARLEQUIN v3.01 [78] with 16000 permutations. Populations with three samples or fewer were excluded from this analysis.

To determine whether or not populations are structured by geography or habitat distribution, three analyses of molecular variance (AMOVAs; [79]) were performed based on pairwise differences using ARLEQUIN with populations treated as (a) a single group to determine the amount of variation partitioned among and within populations, and (b) grouped into east and west of the Isthmus of Tehuantepec, (c) three groups according to habitat distribution or (d) grouped into six groups according to geography and mountain range (Additional files 1 and 2). AMOVAs were performed using the Tamura-Nei model and 16,000 permutations to determine the significance of each AMOVA.

To estimate the relationships among groups of populations (Additional files 1 and 2), we used ITS and trnL-F sequences for P. schiedeanus samples and *BEAST [80] with the multispecies coalescent model implemented in BEAST. *BEAST models the lineage sorting process between units for groups of individuals not connected by gene flow above, at, or below the species level [80]. Psittacanthus mayanus samples were used as outgroups. Nucleotide substitution models, JC for ITS and F81 for trnL-F selected with jMODELTEST were incorporated as the Hasegawa-Kishino-Yano (HKY) model for both markers. Species were assigned based on the identification of each sample using haplotype or ribotype distribution and locality. The simulation was first run with all samples as two lineages separated by the Isthmus of Tehuantepec (WEST and EAST), then with samples as six lineages corresponding to geography and mountain ranges (nSMO, cSMO, sSMO, CHIS, OAX, BREE) or three separate lineages according to habitat (SCHI = cloud forests from San Luis Potosí to Oaxaca and Chiapas, CALY = xeric vegetation in central Oaxaca, BREE = tropical deciduous forests in Chiapas). We ran each BEAST analysis two times for 30 million generations, sampling every 1000 steps, using a Yule speciation tree prior, relaxed clock model with an uncorrelated lognormal distribution, and the mean mutation rates of 4.13 × 10⁻⁹ s/s/y for ITS [73] and 8.24 × 10⁻⁹ s/s/y for the trnL-F [74]. After the analysis in BEAST, log and tree files were combined using LOGCOMBINER and summarized as a maximum clade credibility tree using TREEANNOTATOR with a burn-in of 25 %. We used TRACER to visualize the results of the runs and to check the ESS (cut off values >50) of each parameter.

Demographic history

The demographic history of each P. schiedeanus resulting group was determined by means of neutrality tests and mismatch distributions carried out in ARLEQUIN. To test whether populations evolve under neutrality, Fu’s Fs test [81] and Tajima’s D [82] were calculated, and mismatch distributions [83] were calculated using the sudden expansion model of Schneider and Excoffier [84] with 1000 bootstrap replicates. The validity of the sudden expansion assumption was determined using the sum of squares differences (SSD) and Harpending’s raggedness index (Hri; [85]), both of which are higher in stable, non-expanding populations [86].

Bayesian skyline plots [87] were also used to infer changes in effective population size (N _e) of each P. schiedeanus group (see Results) through time in 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 10 million generations each were run, with trees and parameters sampled every 1000-iterations, with 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 50.

We used the software IMa [88] on P. schiedeanus genetic groups (see Results; see dataset on Dryad: doi:10.5061/dryad.t49h6) to estimate the effective population size of the ancestral (q _a) and the two descendant populations (q ₁ and q ₂), effective number of migrants per generation in both directions (m _1-to-2 and m _2-to-1 ), and time since divergence (t) at which the ancestral population gave rise to the descendant populations. The isolation-with-migration (IM) model [89] 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 runs of 10,000 steps (following 100,000 iterations as burn-in) to assess mixing and to fine-tune the parameter space. We then conducted the simulation for a burn-in of 1 million generations and 30 million steps, under the HKY model of sequence evolution. Two independent runs were performed with different seed numbers to guarantee convergence of samples [89]. 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 two runs and the 90 % highest posterior densities (HPD) intervals of each parameter. We used the geometric mean of 2.49 × 10⁻⁶ of the mutation rates (per year) for both loci (1.33 × 10⁻⁶ for ITS and 4.69 × 10⁻⁶ for the trnL-F). The mutation rates were converted to per locus rate by multiplying by the fragment length in base pairs for conversion to demographic units, as required by IMa [88]. We used an 11-year generation time based on the observation that the age at maturity (seed production) begins ca. 2 years after seed germination and an assumed annual survivorship of 0.9 based on mistletoe seedling survivorship [30] to convert the effective populations size estimates. The approximate average generation time (T) is calculated according to T = a + [s ⁄ (1–s)] [90], where a is the time to maturity and s is the adult annual survival rate. Based on this, 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 mean rate described above.

Evolutionary transition of mistletoe invasion

The evolutionary transition of mistletoe invasion was inferred from an approximate Bayesian computation in which evolutionary scenarios were simulated and compared through posterior probabilities [91] using the DIYABC software [92]. Populations covering the whole species’ distribution were analyzed to infer the history of the genetic structure indicated by phylogenetic and phylogeographic analyses of ITS and trnL-F sequences. Because three haplotype clusters are hypothesized to be associated with individuals/populations separated by habitat (BREE, SCHI, CALY; Additional file 2) with post-glacial population expansions and potentially corresponding to post-glacial invasion types (see Results), five evolutionary scenarios were built and tested considering haplotype network and *BEAST analyses: (a) three populations (Pop1, Pop2 and Pop3) have diverged simultaneously from an ancestral population at t1 (scenario 1, null model), which corresponded to the plant genetic groups, BREE, SCHI and CALY, respectively; (b) isolation split model 1 (scenario 2), in which Pop1 (BREE) merged with Pop2 (SCHI) at t1 and subsequently with Pop3 (CALY) at t2; (c) isolation split model 2 (scenario 3), in which Pop3 merged with Pop2 at t1, then both populations merged with Pop1 at t2; (d) isolation split model 3 (scenario 4), in which Pop3 merged with Pop1 at t1 and then both populations merged with Pop2 at t2; and (e) isolation with admixture model (scenario 5), in which Pop2 was generated by admixture of Pops 1 and 3 at t1, then Pop1 merged with Pop3 at t2. Although there are numerous possible scenarios of divergence, we considered that these five scenarios represent the close relationships among the groups and the complex patterns of genetic divergence and admixture (see Results).

We generated 1 million simulated datasets per scenario considering a generalized stepwise-mutation model, uniform prior distributions for effective population sizes (10–100 000) and splitting events (100–50 000 generations). Comparison of scenarios was implemented in the DIYABC software [92]. The posterior probability of scenarios was assessed using a weighted polychotomous logistic regression on the 1 % of simulated datasets closest to the observed data [93]. For the best-supported scenario, we performed a model checking procedure by applying a PCA on test statistic vectors to visualize the fit between the simulated and the observed datasets. To assess the confidence in scenario choice, we simulated 500 pseudo-observed datasets (PODs) under each scenario to estimate Type I and Type II error rates [85]. Finally, point estimates for demographic and temporal parameters were obtained by local linear regression on the 1 % of simulations closest to the observed data set for the best-supported scenario [91, 92].

Results

Paleodistribution modeling and environmental variation

The SDM yielded a good fit for the current geographic distribution of P. schiedeanus. The AUC value for the replicate runs was 0.992 ± 0.003 (Fig. 1). Although there was some over prediction in certain geographic areas (the Trans-Mexican Volcanic Belt that was predicted by the model with moderate probability), the modeled distribution closely matches the current known distribution of P. schiedeanus. The paleodistribution modeling revealed that suitable habitat for P. schiedeanus during the LIG was more fragmented and restricted to mountain areas along the Pacific slope (Fig. 1a). However, the predictions for the LGM revealed that conditions of suitable habitat potentially expanded the distribution of P. schiedeanus, particularly in the eastern side of the Sierra Madre Oriental, Oaxaca and Chiapas applying both CCSM and MIROC simulations (Fig. 1b–c), and with a low probability into the Yucatan Peninsula under the CCSM scenario (Fig. 1b). Overall, the comparison between present and past climatic conditions predicted by the models suggests that climatic conditions of suitable habitat for P. schiedeanus experienced range shifts from the Pacific slope to east, expansion from LIG to LGM and contraction from LGM to present (Fig. 1a–d).

Fig. 1

Results from the MAXENT analyses showing species distribution models for Psittacanthus schiedeanus (a) Last Inter Glacial (LIG, 140–120 ka), (b) Last Glacial Maximum (LGM, CCSM, 21 ka), (c) Last Glacial Maximum (LGM, MIROC, 21 ka), and (d) at present. The output of MAXENT consists of grid maps with each cell having an index of suitability between 0 and 1. Low values (white to green) indicate that the conditions are unsuitable for the species to occur, whereas high values (orange to red) indicate that the conditions are suitable. There is clear evidence that populations are connected during glacial cycles. Note that for maps b and c projections extend out into the ocean because of changes in sea levels

The PCA of environmental data indicated two niche axes that together explained 77.9 % of the variation in P. schiedeanus (Fig. 2). The first niche axis (41.2 % of variation) was associated with precipitation variables (BIO12–BIO19), and mean temperature diurnal range (BIO2), isothermality (BIO3), temperature seasonality (BIO4) and temperature annual range (BIO7). The second niche axis (36.7 %) was associated with temperature variables (BIO1, BIO5–BIO6, BIO8–BIO11).

Fig. 2

Principal components analysis (PCA) on the 19 BIOCLIM variables with the first principal component (PC1) largely a measure of precipitation conditions, and the second principal component (PC2) mainly determined by temperature measures. Numbers refer to collection sites (see Additional file 2) and pie charts represent trnL-F haplotypes found in each sampling locality (see Fig. 5)

Phylogenetic analysis, divergence time estimation and haplotype relationships

The aligned sequences of ITS region (n = 292) and the trnL-F intergenic spacer (n = 283) for Psittacanthus species were 569 and 321 base pairs (bp) long, respectively. Out of 1018 bp of the concatenated sequences, 588 sites were variable and out of these 464 were parsimony informative. In P. schiedeanus samples, we detected polymorphic sites in nine samples of the ITS region, with four of the 569 aligned sites (0.7 %). These were coded as ambiguous.

The trees obtained from the Bayesian analyses of the individual molecular markers (ITS and trnL-F) and the combined data set all differ in topology and degree of resolution. The trnL-F gene produced the poorest resolved tree (Additional file 4) and also contains the smallest number of nodes with posterior probability (PP) values above 0.90. The combined ITS and trnL-F data set produced the most resolved tree and contains the highest number of nodes with PP values above 0.90. The topology of the resulting BI tree showed monophyly of Psittacanthus samples (PP = 1; Additional file 4), with an early divergence of a clade composed of P. ramiflorus, P. sonorae and P. palmeri, followed by a clade with moderate support containing samples of P. rhynchanthus (PP = 1) and a clade with moderate support (PP = 0.86) containing samples of the Brazilian species (P. robustus, P. biternatus, P. cordatus, P. acinarius) and samples of P. macrantherus, P. mayanus, P. auriculatus, P. calyculatus and P. schiedeanus. Samples of P. calyculatus and samples of P. schiedeanus from the eastern Mexico cloud forests, central Oaxaca and Chiapas turned out to be monophyletic, but little support (PP = 0.67) and genetic structure was observed within this clade. The Bayes factor tests showed extensive incongruence between the chloroplast and the nuclear partitions. Combining ITS and trnL-F data sets produced a higher likelihood score than those for the alternative ITS or trnL-F data sets (2 × ln BF; ITS versus combined = 5031.742; trnL-F versus combined = 15946.436; ITS versus trnL-F = 10914.694).

The BEAST analyses placed the origin of the Psittacanthus crown clade in the Early Pliocene (Fig. 3). Although a split between Psittacanthus and a clade composed of samples of Oryctanthus, Phthirusa, Tripodanthus, Cadocolea and Struthanthus is recovered in the Bayesian analysis with strong support (PP = 1.0), the relationships within this clade were poorly resolved (Fig. 3). Within Psittacanthus the group formed by samples of P. ramiflorus, P. palmeri and P. sonorae split from other Psittacanthus species in the ingroup 5.41 MYA (95 % HPD 6.75–1.11 MYA) assuming constant population size. Samples of P. schiedeanus formed a well-supported monophyletic group (PP = 0.94), with diversification occurring in the early Pleistocene (2.55 MYA, 95 % HPD 1.57–0.23 MYA) and a split between P. schiedeanus and P. auriculatus estimated at 2.80 MYA (95 % HPD 2.08–0.31 MYA) assuming constant population size (Fig. 3).

Fig. 3

a Chronogram based on a Bayesian approach using a coalescent prior under an uncorrelated lognormal relaxed clock model and assuming constant population size of Psittacanthus schiedeanus sequences in BEAST. Pink bars indicate 95 % highest posterior density (HPD) intervals for nodes of particular interest. These nodes all have posterior probabilities above 0.9. b *BEAST model with simulation run with all samples as two lineages separated by the Isthmus of Tehuantepec (WEST and EAST), (c) samples as six lineages corresponding to geography and mountain ranges (nSMO, cSMO, sSMO, CHIS, OAX, BREE) or (d) three separate lineages according to habitat (SCHI cloud forests from San Luis Potosí to Oaxaca and Chiapas, CALY xeric vegetation in central Oaxaca, BREE tropical deciduous forests in Chiapas). Maximum-likelihood scores are indicated for *BEAST analyses (b–d). The likelihood scores under the three species delimitation hypothesis (SCHI, BREE, CALY) are compared with ln Bayes factors (BF) tests. The three-species hypothesis produced a higher likelihood score than those for the alternative hypotheses although the difference was not very strong (2 × ln BF; six-species hypothesis versus three-species hypothesis = 11.228; two-species hypothesis versus three-species hypothesis = 4.046; six-species hypothesis versus two-species hypothesis = 7.182)

The aligned ITS data set for 251 samples of P. schiedeanus (559 bp) yielded 44 ribotypes (Additional file 5). In the ribotype network (Fig. 4), one ribotype (R5) was found in 113 of the 251 samples for 22 of the 29 sampled populations and forms the core of the network, from the northernmost population (Xilitla, San Luis Potosí) to Chiriquí, Panama. The most frequent and widespread ribotype (45 % of the individuals and 75.8 % of the sampled populations), inferred as ancestral based on its frequency and position in the network (Fig. 4), was separated by one or two mutations of most ribotypes in the network and they were singletons; it was not retrieved in two easternmost populations in Oaxaca and in several populations from Chiapas. The second most frequent and most widespread ribotype (R23) was found (31 samples) in three populations of Oaxaca (16 samples) and six populations from central Veracruz (13 samples) and one population from Chiapas (2 samples). Although ribotypes R1, R2, R19, R31, and R35 were not singletons they had certain geographic structure, with R1 and R2 mostly distributed in the northernmost populations (localities 1 and 2), R19 in two populations from Chiapas (5 samples), one from Veracruz (1 sample) and one population from Oaxaca (3 samples), R31 exclusively found in one population from Chiapas (9 samples), and R35 found in one population from Veracruz (1 sample) and three populations from Chiapas (11 samples).

Fig. 4

Geographic distribution and statistical parsimony network of 44 ITS ribotypes of Psittacanthus schiedeanus. Current natural range of cloud forests (indicated by green shading) is overlaid on a relief map of eastern Mexico. Pie charts represent ribotypes found in each sampling locality. Ribotype distributions at central Veracruz are shown separately in the inset

The aligned trnL-F data set for 242 samples of P. schiedeanus (282 bp) yielded 9 haplotypes (Fig. 5; Additional file 5). Statistical parsimony retrieved a well-resolved network, in which three main haplotypes could be distinguished (Fig. 5). Samples from cloud forest localities shared the most frequent (35.5 % of the individuals and 80.6 % of the sampled populations) haplotype (H3) forms the core of the “cloud forest” haplogroup, from Xilitla, San Luis Potosí to Santiago Comaltepec, Oaxaca and populations from Chiapas; it was not retrieved in the populations with more xeric vegetation in central Oaxaca or tropical deciduous forests in Chiapas. Haplotypes connected to H3 by one or more steps consists of haplotypes exclusively found in populations from central Veracruz (H1–H2, H4–H6), H7–H9 located in populations in central Oaxaca and H7 and H9 in populations from Chiapas with tropical deciduous forests (Fig. 5).

Fig. 5

Geographic distribution and statistical parsimony network of 9 trnL-F haplotypes of Psittacanthus schiedeanus. Current natural range of cloud forests (indicated by green shading) is overlaid on a relief map of eastern Mexico. Pie charts represent haplotypes found in each sampling locality. Haplotype distributions at central Veracruz are shown separately in the inset

Population genetic diversity indices, geographic structure and relationships among populations

Differentiation among populations based on ITS variation indicated that P. schiedeanus is genetically subdivided (G _ST = 0.295, SE = 0.0629). Genetic diversity across all populations (h _T = 0.804, SE = 0.0431; v _T = 0.807, SE = 0.1171) was higher than the average within-population (h _S = 0.567, SE = 0.0499; v _S = 0.484, SE = 0.0683). PERMUT analysis showed that N _ST (0.400, SE = 0.0639) was significantly higher (p <0.05) than that for G _ST, indicating phylogeographical structuring. Similarly, phylogeographical structuring was observed in P. schiedeanus based on trnL-F variation (G _ST = 0.692, SE = 0.1058, N _ST = 0.741, SE = 0.1122; N _ST > G _ST, p <0.05), with genetic diversity across all populations (h _T = 0.486, SE = 0.0925; v _T = 0.487, SE = 0.1221) much higher than the average within-population (h _S = 0.150, SE = 0.0497; v _S = 0.126, SE = 0.0533). Pairwise comparisons of F _ST values were not significant for ITS when groups of populations separated by habitat type were compared, whereas pairwise comparisons of F _ST values for trnL-F were high and significant (Table 1).

Table 1 Pairwise comparisons of F _ST values of ITS (above the diagonal) and trnL-F (below the diagonal) among populations of Psittacanthus shiedeanus grouped by habitat type

The AMOVA results showed that 82.5 % of the genetic variation for ITS and 21.1 % for trnL-F was explained by differences within populations and 17.5 % and 78.9 %, respectively, by differences between populations when all locations were treated as a single group (Table 2). The AMOVA for the trnL-F revealed population structure, with highest F _CT value (F _CT = 0.78) obtained when populations are grouped by habitat type but differences between groups were not significant in the ITS data set (Table 2), suggesting that seed-mediated dispersal is more limited than pollen-mediated dispersal. When sampling sites are grouped as separated by the Isthmus of Tehuantepec or grouped by mountain range, significant but smaller proportion of the variation was attributed to differences between groups (Table 2). Nucleotide diversity (π) was low for both data sets (Additional file 6).

Table 2 Results of AMOVA models on Psittacanthus schiedeanus populations with no groups defined a priori (a), and grouped into east and west of the Isthmus of Tehuantepec (b), three groups according to habitat distribution (c) or grouped into six groups according to geography and mountain range (d)

The *BEAST tree (Fig. 3b) of multilocus data for differentiation between clades separated by the Isthmus of Tehuantepec showed strong support (PP = 1; tree inset of Fig. 3), and the split occurred at ca. 0.098 MYA (95 % HPD: 0.16–0.047 MYA). However, this scenario was not the best supported compared with alternative species assignments (Fig. 3c–d).

Demographic history

When populations are grouped as separated by habitat type, neutrality test values were negative and significant for SCHI (except Fu’s Fs in the trnL-F), indicating demographic expansion for SCHI (Additional file 6). In the mismatch distribution, sudden demographic expansion was not rejected for SCHI and CALY (Hri values) in the ITS, and in the trnL-F the Hri value was not rejected for SCHI (Additional file 6). The Bayesian skyline plots suggested that the effective population size increased in the SCHI population and the population as a whole around 80,000–10,000 years ago (Fig. 6). Changes in the effective population size were only observed in the concatenated and ITS data sets.

Fig. 6

Bayesian skyline plots showing historical demographic trends of Psittacanthus schiedeanus for all pooled populations (a–c) and for the SCHI population (d–f) using the concatenated, ITS and trnL-F data sets, respectively. 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. SCHI cloud forests from San Luis Potosí to Oaxaca and Chiapas

IMa results summarized in Table 3 are reported as the average of two runs of mean parameter estimates and the 95 % HPD intervals of each parameter. The ancestral population size (N _E) is estimated to be higher than population sizes of the SCHI and CALY+BREE descendant populations (Table 3). Divergence time between the SCHI and the other populations separated by habitat type was estimated at ca. 1 MYA. When testing for migration during the time since population splitting, migration among populations occurred predominantly in one direction, moving from CALY+BREE to SCHI (Table 3).

Table 3 Results of isolation-with-migration model (IMa) for the split between groups of populations of Psittacanthus schiedeanus

Evolutionary transition of mistletoe invasion

Considering the five scenarios tested for nrDNA and cpDNA sequences, DIYABC analysis indicated that a simultaneous split (scenario 1) is the best-supported scenario (Fig. 7), with a posterior probability value much higher than those for the other scenarios. The 95 % confidence intervals for this model did not overlap with those obtained for the other scenarios (Table 4). Based on model performance, we were able to discriminate scenario 1 from other scenarios. This is corroborated by model checking, which showed a large cloud of data from the prior and observed datasets centered on a small cluster from the posterior predictive distribution, suggesting that the best supported scenario explained the observed data well (Results not shown). Analyses to estimate confidence in scenario choice based on 500 PODs indicate that Type I and Type II errors for the best-supported scenario were low (Table 4). Under scenario 1, posterior mean parameter estimates indicated that the simultaneous divergence occurred 21,780 generations ago. Assuming an 11-year generation time, the simultaneous split was dated to 217,800 years before present (Table 5), and probably pre-dates the Illinoian glacial stage (ca. 191,000–130,000 years ago). In accordance with the isolation and migration model, DIYABC estimated a smaller Ne for the descendant populations compared with that of the ancestral populations (Table 5). Estimated mean mutation rates of nrDNA and cpDNA were estimated to be 3.97 × 10⁻⁸ and 5.13 × 10⁻⁸, respectively (Table 5).

Fig. 7

Competing demographic scenarios of Psittacanthus schiedeanus divergence and admixture. Because three haplotype clusters are hypothesized to be associated with individuals/populations separated by habitat (BREE, SCHI, CALY) with post-glacial population expansions and potentially corresponding to post-glacial invasion types (see Results), five evolutionary scenarios were defined and tested: (a) three populations (Pop1, Pop2 and Pop3) have diverged simultaneously from an ancestral population at t1 (scenario 1, null model), which corresponded to the plant genetic groups, BREE, SCHI and CALY, respectively; (b) isolation split model 1 (scenario 2), in which Pop1 (BREE) merged with Pop2 (SCHI) at t1 and subsequently with Pop3 (CALY) at t2; (c) isolation split model 2 (scenario 3), in which Pop3 merged with Pop2 at t1, then both populations merged with Pop1 at t2; (d) isolation split model 3 (scenario 4), in which Pop3 merged with Pop1 at t1 and then both populations merged with Pop2 at t2; and (e) isolation with admixture model (scenario 5), in which Pop2 was generated by admixture of Pops 1 and 3 at t1, then Pop1 merged with Pop3 at t2. Comparison of the scenarios was implemented in the DIYABC software [92]. Habitat type abbreviations are as follows: SCHI cloud forests from San Luis Potosí to Oaxaca and Chiapas, CALY xeric vegetation in central Oaxaca, BREE tropical deciduous forests in Chiapas (Additional file 2)

Table 4 Posterior probability of each scenario and 95 % confidence intervals (CI) based on the logistic regression approach for approximate Bayesian computation (ABC) analyses considering the groups of populations of CALY, BREE and SCHI mistletoes. Simulations and ABC analyses were performed considering both nrDNA and cpDNA sequences

Table 5 Posterior parameter estimates (median and 95 % confidence intervals) for the best-supported scenario (scenario 1) considering the three mistletoe groups (BREE, SCHI and CALY). 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 cpDNA sequences

Discussion

Mistletoes in the Mesoamerican cloud forests

In northern Mesoamerica, cloud forests encompass an extremely heterogeneous mixture of North American temperate tree species that are present in the region since the Tertiary and tropical species with South American origins during the early Miocene [39–42]. The dynamic Mesoamerican geological landscape from the Miocene to the Pliocene and the repeated cycles of cloud forest contraction and expansion due to Pleistocene climatic conditions [94, 95] affected the distribution and composition of cloud forests in the region, and contributed to both ancient [32, 46, 96, 97] and more recent divergences producing a common phylogeographic break at the Isthmus of Tehuantepec (e.g., [32, 34, 98]). The emerging phylogeographical patterns of temperate tree species that have migrated south from North America have been attributed to isolation by arid conditions during the Pliocene and climate changes during Pleistocene glaciations that promoted the expansion, contraction and divergence of populations [46, 99]. However, plant species that presumably colonized northern Mesoamerica from South America are little studied.

Seven studies have explored the phylogeography of widespread, cloud forest-adapted plant species across northern Mesoamerica. High levels of genetic variation and geographical structuring were inferred for wind-dispersed Alsophila firma (Baker) D.S. Conant (Cyatheaceae) [58] and Liquidambar styraciflua L. (Altingiaceae) [46] or bird-dispersed Podocarpus matudae Lundell (Podocarpaceae), Palicourea padifolia (Humb. & Bonpl. ex Schult.) C.M. Taylor & Lorence (Rubiaceae) and Rhipsalis baccifera (Sol.) Stearn (Cactaceae) [34, 98, 100]. On the other hand, for wind/gravity dispersed Begonia heracleifolia Cham. & Schltdl. (Begoniaceae) [101] and Moussonia deppeana (Schltdl. & Cham.) Hanst. (Gesneriaceae) [33], the genetic variation and genetic structure were inferred to be stronger. The high genetic differentiation of populations with strong spatial structuring provides support for the persistence of populations in situ during dry periods in the Pleistocene, in which populations remained isolated within the current fragmented distribution of the cloud forest [33, 101], suggesting that the vegetation along the Isthmus of Tehuantepec has acted as a barrier to seed dispersal and that the multiple vicariance events observed in various taxa occurred at earlier different times [32, 34, 98]. However, the widely documented genetic break at the Isthmus of Tehuantepec in several taxa [32, 96], though significant for the plastid trnL-F data, only explained 28 % of the total variation according to the AMOVA model, suggesting that low levels of seed flow have occurred across this barrier.

We found population genetic structure in P. schiedeanus, with differentiation between groups of populations separated by habitat type. Despite the genetic differentiation of some P. schiedeanus populations, there are widespread ribotypes and haplotypes in the ITS and trnL-F markers, respectively. In particular, the lower variation found in ITS and the presence of widespread ribotypes in the region are suggestive of effective nuclear gene flow via pollen (e.g., [34, 102]). In contrast, the population differentiation of P. schiedeanus shown by the plastid trnL-F data suggests that seed flow has been more restricted in some cases. Our results from the pairwise F _ST comparisons and AMOVA analysis suggest fairly high levels of genetic diversity among populations and habitat types of P. schiedeanus, especially when considering the cpDNA. However, significant evidence of population expansions was found in cpDNA sequences and the most frequent haplotypes were retrieved in most cloud forest areas of P. schiedeanus. BSPs, negative values of Tajima’s D, Fu’s Fs and the results of mismatch distribution for the plastid marker suggest that the population in the cloud forests (SCHI) has experienced recent demographic expansion, most likely during the LGM.

Predictions of SDMs under past LGM climatic conditions and the significant signal of demographic expansion suggest that P. schiedeanus (SCHI) populations experienced a range expansion tracking the conditions of the cloud forest distribution and shifted to the lowlands with population connectivity during the LGM, as predicted by the moist forests model [57, 58]. In addition, the Bayesian skyline plot for the SCHI population detected a significant increase in the overall effective population size dating to sometime during the LGM, when sudden demographic expansion occurred and descent to lower elevations allowed populations to gene flow. However, patterns of high genetic diversity among different populations do not fully support this hypothesis. The distribution of haplotype diversity suggests that the low genetic structure observed among populations of cloud forest areas is due to population connectivity across low-elevation barriers. Accordingly, we would hypothesize that during glacial periods, populations of P. schiedeanus expanded to the lowlands, gene flow was extensive, and habitats in which populations had gone extinct during the interglacials were probably re-colonized. Throughout interglacial periods, populations fragmented and contracted as they moved up in elevation with the warming climate and genetic divergence commenced across cloud forest areas dissected by mountains and low-elevation barriers.

Despite the appeal of the contraction/expansion model of populations of cloud forest-adapted species throughout glacial-interglacial climate cycles by shifting elevation, there is little evidence for a common demographic pattern of co-distributed species [57]. Ruiz-Sanchez and Ornelas [46] analyzed the phylogeographical structure of Liquidambar styraciflua (Altingiaceae), a deciduous tree commonly parasitized by P. schiedeanus with a disjunct distribution between the deciduous forest of the southeastern USA and the Mesoamerican cloud forest, and found signals of demographic expansion of Mesoamerican populations and habitat suitability connected and expanded to lower coastal areas under LGM climatic conditions. However, these patterns of genetic structure are not shared by P. schiedeanus. The observed inconsistency between demographic histories of P. schiedeanus and its most frequent host tree species in cloud forests can be the result of a differential response of species to climate change.

The role of historical climate changes on mistletoe evolution

The phylogeographical structure found in P. schiedeanus suggests that three lineages occupy geographic regions separated by habitat type, driving isolation and genetic divergence between groups of populations enhanced by environmental differences. Despite the high levels of genetic diversity and genetic structure of cpDNA among groups of populations (significant divergence in F _ST values and AMOVA), the occurrence of widespread plastid haplotypes in most areas (i.e., post-glacial population expansion) of the range of P. schiedeanus is consistent with the hypothesis of effective seed dispersal episodes. Accordingly, the mobility of bird dispersal agents has erased the genetic structure predicted to groups of populations along the cloud forest distribution (e.g., [32, 34]. It is possible that large-scale seed recruitment patterns produced by frugivorous birds that consume and eventually disperse mistletoe seeds might have led to increased gene flow among cloud forest areas.

Several bird species consuming fruits of P. schiedeanus mistletoes have been reported, the most important being cedar waxwings (Bombycilla cedrorum Vieillot, Bombycillidae), grey silky-flycatchers (Ptilogonys cinereus Swainson, Ptilogonatidae), and social flycatchers (Myiozetetes similis Spix, Tyrannidae) [14, 20, 36, 37]. The foraging and flocking behavior and local abundance bird dispersers differ widely [14, 36], and, consequently, affect differently the spatial patterns of mistletoe seed deposition [14]. On their winter migration cedar waxwings usually arrive in groups of 15–80 individuals on a single infected tree, whereas grey silky-flycatchers and social flycatchers occur in pairs or small flocks of 3–6 individuals. The winter migration route of cedar waxwings overlaps the distribution of P. schiedeanus, and the timing of its arrival coincides with the mistletoe fruiting season [36]. Bird migration could produce a large-scale seed rain pattern in a north to south direction, and this current directional gene flow could have erased past phylogeographic patterns.

Because cedar waxwings and grey silky-flycatchers tend to perch on mistletoe infected trees, their residence time on the most abundant host species ranges from 9 to 30 min [14] and defecate P. schiedeanus seeds within a short time of ingestion (ca. 20–30 min; [20]), it is also possible that the dispersal of P. schiedeanus may be geographically limited, leading to its local adaptation restricted to available host trees along cloud forests. Cedar waxwings defecate more mistletoe seeds onto branches of L. styraciflua individuals than onto those of other tree species [14], and more seedlings survive on branches of the most frequently infected host species, L. styraciflua [30]. In contrast with other host tree species, the higher germination rates of bird gut-processed mistletoes collected on Liquidambar trees is consistent with the prevalence and spatial distribution patterns reported by previous studies [14, 30, 36], and cross-infection experiments have shown the potential for local host adaptation in P. schiedeanus [20]. In this way, depending on the heterogeneity of host patches at a local scale, some mistletoe plants can develop more specificity on particular host trees that may lead to gene flow changes and the eventual formation of mistletoe races [5, 15, 18, 20]. However, the large host range suggests that P. schiedeanus is not host-specific, and in areas where Liquidambar is absent the parasite-host interaction is predominantly with other host species. Thus, geologic and climate-driven processes implicated in the fragmentation of the Mesoamerican cloud forests and, consequently the distribution of potential host species at a larger geographical scale, as well as selection pressures imposed by floral visitors and seed dispersers, could have influenced the distribution of the genetic variation among populations of P. schiedeanus. If the historical fragmentation of cloud forest patches or selection reduced gene flow among P. schiedeanus populations, host specificity may lead to genetic divergence and ultimately allopatric speciation.

However, the ABC analyses strongly supported the scenario of simultaneous divergence among the three groups dating back ca. 200 k years ago. This divergence time estimate for the groups of populations could have resulted from the low power of ABC analyses and our set of nuclear and chloroplast DNA sequences to detect subtle genetic differences in divergence times among groups of populations that have been isolated more recently. Under this scenario, climatic fluctuations throughout the Pleistocene would have altered the distribution of suitable habitat for mistletoes throughout Mesoamerica leading to variation in population continuity and isolation. Species that occur in topographically complex regions will frequently move between high and low elevation habitats in response to climatic fluctuations during glacial cycles [57]. Furthermore, the plant density and the clumped spatial distribution of conspecifics may affect hummingbird foraging behavior. Hummingbirds may forage more locally in search of nectar in high-density populations of P. schiedeanus, thereby reducing the dispersal distances of pollen. Therefore, pollen dispersal distance is expected to be lower in a high-density species distributed in a naturally fragmented habitat where individuals are locally aggregated on cloud forest host trees. However, its wide host and geographical ranges and large-scale seed recruitment patterns produced by long-distance bird migrants that consume and eventually disperse mistletoe seeds might have lead to increased gene flow, thus obscuring the geographical genetic structure created by cloud forest isolation before and during the Pleistocene.

Conclusions

Our phylogeographic approach allowed us to clarify relationships among populations of P. schiedeanus mainly distributed among cloud forests (SCHI) and the distribution boundaries with phenotypically similar Psittacanthus populations allopatrically distributed in more xeric environments (BREE in tropical dry forests of the Central Depression in Chiapas and CALY in more xeric plant communities of central Oaxaca). Genetic differentiation among these groups was strong and significant for the plastid trnL-F marker and non-significant for the nuclear ITS. Collectively, our data suggest that geographical, ecological, historical and environmental factors are the main drivers of genetic differentiation of these Mesoamerican parrot-flower mistletoes. According to SDMs, population expansion tests, skyline plots and divergence time estimations, the invasion of P. schiedeanus to cloud forests in eastern Mexico is relatively recent, and the invasion occurred to host tree species with longer histories in the region (e.g., [46]). Nonetheless, to better understand and manage threats of mistletoe populations, including potential transmission pathways of these commercially and ecologically important pathogens [1, 4] and genetic dilution from hybridization between P. schiedeanus and other sympatric Psittacanthus species, a more detailed sampling within the distribution of closely related P. calyculatus and P. breedlovei along the Trans-Mexican Volcanic Belt and tropical dry forests in the Central Depression of Chiapas is needed. Also, the genotyping of individuals using specifically designed microsatellites [103] will be critically important to the understanding of Psittacanthus species boundaries.

Declarations

Ethics approval

No aspect of this study required ethics approval.

Consent for publication

No aspect of this study required written informed consent to participate.

Availability of data and materials

Sequence reads can be accessed through GenBank under the Accession Numbers KU922961–KU923047 (ITS), trnL-F: KU923270–KU923321 (trnL-F). The data sets (alignments of ITS and trnL-F and IMa and MAXENT input files) supporting the results of this article are available in the Dryad Digital Repository (http://dx.doi.org/10.5061/dryad.t49h6) as Ornelas et al. [104].