Nuclear phylogeography of the temperate tree species Chiranthodendron pentadactylon (Malvaceae): Quaternary relicts in Mesoamerican cloud forests
BMC Evolutionary Biology volume 20, Article number: 44 (2020)
The Mexican hand tree or Canac (Chiranthodendron pentadactylon) is a temperate tree species of cloud and pine-oak forests of southern Mexico and Guatemala. Its characteristic hand-shaped flower is used in folk medicine and has constituted the iconic symbol of the Sociedad Botánica de México since 1940. Here, the evolutionary history of this species was estimated through phylogeographic analyses of nuclear DNA sequences obtained through restriction site associated DNA sequencing and ecological niche modeling. Total genomic DNA was extracted from leaf samples obtained from a representative number (5 to 10 per sampling site) of individuals distributed along the species geographic range. In Mexico, population is comprised by spatially isolated individuals which may follow the trends of cloud forest fragmentation. By contrast, in Guatemala Chiranthodendron may constitute a canopy dominant species near the Acatenango volcano. The distributional range of this species encompasses geographic provinces separated by the Isthmus of Tehuantepec.
The objectives of the study were to: (i) estimate its genetic structure to define whether the observed range disjunction exerted by the Isthmus of Tehuantepec translates into separate populations, (ii) link population divergence timing and demographic trends to historical climate change, and (iii) test hypotheses related to Pleistocene refugia.
Patterns of genetic diversity indicated high levels of genetic differentiation between populations separated by the Isthmus. The western and eastern population diverged approximately 0.873 Million years ago (Ma). Demographic analyses supported a simultaneous split from an ancestral population and rapid expansion from a small stock approximately 0.2 Ma corresponding to a glacial period. The populations have remained stable since the LIG (130 Kilo years ago (Ka)). Species distribution modelling (SDM) predicted a decrease in potential distribution in the Last Interglacial (LIG) and an increase during the Last Glacial Maximum (LGM) (22 Ka), Mid-Holocene (6 Ka) and present times.
Divergence time estimations support the hypothesis that populations represent Quaternary relict elements of a species with broader and northernmost distribution. Pleistocene climatic shifts exerted major influence on the distribution of populations allowing dispersion during episodes of suitable climatic conditions and structuring during the first interglacial with a time period length of 100 Kilo years (Kyr) and the vicariant influence of the Isthmus. Limited demographic expansion and population connectivity during the LGM supports the moist forest hypothesis model.
Chiranthodendron pentadactylon Sessé ex. Larreat is an evergreen temperate tree species of tropical montane cloud forests and pine-oak forests of southern Mexico  (Guerrero, Oaxaca and Chiapas) and Guatemala  (El Progreso, Zacapa, Sacatepéquez, Chimaltenango, Sololá, Totonicapán, Quetzaltenango, Quiché, Huehuetenango, and San Marcos). It is distinguished by dark reddish flowers in which the androecium portion that harbors the stamens resembles the fingers of a hand hence the vernacular names of Mexican hand tree or tree of the little hands in Mexico and Canac (hand-flower) in Guatemala. The flowers are currently used in folk medicine to treat heart conditions, high pressure  and gastrointestinal disorders such as diarrhea and dysentery . Furthermore, the flower has constituted the iconic symbol of the Sociedad Botánica de México since the early 40s. In Mexico, the population spatial structure is constituted by spatially isolated individuals which may follow the trends of cloud forest fragmentation. On the other hand, Véliz  observed that C. pentadactylon represents a canopy dominant species of the cloud forest community found near the Acatenango volcano in Guatemala.
In the IUCN Red List of Threatened Species , C. pentadactylon is listed as vulnerable in Mexico and as near threatened in Guatemala . The Norma Oficial Mexicana (NOM-059-SEMARNAT-2010)  however has listed C. pentadactylon as threatened. This implies that populations could become endangered in the short or medium term if the factors that negatively affect their viability continue to operate including deterioration or modification of their habitat or direct reduction of their population sizes .
The monospecific genus Chiranthodendron Larreat. along with Fremontodendron Coville, form the Fremontodendreae tribe within the Bombacoideae subfamily of the Malvaceae family. Fremontodendron bears three species: F. californicum (Torrey) Coville, F. mexicanum Davidson and F. decumbens R. M. Lloyd. They are shrubs and small trees that constitute distinctive elements of the chaparral vegetation in the Sierra Nevada foothills and coastal ranges of California extending to scattered locations in central and western Arizona and northern Baja California in Mexico . Chiranthodendron and Fremontodendron share leaf-opposed hermaphroditic flowers greater than 2 cm in diameter and the production of large amounts of nectar .
The high biotic diversity of Mesoamerican forests [9, 10] is the result of the interchange between Nearctic and Neotropical floristic components [10,11,12,13]. This interchange has taken place during different time periods of climate change and depends on both historical and ecological factors. For instance, pollen records from the Elsinore flora in San Diego, suggest a cooling and drying trend beginning in the middle Eocene (52–40 Ma) and Late Eocene (40–38 Ma)  to early Miocene (23.3–16.3 Ma)  which resulted in the establishment of a winter-dry mixed deciduous hardwood-coniferous forest in the uplands, and the replacement of the Paleocene (65–56.5 Ma) to earliest Lutetian (48.6 Ma) tropical and paratropical lowland vegetation by more seasonally dry communities.
The drying trend is evident from pollen records suggesting the last appearance of Chiranthodendron-Fremontodendron potential predecessors from the early Lutetian. Hence, species with north temperate affinity constitute forms that survived the progressive drying and cooling trends of the middle Eocene and early Miocene. Subsequent colonization events have ensued during the Pliocene epoch (5.33–2.59 Ma) whence the cooling trend subsided gradually culminating with higher temperatures and wetter climatic conditions in the mid-Pliocene (3.6 Ma). Thereupon, during the late Pliocene (3 Ma) temperatures decreased leading to the onset of the Pleistocene (2.59–0.01 Ma) glacial cycles with the formation of the North American Laurentide sheet approximately 2.6 Ma . These have been regarded as more recent drivers of species diversification and plant distribution dynamics [16,17,18,19,20,21,22]. Species responses to the Pleistocene ice ages include extinction over a large part of their range, dispersion to new locations, survival in refugia and postglacial colonization . The extent of change in species ranges depended on latitude and topography. High latitudes were covered with ice or permafrost while temperate and tropical regions were compressed towards the equator .
The genetic and demographic consequences suggested for Mesoamerican cloud forest species have been interpreted under refugia hypotheses, which describe differing precipitation regimes. These include the moist forest hypothesis and the dry refugia model. The moist forest hypothesis states that during glacial periods there was no significant change in precipitation allowing the development of population connectivity, downslope migration and range expansion . The degree of genetic structure homogenization and population growth is positively correlated with the species dispersal abilities . During dry and warm interglacial periods there was upslope migration, range contraction, isolation and genetic differentiation . Conversely, the dry refugia model suggests a significant reduction in precipitation during glacial periods in which there was downslope migration but no population connectivity or gene flow among populations. Populations were compressed into refugia by the opposing effects of aridity and cooling subsequently becoming isolated and genetically differentiated. During humid and warm interglacial periods, these populations expanded and recolonized their former range. The degree of population differentiation depended on the extent of former refugia and on their past levels of gene flow .
Presently, Mesoamerican cloud forests take place in narrow altitudinal zones  and occur in restricted and limited areas . They represent evergreen forests  characterized by a constant, frequent or seasonal cloud cover . The geographic range of Chiranthodendron covers several geographic provinces separated mainly by the Isthmus of Tehuantepec and may represent a potential physical and ecological barrier  to dispersal and gene flow. Therefore, lineages west and east of the Isthmus might constitute separate populations. However, their current distribution may be the result of the influence of past climate changes as opposed to mainly the potential vicariant effect of the Isthmus. Furthermore, the long lifespan of trees hampers their rapid adaptation to environmental changes  which make them particularly vulnerable to climate changes influencing their population dynamics, genetic diversity patterns and demographic history.
To gain insight into the potential historical factors influencing C. pentadactylon present distribution, the evolutionary history is examined through phylogeographic, population genetic and ecological niche modelling approaches. The main objectives of the study were to use C. pentadactylon to: (i) estimate the genetic structure to define whether the observed range disjunction exerted by the Isthmus of Tehuantepec translates into separate populations, (ii) link population divergence timing and demographic trends to historical climate change, and (iii) test hypotheses related to Pleistocene refugia.
Genetic analyses were conducted using nuclear loci obtained through restriction-site associated DNA sequencing (RADSeq). Nuclear genomic DNA may be the source of numerous unlinked loci of enough resolution to detect spatial distributional patterns of gene genealogies . Moreover, nuclear markers are biparentally inherited reflecting both seed and pollen movement  as opposed to the traditionally used chloroplast DNA markers [11, 30] which reflect only seed movement  and the spatial genetic diversity patterns are thus only partially analyzed.
Genetic diversity and genetic differentiation
The overall nucleotide diversity (Pi) was low (0.01463). The STRUCTURE analysis yielded a K value of 2 for the number of groups (Fig. 1) with the highest likelihood (Additional file 5). The first group includes the individuals sampled west of the Isthmus of Tehuantepec (Guerrero and Oaxaca) and the second group includes the individuals sampled east of the Isthmus (Chiapas and Guatemala). This arrangement of groups was further introduced in Arlequin to evaluate the degree of genetic differentiation through a global analysis of molecular variance (AMOVA). The AMOVA (Table 1) showed that most of the variation is attributed to differences among groups (62.87%) with a fixation value (FST) of 0.62879 followed by moderate levels of within population variation (37.12%). The estimated degree of genetic differentiation among groups separated west and east of the Isthmus supports their attribution as separate populations referred to as western and eastern population respectively.
Phylogenetic network and divergence times
The phylogenetic network is shown in Fig. 2 and shows the same set of populations inferred by STRUCTURE. The StarBEAST2 tree (Fig. 3) showed a strong support (PP = 1) for the differentiation of the clades west and east of the Isthmus of Tehuantepec. The western and eastern population diverged approximately 0.873 Ma (HPD 95%: 0.693–1.022 Ma).
The BSPs estimated for the SNPs sequences indicated an ongoing increase in effective population size from 0.5 Ma onwards followed by sudden demographic expansions. The western population (Fig. 4a) underwent one sudden demographic expansion approximately 0.2 Ma while the eastern population (Fig. 4b) two sudden demographic expansions, the first approximately 0.4 Ma and the second approximately 0.15 Ma. After the latest expansions, effective population sizes have remained stable until recent times.
Considering the five scenarios (Fig. 5) tested for the SNPs sequences, DIYABC analysis indicated that a simultaneous split from an ancestral population (scenario 1) is the best-supported scenario with a higher posterior probability value than estimated for other scenarios (Table 2). Confidence in scenario choice evaluation based on 500 PODs yielded low values for the Type I and Type II errors.
The mean area under the receiver operating characteristic curve (ROC) for the replicate runs displays a good model performance (0.970). The species distribution modelling (SDM) for the LIG shows a major reduction in the species ecological niche as compared with the SDM for the LGM and Mid-HLC in which the species ecological niche is increased. The potential geographic distribution is maintained in the LGM through the Mid-HLC but is further increased in recent times (Fig. 6).
The vicariant role of the isthmus of Tehuantepec in shaping patterns of genetic diversity
The distributional range of C. pentadactylon encompasses geographic provinces separated by the Isthmus of Tehuantepec. The Isthmus of Tehuantepec is a savannah-like valley area  of high seismicity that unites the North American continent with Nuclear Central America. At the North American continent junction decreases to 200 Km width . Geological evidence from the late Miocene to late Pliocene suggests a reduction of the highlands and possible marine embayment as consequence of an extensive down dropping of the eastern block of the Tehuantepec fault [32, 33]. The Isthmus of Tehuantepec has been regarded as a biogeographical barrier to many taxa. Studies on Mesoamerican cloud forest species currently codistributed alongside and near the Isthmus of Tehuantepec in Mexico, observed a genetic and distributional range disjunction produced by the occurrence of the Isthmus. For instance, phylogeographic studies of the distylous shrub Palicourea padifolia  and of the mistletoe cactus Rhipsalis baccifera which followed the trends of cloud forest demographic dynamics during periods of climate change  revealed a genetic and population disjunction at the Isthmus. Moreover, a comparative phylogeographic study that included the cloud forest tree species Podocarpus matudae and Liquidambar styraciflua, the shrub P. padifolia, the herb Moussonia deppeana and the epiphytic Rhipsalis baccifera as well as the rodent species Peromyscus aztecus, Reithrodontomys sumichrasti, and Habromys lophurus, the hummingbird species Amazilia cyanocephala, Campylopterus curvipennis and Lampornis amethystinus, the woodcreeper Lepidocolaptes affinis and the passerine species Chlorospingus opthtalmicus, Buarremon brunneinucha and Basileuterus belli  likewise supported a genetic and population disjunction at the Isthmus. Conversely, in the mistletoe Psittachanthus schiedeanus  the vicariant role of the Isthmus was not evident due to effective gene flow mechanisms resulting in shallow levels of genetic structure.
To determine whether the distribution disjunction exerted by the Isthmus of Tehuantepec represents a biogeographical barrier to dispersal and gene flow, population estimates were conducted based on nuclear data derived from RAD sequencing. The degree of population substructure was estimated under the Bayesian approach implemented by the software STRUCTURE. The Bayesian algorithm detected a genetic differentiation among C. pentadactylon lineages distributed west and east of the Isthmus of Tehuantepec and were thenceforth considered as separate populations. The degree of genetic differentiation among populations was further estimated by the global AMOVA implemented by Arlequin which indicated high levels of genetic differentiation among populations and moderate levels of within population variation. The high vagility of the pollination vectors such as birds and bats might induce low levels of among population variation. Nevertheless, the probability of long-distance dispersal events along the Isthmus is reduced given the hermaphroditic floral morphology and therein breeding system. On the other hand, the dichogamous nature of the flowers promotes a mixed-mating mechanism which derives in moderate levels of within population genetic variation. Further, low levels of within population subdivision are expected given the highly vagile pollination vectors such as birds and bats since they are likely to visit numerous plants . Seed dispersal mediated by birds might also result in low levels of within population subdivision.
Finally, the relationships among individuals estimated by the Neighbor-Net algorithm consisted on the two populations detected by previous analyses. Moreover, the StarBEAST2 tree showed a strong support for the differentiation of the populations.
Population divergence timing and demographic history relations to historical climate change
Changes in the distribution of cloud forests throughout time have been determined by factors such as geological events and climate change. In addition to the occurrence of the Isthmus of Tehuantepec, the geographic range of C. pentadactylon might have been influenced by the climatic oscillations of the Quaternary ice ages. To relate genetic divergence to pre-Pleistocene and Pleistocene events, calibrated time trees were calculated under the multispecies coalescent model implemented in StarBEAST2. Chiranthodendron and its sister genus Fremontodendron bring together the Fremontodendreae tribe within the Bombacoideae subfamily. These might have derived from predecessors with Laurasian origin  which migrated to North America and subsequently to Mesoamerica . According to Richardson et al. , the MRCA of the Fremontodendreae tribe dates to approximately 5 Ma in the early Pliocene. Using this date as calibration point and 36 polymorphic nuclear loci, populations divergence was estimated to have a median age of 0.873 Myr. This date approximates to the latest stage of a slow cooling progress developed since the last 50 Myr. At the start of the Pleistocene epoch, alternations between cold glacial periods and warmer intervals took place and between 1.2 and 0.6 Ma weaker cycles with a period of 40 Kyr gave way to stronger cycles of approximately 100 Kyr mostly beginning 800 Ka . The extended warming of this time period may have reduced the potential distribution of Chiranthodendron and contracted its range to the highlands where suitable climatic conditions remained. The DIYABC analysis based on SNPs sequences supported a simultaneous split from an ancestral population scenario. This scenario suggests a once widespread ancestral species that contracted its range during a period of unsuitable climatic conditions as those inferred for the last 800 Kyr leading to variation in population continuity and isolation. Range contraction might have included a northernmost distribution. This is supported by the upheld distribution of its sister genus Fremontodendron in southwestern United States and northwestern Mexico. This genus is comprised by shrubs and small trees of chaparral communities  and dry temperate sclerophyllous floras . Some predecessors might have survived the climatic shifts of the Pleistocene epoch and adapted successfully to these climatic conditions. While some dispersed to northwestern Mexico possibly in the early Pleistocene during glacial intervals from which low-pressure systems allowed a deeper penetration of cyclones from the Pacific Ocean and the Gulf of Mexico . Later, some predecessors may have reached southern Mexico and Guatemala. The distribution might have become restricted to a southernmost distribution after more severe cold periods beginning approximately 1.1 Ma  and the subsequent warming and drying trend of the late Holocene towards modern climatic conditions in northern Mexico .
Changes in population size over time were estimated through BSPs based on SNPs sequences. These analyses detected ongoing increases in effective population size from 0.5 Ma onwards followed by sudden demographic expansions. The latest demographic expansions might coincide with a glacial period dated approximately 0.2 Ma . These expansions support the moist forest hypothesis in which downslope migration, range expansion and population connectivity are predicted . High dispersal abilities might have allowed increases in population sizes  as opposed to the traditionally assumed little or no demographic expansion under this model. On the other hand, the low level of nucleotide diversity observed indicates demographic expansions from small size populations. Conversely, the absence of population size reductions or bottlenecks does not support their survival in refugia as expected under the dry refugia model in which populations were compressed into refugia by the opposing effects of aridity and cooling .
Phylogeographic studies on additional cloud forest species allowed the recognition of signatures of rapid demographic expansion. For instance, haplotype diversity and topology estimations coupled with neutrality tests suggested a rapid increase in population size from a small stock in the distylous shrub P. padifolia . However, the high variation within populations as expected from species with high outbreeding rates, was not attributed to their survival in the small refugia proposed by Toledo  during the coldest periods of the Pleistocene epoch. Further studies on the geographic structure of genetic variation of highland populations of Podocarpus spp., distributed in Mesoamerica along with the reconstruction of potential distributions during the LIG and LGM  allowed the recognition of demographic signatures related to the long-term in situ permanence in multiple refugia during the LGM which comprise the absence of demographic expansion and limited gene flow. The moist forest hypothesis model has also been inferred for L. styraciflua , R. baccifera , P. schiedeanus  and Podocarpus spp.  of which predicted distribution suggests a continuous spatial trend during the LGM with no significant demographic expansion.
In C. pentadactylon, reduction in potential distribution range predicted for the LIG may have prevented further demographic expansions but nevertheless reductions in effective population sizes were not observed. The latter suggests that populations might have survived in interglacial refugia and were able to maintain their effective population sizes. The absence of demographic expansions predicted for the LGM coincides with the demographic trends expected under the moist forest hypothesis.
The populations have remained constant from the LGM to present times even at the onset of an increase in potential distribution, suggesting the influence of additional variables such as ecological and anthropic factors.
The phylogeographic approach conducted allowed the recognition and attribution of the Isthmus of Tehuantepec as a current driver of population differentiation and of the climatic oscillations of the Pleistocene epoch as a determining factor for their current distribution. Demographic analyses and divergence timing estimation suggest that populations derived from an ancestral population possibly of more widespread and northernmost distribution that underwent range contraction around the onset of the first interglacial with a time period length of 100 Kyr (800 Ka). Estimation of changes in population size over time suggests a demographic expansion during the glacial period prior to the LIG. The genetic diversity measures suggest an expansion from a small size population. The SDMs predicted a reduction in potential distribution for the LIG. However, population size reductions were not detected. Populations have not undergone demographic expansions since the unsuitable climatic conditions of the LIG and appear to remain stable. Nevertheless, the low genetic nucleotide diversity, low effective population sizes and restricted distribution to specific biomes limited in spatial range, make C. pentadactylon a species with limited adaptability to recent climate change and global warming.
Furthermore, the reduction of forest areas due to anthropic influence as established by the IUCN and the Norma Oficial Mexicana (NOM-059-SEMARNAT-2010), could render it endangered and consequently at risk of extinction. The disjunct distribution between Chiranthodendron and its sister genus Fremontodendron as well as the recent divergence time estimated for the extant populations, suggest that it represents a relict species and thus bears an additional value for conservation.
Additional ecological and evolutionary studies are required to determine the impact that C. pentadactylon provides to forest ecosystems, and to determine the characters subject to selection. This information may provide useful insights into the reproductive requirements and success of the species to efficiently conduct conservation efforts.
Chiranthodendron pentadactylon is a tree of up to 30 m tall (Fig. 7a) with a minimum generative growth time of more than 3 years . It is resistant to low temperatures (down to 5 °C) and dry conditions . The dichogamous red flowers (Fig. 7b) can be cross-pollinated or geitonogamously pollinated . Pollination is mediated by birds  and bats . The flowering period spans from November to April and the fruiting period from April to May . Fruits are dehiscent capsules (Fig. 7c) and the strophiolated seeds  dispersal is potentially mediated by birds .
Sampling and DNA sequencing
Leaf samples from 5 to 10 individuals per sampling site found along the range of C. pentadactylon in Mexico and Guatemala (Fig. 8) were collected and stored and dried with silica gel at room temperature until further processing. Two additional leaf samples from Fremontodendron californicum and F. mexicanum were collected from the Rancho Santa Botanic Garden in Claremont, California to be incorporated as outgroup species. No special permission was required to collect the samples. Taxon identification and geographic corroboration were performed during fieldwork after reviewing voucher specimens deposited in the herbarium of the Universidad Autónoma de Aguascalientes (Herbario UAA), the Instituto de Ecología AC at Pátzcuaro, the Facultad de Agronomía de la Universidad de San Carlos de Guatemala (AGUAT), the Escuela de Biología, USAC (Herbario BIGU), the Missouri Botanical Garden (MO) and the Chicago Natural History Museum. Each later defined population has an accompanying voucher deposited at the Herbarium UAA (Additional file 1).
Total genomic DNA was extracted following the cetyltrimethylammonium bromide (CTAB) protocol  with an additional incubation step with the enzyme pectinase to remove the excess of carbohydrates. Thereupon, a RAD sequencing library was prepared at Rancho Santa Ana Botanic Garden (Claremont, California) and at the University of California Riverside as described by Etter et al.  and modified by Medina  using SbfI-HF restriction enzyme (New England Biolabs). Samples were pooled into a multiplexed library and sequenced on an Ilumina NextSeq500 platform to generate 150-bp single-end readings.
The raw sequence data were analyzed and assembled de novo in the software pipeline ipyrad v0.7.28 , which identifies two sequences as orthologous as determined by a specified degree of similarity (0.90). Filtering parameters were set to trim all reads to a length of 140 bp. Reads were clustered allowing a maximum number of ambiguous base calls (Q < 20) of 0 and a maximum depth of coverage of 10 k. To date, karyotype information and ploidy level have not been determined, therefore loci containing more than two alleles were discarded to exclude potential paralogs.
When clustered across samples a minimum number of samples per locus were set to 66 to maximize the representativeness of the total samples per locus. The samples had an average of 565 K reads that passed quality filtering. These clustered into an average of 8837 clusters per sample with a mean depth of 4526 producing a mean of 4075 consensus sequences per sample and a maximum total of 70 loci were recovered (Additional file 2).
Loci retrieved from the RADSeq genome-partitioning approach are mainly from the nuclear genome . However, 5 loci were highly similar to mitochondrial loci from closely related species as identified through the basic local alignment search tool (BLAST) option from the National Center for Biotechnology Information (NCBI) database and were discarded from further analysis, resulting in a total of 65 nuclear loci. Mitochondrial DNA substitution rate is at least 5 times slower than nuclear genes  and thence may not provide the enough resolution to detect infraspecific historical processes. Additional data files recovered from the alignment include the variant call format (vcf) file and the single-nucleotide polymorphism (SNP) file.
Genetic diversity and genetic differentiation
The genetic diversity was estimated through the nucleotide diversity (Pi) measure proposed by Nei (1987) and implemented by the software DnaSP v6 . The estimation was based on the sequence variation information stored in the variant call format (vcf) retrieved from the ipyrad alignment.
The underlying degree of genetic structure was determined through the Bayesian clustering approach implemented by the software STRUCTURE v2.3.4 . The substructure was estimated from the set of 65 nuclear loci. Posterior probabilities for the ancestry model and allele frequency model parameters were estimated under the admixture model and the correlated allele frequencies model respectively. The MCMC process estimation was conducted using a burn-in of 1,000,000 followed by 100,000 chains for each of the 10 replicas assigned to a range of 1 to 9 K assumed groups. The most probable K value was estimated using the Evanno method  implemented in STRUCTURE HARVESTER  and the corresponding individual membership coefficients matrix was used as input file in STRUCTURE PLOT  to visualize the degree of genetic differentiation among groups.
Genetic structure was further evaluated through a global analysis of molecular variance (AMOVA) using the software Arlequin v3.5 . The analysis was performed based on a weighted average over all loci and confidence intervals for the fixation index were calculated by bootstrapping with 20,000 replicates.
Phylogenetic network and divergence times
To visualize the phylogenetic relationships among and within populations, a phylogenetic network was constructed using the Neighbor-Net algorithm implemented by the software SplitsTree5 . The relationships were estimated using the linked SNPs file. The Neighbor-Net method consists of the agglomeration of weighted collection of splits or partition of the set of taxa which constitute the building blocks of a phylogenetic tree and provides the visualization of the space of feasible trees. Thereof, constitutes a useful tool for the representation of the relationships of taxa in which the underlying evolutionary history may not be treelike due to processes such as recombination, hybridization, gene conversion and gene transfer .
To relate the genetic divergence among populations to pre-Pleistocene and Pleistocene events, calibrated time trees were calculated under the multispecies coalescent model implemented in StarBEAST2 . StarBEAST2 models the incomplete lineage sorting process between and within species with no recent history of gene flow  and allows variation in substitution rates of different genes and species by using gene tree relaxed clock models . The time trees were calculated from a selection of 36 polymorphic nuclear loci (Additional file 3) of approximately 140 bp. The best DNA substitution model (Table 1 in Additional file 4) and gamma rate heterogeneity for each locus were determined using jModelTest v2.1.10  through the implementation of the Akaike information criterion (AIC). The corresponding site models were set with four rate categories and an uncorrelated lognormal molecular clock was established for the clock model. The constant populations model was selected for the population model parameter and the Calibrated Yule model for the tree prior. The time trees were calibrated using a secondary calibration point due to the absence of paleontological information to constrain a minimum age of divergence for C. pentadactylon. This calibration point derives from a dated phylogeny for the Malvaceae family, which was used to infer divergence dates and diversification rates within the Theobroma genus . The resulting estimated age of divergence between C. pentadactylon and its sister genus Fremontodendron was on the order of 5 Ma. This was set as the median age for the informative lognormal prior with a standard deviation of 0.01. The MCMC chain length was 180,000,000 and logged every 5000. The log files were analyzed with Tracer v1.7.1 for parameter convergence. The burn-in was set to remove 10 % of the tree files before the generation of a maximum clade credibility tree with median heights in TreeAnnotator.
Demographic and evolutionary history
To assess changes in population size over time and determine their correspondence to pre-Pleistocene or Pleistocene events, a Coalescent Bayesian Skyline Plot (BSP) was estimated for each group defined in STRUCTURE and later defined as separate populations, using BEAST v.2.5.2. Linked SNPs data were set as input files. The best substitution model for each set of SNPs was determined using jModelTest through the implementation of the AIC (Table 2 in Additional file 4). The corresponding site models were set with four rate categories and a strict molecular clock was established for the clock model. The time scale was calibrated with the divergence age of 0.873 Ma previously estimated and a standard deviation of 0.2 under a lognormal distribution. The MCMC chain length was set to 10,000,000 and to 15,000,000 for the western and eastern population respectively and logged every 1000. The log files were analyzed in Tracer for parameter convergence.
The population history of C. pentadactylon was inferred from an approximate Bayesian computation (ABC) in which evolutionary scenarios were simulated and compared through posterior probabilities using the DIYABC v2.1.0  software. The linked SNPs data file was used to simulate 1 million datasets per scenario. The evolutionary scenarios were built considering the genetic structure and StarBEAST2 analyses: (i) two populations (Pop1 and Pop2) have diverged simultaneously from an ancestral population at t1 (scenario1), (ii) Pop2 (Eastern population) diverged from Pop1 (Western population) at t1 (scenario 2), (iii) Pop1 diverged from Pop2 at t1 (scenario 3), (iv) Pop2 derived from few western immigrants at t2 and diverged at t1 (scenario 4, and (vi) Pop1 derived from few eastern immigrants at t2 and diverged at t1 (scenario 5).
Posterior probabilities of scenarios were assessed using a polychotomic weighted logistic regression on 1% of the simulated datasets. The fit between the simulated and observed datasets was evaluated through the implementation of a model checking procedure based on a principal component analysis (PCA). Confidence in scenario choice was assessed by means of the simulation of 500 pseudo-observed datasets (PODs) under each scenario to estimate Type I and Type II error rates.
The potential distribution of C. pentadactylon for the present day, Mid-Holocene (Mid-HLC, 6 Kyr), Last Glacial Maximum (LGM, 22 Kyr) and Last Interglacial (LIG, 130 Kyr) was predicted using Maxent 3.4.1 . The georeferenced sampling locations were used as occurrence data in addition to data retrieved (15 March 2017) from the Global Biodiversity Information Facility (GBIF; Chiranthodendron pentadactylon Larreat in GBIF Secretariat (2019). GBIF Backbone Taxonomy. Checklist dataset https://doi.org/10.15468/39omei), TROPICOS (http://legacy.tropicos.org/Name/3900594, Missouri Botanical Garden, St. Louis, MO, USA; available from: http://www.tropicos.org) and REMIB (http://www.conabio.gob.mx/remib/doctos/remib_esp.html) databases and represented 152 unique localities after data cleansing including the removal of duplicate points. The models were built from 9 climate layers selected through the jackknife method. These had a 2.5 arc minute spatial resolution and were acquired from the WorldClim database . Layers for the LGM and the Mid-HLC were obtained from two paleoclimate models (CCSM and MIROC) that simulate differing climatic predictions . Distribution models were built with 10 replicates using the default settings and 25% of the occurrence data was selected for model validation.
Availability of data and materials
The data that support the findings of this study is openly available in the Dryad Digital Repository: https://doi.org/10.5061/dryad.dfn2z34w9
Approximate Bayesian Computation
Akaike information criterion
Analysis of molecular variance
Bayesian skyline plot
Community climate system model
Global biodiversity information facility
Last glacial maximum
Markov Chain Monte Carlo
Model for interdisciplinary research on climate
Most recent common ancestor
National Center for Biotechnology Information
Restriction-site associated DNA sequencing
Receiver operating characteristic curve
Species distribution modeling
González-Espinosa M, Meave JA, Lorea-Hernández FG, Ibarra-Manríquez G, Newton AC. The red list of Mexican cloud Forest trees. Cambridge: Fauna and Flora International; 2011.
Véliz PME. Caracterización de la comunidad de Canac (Chiranthodendron pentadactylon Larreategui) en el volcán de Acatenango [undergraduate dissertation]. Guatemala: Universidad de San Carlos de Guatemala. Facultad de Agronomía; 1989.
Osuna-Fernández R, Laguna-Hernández G, Brechú-Franco A, Orozco-Segovia A. Germinación de Chiranthodendron pentadactylon Larr. (Sterculiaceae), en respuesta a la escarificación, temperatura y luz. Bol Soc Bot México. 1997;60:5–14.
Velázquez C, Calzada F, Esquivel B, Barbosa E, Calzada S. Antisecretory activity from the flowers of Chiranthodendron pentadactylon and its flavonoids on intestinal fluid accumulation induced by Vibrio cholerae toxin in rats. J Ethnopharmacol. 2009;126:455–8.
Vivero JL, Szejner M, Gordon J, Magin G. The red list of trees of Guatemala. Cambridge: Fauna and Flora International; 2006.
SEMARNAT (Secretaría de Medio Ambiente y Recursos Naturales). NORMA Oficial Mexicana NOM-059-SEMARNAT-2010. Protección ambiental-Especies nativas de México de flora y fauna silvestres-Categorías de riesgo y especificaciones para su inclusión, exclusión o cambio-Lista de especies en riesgo. Diario Oficial de la Federación, México; 2010.
Kelman W, Broadhurst L, Brubaker C. Genetic relationships among Fremontodendron (Sterculiaceae) populations of the Central Sierra Nevada foothills of California. Madroño. 2006;53(4):380–7.
Bayer C, Kubitzki K. Malvaceae. In: Kubitzki K, Bayer C, editors. The Families and Genera of Vascular Plants. Springer Editorial; 2003. p. 225–311.
Ruíz-Sánchez E, Ornelas JF. Phylogeography of Liquidambar styraciflua (Altingiaceae) in Mesoamerica: survivors of a Neogene widespread temperate forest (or cloud forest) in North America? Ecol Evol. 2014;4(4):311–28.
De Albuquerque FS, Benito B, Beier P, Assunção-Albuquerque MJ, Cayuela L. Supporting underrepresented forests in Mesoamerica. Nat Conserv. 2015;13:152–8.
Luna-Vega I. Aplicaciones de la biogeografía histórica a la distribución de las plantas mexicanas. Rev Mex Biodivers. 2008;79:217–41.
Morrone JJ. Fundamental biogeographic patterns across the Mexican transition zone: an evolutionary approach. Ecography. 2010;33:355–61.
Corrales L, Bouroncle C, Zamora JC. An overview of forest biomes and ecoregions of Central America. In: Chiabai A, editor. Climate change impacts on tropical forests in Central America. London: Routledge editorial; 2015. p. 17–38.
Frederiksen NO. Pulses of middle Eocene to earliest Oligocene climatic deterioration in Southern California and the Gulf coast. SEPM. 1991;6(6):564–71.
Graham A. Late cretaceous and Cenozoic history of north American vegetation. New York: Oxford University Press; 1999.
Van der Hammen T. The Pleistocene changes of vegetation and climate in tropical South America. J Biogeogr. 1974;1(1):3–26.
Webb T III, Bartlein PJ. Global changes during the last 3 million years: climatic controls and biotic responses. Annu Rev Ecol Syst. 1992;23:141–73.
Flenley JR. Tropical forests under the climates of the last 30,000 years. Clim Chang. 1998;39:177–97.
Hewitt G. The genetic legacy of the quaternary ice ages. Nature. 2000;405:907–13.
Foster P. The potential negative effects of global climate change on tropical montane cloud forests. Earth-Sci Rev. 2001;55:73–106.
Jansson R, Dynesius M. The fate of clades in a world of recurrent climatic change: Milankovitch oscillations and evolution. Annu Rev Ecol Syst. 2002;33:741–77.
Hewitt GM. Genetic consequences of climatic oscillations in the quaternary. Phil Trans R Soc Lond B. 2004;359:183–95.
Ramírez-Barahona S, Eguiarte LE. The role of glacial cycles in promoting genetic diversity in the Neotropics: the case of cloud forests during the last glacial maximum. Ecol Evol. 2013;3(3):725–38.
Ornelas JF, Ortíz-Rodríguez AE, Ruíz-Sánchez E, Sosa V, Pérez-Farrera MA. Ups and downs: genetic differentiation among populations of the Podocarpus (Podocarpaceae) species in Mesoamerica. Mol Phylogenetics Evol. 2019;138:17–30.
Hamilton LS, Juvik JO, Scatena FN. The Puerto Rico tropical cloud Forest symposium: introduction and workshop synthesis. In: Hamilton LS, Juvik JO, Scatena FN, editors. Tropical montane cloud forests. New York: Springer-Verlag; 1995. p. 1–18.
Luna-Vega I, Alcántara AO, Espinosa OD, Morrone JJ. Historical relationships of the Mexican cloud forests: a preliminary vicariance model applying parsimony analysis of Endemicity to vascular plant taxa. J Biogeogr. 1999;26:1299–305.
Aldrich M, Bubb P, Hotstettler S, van de Wiel H. Tropical montane cloud forests. Time for action. Cambridge: WWF International/IUCN. The World Conservation Union; 2000.
Marshall CJ, Liebherr JK. Cladistic biogeography of the Mexican transition zone. J Biogeogr. 2000;27:203–16.
Lindner M, Maroschek M, Netherer S, Kremer A, Barbati A, Garcia-Gonzalo J, Seidl R, Delzon S, Corona P, Kolström M, Lexer MJ, Marchetti M. Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. For Ecol Manag. 2010;259:698–709.
Avise JC. Phylogeography: retrospect and prospect. J Biogeogr. 2009;36:3–15.
Morris AB, Ickert-Bond SM, Brunson DB, Soltis DE, Soltis PS. Phylogeographical structure and temporal complexity in American sweetgum (Liquidambar styraciflua; Altingiaceae). Mol Ecol. 2008;17:3889–900.
Ornelas JF, Sosa V, Soltis DE, Daza JM, González C, Soltis PS, Gutiérrez-Rodríguez C, de los Monteros AE, Castoe TA, Bell C, Ruíz-Sánchez E. Comparative Phylogeographic Analyses Ilustrate the Complex Evolutionary History of Threatened Cloud Forests of Western Mesoamerica. PLoS One. 2013;8(2):e56283.
Barrier E, Velasquillo L, Chavez M, Gaulon R. Neotectonic evolution of the isthmus of Tehuantepec (southeastern Mexico). Tectonophysics. 1998;287:77–96.
Gutiérrez-Rodríguez C, Ornelas JF, Rodríguez-Gómez F. Chloroplast DNA phylogeography of a distylous shrub (Palicourea padifolia, Rubiaceae) reveals past fragmentation and demographic expansion in Mexican cloud forests. Mol Phylogenetics Evol. 2011;61:603–15.
Ornelas JF, Rodríguez-Gómez F. Influence of Pleistocene glacial/interglacial cycles on the genetic structure of the mistletoe Cactus Rhipsalis baccifera (Cactaceae) in Mesoamerica. J Hered. 2015;106(2):196–210.
Ornelas JF, Gándara E, Vásquez-Aguilar AA, Ramírez-Barahona S, Ortíz-Rodríguez AE, González C, Mejía SMT, Ruíz-Sánchez E. A mistletoe tale: postglacial invasion of Psittacanthus schiedeanus (Loranthaceae) to Mesoamerican cloud forests revealed by molecular data and species distribution modeling. BMC Evol Biol. 2016;16(78):1–20.
Loveless MD, Hamrick JL. Ecological determinants of genetic structure in plant populations. Ann Rev Ecol Syst. 1984;15:65–95.
Raven PH, Axelrod DI. Angiosperm biogeography and past continental movements. Ann Missouri Bot Gard. 1974;61:539–673.
González-Medrano F. Algunos aspectos de la evolución de la vegetación de México. Bol Soc Bot México. 1996;58:129–36.
Richardson JE, Whitlock BA, Meerow AW, Madriñán S. The age of chocolate: a diversification history of Theobroma and Malvaceae. Front Ecol Evol. 2015;3(120):1–14.
Past Interglacials Working Group of PAGES. Interglacials of the last 800,000 years. Rev. Geophys. 2016;54:162–219.
Pavek DS. Fremontodendron californicum. In: Fire Effects Information System. U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station, Fire Sciences Laboratory. 1993. https://www.fs.fed.us/database/feis/plants/shrub/frecal/all.html Accessed 30 Oct 2019.
Metcalfe SE, O’Hara SL, Caballero M, Davies SJ. Records of Late Pleistocene-Holocene climatic change in Mexico – a review. Quat Sci. 2000;19:699–721.
Toledo VM. Pleistocene changes of vegetation in tropical Mexico. In: Prance GT, editor. Biological diversification in the tropics. New York: Columbia University Press; 1982. p. 93–111.
Chimera C. TAXON; Chiranthodendron pentadactylon Larreat. Plantpono. 2016; https://plantpono.org/wp-content/uploads/Chiranthodendron-pentadactylon.pdf. Accessed 17 May 2018.
Toledo VM. Chiranthodendron pentadactylon Larreategui (Sterculiaceae): una especie polinizada por aves percheras. Bol Soc Bot México. 1975;35:59–67.
Fleming TH, Geiselman C, Kress WJ. The evolution of bat pollination: a phylogenetic perspective. Ann Bot. 2009;104:1017–43.
Doyle JJ, Doyle JL. Isolation of plant DNA from fresh tissue. Focus. 1990;12:13–5.
Etter PD, Bassham S, Hohenlohe PA, Johnson EA, Cresko WA. SNP discovery and genotyping for evolutionary genetics using RAD sequencing. Methods Mol Biol. 2011;772:157–78.
Medina N. 2017. RADSeq library preparation protocol. Modified from: Etter PD, Bassham S, Hohenlohe PA, Johnson EA, Cresko WA. SNP discovery and genotyping for evolutionary genetics using RAD sequencing. Methods Mol Biol. 2011;772:157–78.
Eaton DAR. PyRAD: assembly of de novo RADSeq loci for phylogenetic analyses. Bioinformatics. 2014;30(13):1844–9.
Yu X, Yang D, Guo C, Gao L. Plant phylogenomics based on genome-partitioning strategies: Progress and prospects. Plant Diver. 2018;40:158–64.
Wolfe KH, Li W-H, Sharp PM. Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear DNAs. Proc Natl Acad Sci. 1987;84:9054–8.
Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34(12):3299–302.
Pritchard JK, Wen X, Falush D. Documentation for structure software: Version 2.3. 2010. http://pritch.bsd.uchicago.edu/structure.html Accessed 2 Feb 2020.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.
Earl DA. vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genet Resour. 2012;4:359–61.
Ramasamy RK, Ramasamy S, Bindroo BB, Naik VG. STRUCTURE PLOT: a program for drawing elegant SRUCTURE bar plots in user friendly interface. Springerplus. 2014;3(431):1–3.
Excoffier L, Lischer H. Arlequin ver 3.5. An integrated software package for population genetics data analysis: Swiss Institute of Bioinformatics; 2015. http://www.cmpg.unibe.ch/software/arlequin3/. Accessed 27 Feb 2020.
Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23(2):254–67.
Bryant D, Moulton V. Neighbor-net: an agglomerative method for the construction of phylogenetic networks. Mol Biol Evol. 2004;21(2):255–65.
Ogilvie HA, Bouckaert RR, Drummond AJ. StarBEAST2 brings faster species tree inference and accurate estimates of substitution rates. Mol Biol Evol. 2017;34(8):2101–14.
Heled J, Drummond AJ. Bayesian inference of species trees from multilocus data. Mol Biol Evol. 2010;27(3):570–80.
Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.
Cornuet JM, Pudlo P, Veyssier J, Dehne-García A, Estoup A. DIYABC version 2.0. A user-friendly software for inferring population history through Approximate Bayesian Computations using microsatellite, DNA sequence and SNP data. Institut National de la Recherche Agronomique; 2014. http://www1.montpellier.inra.fr/CBGP/diyabc/. Accessed 19 Feb 2020.
Phillips SJ, Dubik M, Schapire RE. Maxent software for modeling species niches and distributions (Version 3.4.1). 2006. http://biodiversityinformatics.amnh.org/open_source/maxent/ Accessed 23 Mar 2017.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.
Varela S, Lima-Ribeiro MS, Terribile LC. A short guide to the climatic variables of the last glacial maximum for biogeographers. PLoS One. 2015;10(6):e0129037.
We give special thanks to Mario Esteban Véliz Pérez and Luis Velázquez (Universidad de San Carlos de Guatemala), Ernesto López, Prócoro Almazán Rodríguez, Andrés Ernesto Ortíz Rodríguez (Universidad Nacional Autónoma de México) and Gabriel González Adame (Universidad de la Sierra Juárez) for field assistance; Matthew Collin and Clay Clark (University of California, Riverside Genomics Core), Loraine Washburn (Rancho Santa Ana Botanic Garden, Claremont, California) and Rebeca Hernández Gutiérrez (Universidad Nacional Autónoma de México) for laboratory assistance; and Pablo César Hernández Romero (Universidad Autónoma de Aguascalientes) for his assistance in ecological niche modeling.
Funds were provided by CONACyT grant (358420).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Taxon sampling, geographic location and voucher information. The two specimens collected as vouchers are deposited in the herbarium of the Universidad Autónoma de Aguascalientes (HUAA).
RADSeq de novo assembly stats of raw sequences.
Sequence matrix of 36 polymorphic nuclear loci.
Nuclear loci and SNPs data sets substitution models.
Evanno method results from STRUCTURE HARVESTER.
About this article
Cite this article
Hernández-Langford, D.G., Siqueiros-Delgado, M.E. & Ruíz-Sánchez, E. Nuclear phylogeography of the temperate tree species Chiranthodendron pentadactylon (Malvaceae): Quaternary relicts in Mesoamerican cloud forests. BMC Evol Biol 20, 44 (2020). https://doi.org/10.1186/s12862-020-01605-8