Population Sampling

Leaf tissue samples were collected from 53 individuals in six locations (5–9 individuals per location) throughout the species range in Mexico (Table 1, Figure 1). Target sampling localities were chosen based on localities of occurrence data acquired from Mexican herbarium records and Kuijt (2009). Only one P. auriculatus plant was sampled per individual host tree. Leaf tissue samples were preserved in silica gel desiccant until DNA extractions were performed.

Table 1.

Geographic Information and Sample Sizes of the Psittacanthus auriculatus Localities of Oaxaca Sampled in This Study.

Figure 1.

Collection localities (population numbers as in Table 1) and geographical distributions of 8 chloroplast haplotypes (atpB-rbcL/trnL-F) found in six populations of Psittacanthus auriculatus in Oaxaca. The blue haplotypes correspond to the northern populations group and the green haplotypes to the southern populations group identified by AMOVA. Small black dots are nonsampled haplotypes.

Chloroplast Sequencing

Total genomic DNA was extracted from silica-dried material using a modified 2×cetyl trimethyl ammonium bromide protocol (Doyle & Doyle, 1987) or the DNeasy Plant Mini kit (Qiagen, Valencia, CA, USA) using the manufacturer protocol. Sequences of noncoding spacer region that lies between the large subunit of ribulose-1, 5-bisphosphate-carboxyalse (rbcL) and the beta-subunit of the chloroplast ATP-synthase (atpB; hereafter atpB-rbcL) in combination with sequences of the plastid trnL-trnF intergenic spacer region (hereafter trnL-F) from the chloroplast genome (cpDNA) were amplified by PCR and sequenced. These markers show the appropriate level of variation for mistletoe studies at the species and population levels (e.g., Amico, Vidal-Russell, García, & Nickrent, 2012; Ornelas et al., 2016; Pérez-Crespo et al., 2017). Amplification of the genes atpB and rbcL was conducted with the primers atpB-1 reverse and rbcL-1 reverse described in Chiang, Schaal, and Peng (1998), whereas for trnL-F region we used the universal primers trnL ‘c’ and ‘f’ described in Taberlet, Gielly, Pautou, and Bouvet (1991). The 25-µL PCR mix for atpB-rbcL contained 5 µL of 5×buffer (1×; Promega, Madison, WI, USA), 2.5 µL MgCl₂ (2.50 mM), 2.5 µL dNTPs mix (0.8 mM), 1 µL of each primer (0.4 µM), 0.4 µL Taq polymerase (0.08U/µL; Promega), 3 µL of template DNA, and finally dH₂O added to bring to volume. Amplification of atpB-rbcL used the following profile: initial denaturation at 94℃ for 2 min, followed by 35 cycles of 94℃ for 45 s, 50℃ for 1.15 min, 72℃ for 1.15 min, and a final step of 72℃ for 10 min. For the trnL-F, the 25-µL PCR mix for contained 3.55 µL of 5×buffer (0.71×; Promega, Madison, WI, USA), 3.6 µL MgCl₂ (3.60 mM), 1.80 µL dNTPs mix (0.58 mM), 0.45 µL of each primer (0.18 µM), 0.20 µL Taq polymerase (0.04U/µL; Promega), 1.5 to 3 µL of template DNA, and finally dH₂O added to bring to volume. Amplification of trnL-F used the following profile: initial denaturation at 95℃ for 5 min, followed by 40 cycles of 95℃ for 10 s, 55℃ for 1 min, 72℃ for 20 s, and a final step of 72℃ 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, CA, USA). The products were analyzed on a 310 automated DNA sequencer (Applied Biosystems) at University of Washington High Throughput Genomics Unit, Seattle, WA. Finally, edited sequences were imported into SE-AL 2.0a111 ( http://tree.bio.ed.ac.uk/software/seal) and aligned manually. For atpB-rbcL and trnL-F, indels were introduced during alignment but the resulting gaps were not coded. All unique sequences of P. auriculatus have been submitted to GenBank (Accession nos. atpB-rbcL: MF983856–MF983863, trnL-F: MF983864–MF983871).

Haplotype Relationships

To infer genealogical relationships among haplotypes, a statistical parsimony network for the combined atpB-rbcL/trnL-F data set was constructed as implemented in TCS 1.2.1 (Clement, Posada, & Crandall, 2000), with the 95% connection probability limit and treating gaps as single evolutionary events. Loops were resolved following the criteria given by Pfenninger and Posada (2002). Since the chloroplast is inherited as a single unit, the results we report herein are for the combined atpB-rbcL/trnL-F data set.

Genetic Diversity and Population Structure

Average gene diversity within populations (h_S, v_S), total gene diversity (h_T, v_T), and the differentiation for unordered allelles (G_ST) and for ordered alleles (N_ST) based on 10,000 permutations were estimated using PERMUT 2.0 (Pons & Petit, 1996). If N_ST is significantly larger than G_ST, this means that two alleles sampled from the same region are phylogenetically more closely related than two alleles sampled from different regions, evidencing the presence of a significant phylogeographical signal in the data (Pons & Petit, 1996). Molecular diversity indices, including the number of segregating sites (S), the number of haplotypes (N_H), nucleotide diversity (π), and pairwise comparisons of F_ST values between populations were calculated using ARLEQUIN 3.1.1 (Excoffier, Laval, & Schneider, 2005) with 1,000 permutations. Populations with three samples or fewer were excluded as required for this analysis.

To determine whether or not populations are structured, two analyses of molecular variance (AMOVAs; Excoffier, Smouse, & Quattro, 1992) 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 northern or southern according to parsimony network results. Locus-by-locus AMOVAs were performed using the Jukes and Cantor model for atpB-rbcL/trnL-F sequences and 16,000 permutations to determine the significance of each AMOVA.

Demographic History

Signatures of demographic changes or selection in the recent history of P. auriculatus were addressed using analyses based on different coalescent approaches. Data were first tested against a neutral Wright–Fisher model by means of neutrality tests and mismatch distributions carried out in ARLEQUIN. To test whether populations evolve under neutrality, Fu's Fs test (Fu, 1997) and Tajima's D (Tajima, 1989) neutrality test statistics were calculated, and to assess possible expansions, mismatch distributions (Harpending, 1994) were calculated using the sudden expansion model of Schneider and Excoffier (1999) with 1,000 bootstrap replicates. The validity of the sudden expansion assumption was determined using the sum of squares differences (SSD) and Harpending's raggedness index (Hri; Harpending, 1994), both of which are higher in stable, nonexpanding populations (Rogers & Harpending, 1992).

Ecological Niche Modeling

The observed demographic patterns of arid land plant populations were tested for P. auriculatus using ENM (Elith et al., 2011). ENM was used predict the distribution of P. auriculatus at the Mid-Holocene (MH; ca. 6,000 years ago), LGM (ca. 20,000 years ago), and at the LIG (ca. 120,000–140,000 years ago). Mistletoe occurrence data were assembled mainly from herbarium records from the Global Biodiversity Information Facility ( http://www.gbif.org/species/4003073) and supplemented with records from Kuijt (2009) and our geo-referenced records from field collection. After careful verification of every data location, excluding duplicate occurrence records or in close proximity to each other (ca. 1 km²) to reduce the effects of spatial autocorrelation, we restricted the data set to 49 unique presence records for the analysis. Distributional records were input and analyzed to infer an ENM with the maximum entropy algorithm in MAXENT 3.4.1 (Phillips, Anderson, & Schapire, 2006) and using the raster package 2.5-8 (Hijmans et al., 2016) in R 3.4.1 (R Development Core Team, 2017) to extract the data. Present-day and LIG temperature and precipitation data (BIO 1–19 variables) were drawn as climate layers from the WorldClim database at a spatial resolution of ca. 1 km² in each raster (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005;  http://www.worldclim.org). We also obtained data for past climate layers for two MH (3 arc-minutes, pixel size ca. 21.62 km²) and LGM (3 arc-minutes) past conditions based on the Community Climate System Model (CCSM) and Model Interdisciplinary Research on Climate (MIROC) global models (Braconnot et al., 2007; Hijmans et al., 2005).

We first constructed a model to the present using all climatic variables with 70% of the locality records as training data and the other 30% as testing data in addition to response curves, jackknife tests, logistic output format, and random seed parameters. We ran 1,000 iterations with no extrapolation in order to avoid artificial extrapolations from extreme values of the ecological variables; as such parameters are biased toward the environmental envelope of background points and occurrence data (Phillips et al., 2006). All other parameters in MAXENT were maintained at default settings. The logistic format was used to obtain the values of habitat suitability and the variable importance was determined comparing percentage contribution values and jackknife plots. Then, correlation coefficients between present variables were calculated using PAST 2.12 (Hammer, 2011) with the objective to identify and remove the variables that were highly correlated (correlation values ≥ 0.8) and less explanatory, having as a rule of decision, the relative contributions of each of them in the MAXENT model. After removing the highly correlated variables, five variables were used in the final analysis for the mistletoe (BIO6 = Minimum Temperature of Coldest Month, BIO14 = Precipitation of Driest Month, BIO15 = Precipitation Seasonality, BIO16 = Precipitation of Wettest Quarter, BIO19 = Precipitation of Coldest Quarter). The models were calibrated using the biogeographic provinces (Morrone, 2005, 2014) as hypothesis of the accessibility area of P. auriculatus, and vegetation range limits (Kuijt, 2009), which represent the species historical barriers to dispersal (Kuijt, 2009). Resulting species distribution under current climate conditions were then projected onto the LIG and onto past climate reconstructions for the MH and LGM (Braconnot et al., 2007; Hijmans et al., 2005). For future climate conditions, we projected the current conditions to the layers corresponding to two concentration pathways RCP 4.5 (intermediate emission scenarios achieving an impact of 4.5 watts per square meter by 2,100) and RCP 8.5 (hard emission scenario accomplishing an increase of 8.5 watts per square meter by 2,100) (Intergovernmental Panel on Climate Change, 2007; van Vuuren et al., 2011), from coupled atmospheric-ocean general circulation models based on the CCSM (Gent et al., 2011) and MIROC; (Watanabe et al., 2010). Changes were investigated separately for two temporal frameworks representing the 21st century: 2041 to 2060 (average in 2050) and 2061 to 2080 (average in 2070).

Final models were constructed with 10 cross-validation replicates without clamping and extrapolation, and considering the average output grids as the final predictive models. To convert initial model outputs to binary predictions, we used the equal training sensitivity and specificity threshold over the average predicted suitability across replicates and we calculated the sensitivity (the proportion of observed presences in relation to those that were predicted, which quantifies the omission errors), the specificity (the proportion of observed absences compared to those that were predicted, which quantifies the commission errors), the overall accuracy, and the TSS (True Skill Statistic; Allouche, Tsoar, & Kadmon, 2006). The TSS test corrects the overall accuracy of the model prediction by the accuracy expected by chance. This test provides a score between −1 and +1, with values > 0.6 considered good, 0.2 to 0.6 fair to moderate, and <0.2 poor (Jones, Acker, & Halpern, 2010). We also used a threshold-independent method of model validation, the receiver operating characteristic (ROC) curve analysis. The ROC analysis characterizes the predictive performance of a model at all possible thresholds by a single number, the area under the curve (AUC; Phillips et al., 2006). The value of the AUC is between 0.5 and 1.0. If the value is 0.5, the model is no better than random; an area with a value equal or close to 1 indicates a good model (Fielding & Bell, 1997).

Environmental Variation

Measurements of temperature and precipitation (BIO 1–19 variables) were extracted for each location with occurrence of P. auriculatus from WorldClim (Hijmans et al., 2005) as explained above. We performed principal components analysis (PCA) using SPSS Statistics for Mac 17 (SPSS Inc., Chicago, IL, 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. Differences in PC scores between genetic groups were tested with one-way analyses of variance in SPSS.

Haplotype Relationships

We obtained sequences of the chloroplast intergenic spacers atpB-rbcL (n = 43) and trnL-F (n = 53) for the 53 collected individuals from six populations. The aligned concatenated sequences of the chloroplast intergenic spacers (n = 43) yielded a 1,393 bp fragment and eight haplotypes. The results we report herein are for the combined data set.

Statistical parsimony retrieved a well-resolved network in which two haplogroups could be distinguished for P. auriculatus (Figure 1). The most frequent (46.5% of the individuals and 66.6% of the sampled populations) haplotype (H4) forms the core of the southern haplogroup; it was not retrieved in the southern populations. Haplotypes connected to H4 by one step consist of haplotypes exclusively found in the southern population (H1–H3, H5–H6). The second most frequent (34.9% of the individuals and 33.3% of the sampled populations) haplotype (H7) forms the core of the second group of haplotypes distributed in northern populations connected to singleton H8. Haplotypes in the northern haplogroup are connected to those in the southern haplogroup by three or more steps (Table 1, Figure 1).

Genetic Diversity and Population Structure

Differentiation among populations based on atpB-rbcL/trnL-F variation (G_ST = 0.561, SE = 0.1474) indicated that P. auriculatus is genetically subdivided. Genetic diversity across all populations (h_T = 0.793, SE = 0.0573; v_T = 0.821, SE = 0.1097) was higher than the average within-population value (h_S = 0.348, SE = 0.1211; v_S = 0.186, SE = 0.0856). PERMUT analysis showed that N_ST (0.773, SE = 0.1184) was significantly higher than the G_ST (p < .05), indicating phylogeographical structuring. When groups of populations (northern vs. southern) were compared, pairwise comparisons of F_ST value was high and significant for the atpB-rbcL/trnL-F data set (F_ST = 0.838, p < .001).

The AMOVAs showed that 85.5% for the atpB-rbcL/trnL-F was explained by differences within populations and 14.5% by differences between populations when all locations were treated as a single group (F_ST = 0.89, p < .0001). The AMOVA for the atpB-rbcL/trnL-F revealed population structure, with a high F_CT value (F_CT = 0.89, p < .0001) obtained when populations are grouped as northern and southern, suggesting that seed-mediated dispersal is limited (Table 2). Genetic diversity (h) and nucleotide diversity (π) were low for both data sets (Table 3).

Table 2.

Results of Locus by Locus AMOVA Models on Psittacanthus auriculatus Populations With no Groups Defined a priori (a), and Grouped into Two Lineages Corresponding to P. auriculatus Populations with Northern and Southern Distribution in Oaxaca (b).

Table 3.

Summary Statistics of Neutrality Tests and Demographic Analysis of Psittacanthus auriculatus Samples by Genetic Groups Using the trnL-F/atpB-rbcL Data set to Infer Demographic Range Expansion.

Demographic History

When populations are treated as two groups (northern and southern), the Tajima's D and Fu's Fs values were negative for both the northern and southern groups and for the species (considering all populations as one group) but only significant in the atpB-rbcL/trnL-F data set for the southern group, indicating deviation from the model of neutral evolution for the southern P. auriculatus population (Table 3). In the mismatch distribution, SSD and Hri results were low and nonsignificant (p > .05) in the atpB-rbcL/trnL-F data set, providing evidence for sudden demographic expansion in the past for all groups of populations. Low and nonsignificant values of Hri and SDD indicated a good fit between the observed and the expected values of the sudden expansion model (Table 3).

Ecological Niche Modeling

The ENM yielded a good fit for the current geographic distribution of P. auriculatus (AUC value = 0.884 ± 0.062) and threshold probability maps corresponded well to present species range (TSS = 0.650 ± 0.112; Figure 2). ENMs revealed that suitable habitat for P. auriculatus was more expanded and continuous under LGM conditions than the more contracted and fragmented distribution under LIG and present conditions (Figure 2). Predictions for the LGM applying both CCSM and MIROC simulations revealed that conditions of suitable habitat potentially contracted and fragmented the distribution of P. auriculatus, particularly in the western and central valleys of Oaxaca according to both the CCSM and MIROC models (Figure 2). Interestingly, conditions of suitable habitat were contracted and fragmented during the MH, with suitable habitat for P. auriculatus restricted to small areas in the western and central valleys, Tehuacán-Cuicatlán Valley, and the central portion of Puebla State. Overall, the comparison between present and past climatic conditions predicted by the models suggests that climatic conditions for P. auriculatus suitable habitat experienced range expansions/contractions in Oaxaca, with contraction and fragmentation from LIG to LGM, contraction and fragmentation from LGM to MH, and expansion from MH to present.

Figure 2.

Results from the MAXENT analyses showing binary transformation of ecological niche models for Psittacanthus auriculatus at LIG (140–120 ka), Last Glacial Maximum (LGM, CCSM and MIROC scenarios, 21 ka), Mid-Holocene (MH, CCSM, and MIROC scenarios, 6 ka), and at present. Light blue indicates binary transformation of occurrence using equal training sensitivity and specificity threshold that maximized the true skill statistics (TSS). Black dots shown on the left map below indicate occurrence records. Vegetation and land use on right below map (scale 1:250000; INEGI 2007), with white dots indicating occurrence records. Psittacanthus auriculatus in a mesquite-grassland at Santiago Matatlán, Oaxaca, Mexico.

Neither CSSM nor MIROC scenarios for 2050 and 2070 predicted major future distributional change for P. auriculatus in Oaxaca (Figure 3). Only minor changes were predicted under extreme conditions in the northern portion of its present distribution, close to northern limit of its potential range, and the region between the western and central valleys of Oaxaca (Figure 3).

Figure 3.

Predicted distribution of Psittacanthus auriculatus under two climate changes scenarios, from conservative RCP 4.5 (blue) to extreme RCP 8.5 (dark blue) conditions, for (a) 2050 and (b) 2070 with CCSM, and (c) 2050 and (d) 2070 with MIROC. Yellow represents the combined predicted areas under present (light blue) and the two concentration pathways (RCP).

Environmental Variation

The PCA indicated four niche axes that together explained 93.6% of the environmental variation in P. auriculatus (Table 4). The first niche axis (48.9% of variation) was positively associated with precipitation variables (Table 4). The second niche axis (24%) was positively associated with mean diurnal range and precipitation seasonality and negatively associated with precipitation of driest month, precipitation of driest quarter, and precipitation of coldest month (Table 4). The third niche axis (14.2%) was positively associated with temperature seasonality and temperature annual range and negatively associated with isothermality. Environmental variation between groups was only statistically different in PC3 scores (Table 4). Univariate analysis of variance of PC3 using genetic group as a fixed factor showed significant differences between P. auriculatus populations (northern and southern; PC3: F_1,52 = 14.2, p < .0001; Table 4). On average, temperature seasonality (BIO4) and temperature annual range (BIO7) values from northern locations are significantly higher than those at southern locations.

Table 4.

Factor Loadings from the Principal Components Analysis of Psittacanthus auriculatus in Oaxaca on Temperature and Precipitation Variables From WorldClim.

Genetic Diversity and Population Structure

Our survey of chloroplast sequence data in populations of P. auriculatus identified two haplogroups in Oaxaca. Phylogeographic and distributional projections derived from ENM are consistent with the hypothesis that isolation by the arid pine-oak forests along the Nudo Mixteco mountain range (Nudo Mixteco, Volcan Prieto, Cerro Piedra de Olla) may have restricted northward gene flow between the western and central valleys (northern and southern locations). Together, our results suggest a recent geographical northern expansion into the xeric vegetation of the western valleys (arid subtropical scrub, mesquite scrub). These are consistent with a model of population expansion/contraction for species responding to climate changes during the Pleistocene and isolation and reduced gene flow by both physical and environmental barriers. Given the few phylogeographic studies in the area, the phylogeography of P. auriculatus in combination with a niche-modeling approach is particularly important to understand the evolution of the Oaxaca biodiversity.

An alternative explanation is that the observed haplotype distribution in P. auriculatus may have arisen from habitat and population fragmentation. Two haplotypes of the combined atpB-rbcL/trnL-F (H7–H8) were found to be exclusive to the northern locations and six haplotypes separated by three mutational steps to the southern location (H1–H6). The lack of haplotype sharing between groups favors the hypothesis of restricted bird seed-mediated gene flow between areas. Significant differences in current environmental conditions between groups of locations are interpretable as an index of temperature variability, with temperature seasonality values (BIO4, BIO7) from northern locations significantly higher than those at southern locations. However, the fact that there are some environmental differences between areas does not necessarily point to adaptation as the driver of genetic differentiation, and thus these analyses are correlative. It is possible that for these mistletoes, the arid pine-oak forests may represent an impediment to gene flow between the western and central valleys in which arid tropical scrub and much of the central valleys are covered with scattered pockets of arid scrub and intermingled with scattered trees or shrubs often related to the surrounding habitat and mescal cultivars (Binford, 1989; Gómez-Mendoza et al., 2008). With seed dispersal mediated by highly vagile bird species (Díaz Infante et al., 2016), this hypothesis is unlikely. The genetic diversity of P. auriculatus was highly relative to that observed in more widespread Psittacanthus congeners (Ornelas et al., 2016; Pérez-Crespo et al., 2017). Because the fruits of P. auriculatus are consumed and dispersed by several bird species, mainly the widespread Tyrannus vociferans followed by the altitudinal migrant Ptilogonys cinereus (Díaz Infante et al., 2016), it is possible that the observed genetic structure and genetic differentiation has not been eroded by large-scale seed recruitment or historical rates of gene flow (Ornelas et al., 2016). Distinguishing between these possibilities will require multilocus data sets not currently available. More intensive fine-scale sampling and genotyping more loci (e.g., microsatellites) could reveal further genetic structuring linked to host usage.

Population Glacial Refugia

The levels of genetic differentiation found in P. auriculatus reveal the presence of two genetic groups, and the pattern of differentiation provides support for the predominant role of isolation and environmental factors in driving genetic differentiation between populations of P. auriculatus. The ENMs revealed a demographic history for P. auriculatus contracted its range from the LIG to the LGM, which was consistent with neutrality and population expansion tests. ENM results indicate that the range of P. auriculatus contracted and fragmented from the LIG to the LGM, in both CCSM and MIROC scenarios, and expanded from the MH to present (Figure 2). However, the ENM predicted suitable habitat and population fragmentation of P. auriculatus between western and central valleys since the LGM 21,000 years ago (Figure 2). The LIG to LGM contraction, LGM to MH contraction, and expansion from MH to present scenarios are generally consistent with the glacial refugia hypothesis (see also Cornejo-Romero et al., 2017). Based on the phylogenetic evidence and divergence time estimation suggesting that the species originated before the LIG with a split from the sister species clade composed of P. schiedeanus and P. calyculatus estimated at 2.8 to 1.77 Mya (Ornelas et al., 2016; Pérez-Crespo et al., 2017), we do interpret the genetic differentiation between P. auriculatus populations in the western and central valleys of Oaxaca as a hypothesis of population and habitat fragmentation during glacial periods.

Palynological data from the Eocene Mequitongo Formation (Ramírez-Arriaga, Prámparo, Nieto-Samaniego, & Valiente-Banuet, 2017) and the Tehuacán Formation geochronologically dated as Middle Miocene (Ramírez-Arriaga et al., 2014) suggest that conditions were more humid than those that exist today in the Tehuacán-Cuicatlán Valley, a xeric region north of the central valleys of Oaxaca. The overall palynoflora suggests that climatic conditions led to the development of temperate Pinus forest, Pinus-Quercus forest, deciduous forest, and cloud forest, which grew in mountain ranges and components of tropical deciduous forests, such as Burseraceae, Leguminosae, and Cactaceae, were also present, indicating semiarid conditions. Palynological records of the Middle Miocene Tehuacan Formation suggest a change from a relatively humid to a semiarid climate, which favored the expansion of groups currently abundant in the xerophilous scrub and low-growth seasonally dry tropical forest (Cornejo-Romero et al., 2017). This study suggests that during Late Quaternary, arid conditions persisted at LGM in both the western and central valleys of Oaxaca. Therefore, the ENM results support the xerophilous refugia hypothesis in which the paleodistribution of P. auriculatus was contracted and fragmented into multiple refugia during the LGM, while during the interglacials it was continuous. Future phylogeographic studies are needed on other ecologically important and codistributed species to determine whether the responses of those species to isolation and Pleistocene climate change are congruent with that observed in P. auriculatus.