Genetic Differentiation and Fragmentation in Response to Climate Change of the Narrow Endemic Psittacanthus auriculatus

Abstract

The state of Oaxaca is positioned in a rather unique biogeographical position with the highest diversity of vascular plants in Mexico. The isolation of xeric valleys surrounded by complex mountain ranges in Oaxaca supplies an excellent opportunity to investigate the influence of the Pleistocene events on xeric species. To test for the alternative hypotheses of Pleistocene glacial refugia, we used sequences of two chloroplast markers to examine the phylogeographic patterns of the endemic mistletoe species Psittacanthus auriculatus (Loranthaceae) across its known range in Oaxaca and conducted ecological niche modeling (ENM) to explore changes of its distribution range through present, future, and palaeo periods. Our results revealed two groups corresponding to the distribution of individuals/populations from the northern locations (western valleys), and those from southern localities at central valleys of Oaxaca. A significant genetic signal of differentiation, demographic stability, and contraction of suitable habitat during the Last Glacial Maximum predicted by ENMs strongly supported a scenario of habitat fragmentation during the Late Pleistocene. We conclude that the genetic differentiation of P. auriculatus is consistent with a model of range contraction during glacial cycles and expansion during interglacials with no major range changes under future scenarios of climate change. The findings verified the profound effects of Pleistocene climatic fluctuations on this endemic mistletoe species, and the low genetic diversity within populations highlights, paradoxically, the urgency of preserving vulnerable populations of endemic yet parasitic mistletoes.

Introduction

The Mexican highlands with a complex biogeographical, climatic, and geological history encompass one of the 25 world biodiversity hotspots (Mittermeier et al., 2005; Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000). Within these complex mountain systems, the state of Oaxaca is located in a rather unique biogeographical position with high diversity of endemic vascular plants (8.3% of the vascular flora of Oaxaca and 21% of the vascular flora of Mexico; García-Mendoza, Ordóñez, & Briones-Salas, 2004). The Valles Centrales, a region in the heart of the state of Oaxaca, is a unique geographical and cultural region consisting of several interconnected valleys surrounded by large rivers and high mountain ranges. These valleys and the adjacent Tehuacán-Cuicatlán Valley are arid and hot, with a dry-cold climate and a mean annual precipitation of 500 to 800 mm and an annual average temperature of 18℃ to 22℃ (García-Mendoza et al., 2004), and are especially rich in endemic plant species with very restricted distributions (Dávila et al., 1993; García-Mendoza et al., 2004; Sosa & De-Nova, 2012). The geological history and the glacial oscillations that occurred during the Pliocene and Pleistocene probably caused the fragmentation and long-term isolation of populations, subsequently leading to highly divergent lineages and endemism in this region (Bryson, García-Vázquez, & Riddle, 2011; Rojas-Soto, Espinosa de los Monteros, & Zink, 2007; Ruiz-Sanchez & Specht, 2013).

Integrative application of phylogeography and ENM has become a key approach for explaining intraspecific evolutionary patterns associated with climate oscillations (e.g., Bonatelli et al., 2014; Carstens & Richards, 2007; Ornelas et al., 2015; Perktas, Gür, & Ada, 2015; Peterson & Nyári, 2007; Waltari et al., 2007). Numerous studies using this approach indicate that the genetic differentiation of many extant taxa and species ranges were affected by climatic and environmental discontinuities as barriers and Pleistocene glacial/interglacial cycles, most notably in the temperate regions of northern Europe and North America (Carnaval, Hickerson, Haddad, Rodrigues, & Craig Moritz, 2009; Hewitt, 2000). Yet, Pleistocene climatic oscillations have had a profound effect on the distribution, demographic dynamics, and patterns of genetic variation of arid adapted species (Allal et al., 2011; Byrne, 2008; Meng, Gao, Huang, & Zhang, 2015; Riddle & Hafner, 2006; Turchetto-Zolet, Pinheiro, Salgueiro, & Palma-Silva, 2013). In the New World, cyclic range contractions and expansions shaped the current distribution of arid plants, their concomitant population dynamics and genetic differentiation, particularly during the Late Quaternary climatic changes (Collevatti, Rabelo, & Vieira, 2009; Rebernig et al., 2010). However, the influence of these climatic fluctuations on phylogeographic patterns of taxa occurring at the intertropical drylands remains largely unknown.

