Quaternary climate change drives allo-peripatric speciation and refugial divergence in the Dysosma versipellis-pleiantha complex from different forest types in China

Subtropical China harbours the world’s most diverse temperate flora, but little is known about the roles of geographical and eco-climatic factors underlying the region’s exceptionally high levels of species diversity and endemism. Here we address this key question by investigating the spatio-temporal and ecological processes of divergence within the Dysosma versipellis-pleiantha species complex, endemic to subtropical China. Our cpDNA phylogeny showed that this monophyletic group of understory herbs is derived from a Late Pliocene ancestor of the Qinghai-Tibetan Plateau (QTP)/Southwest China. Genetic and ENM data in conjunction with niche differentiation analyses support that the early divergence of D. versipellis and D. pleiantha proceeded through allo-peripatric speciation, possibly triggered by Early Pleistocene climate change, while subsequent climate-induced cycles of range contractions/expansions enhanced the eco-geographical isolation of both taxa. Furthermore, modelling of population-genetic data indicated that major lineage divergences within D. versipellis likely resulted from long-term allopatric population isolation in multiple localized refugia over the last glacial/interglacial periods, and which in turn fostered endemic species formation (D. difformis, D. majoensis) from within D. versipellis in Southwest China. These findings point to an overriding role of Quaternary climate change in triggering essentially allopatric (incipient) speciation in this group of forest-restricted plant species in subtropical China.

. Although the 'western' and 'eastern' lineages appeared to be sisters, this relationship received very low support (PP = 0.94, ML = 56%). Individual mean age estimates for relevant nodes and their 95% highest posterior densities (HPDs) are provided for the different coalescent-type models in Supplementary Table S3 Table S3).
In the rooted tcs network (Fig. 1c), all haplotypes of D. pleiantha were recovered as sister clade to the remainder, separated by five mutational steps. Within this latter clade (hereafter named 'D. versipellis s. lat.'), haplotypes of the 'western' vs. 'central-eastern' lineages formed interior (ancestral) vs. tip (derived) subclades, again, separated by five steps. Notably, the low frequency haplotypes of both D. difformis (H52-H54) and D. majoensis (H55) appeared to be mutational derivatives of the most common and widespread haplotype of D. versipellis (H1) as part of the 'central-east' subclade.
Demographic history. The mismatch distribution for cpDNA haplotypes of the 'western' lineage of D. versipellis s. lat.
was bimodal, and thus differed strongly from that predicted under the spatial and demographic expansion models ( Supplementary Fig. S1). This difference was also registered by non-significant values of both Fu's F S and Tajima's D (Table 2). By contrast, both parameter values were significantly negative for the 'central-east' lineage of D. versipellis s. lat., indicative of population expansion, and the same applied to D. pleiantha in terms of F S (Table 2). Moreover, for either lineage, the MDAs gave unimodal graphs as well as non-significant SSD and H Rag values under both the spatial and demographic models (Supplementary Fig. S1; Table 2). Based on the corresponding τ values and the assumed substitution rate of 1.3 × 10 −9 s/s/y, we calculated a mid-Pleistocene expansion time for the 'central-east '    In the structure analysis of the entire species complex using all 15 loci, ln P(D) progressively increased with K, whereas the ad hoc statistic Δ K showed the highest likelihood at K = 4 ( Supplementary Fig. S2). Based on the latter model, all samples of D. pleiantha formed a distinct cluster (red in Fig. 3a,b), while those of D. versipellis segregated into three regional clusters, with populations originating from (i) both the west and the north (dark blue); (ii) the south (light blue); and (iii) the east (green). All samples of D. difformis and D. majoensis were assigned to the southern cluster of D. versipellis (Fig. 3a,b). Similarly, in the corresponding PCoA, the first axis distinguished individuals of D. pleiantha from all others, whereas individuals of D. difformis and D. majoensis broadly overlapped with the southern cluster of D. versipellis in genetic space (Fig. 3c). After removing the six outlier loci, Δ K favoured a model with K = 3. Running structure and using K = 3, D. pleiantha and eastern populations of D. versipellis were assigned to two disctinct gene pools, whereas the reminder made up one single gene pool (see Supplementary Fig. S3). However, with K = 4, results were consistent with the clustering patterns inferred using all 15 loci ( Supplementary Fig. S3).
The migrate analysis using all 15 loci revealed all pairwise estimates of migration/gene flow (M) among D. pleiantha, D. difformis, D. majoensis and the three regional clusters of D. versipellis to be low, ranging from 0.351 (D. difformis to west-northern cluster of D. versipellis) to 2.160 (eastern to west-northern cluster of D. versipellis) ( Table 3). Basically the same results were obtained after removing the outlier loci (Table 3). Overall, these structure, PCoA and migrate analyses indicate only relatively minor effects of outlier loci on the overall population structure and gene flow of this species complex.

