Idiosyncratic responses of evergreen broad-leaved forest constituents in China to the late Quaternary climate changes

Subtropical evergreen broad-leaved forest (EBLF) is one of the most important vegetation types in China. Inferences from palaeo-biome reconstruction (PBR) and phylogeography regarding range shift history of EBLF during the late Quaternary are controversial and should be reconciled. We compared phylogeographic patterns of three EBLF constituents in China, Castanopsis tibetana, Machilus thunbergii and Schima superba. Contrary to a chorus of previous phylogeographic studies and the results of species distribution modelling (SDM) of this study (in situ survival during the LGM), the three species displayed three different phylogeographic patterns that conform to either an in situ survival model or an expansion-contraction model. These results are partially congruent with the inference of PBR that EBLF was absent to the north of 24° N at the LGM. This study suggests that the constituents of EBLF could have responded idiosyncratically to climate changes during the Late Quaternary. The community assemblages of EBLF could have been changing over time, resulting in no palaeo-analogs to modern-day EBLF, which may be the main reason responsible for the failure of PBR to detect the occurrence of EBLF north of 24° N at the LGM.

is preffixed by a letter C in (b), M in (c) and S in (d) that particularly denotes the sampled populations of Castanopsis tibetana, Machilus thunbergii and Schima superba, respectively. One population (M69) remote from others in Machilus thunbergii is not shown in (c). The red dashed lines show the approximate northern borders of EBLF during the LGM based on PBR (modified after Harrison et al., 2001). Maps were generated using ArcGIS version 9.3 (http://www.esri.com/software/arcgis/arcgis-for-desktop) and Adobe Illustrator CS3 13.0 and modified using Adobe Photoshop CS 8.0.
whether changes in distribution and population size represent the influence of extrinsic factors that affected whole communities, or whether they can be ascribed to stochastic variation and other chance factors particular to each organism 23 . Some comparative studies have succeeded in identifying patterns common across many groups of organisms [24][25][26] . However, idiosyncratic phylogeographic patterns among co-distributed species are commonly reported [27][28][29] , suggesting that community composition is highly dynamic 29 . If the latter model applies to EBLF in subtropical China, it is conceivable that the absence of EBLF in subtropical China at the LGM might have been caused by alterations in community composition and the changes in relative abundance of EBLF components, rather than by the complete retreat to the south, thus settling the dispute between PBR and phylogeography. However, such a hypothesis has never been proposed and explicitly tested by phylogeographic researchers in China, possibly due to the difficulty of gathering large data sets from multiple co-distributed species.
To reconcile the two scenarios suggested by PBR and phylogeography, this study used the CP approach to test for congruent versus idiosyncratic phylogeographic patterns of multiple species from EBLF in subtropical China. The EBLF zone within China is divided into two subregions (Fig. 1a): the moist EBLF in eastern China influenced by the Pacific monsoon and the semi-moist EBLF in the west determined by the Indian monsoon 2 . Because the climate in the two subregions is quite different (e.g., the annual rainfall is 1000-2000 mm in the moist east and 900-1200 mm in the semi-moist west 1 ), very few dominant tree species are shared across the two subregions. Therefore, we only focused on three tree species in the eastern subregion of EBLF, i.e., Castanopsis tibetana Hance  Table S2; Fig. 1). These species satisfy three criteria: 1) they are all dominant broad-leaved tree species in EBLF and thus play a fundamental role in the ecological functions of EBLF; 2) they belong to different families and genera that avoid the influence of conserved niche within lineages (and thus similar responses to climate change); and 3) they are components of the same forest communities. We surveyed chloroplast DNA (cpDNA) sequence variation in these three tree species and determined their phylogeographic and historical demographic patterns. We also used SDM to hindcast the geographical distribution of the three tree species at the time of the LGM. By providing quantitative and high-resolution predictions of individual palaeodistributions, this method can complement fossil and genetic evidence concerning the past distribution of species and forest communities 30 . By using molecular and SDM approaches, we address the following issues: i) whether the focal species conform to the ISS or the EI model during the late Quaternary; ii) whether they display congruent phylogeographic patterns and demographic histories across a shared landscape; and iii) whether CP is able to reconcile the dispute between PBR and previous phylogeographic studies.
Estimates of haplotype diversity (h) and nucleotide diversity (π ) for each population of each of the three species are summarized in Supplementary Table S2. The overall haplotype diversity was equivalent in C. tibetana and S. superba (h T = 0.797 and 0.765, respectively), but higher than that in M. thunbergii (h T = 0.644). Machilus thunbergii displays the lowest nucleotide diversity, while C. tibetana has the highest. The average population genetic diversity values are all considerably low in C. tibetana, M. thunbergii, and S. superba (h S = 0.151, 0.165, and 0.136, respectively). Phylogeographic structure. In C. tibetana, samova failed to uncover any reliable population genetic grouping, and there was no apparent structure in the haplotype network (Fig. 1b). Except for two relatively widespread haplotypes, CH1 and CH6, which were found in 22 (49%) and 11 (24%) populations, respectively, most haplotypes were confined to a specific region. Fifteen polymorphic populations and 13 haplotypes (CH2-5, CH8-10, CH13-16 and CH18-19) private to a single population were scattered across the species' range.
In the samova analysis of S. superba, the value of F CT reached a plateau at 0.894 when the number of groups was two (K = 2). As shown in Fig. 1d, the two groups largely correspond to western and eastern parts of the distribution range. The western group was dominated by SH8. The eastern group had two subgroups, the southern one nearly fixed for SH2 and the northern one for SH5. The haplotype network recovered all western haplotypes (SH6-9) as a clade, which was separated by five mutational steps from an eastern clade (containing SH1-5). Several populations in the contact area of the two groups showed a mixture of two clades (S21, S28, S33, and S57).
In the multiple regression analyses of genetic variation, haplotype diversity was correlated negatively with latitude in M. thunbergii (P < 0.001), but not in C. tibetana or S. superba (Fig. 2). Nonetheless, a U-test showed that N ST was significantly larger than G ST Table 1). Despite a slight recent decrease in population size detected by our BSP analysis ( Supplementary Fig. S1a), the data suggest [we prefer] that C. tibetana likely experienced past expansion. Based on the corresponding τ value, the demographic expansion event was estimated to have occurred at 9,600 (95% highest posterior density, HPD: 0-53,000) years bp, although this estimate should be treated with caution given the wide confidence interval ( Table 1).
For M. thunbergii, the negative Tajima Table 1), the mismatch distribution for S. superba was bimodal, which may reflect the regional population structure of the species due to the west-east divergence (Fig. 3c). In support of this explanation, when analyzed separately, the western populations showed a unimodal distribution (Fig. 3b) that was statistically consistent with the expansion model (P values > 0.05 for SSD and H Rag ; Table 1). This expansion event was dated at 6,100 (95% HPD: 0-115,100). When examined separately, the eastern populations showed a bimodal distribution (Fig. 3c), because this group was further structured into two subgroups (Fig. 3c). A BSP analysis detected a slight population increase in both the western and eastern groups ( Supplementary Fig. S1c, S1d). Taken together, the data suggest that S. superba may have experienced regional post-glacial expansion from localized refugia.