Phylogeographic studies of taxa from arid regions have revealed a pattern of lineage divergence associated to contraction to and expansion from xerophilous refugia (Bonatelli et al., 2014; Collevatti, Castro, Lima, & Telles, 2012; Collevatti, Grattapaglia, & Hay, 2003; Garrick, Nason, Meadows, & Dyer, 2009; Nason, Hamrick, & Fleming, 2002). The glacial refugia hypothesis (Haffer, 1969), which was originally conceived for tropical biogeographical realms, posits that species contracted to one or more southerly refugia during cold/dry glacial periods and expanded out from them in warm/humidity interglacials (Avise, 2000; Hewitt, 2000). The glacial refugia hypothesis is based on the assumption that populations were restricted to refugia during glacial maxima within the bounds of the current observed distribution of a species. Contrary to predictions of the glacial refugia hypothesis, recent studies using species distribution models with coalescent simulations based on DNA sequences indicate more fragmentation of suitable habitats during the Last Interglacial (LIG; 140,000–120,000 years BP) and the present than in the Last Glacial Maximum (LGM; ca. 20,000 years BP), probably due to topography, and more expansion of suitable climatic conditions during the LGM (e.g., Leite et al., 2016). Two contrasting responses of arid-adapted species to climatic fluctuation have been observed: (a) genetic divergence linked to population contractions into patches of xerophilous vegetation with warm temperatures during glacial maxima and expansions into different directions during the large-scale aridification phase of the Holocene (Clark-Tapia & Molina-Freaner, 2003; Fehlberg & Ranker, 2009; Gutiérrez-Flores, García-De León, León-de la Luz, & Cota-Sánchez, 2016; Sosa, Ruiz-Sánchez, & Rodríguez-Gómez, 2009). This response is consistent with predictions of the glacial refugia hypothesis (Haffer, 1969), which is widely accepted for temperate species in the mid latitudes of the Northern Hemisphere and (b) populations of species contracted to warm and humid refugia during the interglacials, and expanded out under the cold/dry climate of the LGM (Bonatelli et al., 2014), according to the interglacial refugia hypothesis (Haffer, 1969; Hewitt, 2004).

