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 MgCl2 (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 dH2O 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 MgCl2 (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 dH2O 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 km2 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−9 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−9 s/s/y for ITS of herbaceous annual/perennial plants [73] and 8.24 × 10−9 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−9 s/s/y for ITS [73] and 8.24 × 10−9 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
1 and q
2), 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−6 of the mutation rates (per year) for both loci (1.33 × 10−6 for ITS and 4.69 × 10−6 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].