Species distribution modeling. For each species at present and LGM, the areas under the 'Receiver
Operating Characteristic (ROC) Curve' (AUC) all have values > 0.95, indicating good predictive model performance. In general, the predicted distributions of the three species under present conditions are similar to their actual distributions in mainland China, each with a potentially patchy range in mountainous areas of subtropical China and a slightly more continuous range for C. tibetana. At the LGM, the potential ranges of the three species contracted a little to the south in the eastern part of their ranges. However, most of the current occurrences are assigned high suitability at the LGM, suggesting overall range stability (Fig. 4a). Therefore, the three species did not retreat entirely to the tropical south during the LGM. SDM analyses found that the niches for the three species were very similar (Fig. 4b). ENMTOOLS showed that empirically observed values for I and D were comparable to those expected from pseudoreplicated data sets in all paired analyses (C. tibetana vs. M. thunbergii, C. tibetana vs. S. superba and M. thunbergii vs. S. superba) (Fig. 4b).

Idiosyncratic responses of three dominant trees in EBLF to the late Quaternary climate changes.
Recently, phylogeographic patterns of subtropical EBLF constituents in China have received increasing attention (Supplementary Table S1). These studies consistently revealed limited or regional post/inter-glacial expansion from northern localized refugia (ISS model), rather than extensive northward spread from more southerly located refugia (EI model) as inferred from palaeo-biome reconstruction (PBR) based on the fossil pollen record [9][10][11] . However, we found that three dominant tree species in EBLF, C. tibetana, M. thunbergii and S. superba, have different range shift trajectories during the late Quaternary, which cannot be predicted purely from either the ISS or EI model.
Castanopsis tibetana (Fagaceae), displays a pattern of multiple northern refugia (private haplotypes and polymorphic populations occurring in the north, Fig. 1b) and regional postglacial expansion, which is consistent with the ISS model. It exhibits the highest genetic diversity among the three focal species, suggesting that more glacial refugia could have maintained a larger historical population size. All members of Fagaceae in EBLF examined so far show similar phylogeographic patterns (i.e., Castanopsis eyrei 13 , Castanopsis fargesii 14 and Quercus glauca 15 ). The phylogeographic congruence among members of Fagaceae may be explained by shared population histories because all these species are common in EBLF and always constitute the canopy of the forest 31 . However, an alternative explanation for the congruence may be niche conservatism among closely related species 17 . Based on a meta-analysis of phylogenies from many studies in the Southern Hemisphere, Crisp et al. 32 concluded that closely related species are more ecologically similar than would be expected by chance. Therefore, similar ecological requirements due to phylogenetic relatedness among members of Fagaceae in EBLF may account for their phylogeographic congruence, at least in part.
In contrast to C. tibetana, Machilus thunbergii has the lowest genetic diversity of the three focal species and shows a significant negative correlation between genetic diversity and latitude (Fig. 2). The results suggest that it conforms to a typical pattern of pioneer or leading-edge expansion, which is expected by the classical EI model 33 . This is the first study that reports the EI model for constituents of EBLF in subtropical mainland China, although evidence of extensive range expansion has been implicated in two case studies of Taiwanese plants. Wu et al. 34 found that Machilus kusanoi in Taiwan (most of which is covered by EBLF, Fig. 1a) experienced strong post-glacial range expansion. Machilus thunbergii was also investigated in Wu et al. 34 : range expansion is also evident in this species because of low genetic diversity among populations and the presence of a widespread haplotype, similar to the pattern observed in mainland China in our study. For another member of Lauraceae, Cinnamomum kanehirae, Liao et al. 35 also detected significant spatial range expansion in Taiwan. Together, these studies suggest that some constituents of EBLF in subtropical China (e.g., several species of Lauraceae), may have retreated to the south during the LGM and recolonized the present EBLF region during the Holocene. This conclusion is congruent with recent findings of warm-temperate deciduous forest species in subtropical China. While most species of this type of vegetation exhibit a pattern of multiple refugia and limited range expansion (reviewed by Liu et al. 36 ), some species show evidence of extensive range expansion 37,38 .
All populations of Schima superba were divided into two groups by the samova analysis (Fig. 1d). Furthermore, the eastern group can be subdivided into two subgroups. This pattern is an expected consequence of long-term isolation among multiple refugia and subsequent localized range expansion. However, compared to C. tibetana, S. superba shows significant phylogeographic structure because the major groups correspond to two highly differentiated clades (Fig. 1d). The east-west divergence in eastern subtropical China has also been observed in Rhododendron simsii 19 and Tsuga 39 and may reflect the strong landscape effects of the Wuyi Mountains and the Poyang Lake valley as natural barriers to dispersal. However, such barriers may be semi-permeable and not applicable to all species, because similar divergence has not been observed in C. tibetana and M. thunbergii, which occupy the same region.
Taken together, in contrast to phylogeographic congruence of EBLF constituents in subtropical China in previous reports [13][14][15] , this study revealed diverse phylogeographic patterns in three dominant trees, i.e., multiple refugia (across the range) without significant phylogeographic structure in C. tibetana, glacial retreat and post-glacial recolonization in M. thunbergii, and multiple refugia with significant phylogeographic structure in S. superba. This observation agrees with the suggestion of many CP studies that congruent phylogeographic patterns are not necessarily the rule for co-distributed species because of idiosyncrasies in their biological and ecological characteristics, short duration of sharing geographical distribution, and/or different evolutionary histories of genes in the species studied (reviewed by Gutiérrez-García & Vázquez-Domínguez 40 ). The results of this study also support the view that communities are not cohesive units that may be broken up and reformed in different configurations repeatedly and regularly on time scales of a few thousand years 41 .
Reconciling palaeo-biome reconstruction and phylogeography. Although there are contradictory results in a few cases 42,43 , the fossil record and phylogeography are largely synergistic in the inferences of past distribution of organisms 5 . Therefore, the distinct scenarios of the past distribution of EBLF in subtropical China inferred from PBR and phylogeography are quite unusual and need to be reconciled. In this study, we found that the phylogeographic patterns of C. tibetana and S. superba are congruent with the results of previous phylogeographic studies of EBLF constituents. However, M. thunbergii displays a typical pioneer or leading-edge expansion pattern, which conforms to the inference of PBR 9-11 , although there is a small difference in the northern limit at the LGM between M. thunbergii phylogeography (26°N) and PBR (24°N). These results suggest that the inference of previous phylogeographic studies could be biased by stochastic variance among species and phylogenetic niche conservatism within Fagaceae. More phylogeographic surveys, particularly CP studies with the same sampling strategy and molecular markers, may provide an opportunity to reconcile the discrepancy between PBR and phylogeography.
Despite the fact that some species in EBLF such as M. thunbergii retreated to the south during the LGM, the complete absence of EBLF to the north of 24°N at that time suggested by PBR is still contrary to the phylogeographic patterns of most plant species characteristic of this vegetation. Although additional paleopalynological investigations are needed to verify the conclusion of previous PBR, the absence of EBLF to the north of 24°N during the LGM cannot be simply ascribed to the sparse pollen record. In fact, by comparing fossil records of individual forest plants and palaeo-biome maps in North America, Williams et al. 41 also found that the results of PBR sometimes did not match with the original fossil records of individual plants. This observation was attributed to information loss in PBR about the distribution and abundance of individual taxa during the categorization of pollen taxa to one or more plant functional types (PFTs) and biomes 44 and the violation of the assumption of community stability. According to our CP results, EBLF may have survived in the north during the LGM, but with reduced distribution and abundance (as suggested by the strong evidence for regional postglacial expansion in both C. tibetana and S. superba) as well as with a different species composition (as suggested by the evidence that M. thunbergii retreated to the south). This kind of community dynamics makes it likely that PBR would omit some EBLF components in the fossil record or misidentify the "EBLF" of the LGM as other types of vegetation because the LGM "EBLF" has no counterpart in the present vegetation 41 . This conclusion suggests that CP at the community level is capable of providing a full understanding of historical community dynamics 22 , avoiding the stochastic variance inherent in single-species phylogeographic studies and identifying the artifacts specific to PBR. Therefore, CP offers an insightful perspective for testing competing biogeographic hypotheses evoked by different disciplines.
The difference between comparative phylogeography and species distribution modeling. Our predicted distributions of three dominant tree species showed that three species survived in situ during the LGM (Fig. 4a), consistent with the SDM results of a few EBLF constituents [13][14][15] . However, these predictions are only partially congruent with the patterns observed in the CP study conducted here, because M. thunbergii exhibits an expansion-contraction pattern. As with any model, SDMs rely on many underlying assumptions and thus incorporate uncertainties 45 . One of the most important assumptions in SDMs is niche stability over the study period 30 . However, niche stability may be affected by either changes in biotic interactions in the community or evolutionary adaptation to the biotic and abiotic environment (i.e. changes in realized niches and fundamental niches) 4 . While fundamental niches are less likely to evolve over relatively short periods, such as the LGM to the present, shifts in the realized niche may occur due to changes in biotic interactions 30 . For example, Stewart et al. 46 suggested that extinction of competitors owing to environmental change and small patch size in glacial refugia could change community structure and thus alter a species' realized niche. The three species studied here are canopy trees that occupy similar niches (Fig. 4b), implying that competition may exist among them. Such competition could have become much more severe when the EBLF retreated to small glacial refugia because of the low carrying capacity within refugia 47 . It is well known that Fagaceae-dominated vegetation types are climax communities in the EBLF region 48 , indicating that these taxa may have the competitive advantage over other EBLF associates 49 . Therefore, it is reasonable to postulate that M. thunbergii might have been a weak competitor that could have been eliminated from glacial refugia by stronger competitors such as C. tibetana. Thus, the different results from phylogeography (Fig. 1c) and SDMs (Fig. 4a) for M. thunbergii may reflect the fact that it had fewer opportunities to occupy the full extent of its fundamental niche during the LGM because of its poor competitive ability. Future studies using SDMs that account for biotic interactions 50 might provide projected distributions with higher biological realism that might be better reconciled with phylogeographic inferences.