The distributions of organisms in xeric environments of Mexico underwent dramatic changes linked to Pleistocene climate cycles (Metcalfe, O'Hara, Caballero, & Davies, 2000; Shreve, 1942), and plant populations could have experienced repeated geographic isolation during the Pleistocene glacial/interglacial cycles approximately 2.4 to 0.01 million years ago (e.g., Angulo, Ruiz-Sanchez, & Sosa, 2012; Gándara & Sosa, 2014; Hernández-Hernández, Colorado, & Sosa, 2013; Rodríguez-Gómez & Ornelas, 2015; Ruiz-Sanchez, Rodríguez-Gómez, & Sosa, 2012; Vásquez-Cruz & Sosa, 2016). The effects of Pleistocene climate cycles between allopatric populations currently distributed in the southern portion of the Chihuahuan Desert and those distributed in the arid regions of the Valles Centrales and Tehuacan-Cuicatlán Valley of Oaxaca has been sparsely tested, with population divergence between these areas dated during the Pliocene (Gándara & Sosa, 2014) or during the Pleistocene, with demographic expansion during the LGM (Scheinvar, Gámez, Castellanos-Morales, Aguirre-Planter, & Eguiarte, 2017). A pattern of glacial survival into southern refugia followed by postglacial recolonization of northern regions seems to be ubiquitous for a variety of temperate taxa (e.g., Hewitt, 2000; Stewart, Lister, Barnes, & Dalén, 2010). The phylogeographic pattern of southern refugia for temperate species during glacial phases, which generally are confined to the southern portion of the species' distribution during warm climatic phases such as the current interglacial, is not followed by desert-adapted species (Angulo, Amarilla, Anton, & Sosa, 2017; Castellanos-Morales, Gámez, Castillo-Gámez & Eguiarte, 2016; Ruiz-Sanchez et al., 2012; Scheinvar et al., 2017). Phylogeographic patterns for studied species of the Chihuahuan Desert biota suggest the opposite pattern, a greatest range contraction during the LIG, with a later expansion to northern areas over the LGM. This pattern is accompanied with demographic expansion and a gradient of genetic diversity, with higher diversity values at southern latitudes (Angulo et al., 2017; Castellanos-Morales et al., 2016; Rebernig et al., 2010; Ruiz-Sanchez et al., 2012). Overall, results of current ENM and palaeomodeling show that the area currently occupied by studied species is substantially larger than it was during the LIG period and the LGM, suggesting that it has been gradually expanding northward since the LIG period. Previous phylogeographic studies that included plant populations of the Valles Centrales and the Tehuacán-Cuicatlán Valley have shown that their ranges contracted during interglacials and then expanded to colonize northern xeric regions during the glacial periods (the Meztitlán Valley and the Chihuahuan Desert; e.g., Loera, Ickert-Bond, & Sosa, 2017; Ruiz-Sanchez et al., 2012; Scheinvar et al., 2017) or expanded during the warm/humid LIG period and contracted to one or more southerly refugia during cold/dry glacial maxima (Cornejo-Romero et al., 2017).

Mistletoes are considered ecologically important key resources as they provide high-quality nutritional resources for a variety of obligate and partially dependent species, especially during periods of seasonal scarcity (Watson, 2001), and influence indirectly community structure in low productivity systems on bottom-up processes (e.g., Watson, 2009; Watson & Herring, 2014). Most of these parasitic plants depend on avian vectors for pollen and seed dispersal and obtain all of their water and minerals from the host (obligate hemiparasites) through a vascular connection termed haustorium (Calder & Bernhardt, 1983). Given their close connections with and dependence on their hosts (Cocoletzi, Angeles, Ceccantini, Patrón, & Ornelas, 2016; Kuijt, 2009), the geographical ranges of mistletoes are directly related to the availability of suitable host trees, and the geographic structuring of host populations might influence the genetic structuring of mistletoe populations (Norton & Carpenter, 1998; Ramírez-Barahona, González, González-Rodríguez, & Ornelas, 2017). Endemic species of mistletoes might be highly vulnerable to local extinctions, and specific conservation and management planning are needed to further develop strategies to support the persistence of multiple, interacting processes (climatic, hydrological, ecological, demographic, and genetic) of evolutionary refugia and ecological multiple interactions (e.g., host/parasite interaction) in an era of unprecedented environmental change (Davis, Pavlova, Thompson, & Sunnucks, 2013). Mistletoes infect a large number of tree species including commercially important coniferous and other hardwood timber stands (Mathiasen, Nickrent, Shaw, & Watson, 2008). Because of this, the current perception of mistletoes as invasive pests, damaging to individual trees and detrimental to forest health, and noxious weeds needs to be challenged and rather understood as an indicator of habitat health, or in superabundance as a signal of landscape perturbation (Norton & Reid, 1997; Watson, 2001).

Psittacanthus auriculatus (Oliv.) Eichler, a narrow endemic to the Mexican States of Oaxaca and Puebla (Kuijt, 2009), is the only Mesoamerican species in the genus with cordate leaves (Kuijt, 2009). In this region, the species commonly grows on Acacia schaffneri (S. Watson) F.J. Herm., and less frequently on Acacia pennatula (Schltdl. & Cham.) Benth., Acacia cymbispina Sprague & Riley and Prosopis laevigata (Humb. & Bonpl. ex Willd.) M.C. Johnst. trees in dry rocky valleys at an altitude between 1,300 and 2,000 m (Díaz Infante, Lara, Arizmendi, Eguiarte, & Ornelas, 2016; Kuijt, 2009; Pérez-Crespo, Ornelas, Martén-Rodríguez, González-Rodríguez, & Lara, 2016). P. auriculatus plants produce ca. 100 to 200 flowers each that open asynchronously from August to October (Díaz Infante et al., 2016; Pérez-Crespo et al., 2016). The self-compatible, partially protandrous red-pinked flowers (ca. 20-mm long) produce large amounts of nectar (13–30 µL/flower/day), last 2 days, and are mainly pollinated by hummingbirds (Pérez-Crespo et al., 2016). Its ripe purplish-black, fleshy (one-seed) fruits (9 × 13 mm) are consumed by several bird species (Díaz Infante et al., 2016). Local people harvest the mistletoe plants to feed goats. Currently, the ecosystems in the Valles Centrales have become seriously degraded due to natural, historical, and anthropogenic factors, posing significant challenges in the area for vegetation restoration and conservation of mostly endemic plant species. Extensive plant surveys were conducted covering the entire distribution of P. auriculatus and performed plant material collection from several populations, all of which contained a small number of individuals. Therefore, the narrow distribution of this endemic mistletoe species, small population sizes, host specificity, isolation, and exploitation by the local people constitute pressures that can increase its sensitivity to extinction.

In this study, we examine chloroplast DNA (cpDNA) sequence variation to explain genetic diversity in P. auriculatus. We also used ENM to give insight into the effects of past, present, and future climatic changes on the distribution of P. auriculatus. Specifically, we test the following hypotheses: (1) the paleodistribution of P. auriculatus populations is expected to be contracted during the LGM, while during the interglacials they were expanded northwards (glacial refugia hypothesis), (2) fragmented during the interglacials, while during the LGM they are expected to be continuous (interglacial refugia hypothesis), and (3) the paleodistribution of P. auriculatus populations is expected to be contracted and fragmented into multiple refugia during the LGM, while during the interglacials they were continuous (xerophilous refugia hypothesis).

Methods

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.

Pop. no.	Location	Altitude (m)	Latitude (N)	Longitude (W)	n	Haplotype
1	Santiago Chazumba	1709	18.09511	97.68647	7	H7(6), H8(1)
2	Huajuapan de León, El Molino	1695	17.69371	97.60154	9	H7(9)
3	Teotitlán del Valle	1680	17.08060	96.51188	5	H2(3), H4(2)
4	Santiago Matatlán, Guelavía	1782	16.84827	96.37170	6	H3(1), H4(5)
5	San Dionisio Ocotepec	1620	16.79300	96.33406	7	H4(7)
6	Sto. Domingo Tehuantepec	1400	16.52000	96.23000	9	H1(1), H4(6), H5(1), H6(1)

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.

Results

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).

	trnL-trnF/atpB-rbcL
	Sum of squares	Estimated variance	%	Fixation indices