ABC-based inferences of species/cluster divergence.
In the diyabc analysis of the entire species complex based on the nine neutral EST-SSR loci, the highest posterior probabilities derived from the direct estimate approach were found for scenarios 1 (0.31, 95% CI: 0.00-0.71) and 4 (0.29, 0.00-0.69), whereas the values derived from the logistic regression were higher for the former scenario with no overlap of confidence intervals [0.46 (0.45-0.47) vs. 0.27 (0.26-0.28)] ( Supplementary Fig. S4). In consequence, we favoured scenario 1, according to which the three regional structure clusters of D. versipellis s. lat. (west-north, east, south) diverged simultanously from a common ancestor, rather than any alternative hypothesis of 'ordered' relationships among these groups ( Supplementary Fig. S5a, Table 4  Ecological niche modelling and niche identity tests. The maxent models for D. versipellis s. lat. (including D. difformis/D. majoensis) and D. pleiantha had high predictive power and did not overfit the presence data (AUC values = 0.906 ± 0.118 and 0.972 ± 0.041, respectively). The current potential ranges of either group (defined as modelled suitability ≥ 0.75; Fig. 4a,b) were largely consistent with their actual distributions, except for predicted but unsupported occurrences of D. versipellis s. lat. in the far east of China, Taiwan, and South Korea (Fig. 4a). During the LGM (c. 21,000 yr BP), D. versipellis s. lat. experienced a drastic restriction of its current range in Central-East China (with only scattered areas persisting along the middle-Yangtze River), whereas suitable habitat apparently increased around the Sichuan Basin and areas further south (e.g. Southwest China, Yungui Plateau) (Fig. 4c). By contrast, for D. pleiantha, the LGM modelling predicted a more contiguous than currently scattered distribution in East China, extending southward into Taiwan (then connected with Mainland China), plus a higher suitability of areas both west-and eastward, i.e. south of the lower/middle Yangtze and on the then exposed East China Sea basin, respectively (Fig. 4d).
The results of the niche identity tests supported the existence of niche differentiation between the WTD forest-dweller D. versipellis s. lat. and its WTE forest counterpart D. pleiantha. Enmtools showed that the observed similarity values for I and D (0.818 and 0.582, respectively) were both significantly lower (P < 0.01) than the corresponding values expected from the pseudoreplicated datasets in the paired analyses ( Supplementary Fig. S6).

Discussion
Our cpDNA phylogeny (Fig. 2) confirms major relationships within Dysosma as hypothesized previously on the basis of cpDNA/nrDNA data 17 . In both analyses, the D. versipellis-pleiantha complex forms a strongly supported monophyletic group, which, in the present analysis, emerges from a paraphyletic and likewise supported grade of the remainder of the genus, comprised of D. delavayi (as sister to the complex), D. aurantiocaulis, and D. tsayuensis. Even though the species complex and the three latter, high-elevation species are not reciprocally monophyletic, they exhibit different morphologies 16,18 and eco-geographies: while the former group is native to the lowland WTE or montane WTD forests of subtropical China, the latter is endemic to the cool-temperate/subalpine forests of Yunnan and the south-eastern QTP region. Given that the closest relative of Dysosma, Sinopodophyllum hexandrum, is native to the Himalayan regions 16,17,19 , we therefore assume that the most recent common ancestor (MRCA) of the D. versipellis-pleiantha complex originated in the high-altitude habitats of the QTP region/ Southwest China, and then expanded its geographic range across subtropical China. In general, such a niche shift from high-to-lower elevation habitats is thought to be rare 20 , but is also known from other taxa in our study area (e.g.    By contrast, what emerges from our rooted cpDNA haplotype network (Fig. 1c) is striking evidence that D. pleiantha is sister to the remainder of the species complex (D. versipellis s. lat.). This inference of eco-geographically driven, and possibly adaptive divergence is further strengthened by our niche modeling, which suggests that these taxa occupy significantly different climatic environments ( Supplementary Fig. S6). Moreover, in the cpDNA network, the central-eastern clade of D. versipellis s. lat. is placed in a derived (tip) position relative to the western (interior) one, and thus most distantly related to D. pleiantha (Fig. 1c). Hence, the present-day abutting ranges of D. versipellis and D. pleiantha in East China are secondary, rather than reflecting primary divergence between adjacent populations of a previously homogeneous species, i.e. a parapatric mode of speciation. We therefore conclude that the main mechanism, which caused the Early Pleistocene ]. Yet, apparently, these expansions induced no secondary contact, as both taxa are identified here as genetically distinct units in terms of cpDNA haplotype (Fig. 1a,c) and EST-SSR variation (Fig. 3). Moreover, the migrate analyses of EST-SSRs uncovered extremely low levels of post-divergence gene flow between D. versipellis and D. pleiantha ( Table 3), suggesting that their propensity for hybridization has been overstated in the taxonomic literature 16 . Rather, we suspect that both taxa, and their associated forest biomes, largely remained isolated from each other throughout the glacial and inter-/ postglacial periods of the Quaternary. In fact, our ENM results indicate that, during the LGM (and possibly earlier cold periods), D. versipellis s. lat. and D. pleiantha occupied distinct refugia in West/Southwest China and East China, respectively, with only little potential range overlap in Central-East China (Fig. 4c,d). That glacial periods caused the range of D. versipellis s. lat. to retract south-westward (i.e. around the Sichuan Basin, Yungui Plateau), and/or to become fragmented, is similar with findings in other WTD forest species from subtropical China (e.g. Cercidiphyllum japonicum 6 ; Kalopanax septemlobus 5   with fossil-pollen data 3,32 , suggesting that WTE forest largely retracted to the tropical South (≤ 24°N). Instead, our ENM data suggest that this drought-sensitive herb 16 benefitted from relative climatic stability in East China over the last glacial cycle(s), possibly due to warm and moist air masses of monsoons from the ocean 33,34 and/or topographic heterogeneity (e.g. Huang/Tianmu Mts.) that buffered against regional climate change 35 . The almost complete fixation of cpDNA haplotype H48 in D. pleiantha (Fig. 1a), together with its homogenous EST-SSR gene pool (red in Fig. 3), precludes more detailed inferences about the glacial refugia and (re-)colonization routes of this species. Nonetheless, when combined with the MDA evidence, this near lack of genetic structure most likely results from a founder effect in the wake of a leading edge expansion during the Late Pleistocene/Holocene 6 . Assuming this expansion signals recovery from a bottleneck 36,37 , any such earlier reduction in population size in D. pleiantha could have further enhanced its geographic isolation from D. versipellis s. lat. Moreover, any climate amelioration should have likewise promoted the isolation of these taxa. This hypothesis is based on the expectation that inter-/postglacial periods generally fostered the displacement of WTD forest into higher elevations by the expansion of WTE forest at lower elevations, as supported by phylogeographic studies 31,38 and palaeo-data 32 .
In sum, these findings indicate that the divergence between D. versipellis s. lat. and D. pleiantha has been triggered and maintained by geographical and ecological-adaptive processes in response to Pleistocene climate change. Our survey of cpDNA sequence variation at 31 locations of D. versipellis s. lat. revealed substantial population differentiation (Φ ST = 0.84) and significant phylogeographic structure (N ST = 0.818 > G ST = 0.611; P < 0.05). These patterns largely result from the major phylogeographic break across the Sichuan Basin between ancestral 'western' and derived ('founded') central-eastern populations, as well as the marked differences within the latter (Fig. 1). These results are similar to those of two earlier studies in D. versipellis based on cpDNA 31 and AFLPs 39 , which also registered a major west vs. central-east split in this species, despite more limited samples (10 and nine populations, respectively). By contrast, based on our EST-SSR data, structure identified three geographical clusters in D. versipellis s. lat. (west-north, south, east). Hence, all populations of the western cpDNA lineage (EM, DJ) were found to share the same gene pool with five northern populations (ZF, SN, MHP, TP, TT) of the central-eastern lineage (Fig. 3a,b). This failure of the EST-SSRs to register the phylogeographic cpDNA split is probably best explained by extended times for lineage sorting to occur at nuclear compared to cytoplasmic loci 40 rather than contemporary pollen-mediated introgression from the north into the west. The latter appears less likely if one considers not only the limited pollen-dispersal capacity of Dysosma 16,39 but also the migrate results (Table 3)  cpDNA analysis of D. versipellis 31 dated the separation between the western and central-eastern lineages to the mid-Pleistocene, about 0.48 Ma (95% CI: 0.18-0.86 Ma). In the present study, ABC analyses revealed that the divergence of the three regional EST-SSR clusters of D. versipellis s. lat. likewise dates back to the mid-Pleistocene (c. 0.59 Ma, 95% CI: 0.15-1.29 Ma; t 2 in Supplementary Fig. S5a, Table 4). Even though difficult to compare across studies, and despite their broad confidence intervals, both time estimates strongly suggest that major lineage divergences within D. versipellis s. lat. long pre-date the LGM. Moreover, the present cpDNA and EST-SSR data provide a congruent picture of the phylogeographic history of D. versipellis s. lat., as both support the hypothesis that populations of this taxon survived the (Late) Pleistocene glacial maxima in multiple refugia in both West/ Southwest China and Central-East China. In the former region, there were probably three refugia: one, as identified by cpDNA, at the eastern slopes of the QTP (comprising populations DJ, EM; Fig. 1a); and two others, as revealed by the EST-SSRs, extending south vs. north-east of the Sichuan Basin (light vs. dark blue clusters; Fig. 3b,c). In Central-East China, there was probably a major refuge located south of the mid-/lower Yangtze, as demonstrated by the nuclear cohesion of all eastern populations (green cluster; Fig. 3b,c), and possibly a more localized one further north in the Dabie Mts., as might be inferred from the high number of unique cpDNA haplotypes in the 'north-cluster' population TT (Fig. 1a, Supplementary Table S2). Notably, these latter two refugia are also predicted by the ENM of D. versipellis s. lat. for the LGM (Fig. 4c); on the other hand, the genetic data better serve to illustrate that West/Southwest China did not function as a single refugium for this WTD forest species during glacial periods but instead hosted different nuclear and/or cytoplasmic lineages in separate refugia. The occurrence of these 'refugia-within-refugia' 41 not only caused the allopatric divergence of these lineages but probably also fostered the recent speciation of D. difformis and D. majoensis from within D. versipellis.
Dysosma difformis and D. majoensis have proximate but non-overlapping eastern vs. western distributions along the southern fringe of the Sichuan Basin (Fig. 1a,b). These two narrow endemics share the same WTD forest habitat with D. versipellis (s. str.) and are sometimes difficult to distinguish from the latter based on morphology, as most of their diagnosable features relate to aspects of leaf lobing 16,18 . Our genetic data clearly suggest a close relationship among these three taxa, without traces of other Dysosma species (e.g. D. delavayi) as suggested previously for D. majoensis 16 . Rather, all cpDNA haplotypes unique to D. difformis (H52-H54; only population SH) and the single one fixed in D. majoensis (H55) are independent mutational derivatives of the most common and widespread haplotype of D. versipellis (H1; see Fig. 1a,c). Moreover, the EST-SSR data show that D. difformis and D. majoensis, despite sharing the same gene pool (light blue) with all southern populations of D. versipellis (Fig. 3), have likely exchanged only few migrants with each other and the latter species in the distant past (Table 3). Finally, the corresponding ABC analysis favoured a scenario of unresolved relationships, according to which D. difformis, D. majoensis and D. versipellis (southern cluster) evolved independently from a common ancestor, possibly during the mid-Pleistocene [c. 0.43 (0.07-1.47) Ma; t 4 in Supplementary Fig. S5b]. Taken together, these data seem to support a scenario in which climate deterioration following China's penultimate interglacial period (c. 0.33-0.46 Ma; see above) promoted the independent origin of D. difformis and D. majoensis from within D. versipellis through population isolation in two separate refugia along the southern Sichuan Basin. If so, the presence of cpDNA haplotypes H1 and H55 in D. difformis (Fig. 1a,c) most likely results from, respectively, incomplete lineage sorting (ILS; i.e. recent descent from D. versipellis) and past chloroplast introgression from D. majoensis, perhaps facilitated by massive expansion of suitable habitat in the area during the LGM (Fig. 4c). Clearly, the present data cannot fully disentangle these hypotheses of recent speciation, ILS and/or past admixture; however, they do suggest that the combined effects of these historical processes, rather than contemporary gene flow, are responsible for the limited evidence of genetic (cytoplasmic/nuclear) divergence of D. difformis and D. majoensis. These data, therefore, should not be interpreted to reject the hypothesis of geographic-reproductive isolation of these two endemics and their evolutionary distinctiveness as taxonomic units of interest 42,43 ; rather, they underscore the mountain regions of Southwest China of central importance for (micro-)refugial allopatric speciation 31,44 .
In summary, the data presented here for the D. versipellis-pleiantha complex suggest an overriding role of Quaternary climate change in triggering essentially allopatric (incipient) speciation in this group of forest-restricted plant species in subtropical China.
Fifteen EST-SSR primer pairs, developed from the transcriptome of D. versipellis (termed 'EDV' primers), were employed for genotyping all 577 DNA samples. Ten of those primers 37,40,46,52,53,54,59,60,67) were described in Guo et al. 51 (GenBank accession numbers: KJ000290-KJ000299). Respective primers for the other five loci (EDV-81, 82, 102, 118, 119) were developed in this study (Supplementary Table S6). All loci were proved to be consistently variable and transferable among the four species of the D. versipellis-pleiantha complex after a pilot study. Amplification of the EST-SSR loci followed the protocol of Guo et al. 51 . Fluorescently labelled PCR products were supplemented with the internal size standard GS-500 and separated on a MegaBACE 1000 (GE Healthcare Biosciences, Sunnyvale, CA, USA). Alleles were scored manually in genemarker v2.20 (SoftGenetics, State College, PA, USA).
Chloroplast hapotype diversity and population structure. The number of haplotypes, haplotype diversity (h), and nucleotide diversity (π) were calculated at the levels of populations (h S , π S ) and species (h T , π T ) using dnasp v5.0 52 . Non-hierarchical and hierarchical analyses of molecular variance (AMOVAs) were carried out in arlequin v3.5 53 with significance of Φ-statistics tested by 1000 non-parametric random permutations. The presence of phylogeographic structure was tested across all populations of the species complex and for each species and/or lineage separately using permut v1.0 54 . All haplotype sequences identified in the present study were deposited in GenBank with accession numbers (see Supplementary Table S7 for details).
Phylogenetic haplotype relationships and molecular dating. Phylogenetic cpDNA haplotype trees of the species complex were constructed using maximum likelihood (ML) and Bayesian inference (BI) methods, with the remaining Dysosma species treated as part of the ingroup and S. hexandrum and P. peltatum as outgroup (see Supplementary Methods for details). For the phylogeographic cpDNA analyses, a haplotype network was constructed in tcs v1.21 55 under the 95% statistical parsimony criterion with gaps (indels) coded as substitutions (A or T), and rooted with D. delavayi.
Divergence date analyses were conducted on the cpDNA data set in beast v1.80 56,57 using the GTR + G substitution model (see Supplementary Methods for details) under the assumption of an uncorrelated lognormal relaxed clock 56 . As there are no known fossil records in Podophylloideae, we employed a substitution rate of 1.3 × 10 −9 substitutions per site per year (s/s/y) as estimated for non-coding chloroplast regions of D. versipellis 31 . We implemented two coalescent-type models, assuming either a constant population size or population expansion. For each model, four independent Markov chain Monte Carlo (MCMC) runs of 50 million generations each were performed with sampling every 5000 generations, following a burn-in of the initial 10% cycles. We used TreeAnnotator v1.6.1 57 to construct a maximum clade credibility tree.

Demographic analyses.
To infer the demographic history of major cpDNA lineages within the species complex, we tested the null hypotheses of a spatial expansion and a pure demographic expansion using mismatch distribution analysis (MDA) in arlequin. For two expanding lineages identified, the MDA-derived expansion parameter (τ ) and its 95% confidence interval (CI) were converted to absolute estimates of time since expansion (see Supplementary Methods for more details). In addition, we used tests of selective neutrality [Tajima's D 58 ; Fu's F s 59 ] to infer potential population growth and expansion.

Population structure and migration/gene flow inferred from EST-SSR markers. For the EST-SSR
dataset, all 15 genotyped loci were checked for frequencies of null alleles, deviations from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) (see Supplementary Methods for details). For each locus, we also tested potential signatures of selection in lositan 60 , arlequin, and bayescan v2.1 61 . Following De Mita et al. 62 , we considered a locus to be under selection if at least two outlier tests were significant for that particular locus. Using fstat, the following diversity and inbreeding parameters were calculated for each population and across all 15 loci: N A , numbers of alleles; R S , allelic richness 63 ; H E ; and the average inbreeding coefficient (F IS ). Hierarchical and non-hierarchical AMOVAs were carried out in arlequin using R-statistics and significance tests as described above for cpDNA.
Genetic subgroups of the species complex were identified by Bayesian analyses in structure v2.3.4 64 , using the admixture model and assuming independent allele frequencies among populations (see Supplementary Methods for details). Clustering analyses for the EST-SSRs were complemented with an individual-based principal coordinates analysis (PCoA) using software genalex v6.5 65 .
We also used coalescent-based methods implemented in migrate-n v3.6 66 to obtain pairwise estimates of mutation-scaled migration rate (M) over c. 4N e generations in the past 67 between species and/or intraspecific regional clusters (only D. versipellis; see structure results). All analyses were carried out with and without outlier loci included to assess their effects on the inference of population structure and gene flow 68 . More software running details can be found in the Supplementary Materials.
Testing for the times and orders of species/cluster divergence using ABC. We used ABC simulations in diyabc v2.0 69,70 to gain further insights into the times and orders of divergence among clusters identified by structure based on the nine neutral EST-SSR loci. In a first analysis ('a'), we tested the simultanous divergence of the three clusters of 'D. versipellis s. lat.' (i.e. west-north, east, south, with the latter including D. difformis/D. majoensis; see Results) from a common ancestor against three alternative models, reflecting all possible relationships among these clusters (see scenarios 1-4 in Supplementary Fig. S5a and Table S5). For each scenario, D. pleiantha was assumed as sister, based on the cpDNA haplotype network analysis (see Results). In Scientific RepoRts | 7:40261 | DOI: 10.1038/srep40261 a second analysis ('b'), we tested four similar divergence order hypotheses between the taxonomic units of the southern cluster, i.e. D. versipellis (south), D. difformis and D. majoensis (scenarios 1-4 in Supplementary Fig. S5b and Table S5). Additional details in data simulation, scenario comparison, parameter estimation, etc. can be found in the Supplementary Materials. Ecological niche modelling and niche identity tests. Ecological niche modelling (ENM) was carried out in maxent v3.3.1 71,72 to predict suitable past and present climate envelopes for, respectively, the WTD forest species (i.e. D. versipellis s. lat., including D. versipellis, D. difformis, D. majoensis) and the WTE forest species D. pleiantha. Information on the geographic distribution of the three former species was based on the 31 populations included in this study, combined with 102 presence records obtained from the Chinese Virtual Herbarium (CVH) (www.cvh.org.cn) and the National Specimen Information Infrastructure of China (NSII) (www.nsii.org. cn). In addition, 44 records were available for D. pleiantha, i.e. 11 localities of this study plus 33 obtained from CVH and NSII.
Current distributions of the two sets of species were modelled using six bioclimatic data layers (annual mean temperature, annual precipitation, precipitation of wettest, driest, warmest and coldest quarter) available from the WorldClim database (http://www.worldclim.org) 73 at 2.5 arc-min resolution for the present . This restricted bioclimatic dataset avoided including highly correlated variables, and thus prevented potential over-fitting 74 . This model was then projected onto the palaeoclimate dataset simulated by the community climate system model v3.0 75 to infer the extent of suitable habitat during the Last Glacial Maximum (LGM; c. 21 kya BP). Model performance was evaluated using the area under the 'Receiver Operating Characteristic (ROC) Curve' (AUC) calculated by maxent. Values between 0.7 and 0.9 indicate good discrimination 76 .
We also performed niche identity tests in enmtools v1.4.3 77 based on all the 19 BIOCLIM variables from the WorldClim dataset for testing the null hypothesis that the WTD vs. WTE forest species (as defined above) are occupying identical climatic environments ('niches'). However, we did not further evaluate the niche similarity between D. difformis and D. majoensis because of the very low number of known occurrence records for either species. Niche overlap was quantified using the standardized Hellinger distance (I) and Schoener's D 78 .