The limitations of this study.
Comparative phylogeography is a powerful tool for addressing community dynamics in distribution and abundance at the regional scale. However, the conclusions of this study might be subject to two limitations. First, because of financial constraints, this study was limited to two cpDNA markers to allow for a large number of individuals and populations to be surveyed. However, a given spatial genetic structure inferred by cpDNA reflects the history of maternal lineages, as recorded by seed movement, whereas pollen flow, as measured by the nuclear genome, may be more extensive than seed movement, and stochastic processes may lead to a wide range of genetic patterns 27 . Therefore, the patterns observed here should ideally be confirmed using nuclear molecular markers. Second, the three species studied here have a variety of congeners in the EBLF of subtropical China 31 , and hybridization among them is a point of concern for phylogeographic and population genetic studies 51 . However, this situation is difficult to avoid because almost all dominant tree genera in EBLF are species-rich taxa 31 . To limit the potential influence of hybridization on our work, we only included samples from individuals exhibiting typical morphological traits identified by experienced botanists, and we abandoned individuals or populations where morphologically similar species co-occurred. Further investigations on the influence of hybridization by comparing phylogeographic structures of congeneric species (e.g., Saeki et al.) 51 in EBLF may complement the conclusions of this study. Despite these caveats, we believe that the basic conclusions of this study are robust because the phylogeographic pattern of multiple refugia in subtropical China observed in C. tibetana and S. superba has been repeatedly revealed in the past (see review in Liu et al.) 36 . Furthermore, the pattern of glacial retreat and post-glacial recolonization, although rarely reported in the past, has been stressed in a few recent studies 37,38 . The main finding of idiosyncratic responses of EBLF constituents in China to the late Quaternary climate changes suggests that the community assemblages of EBLF have changed over time and there is no palaeo-analog to the modern-day EBLF, which explains the discrepancy between PBR and phylogeography.