(a) No groups defined
Among populations	31.940	0.87792	85.516
Within populations	5.502	0.14869	14.484	FST = 0.85***
Total	37.442	1.02661
(b) Two lineages
Among groups	30.504	1.49911	89.336	FCT = 0.89***
Among populations within groups	1.436	0.03026	1.803	FSC = 0.16*
Within populations	5.502	0.14869	8.861	FST = 0.91***

*p < .05, **p < .001, ***p < .0001.

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.

	Two groups
Parameter	Northern	Southern	Global
N	16	27	43
N H	2	6	8
h, gene diversity	0.1250 ± 0.1064	0.3989 ± 0.1172	0.6501 ± 0.0510
π, nucleotide diversity	0.0001 ± 0.0002	0.0005 ± 0.0004	0.0017 ± 0.0011
Tajima's D	−1.1622	−1.5975*	−0.4348
Fu's Fs	−0.7001	−3.6993**	−0.9491
SDD	0.0423	0.0019	0.0796
Hri	0.5781	0.1516	0.2225

Note. D_T and F _S positive values are indicative of mutation-drift-equilibrium, which is typical of stable populations, and negative values that result from an excess of rare haplotypes indicate that populations have undergone recent expansions, often preceded by a bottleneck. Significantly negative values (at the 0.05 level) reveal in both tests historic demographic expansion events. Significant (p ≤ .05) SSD and Hri values indicate deviations from the sudden expansion model. *p < .05; **p < .01. N = number of individuals; N _H  = number of haplotypes; SSD = differences in the sum of squares or mismatch distribution; Hri = Harpending's raggedness index.

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.

	Factor loading