Materials and Methods
Population sampling. Our sampling strategy for phylogeographic analyses was to collect a few individuals Phylogeographic and population demography analyses. All sequences for each region were edited with sequencher (GeneCodes Corporation, Ann Arbor, MI, USA) and were aligned using clustal_x 1.81 53 . Indels and inversions were treated as single mutations. Haplotype relationships were assessed using the median-joining network method 54 in network 4.1.0.8 (http://www.fluxus-engineering.com).
We used the program haplonst 55 to calculate the total and within-population haplotype diversity (h T and h S ), as well as population differentiation based on ordered (N ST ) and unordered (G ST ) haplotypes. The values of N ST and G ST were then compared using U-statistics to test for the presence of phylogeographic structure. To test the EI model for each species, we calculated the correlation of genetic variation with latitude by performing a multiple regression analysis in spss version 13.0 (SPSS Inc, Chicago, IL). Spatial genetic structure of each species was analyzed by the program samova (Spatial Analysis of Molecular Variation) 56 , which seeks for the best grouping (K) of populations that are geographically homogeneous and maximally differentiated from each other.
Tajima's D 57 and Fu's F S 58 were tested for each species to infer potential population growth and expansion. Mismatch distribution analyses 59 were implemented in arlequin (ver. 3.1) 60 under the model of demographic expansion. The goodness-of-fit was tested with the sum of squared deviations (SDD) between observed and expected mismatch distributions, and raggedness index (H Rag ), using 1000 parametric bootstrap replicates. The expansion time (in generations) for expanding species or groups was calculated using the formula T = τ/2 μkg 59 , where μ is the substitution rate in substitution/site/year, k is the average sequence length used for analysis (C. tibetana: 1826 bp; M. thunbergii: 1450 bp; S. superba: 1399 bp), and g is the generation time in years. The μ of C. tibetana was set as 0.71 × 10 −9 substitutions per site year −1 (s/s/y), which is the mean substitution rate of two chloroplast intergenic spacers (ndhJ-trnF and atpI-atpH) in a closely related genus Fagus we previously studied 61 . This rate is similar to μ of three chloroplast intergenic spacers (trnH-psbA, trnT-trnL, and atpI-atpH) in another Fagaceae member, Quercus glauca (0.96 × 10 −9 or 0.69 × 10 −9 s/s/y) 15 . To approximate the substitution rate of Scientific RepoRts | 6:31044 | DOI: 10.1038/srep31044 chloroplast genomic non-coding regions in M. thunbergii,a phylogenetic tree of Lauraceae was reconstructed based on the chloroplast sequences (rpl16, trnL-F, and psbA-trnH) of species in Fig. 3 of Nie et al. 62 , and the substitution rates across the trees were computed by beast v.1.7.1 using the same settings as that study used. We chose the local substitution rate (4 × 10 −9 s/s/y) of clade D of Fig. 3 in Nie et al. 62 to represent the substitution rate of psbA-trnH and rpl32-trnL (UAG) in M. thunbergii, because substitution rates varied among clades and clade D is a minimally monophyletic group that includes M. thunbergii when only chloroplast sequences were used for phylogenetic reconstruction. Similarly, we reconstructed the phylogenetic tree of Theaceae in Fig. 2 of Li et al. 63 using beast v.1.7.1 to approximate the substitution rate of chloroplast genome in S. superba. The local substitution rate (0.6 × 10 −9 s/s/y) of the clade C1 [see Fig. 2 in Li et al. 63 ] in which the genus Schima was embedded, was adopted in this study. There is no accurate record for the first reproduction age of the three species. For C. tibetana, 25-year was used based on the generation time reported in another Castanopsis species 64 . In Lauraceae, Ocotea tenera takes at least 5 years to reach its reproductive maturity 65 , but the age at maturity of Sassafras was reported as 10 years 66 . Natural populations are likely to take a longer time to reach reproductive maturity. Therefore, we used 10 years as the generation time for M. thunbergii. We used 8 years to approximate the generation time of S. superba based on observations on age of first reproduction of cultivated trees at the arboretum of Jiangxi Agricultural University (Z.Y. Zhang, personal observation).
We used the Bayesian skyline plots (BSPs) method 67 as implemented in beast v. 1.7.1 for estimating fluctuations in the effective population size of each species using the above substitution rates accordingly. MCMC chains were run for 20,000,000 generations for C. tibetana, M. thunbergii and S. superba under the GTR, GTR + G and HKY model chosen by jModelTest 2.1.5 respectively 68 .
Species distribution modeling. We employed the maximum entropy approach (maxent) 69 to predict the distribution of the three species at the present time and at the time of the LGM (0.021-0.018 Ma). Nineteen environmental variables for present and the time of the LGM (MIROC 3.2 scenario) were compiled from the WORLDCLIM database with a resolution of 2.5 arc-min (http://www.worldclim.org) 70 for each environmental layer. To avoid over-fitting of niche models, seven variables with pairwise Pearson correlation coefficients r ≤ 0.70 (annual mean temperature, BIO 1 , mean diurnal range, BIO 2 , isothermality, BIO 3 , temperature seasonality, BIO 4 , annual precipitation, BIO 12 , precipitation seasonality, BIO 15 and precipitation of coldest quarter, BIO 19 ) were retained for subsequent analyses.
SDMs were constructed using 205/127/169 presence records of Castanopsis tibetana, Machilus thunbergii, and Schima superba, respectively, including all sites recorded in the field work of this study and for which specimen records with GPS (Chinese Virtual Herbarium) data are available. Each model was run 10 times using the default parameters (convergence threshold of 10 −5 , maximum iterations of 500 and regularization multiplier of 1) and the following user-selected features: application of a random seed, duplicate presence records removal and logistic probabilities used for the output. 80% of species records were used to train the model and 20% to test the model. Model accuracy was measured by the area under the 'Receiver Operating Characteristic (ROC) Curve' (AUC) 71,72 . A score between 0.7 and 1.0 indicates that the model performs better than random and was considered acceptable discrimination 71 .
We then measured niche differences between each species' ENM following the suggestions of Warren et al. 73 . enmtools, version 1.3, was used to calculate Schoener's D 74 and a standardized version of Hellinger distance (calculated as I) 73 . Both D and I ranged from 0 (no niche overlap) to 1 (identical niches). Next, we conducted an identity test to build niche models based on a set of pseudoreplicates generated from a random sampling of data points pooled for each pair of species. A total of 100 replicates were used to generate the pseudoreplicated data sets. The observed measures of niche similarity between species were compared with the null distribution.