Variables	PC1	PC2	PC3	PC4
Annual mean temperature (BIO1)	0.932	0.301	0.143	−0.111
Mean diurnal range (BIO2)	0.512	0.566	0.389	0.459
Isothermality (BIO3)	−0.072	0.533	−0.705	0.301
Temperature seasonality (BIO4)	−0.253	−0.095	0.931	0.026
Max. temperature of warmest month (BIO5)	0.780	0.413	0.448	−0.024
Min. temperature of coldest month (BIO6)	0.945	0.082	−0.172	−0.257
Temperature annual range (BIO7)	−0.466	0.340	0.734	0.334
Mean temperature of wettest quarter (BIO8)	0.910	0.321	0.182	−0.103
Mean temperature of driest quarter (BIO9)	0.966	0.208	0.027	−0.125
Mean temperature of warmest quarter (BIO10)	0.910	0.289	0.255	−0.119
Mean temperature of coldest quarter (BIO11)	0.941	0.306	−0.007	−0.113
Annual precipitation (BIO12)	0.813	−0.428	0.160	0.249
Precipitation of wettest month (BIO13)	0.832	−0.292	0.105	0.339
Precipitation of driest month (BIO14)	0.198	−0.943	0.087	−0.011
Precipitation seasonality (BIO15)	0.339	0.851	−0.111	0.088
Precipitation of wettest quarter (BIO16)	0.864	−0.334	0.115	0.257
Precipitation of driest quarter (BIO17)	0.245	−0.910	0.109	0.334
Precipitation of warmest quarter (BIO18)	0.595	0.088	−0.441	−0.103
Precipitation of coldest quarter (BIO19)	0.484	−0.645	−0.149	−0.125
Proportion of variance explained	48.94%	24.03%	14.21%	6.40%
ANOVA	F(1,52) = 0.4; p = .501	F(1,53) = 3.3; p = .077	F(1,53) = 27.4; p < .0001	F(1,53) = 2.3; p = .133

Note. PC scores above 0.5 indicated in bold.

Discussion

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.

Implications for Conservation

Despite their positive role as ecological keystones in many forests and woodlands (Watson, 2001; Watson & Herring, 2014) and facilitative effects of hemiparasitic mistletoes on community structure (March & Watson, 2010; Ndagurwa, Dube, & Mlambo, 2014; Watson, 2015; Watson & Herring, 2012), the loss of mistletoe species should not represent intuitively a threat to associated hosts and could even provide benefit (Strona, 2015; Strona, Galli, & Fattorini, 2013). For instance, the extinction risk of P. auriculatus is low parasitizing widely distributed host species with low vulnerability to extinction like A. schaffneri (86% prevalence; Díaz Infante et al., 2016). Also, our future projections suggest that the distribution of suitable habitat for P. auriculatus would remain similar to current conditions and might even expand its range in which given the key resources provided to pollinators and seed dispersers and its role as soil facilitator, even in a situation of extreme drought.

The results of the study are interesting from a conservation standpoint, as we provided the first estimate of genetic variability and structure of the short-range endemic P. auriculatus. The northern and southern divergent populations are in good congruence with the hypothetical refugia outlined by ENM results. The results exhibit a strong genetic differentiation, which is obviously due to genetic drift in the isolated populations. Our study therefore gives evidence for glacial relict endemism of P. auriculatus in the aridlands of Oaxaca. The narrow distribution of this endemic mistletoe species, small population sizes, isolation, and exploitation by the local people constitute pressures that can increase its sensitivity to extinction. Despite the expanded suitable conditions predicted under future climate change scenarios, the Valles Centrales have become seriously degraded due to natural, historical, and anthropogenic factors, posing significant challenges in the area for vegetation conservation of glacial relict endemism. Our results, while preliminary, illustrate the utility of an integrative approach for identifying geographical regions that harbor evolutionarily distinct populations and caution against using only the distributional patterns of endemics to prioritize regions for conservation. While incorporation of additional species from unrelated taxa will be necessary to draw definitive conclusions about evolutionarily distinct regions, our results suggest a conservation approach for the xeric areas of Oaxaca that would emphasize protection of both northern and southern valleys, thereby maximizing preservation of intraspecific genetic diversity and protection of endemic species.