Skip to main content

More opportunities more species: Pleistocene differentiation and northward expansion of an evergreen broad-leaved tree species Machilus thunbergii (Lauraceae) in Southeast China

Abstract

Background

The broad continuum between tropical and temperate floras in Eastern Asia (EAS) are thought to be one of the main factors responsible for a prominent species diversity anomaly of temperate plants between EAS and eastern North America (ENS). However, how the broad continuum and niche evolution between tropical and temperate floras in EAS contributes to lineage divergence and species diversity remains largely unknown.

Results

Population genetic structure, demography, and determinants of genetic structure [i.e., isolation-by-distance (IBD), isolation-by-resistance (IBR), and isolation-by-environment (IBE)] of Machilus thunbergii Sieb. et Zucc. (Lauraceae) were evaluated by examining sequence variation of ten low-copy nuclear genes across 43 populations in southeast China. Climatic niche difference and potential distributions across four periods (Current, mid-Holocene, the last glacial maximum, the last interglacial) of two genetic clusters were determined by niche modelling. North and south clusters of populations in M. thunbergii were revealed and their demarcation line corresponds well with the northern boundary of tropical zone in China of Zhu & Wan. The divergence time between the clusters and demographic expansion of M. thunbergii occurred after the mid-Pleistocene climate transition (MPT, 0.8–1.2 Ma). Migration rates between clusters were asymmetrical, being much greater from north to south than the reverse. Significant effects of IBE, but non-significant effects of IBD and IBR on population genetic divergence were detected. The two clusters have different ecological niches and require different temperature regimes.

Conclusions

The north-south genetic differentiation may be common across the temperate-tropical boundary in southeast China. Divergent selection under different temperature regimes (possibly above and below freezing temperature in winter) could account for this divergence pattern. The broad continuum between tropical and temperate floras in EAS may have provided ample opportunities for tropical plant lineages to acquire freezing tolerance and to colonize the temperate regions during the late-Cenozoic global cooling. Our findings shed deeper insights into the high temperate plant species diversity in EAS.

Peer Review reports

Background

Owing to the general similarities in climate, vegetation types, and floristic composition between Eastern Asia (EAS) and eastern North America (ENA), a prominent species diversity anomaly of temperate plants between them has been bewildering generations of evolutionary biologists and ecologists [1,2,3,4,5,6,7]. This phenomenon is even more intricate in view of the fact that the contemporary floras of the two regions are largely derived from a single paleoflora, the “Boreotropical flora”, that was broadly distributed across the Northern Hemisphere during the early Cenozoic [8, 9]. The anomaly in species diversity has been associated with a few factors in EAS, such as high climatic and topographic heterogeneity, a stronger monsoon climate, and the broad continuum between tropical and temperate floras [3, 4, 7], which have stimulated many hypothesis-based phylogenetic and phylogeographic studies to elucidate the mechanisms underlying this prominent pattern [10,11,12,13].

The India-Eurasia collision since 55–40 million years ago (Ma) leads to the formation of the highest plateau of the world (Qinghai-Tibetan Plateau, QTP) and creates huge topographic and climatic heterogeneity in Eastern Asia [14]. Higher diversification rate has been demonstrated in a wide array of organisms on the QTP and adjacent regions [10, 15, 16]. The India-Eurasia collision also dramatically modified climates in EAS, a major climatic reorganization (from a planetary climate system to a monsoon-dominated pattern) occurred since the Early (23.0–16.0 Ma) to Middle Miocene (16.0–11.6 Ma) [17]. The mild monsoon climate, coupled with the complex topography of EAS, provided refugia for many ancient lineages, which has added species diversity as well as phylogenetic diversity to EAS [3, 12, 18]. Inexplicably, how the broad continuum between tropical and temperate floras contributes to the anomaly between EAS and ENA has received much less attention and seldom been tested explicitly.

According to the theory of niche conservatism, species will always inhabit environments that exhibit some similarities to those of their close relatives, for example, half of flowering plant families are strictly tropical and have not spread out of the tropics (tropical niche conservatism) [19, 20]. For a given tropical lineage, the only way to invade temperate climate regime is niche evolution, i.e., the expansion of niche breadth or specialization for new environment [21] because cold and highly seasonal climate especially freezing temperature in winter can cause lethal injuries in living plant tissues [22, 23] and represent a major factor limiting the dispersal of tropical plants into temperate regions (i.e., ecophysiological barrier) [24, 25]. However, species are not able to adapt to ecological conditions that they are never exposed to, the opportunities being afforded to individual species by their geographical location are essential for niche evolution [21].

It is widely known that a relatively continuous, homogenous flora with many tropical and subtropical taxa (the boreotropical flora) [8] covered almost the entire breadth of Eurasia and North America up to Arctic area during the early Cenozoic [2]. As climate cooled down during the middle and late Cenozoic towards the Pleistocene [26], most plants of the boreotropical flora at higher latitudes either migrated to lower latitudes or went extinct, only a subset of ancestrally tropical plant lineages succeeded to adapt to cold and seasonal climate [25]. Because EAS has a continuous latitudinal gradient of the forest vegetation even during the Pleistocene glaciations, whereas organisms in other parts of the Northern Hemisphere such as Europe and North America were heavily extirpated by thick ice sheets [27]. Therefore, the broad continuum between tropical and temperate floras throughout the late Cenozoic may have provided many more opportunities for the niche evolution of tropical plant lineages [25]. This hypothesis, we coin as “more opportunities more species”, offers a plausible explanation for the contribution of the continuum between tropical and temperate floras to the species richness anomaly of temperate plants between EAS and ENA, but has seldom been tested in any particular lineages.

Machilus Nees is a species-rich genus in Lauraceae containing ca. 100 evergreen tree or shrub species distributed in tropical and warm-temperate (subtropical) south Asia, with most species (ca. 82) occurring in south and southeast China [28]. Many species of this genus are dominant components of evergreen broad-leaved forests in tropical and warm-temperate regions of China [29], covering a latitudinal line that is called tropical-temperate divide in southwest America to Central America and straddling the northern boundary of the tropical zone in China recommended by Chinese meteorologists [30]. Across this divide, a shift in the tolerance of freezing temperature (niche evolution) has been reported in a few plant lineages in southwest America to Central America [31]. Recently, a few phylogeographic studies found a north-south break within the temperate-tropical divide of southeast China [32,33,34,35], implying niche evolution is needed for plants to cross this ecophysiological barrier. In this study, we selected Machilus thunbergii Sieb. et Zucc., a widespread and dominant components in tropical monsoon and warm-temperate evergreen broad-leaved forests [28] as a study system. Our previous phylogeographic study using chloroplast DNA sequences indicates that this species may have experienced extensive post-glacial range expansion [33]. Because a single maternally-inherited gene can hardly capture all major events in a species history, more insightful understandings of its evolutionary history should be gained by using multi-locus sequence data. In this study, we used 10 nuclear loci to analyze the phylogeographic structure and determined the potential drivers of spatial genetic variation of M. thunbergii. In particular, we addressed the question of whether or not a clear north-south differentiation can be observed within this species. Specifically, based on multi-locus sequence data, we applied an approximate Bayesian computation (ABC) procedure to test competing scenarios of population demographic history that has potentially contributed to the contemporary distribution of the species. In addition, we used landscape genetic approaches to examine the relative contributions of geography and ecology to genetic divergence. Together, these analyses would contribute to a better understanding of how the broad continuum and niche evolution between tropical and temperate floras in EAS contribute to lineage divergence and thus species diversity.

Results

Sequence variation and neutrality test of nuclear loci

Ten low-copy nuclear genes were sequenced for all 211 individuals from 43 populations. The total aligned length was 3130 bp, with loci ranging from 221 to 389 bp (Table S1). A 6-bp indel locating in a non-coding region of MT159 was detected and excluded from subsequent analyses. A total of 186 segregating sites, including 52 singletons, were found from 10 nuclear genes.

Average genetic diversity parameters, such as total nucleotide diversity, silent nucleotide diversity, non-synonymous nucleotide diversity, and haplotype diversity, were relatively low at species level (πt = 0.00356, πs = 0.00663, πa = 0.00410, Hd = 0.565). As for the two population clusters (see below), average genetic diversity was much higher in south cluster (πt = 0.00439, πs = 0.00751, πa = 0.00336, Hd = 0.655) than in north cluster (πt = 0.00230, πs = 0.00456, πa = 0.00162, Hd = 0.418) (Table S2). In each cluster, negative and mostly non-significant Tajima’s D values and Fu and Li’s D* and F* values were detected at the majority of 10 nuclear loci (Table S2). There was non-significant deviation from neutrality in each cluster using multilocus Hudson-Kreitman-Aguadé tests. In addition, at both cluster and species levels, MFDM tests rejected the likelihood of natural selection acting on individual loci (results not shown).

Population genetic structure and demographic history

The LnP(D) and ΔK statistics inferred from the Bayesian clustering algorithm (STRUCTURE) showed the most likely number of clusters for the entire dataset was K = 2 (Fig. S1). The two clusters clearly correspond to the geographic distribution: the north cluster (N) comprised of populations 1–32, and the south cluster (S) contained populations 33–43 (Fig. 1a and b). Two population clusters was also found in PCA (Fig. 1c). The demarcation line of the two clusters corresponds well with the northern boundary of tropical zone in China of Zhu & Wan [30] with a few exceptions possibly due to microhabitat differences.

Fig. 1
figure 1

a Colour-coded grouping of the 43 Machilus thunbergii populations according to STRUCTURE with the most likely group number K = 2. The colored pies were drawn on the China map which was downloaded from https://www.webmap.cn/;b histogram of the STRUCTURE assignment for all individuals at K = 2; c principle component analysis on pairwise genetic differentiation of populations. The red dotted line shows the northern boundary of tropical zone in broad sense [30, 36]

The divergence time between clusters N and S was estimated to 0.58 million years ago [95% highest posterior density (HPD) intervals: 0.41–0.71 Ma] using IMa2 (Table S3). This time is slightly older than the divergence time estimated in ABC modeling [0.46 Ma, 95%CI (0.085–0.9 Ma)]. Migration rates (m) estimated under the IM model using IMa2 were asymmetrical, with much higher gene flow from N to S (1.468) than in the reverse direction (0.456) (Fig. 2, Table S3).

Fig. 2
figure 2

Posterior probability distributions of (a) the migration rate (m) between two clusters (N and S) from both directions, (b) effective population size (θ) for each cluster (N, S and all populations), and (c) divergence time (t) between two clusters (N and S) estimated by IM model

Among eight potential demographic scenarios in ABC modeling (Fig. 3), M. thunbergii as a whole, clusters N and S all conformed to scenario 2 (91.03, 61.58, 65.46% posterior probability, respectively; Table S4, Fig. S2), which indicates stepwise population expansions. One expansion time (T3) was estimated at 0.68 Ma (95% CI: 0.25–0.98 Ma), 0.52 Ma (95% CI: 0.12–0.94 Ma) and 0.46 Ma (95%CI: 0.09–0.93 Ma), respectively, and the other (T1) was at 0.10 Ma (95% CI: 0.02–0.19 Ma), 0.15 Ma (95% CI: 0.02–0.28 Ma) and 0.07 Ma (95% CI: 0.01–0.14 Ma), respectively (Table 1, Fig. S3). The LAMARC analysis also suggested that the two clusters have undergone changes in population size. For clusters N and S, the estimated g values were 3005.49 and 7234.10, respectively, suggestive of population growth. The EBSP simulations indicated the two clusters and M. thunbergii as a whole all experienced continuous recent expansions (Fig. S4).

Fig. 3
figure 3

Eight assumed demographic scenarios tested for clusters N, S, and A in DIYABC. N1 and N2 represent current population sizes, and NA represents the ancestral population size. N1a, N1b, N1c, N2a, N2b and N2c represent population sizes between the ancestral population and the current population. t1, t2 and t3 represent times of population changes

Table 1 Posterior estimates of demographic parameters for the best models of population demography in Fig. 3 revealed by Approximate Bayesian Computation

IBD, IBE and IBR

The multiple matrix regression with randomization (MMRR) analysis revealed significant effects of isolation by environment (IBE) on population genetic divergence (FST) (βIBE = 0.482, P = 0.001) but non-significant effects of isolation by distance (IBD) (βIBD = − 0.037, P = 0.636) and isolation by resistance (IBR) (βIBR = − 0.083, P = 0.405). The principal component analysis found clear climatic differentiation among two genetic clusters (N and S) in the first two principal components (PC1: 51.63%, PC2: 28.13%) (Fig. S5). Frequency distributions of N and S clusters for each environmental variable are shown in kernel density plots (Fig. 4). Nonparametric Kruskal-Wallis multiple-range tests suggested that 12 of the 19 bioclimatic variables, except for Bio2 (mean diurnal range), Bio5 (max temperature of warmest month), Bio10 (mean temperature of warmest quarter), Bio12 (annual precipitation), Bio13 (precipitation of wettest month), Bio16 (precipitation of wettest quarter), and Bio18 (precipitation of warmest quarter), distinguished significant differences between two clusters (Fig. 4).

Fig. 4
figure 4

Kernel density plots for 19 environmental parameters of the northern (blue curves) and southern (pink curves) clusters of M. thunbergii. Differentiation among groups was evaluated by a nonparametric Kruskal-Wallis multiple-range test, and the results are indicated in each plot. A lack of significant difference (at the P < 0.05 level) is indicated by an equal sign, while significant differences are coded as either higher or lower

Ecological niche modeling and niche identity tests

ENMs were constructed for the two M. thunbergii clusters to predict their current potential distributions and the model was then projected to past scenarios (Fig. 5). Nine bioclimatic variables (Bio2: mean diurnal range, Bio5: max temperature of warmest month, Bio6: min temperature of coldest month, Bio7: temperature annual range, Bio8: mean temperature of wettest quarter, Bio13: precipitation of wettest month, Bio15: precipitation seasonality, Bio18: precipitation of warmest quarter, Bio19: precipitation of coldest quarter) were retained with r < 0.7 in each pair. In the MAXENT climate models, both groups had high AUC values, indicating that the climate model predicted occurrence well (northern cluster test data AUC = 0.986, SD = 0.007 and southern cluster test data AUC = 0.989, SD = 0.899). Observed measures of niche similarity (D and I) were higher than null distributions for N vs. S, suggesting high niche differentiation between N and S (Fig. 5). A jackknife test of the AUC values for both groups indicated the most effective climate variables for predicting occurrence were minimal temperature of coldest month (northern cluster AUC = 0.952; southern cluster AUC = 0.940) and precipitation of coldest quarter (northern cluster AUC = 0.968; southern cluster AUC = 0.955) (Table S5).

Fig. 5
figure 5

Current potential distributions of two Machilus thunbergii clusters (N and S) predicted by MAXENT, and identity tests results between two groups. The black bars indicate the null distributions of D, while the grey bars indicate I. Arrows indicate values of D (black) and I (grey) in actual MAXENT runs

The predicted present-day, LIG, LGM and MH distributions were nearly identical for the northern group (Fig. S6). For southern group, clear range expansions were predicted for the LGM and MH model relative to the present day (Fig. S6).

Discussion

Late Quaternary divergence history of Machilus thunbergii

In this study, we found a clear north-south population differentiation in an evergreen broadleaved tree species Machilus thunbergii belonging to Lauraceae in southeast China. This pattern has been repeatedly observed in other plant species (e.g., Ficus pumila (Moraceae), Schima superba (Theaceae), Lindera aggregata (Lauraceae), Oreocharis auricula (Gesneriaceae)) [32,33,34,35]. Similarly, in North and Middle America (Mexico and Central America), a north-south phylogeographic break across the so-called tropical-temperate divide was also reported at the similar latitudes (ca. 15–33°N) in live oaks (Quercus series Virentes) [31]. Although the number of study cases is still limited, a common north-south differentiation of plant species straddling the tropical-temperate divide is beginning to emerge in southeast China as well as in North and Middle America. Interestingly, such a steep transition across the tropical-temperate divide has also been observed in several functional traits such as seed size and plant height [37, 38], indicating that there might be a genetically-grounded switch in plant strategy between temperate and tropical zones [38, 39].

The tropical niche conservatism theory predicts that tropical plant lineages may have spread to temperate regions recently only after a shift from tropical to temperate niche had happened during the late-Cenozoic global cooling [19, 20]. In line with the prediction, we found that the north-south separation of M. thunbergii happened between 0.58–0.46 Ma and demographic expansions after 0.52 Ma for each population cluster. In addition, we also detected that much lower genetic diversity in north cluster than that in south cluster, conforming the prediction that temperate lineages are often recently derived from lineages in tropical regions, leading to shallower phylogenetic divergences among temperate lineages than among tropical lineages [21].

Although conforming to theoretical prediction, the separation time of north-south groups is substantially younger than the ages (~ 30–40 Ma) of many temperate plant clades that originated in the tropics [24] when temperate regions expanded [40, 41]. This difference may be simply explained by that intraspecific divergence is often much shallower than interspecific diversification. However, only after the climate changes since the late Pliocene, is it possible for plant species in southeast China to adapt to temperate climate, because the temperate region was confined to relatively high latitudes for most periods of Cenozoic, with tropical climate dominating central to south China until late Miocene [42]. However, global temperatures began to drop sharply during the Pliocene, about 3 Mya, and then, in the Pleistocene, underwent violent fluctuations at intervals of about 100,000 years [26]. During the Last Maximum Glaciation (LGM), for instance, the mean annual temperature dropped by c. 4–6 °C [43] and the evergreen broad-leaved forests was forced to the south of 24°N in China [44].

Note that the time of population differentiation and demographic expansion of M. thunbergii is well after the mid-Pleistocene climate transition (MPT, 0.8–1.2 Ma) [45], a period marked by an increase in the severity of glaciations and the emergence of the ~ 100-kyr glacial cycles [46]. It was reported that the earliest and largest glaciation in China (Wangkun glaciation, 700–500 ka) [47] occurred shortly after the MPT. The ice sheet during Wangkun glaciation was 18 times larger than the present ice sheet on the Qinghai-Tibetan Plateau [48]. Demographic expansions have been reported in several plant species of central to southeast China shortly after the MPT [12, 49]. It is possible that the largest glaciation in China might have exposed the ancestral tropical populations of M. thunbergii to cold climate and the tolerance of freezing temperature might be acquired by some descendants during 0.58–0.46 Ma, producing genetically and possibly geographically distinct populations. In the subsequent warm and humid interglacial (0.50–0.46 Ma) [47], the north populations which had adapted to temperate climate could have experienced strong demographic expansion as suggested in our previous study [33]. Surprisingly, in contrast to the demographic stability of the south populations of another Lauraceae species (Lindera aggregata) [34] in the same region, southern populations of M. thunbergii could have also experienced demographic expansion posterior to the expansion of north populations (0.46 Ma vs 0.52 Ma). This event led to secondary contacts of the diverging lineages and a pattern of asymmetrical gene flow towards south populations, a phenomenon that frequently occurs from the local to the invading lineages [50, 51].

Quantifying the effects of landscape, geography and ecology in spatial genetic divergence of Machilus thunbergii

Physical barriers such as mountains, rivers, inhospitable habitats or even man-made buildings are the most common and straightforward factors responsible for pronounced population differentiation [10, 52,53,54,55,56]. In southeast China, there are several major mountains stretching from west to east or from southeast to northeast (e.g., Nanling Mountains, Wuyi-Xianxia Mountains) which might serve as physical barriers to gene flow. For example, the uplifting of Wuyi-Xianxia Mountains during the Pliocene has been inferred as the causes for the separation of Ficus pumila populations to the north and south of the mountains [32]. However, we did not find significant isolation-by-resistance in M. thunbergii that account for landscape heterogeneity, suggesting the major mountains in southeast China can not impede gene flow of M. thunbergii effectively. Indeed, the Nanling Mountains and Wuyi-Xianxia Mountains, with peaks rarely exceeding 2000 m above sea level, have been frequently proved to be glacial refugia and/or dispersal corridors, rather than to be physical barriers in previous phylogeographic studies [49, 57, 58]. In addition, we did not find significant isolation-by-distance, indicating that genetic drift alone can not produce the north-south differentiation of M. thunbergii, either.

On the contrary, significant effect of isolation-by-environment was detected by the multiple matrix regression with randomization (MMRR) analysis (βIBE = 0.482, P = 0.001). This result indicates that divergent selection associated with environmental heterogeneity, especially changes in climatic regimes, could be the driving force responsible for the population differentiation of M. thunbergii. Many mechanisms such as (natural selection against immigrants, sexual selection against immigrants, reduced hybrid fitness and biased dispersal) can generate IBE [56], direct evidence for local adaptation to different climatic conditions should be gained through common garden experiments [59] or through genomic scanning of the loci underlying natural selection [60]. Although this study did not compare the performance (such as freezing tolerance) of the two provenances of M. thunbergii in common gardens and the genomic loci responsible for local adaptation can not be deciphered due to the paucity of nuclear loci, several lines of evidence suggest that the north-south population differentiation could have been mainly caused by local adaptation to different climatic conditions.

First, it is generally assumed that loci under divergent selection and those tightly physically linked to them may exhibit stronger differentiation than neutral regions with weak or no linkage to such loci. However, divergent selection can also increase genome-wide neutral differentiation by reducing gene flow, thus promoting divergence via the stochastic effects of genetic drift [61]. In addition, neutral loci (i.e., loosely linked loci) that are far along a chromosome from a selected site may be influenced by selection because hitchhiking effects can extend a considerable distance from the selected locus [62,63,64]. Therefore, although the loci we examined may not be the targets of divergent selection, the elevated genetic differentiation between the north and south clusters suggests that population divergence could be caused by hitchhiking effects on loosely linked loci.

Second, climatic adaptation is a highly important component in the evolution of plants with temperature potentially the most important selective agent over an altitudinal/latitudinal cline [65]. Freezing temperatures can cause lethal injuries in living plant tissues [22] and have long been considered a principal restriction on the spreading out of tropics for many tropical pant lineages [25]. This implies that migration from tropical to temperate regions is not easy which may greatly restructure metabolic, morphological, defensive and phenological traits. For a particular tropical plant lineage, an opportunity opposing to divergent selective regimes of temperature is required for generating and maintaining genetic adaptation to temperate climate. The distribution range of M. thunbergii spans the tropical-temperate boundary suggested by Zhu & Wan [30] (Fig. 1a), south of which there is no frost damage in winter [30, 36]. Although the exact freezing temperature of M. thunbergii remains to be determined, it is very likely that the north and south populations of M. thunbergii could have been experiencing different freezing regimes in winter. In addition, the importance of minimal temperature of coldest month in predicting the distribution of M. thunbergii in niche modeling (Table S4) and significant differences of many temperature variables (Fig. 4) in nonparametric Kruskal-Wallis tests also support the strong role of different temperature regimes in driving population divergence of M. thunbergii.

Third, to the south of Zhu & Wan’s Line, the frequency of tropical genera can reach 60% of the regional flora, however, the frequency decrease sharply to the north of the line [66]. For example, Ficus microcarpa, a dominated tree species in tropical forests of southeast China, reach south Jiangxi and central Fujian (25–26°N) [67]. However, transplanted individuals to the north of 25°N always suffer from freezing damage and can not survive the winter outdoors (personal observation). These facts strongly suggest that there is an ecophysiological barrier in southeast China that is associated with the freezing temperature of plants to prevent tropical plants from invasion into temperate zones. Successful colonizers such as the north cluster of M. thunbergii should have gained a new ecological niche via genetic modifications in their evolutionary history.

Implications for understanding plant species diversity anomaly between EAS and ENA

Although population differentiation can not be translated into speciation, the genetic divergence among populations within a species will ultimately lead to speciation or incomplete speciation [68, 69]. Thus, the common phenomenon of north-south differentiation observed in M. thunbergii and other plant species in southeast China may have important implications for understanding plant species diversity anomaly between EAS and ENA.

First, while there are a wealth of tropical species in southeast Asia, they per se may contribute little to the temperate species pool of East Asia, because the freezing temperature represents an ecophysiological barrier for them to colonize the temperate zone. For instance, the distribution range of many tropical lineages of southeast Asia, such as fig trees, seldom extends to the north of 25° N [66]. Second, EAS is characterized by a continuous landmass without being interrupted by sea. This physiographical setting may translate into an amplitude of opportunities for tropical plant lineages to adapt to temperate climate in EAS during the late-Cenozoic global climate cooling. However, such continuum is interrupted in North America by the Gulf of Mexico (and in Europe by the Mediterranean and the Sahara Desert) [3]. Because of the narrow continuum between temperate and tropical floras, north-south genetic differentiation across the tropical-temperate divide of North America has rarely been documented [31], suggesting far less opportunities for the niche evolution and speciation of tropical plant lineages in North America. Third, the commonality and recency of north-south population differentiation imply that the recruitment of tropical plant species into species pool of EAS through adaptive evolution is an ongoing process which might further contribute to the species diversity anomaly between EAS and ENA in the future. Therefore, southeast China is a key area deserving a high conservation priority not only for its high plant species diversity [70, 71], but also for the evolutionary processes that generate the diversity.

Conclusions

Using multi-locus phylogeographic and landscape genetic approaches, this study revealed a north-south genetic differentiation in M. thunbergii and the demarcation line of north-south population clusters corresponds well with the northern boundary of tropical zone in China of Zhu & Wan [30]. Such a divergence pattern may be common across the temperate-tropical boundary in southeast China that could be associated with divergent selection under different temperature regimes (possibly above and below freezing temperature in winter). The results have profound implications for understanding the prominent species diversity anomaly of temperate plants between EAS and ENS in that the broad continuum between tropical and temperate floras in EAS may have provided ample opportunities for tropical plant lineages to acquire freezing tolerance and to colonize the temperate regions during the late-Cenozoic global cooling. Future studies are needed to gain direct evidence for local adaptation to different climatic conditions through common garden experiments or through genomic scanning of the loci underlying natural selection.

Methods

Plant material

In total, leaf tissues from 211 individuals of Machilus thunbergii were sampled from 43 wild populations in Southeast China (Table S6, Fig. 1) and were immediately dried with silica gel. Sampled individuals were spaced at least 50 m for each population. All voucher specimens were deposited in the Herbarium of Lushan Botanical Garden of the Chinese Academy of Sciences (Table S6).

DNA extraction, amplification and sequencing

Total genomic DNA was extracted from silica-dried leaves using a modified cetyltrimethylammonium bromide (CTAB) method [72]. Ten low-copy nuclear genes (MT02, MT04, MT15, MT33, MT55, MT57, MT96, MT115, MT159 and MT164) were developed from transcriptome sequences of M. thunbergii (Table S1) by the method of our previous studies [73, 74]. Transcriptome sequencing was performed on the Illumina NextSeq 500 platform (Illumina, San Diego, CA, USA) at Personal Biotechnology (Shanghai, China). The sequences of transcriptome unigenes were deposited at figshare (https://figshare.com/articles/dataset/Unigene_fa/15090174).

The PCR reaction systems and procedures were the same as those in our previous study [33]. Different annealing temperatures were set for each locus (Table S1). Sequencing reactions were carried out with the corresponding primers in both directions commercially by Sangon Biotech Co., Ltd. (Shanghai, China). All sequences were edited with Sequencher 4.14 (GeneCodes Corporation, Ann Arbor, MI, USA), aligned using Bioedit 7.2 [75] and refined manually in Mega 5.05 [76]. All DNA sequences were deposited in the NCBI GenBank database (accession Nos: MW856287–MW856301, MW872423–MW872637).

Nucleotide diversity and neutrality tests

For each nuclear locus, the sequences were phased using the phase algorithm in Dnasp 5.10 [77]. Then we used Dnasp 5.10 [77] to estimate the basic population genetic parameters, i.e., the number of segregating sites (S), nucleotide diversity (π), Watterson’s parameter (θw) [78], the number of haplotypes (Nh), haplotype diversity (Hd) and the minimum number of recombinant events (Rm). We also calculated Tajima’s D [79] and Fu and Li’s D and F [80] to test the neutral evolution of each locus. The multilocus Hudson–Kreitman–Aguadé test (HKA) [81] was further performed in Dnasp 5.10 to assess the fit of data to neutral equilibrium using Machilus villosa (Roxb.) Hook. f. as an outgroup. Finally, we used the maximum frequency of derived mutations (MFDM) test [82] to examine the likelihood of natural selection acting on individual loci.

Population genetic structure

Population structure was assessed using Structure 2.3.4 [83] with the admixture model. Twenty independent runs were conducted for each number of clusters (K) from 1 to 10 with burn-in setting to 20,000 and Markov chain Monte Carlo (MCMC) run length to 200,000. The LnP(D) [84] and ΔK statistics [85] were used to estimate the most likely number of clusters. Graphics were visualized by Distruct 1.1 [86].

Population divergence and demography

We used the isolation by migration (IM) model to estimate gene flow and divergence time between genetic clusters in the IMa2 software package [87, 88]. The run sets were the same with our previous study [73]. Parameters were estimated based on a mean mutation rate of 3.4×10− 9 per site per year, which is estimated in Lindera aggregata (Lauraceae) [34], another species in Lauraceae. Average generation time was set to 10 years as in previous study [33].

The genetic divergence time between genetic clusters was also tested using ABC approach implemented in Diyabc 2.1.0 [89]. In addition, eight demographic scenarios (Fig. 3) were tested using ABC for cluster N, cluster S and all populations (A), respectively: (1) one expansion (t1), (2) two expansions (t3, t1), (3) one contraction (t3) and one expansion (t1), (4) one expansion (t3), one contraction (t2), and then another expansion (t1), (5) one contraction (t1), (6) two contractions (t3, t1), (7) one expansion (t3) and one contraction (t1), and (8) one contraction (t3), one expansion (t2), and another contraction (t1). Another Bayesian method was also implemented by Lamarc 2.1.8 to estimate the population size changes [90]. The performance settings for DIYABC and Lamarc followed our previous study [73]. DIYABC selects the optimized rate from a given mutation range between 10− 7 and 10− 9. All summary statistics were chosen to test the posterior probability parameters. To ensure statistically robust results, at least 8, 000, 000 simulated datasets were generated for each scenario. In order to identify the best-supported scenario, we estimated posterior probabilities [with 95% confidence intervals (CIs)] of each scenario using logistic regression of 1% of simulated data sets closest to the observed data.

Furthermore, we used extended Bayesian skyline plots (EBSP) analysis in Beast 2.5.1 [91] to infer the temporal changes in the effective population size (Ne) changes of the Structure clusters. We used the default substitution model and the substitution rate of 3.4 × 10− 9 per site per year [34]. We set the MCMC to 20,000,000 steps with sampling every 20,000 steps, and set the burn-in at 20%.

Determinants of population genetic structure

Isolation by distance (IBD), isolation by environment (IBE) and isolation by resistance (IBR) has been considered as the main forces for driving the population structure across a species’ range [92]. Isolation by distance correlates the genetic differences and geographical distances between populations or subgroups [93]. IBE describes the accumulation of genetic divergence when environmental differences between sampled populations increase [56]. Resistance distances are often the outcome of physically isolating barriers (such as mountains or rivers) that can structure spatial genetic divergence [94]. To understand the determinants of population genetic structure, we firstly generated four distance matrices with all populations included. Genetic distances (FST) among each pair of populations were calculated using Arlequin 3.01 [95]; A geographical distance matrix was calculated using the GPS coordinates of populations. We used the sum option in the Ressistance And Habitat Calculator tool in Arcmap 10.8.1 to calculate resistance values from raster layers of elevation, roads, water bodies and rivers. Then the pairwise resistance distance matrix were generated using Circuitscape 4.0 [94]. Nineteen bioclimatic variables were extracted from current climatic data with 30 arc-seconds resolution (http://www.worldclim.org) and used as proxies of environmental gradients. A climatic distance matrix was constructed using population pairwise Euclidean distances across all 19 bioclimatic variables. Then, after standardizing the four distance matrices by Z scores, we quantified the contributions of geographical distance, resistance distance and environmental distance to genetic distance using the ‘Multiple Matrix Regression with Randomization (MMRR)’ function in R 3.6.3 [96].

In this study, we detected significant isolation-by-environment (see Results). To further detect which climatic variable is associated with the genetic differentiation between the two clusters, we conducted a principal component analysis (PCA) on climate data (19 climate variables) obtained from Worldclim (http://www.worldclim.org) located within the sampling areas of the two groups. In addition, we used nonparametric Kruskal–Wallis multiple-range test [97] to examine the difference of each environmental variable between the two clusters. The distribution of each environmental variable was displayed in kernel density plots. We performed these analyses and drew graphical illustrations in R 3.6.3.

Ecological niche modeling

Current potential distributions for northern and southern clusters were predicted using ecological niche modeling (ENM) in Maxent 3.3.4 [98] with the settings as follows: replicates = 20; type = subsample; maximum iterations = 5000; random test points = 25. Nineteen bioclimatic layers with a 2.5 arc minute resolution were downloaded from the WorldClim database. For the mid-Holocene (MH) and last glacial maximum (LGM), the last interglacial (LIG), we employed paleoclimate data based on MIROC-ESM model. The geographical coordinates of 32 sampling sites from northern cluster and 11 from southern cluster were used for developing niche models. The pairwise Pearson correlations (r) among the bioclimatic variables were examined to avoid multicollinearity, and we excluded one of the variables in each pair with r > 0.7. The AUC (area under the receiver operating characteristic curve) values were used to measure model accuracy, and the model with values above 0.9 was considered better discrimination. We also conducted the Jackknife tests to identify the contribution of each climatic variable. The niche differences between the two groups were measured by niche overlap and identity tests in Enmtools 1.4.4 using statistics of Schoener’s D and standardized Hellinger distance (I) [99].

Availability of data and materials

The sequences of transcriptome unigenes were deposited at figshare (https://figshare.com/articles/dataset/Unigene_fa/15090174). All DNA sequences of 10 nuclear loci have been deposited in the NCBI GenBank database under accession numbers MW856287–MW856301, MW872423–MW872637. The datasets supporting the conclusions of this article are included within the article and its additional files. The datasets used and/or analyzed during the current study are available from the authors on reasonable request (Dengmei Fan, dmf.625@163.com).

Abbreviations

EAS:

Eastern Asia

ENA:

Eastern North America

QTP:

Qinghai-Tibetan Plateau

ABC:

Approximate Bayesian computation

MPT:

Mid-Pleistocene climate transition

References

  1. White PS. Eastern Asian-eastern north American floristic relations: the plant community level. Ann Mo Bot Gard. 1983;70:734–47.

    Article  Google Scholar 

  2. Latham RE, Ricklefs RE. Global patterns of tree species richness in moist forests: energy-diversity theory does not account for variation in species richness. Oikos. 1993;67:325–33.

    Article  Google Scholar 

  3. Axelrod DI, Al-Shehbaz I, Raven PH. History of the modern flora of China. In: Zhang A, Wu S, editors. Floristic characteristics and diversity of east Asian plants. Beijing: Higher Education Press; Berlin, Germany: Springer-Verlag; 1996. p. 43–55.

    Google Scholar 

  4. Qian H, Ricklefs RE. Large-scale processes and the Asian bias in species diversity of temperate plants. Nature. 2000;407:180–2.

    Article  CAS  PubMed  Google Scholar 

  5. Qian H, Fridley JD, Palmer MW. The latitudinal gradient of species-area relationships for vascular plants of North America. Am Nat. 2007;170:690–701.

    Article  PubMed  Google Scholar 

  6. Qian H, White PS, Song JS. Effects of regional vs. ecological factors on plant species richness: an intercontinental analysis. Ecology. 2007b;88:1440–53.

    Article  PubMed  Google Scholar 

  7. Qian H, Jin Y, Ricklefs RE. Phylogenetic diversity anomaly in angiosperms between eastern Asia and eastern North America. Proc Natl Acad Sci U S A. 2017;114:11452–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Wolfe JA. Some aspects of plant geography of the northern hemisphere during the late cretaceous and tertiary. Ann Mo Bot Gard. 1975;62:264–79.

    Article  Google Scholar 

  9. Tiffney BH. Perspectives on the origin of the floristic similarity between eastern Asia and eastern North America. J Arnold Arbor. 1985;66:73–94.

    Article  Google Scholar 

  10. Yuan QJ, Zhang ZY, Peng H, Ge S. Chloroplast phylogeography of Dipentodon (Dipentodontaceae) in Southwest China and northern Vietnam. Mol Ecol. 2008;17:1054–65.

    Article  CAS  PubMed  Google Scholar 

  11. Guan BC, Fu CX, Qiu YX, Zhou SL, Comes HP. Genetic structure and breeding system of a rare understory herb, Dysosma versipellis (Berberidaceae), from temperate deciduous forests in China. Am J Bot. 2010;97:111–22.

    Article  CAS  PubMed  Google Scholar 

  12. Kou YX, Cheng SM, Tian S, Li B, Fan DM, Chen YJ, et al. The antiquity of Cyclocarya paliurus (Juglandaceae) provides new insights into the evolution of relict plants in subtropical China since the late early-Miocene. J Biogeogr. 2016;43:351–60.

    Article  Google Scholar 

  13. Yu XQ, Gao LM, Soltis DE, Soltis PS, Yang JB, Fang L, et al. Insights into the historical assembly of east Asian subtropical evergreen broadleaved forests revealed by the temporal history of the tea family. New Phytol. 2017;215:1235–48.

    Article  CAS  PubMed  Google Scholar 

  14. Harris N. The elevation history of the Tibetan plateau and its implications for the Asian monsoon. Palaeogeogr Palaeoclimatol Palaeoecol. 2006;241:4–15.

    Article  Google Scholar 

  15. Xing Y, Ree RH. Uplift-driven diversification in the Hengduan Mountains, a temperate biodiversity hotspot. Proc Natl Acad Sci U S A. 2017;114:E3444–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Favre A, Päckert M, Pauls SU, Jähnig SC, Uhl D, Michalak I, et al. The role of the uplift of the Qinghai-Tibetan plateau for the evolution of Tibetan biotas. Biol Rev. 2015;90:236–53.

    Article  PubMed  Google Scholar 

  17. Guo ZT, Sun B, Zhang ZS, Peng SZ, Xiao GQ, Ge JY, et al. A major reorganization of Asian climate by the early Miocene. Clim Past. 2008;4:153–74.

    Article  Google Scholar 

  18. Chou YW, Thomas PI, Ge XJ, LePage BA, Wang CN. Refugia and phylogeography of Taiwania in East Asia. J Biogeogr. 2011;38:1992–2005.

    Article  Google Scholar 

  19. Wiens JJ, Graham CH. Niche conservatism: integrating evolution, ecology, and conservation biology. Ann Rev Ecol Evol S. 2005;36:519–39.

    Article  Google Scholar 

  20. Wiens JJ, Brandley MC, Reeder TW. Why does a trait evolve multiple times within a clade? Repeated evolution of snake-like body form in squamate reptiles. Evolution. 2006;60:123–41.

    PubMed  Google Scholar 

  21. Wiens JJ, Donoghue MJ. Historical biogeography, ecology, and species richness. Trends Ecol Evol. 2004;19:639–44.

    Article  PubMed  Google Scholar 

  22. Sakai A, Larcher W. Frost survival of plants: responses and adaptation to freezing stress. Berlin: Springer-Verlag; 1987.

    Book  Google Scholar 

  23. Woodward FI, Williams BG. Climate and plant distribution at global and local scales. Vegetatio. 1987;69:189–97.

    Article  Google Scholar 

  24. Mittelbach GG, Schemske DW, Cornell HV, Allen AP, Brown JM, Bush MB, et al. Evolution and the latitudinal diversity gradient: speciation, extinction and biogeography. Ecol Lett. 2007;10:315331.

    Article  Google Scholar 

  25. Donoghue MJ. A phylogenetic perspective on the distribution of plant diversity. Proc Natl Acad Sci U S A. 2008;105:11549–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Zachos J, Pagani M, Sloan L, Thomas E, Billups K. Trends, rhythms, and aberrations in global climate 65 Ma to present. Science. 2001;292:686–93.

    Article  CAS  PubMed  Google Scholar 

  27. Pielou EC. After the ice age. Chicago: University of Chicago Press; 1991.

    Book  Google Scholar 

  28. Li XW, Li J, Huang PH, Wei FN, Lauraceae. In: Wu ZY, Raven P, editors. Flora of China. Beijing: Science Press; St. Louis: Missouri Botanical Garden Press; 2008. p. 102–254.

    Google Scholar 

  29. Wu CY. Vegetation of China. Beijing: Science Press; 1980.

    Google Scholar 

  30. Zhu KZ, Wan MW. Phenology. Beijing: Scientific Education Press; 1963.

    Google Scholar 

  31. Cavender-bares J, Gonzalez-Rodriguez A, Pahlich A, Koehler K, Deacon N. Phylogeography and climatic niche evolution in live oaks (Quercus series Vitentes) from the tropics to the temperate zone. J Biogeogr. 2011;38:962–81.

    Article  Google Scholar 

  32. Chen Y, Compton SG, Liu M, Chen XY. Fig trees at the northern limit of their range: the distributions of cryptic pollinators indicate multiple glacial refugia. Mol Ecol. 2012;21:1687–701.

    Article  PubMed  Google Scholar 

  33. Fan DM, Hu W, Li B, Morris AB, Zheng M, Soltis DE, et al. Idiosyncratic responses of evergreen broad-leaved forest constituents in China to the late Quaternary climate changes. Sci Rep. 2016;6:31044.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Ye JW, Li DZ, Hampe A. Differential quaternary dynamics of evergreen broadleaved forests in subtropical China revealed by phylogeography of Lindera aggregata (Lauraceae). J Biogeogr. 2019;46:1112–23.

    Article  Google Scholar 

  35. Sun ZX. Comparative phylogeography of Oreocharis auricula and O. henryana. CAS Thesis & Dissertation. Database. 2020; http://sciencechina.cn.

  36. Tang YL. An analysis of boundaries of tropical zone. Acta Phytoecol Geobot Sin. 1964;2:135–7.

    Google Scholar 

  37. Moles AT, Ackerly DD, Tweddle JC, Dickie JB, Smith R, Leishman MR, et al. Global patterns in seed size. Glob Ecol Biogeogr. 2007;16:109–16.

    Article  Google Scholar 

  38. Moles AT, Warton DI, Warman L, Swenson NG, Laffan SW, Zanne AE, et al. Global patterns in plant height. J Ecol. 2009;97:923–32.

    Article  Google Scholar 

  39. Swarts K, Gutaker RM, Benz B, Blake M, Bukowski R, Holland J, et al. Genomic estimation of complex traits reveals ancient maize adaptation to temperate North America. Science. 2017;357:512–5.

    Article  CAS  PubMed  Google Scholar 

  40. Behrensmeyer AK, Damuth JD, DiMichele WA, Potts R, Sues HD, Wing SL. Terrestrial ecosystems through time: evolutionary paleoecology of terrestrial plants and animals. Chicago: University of Chicago Press; 1992.

    Google Scholar 

  41. Morley RJ. Origin and evolution of tropical rain forests. Chichester: John Wiley and Sons; 2000.

    Google Scholar 

  42. Fine PVA, Ree RH. Evidence for a time-integrated species-area effect on the latitudinal gradient in tree diversity. Am Nat. 2006;168:796–805.

    Article  PubMed  Google Scholar 

  43. Zheng Z, Yuan BY, Petit-Maire N. Palaeoenvironments in China during the last glacial maximum and the Holocene optimum. Episodes. 1998;21:152–8.

    Article  Google Scholar 

  44. Harrison SP, Yu G, Takahara H, Prentice IC. Palaeovegetation diversity of temperate plants in East Asia. Nature. 2001;413:129–30.

    Article  CAS  PubMed  Google Scholar 

  45. Tzedakis PC, Raynaud D, Mcmanus JF, Berger A, Brovkin V, Kiefer T. Interglacial diversity. Nat Geosci. 2009;2:751–5.

    Article  CAS  Google Scholar 

  46. Clark PU, Archer D, Pollard D, Blum JD, Rial JA, Brovkin V, et al. The middle Pleistocene transition: characteristics, mechanisms, and implications for long-term changes in atmospheric pCO2. Quat Sci Rev. 2006;25:3150–84.

    Article  Google Scholar 

  47. Cui ZJ, Chen YX, Zhang W, Zhou SZ, Zhou LP, Zhang M, et al. Research history, glacial chronology and origins of quaternary glaciations in China. J Quat Sci. 2011;31:749–64.

    Google Scholar 

  48. Shi YF, Li JJ, Li BY. Uplift and environmental changes of Qinghai-Tibetan plateau in the late Cenozoic. Guangzhou: Guangdong Science and Technology Press; 1998.

    Google Scholar 

  49. Tian S, Lei SQ, Hu W, Deng LL, Meng QL, Li B, et al. Repeated range expansions and inter−/postglacial recolonization routes of Sargentodoxa cuneata (Oliv. ) Rehd. Et Wils. (Lardizabalaceae) in subtropical China. Mol Phylogenet Evol. 2015;85:238–46.

    Article  PubMed  Google Scholar 

  50. Currat M, Ruedi M, Petit RJ, Excoffier L. The hidden side of invasions: massive introgression by local genes. Evolution. 2008;62:1908–20.

    PubMed  Google Scholar 

  51. Petit RJ, Excoffier L. Gene flow and species delimitation. Trends Ecol Evol. 2009;24:386–93.

    Article  PubMed  Google Scholar 

  52. Su H, Qu LJ, He K, Zhang Z, Wang J, Chen Z, et al. The Great Wall of China: a physical barrier to gene flow? Heredity. 2003;90:212–9.

    Article  CAS  PubMed  Google Scholar 

  53. Zhang ZY, Zheng XM, Ge S. Population genetic structure of Vitex negundo (Verbenaceae) in three-gorge area of the Yangtze River: the river barriers to seed dispersal in plants. Biochem Syst Ecol. 2007;35:506–16.

    Article  CAS  Google Scholar 

  54. Bai WN, Liao WJ, Zhang DY. Nuclear and chloroplast DNA phylogeography reveal two refuge areas with asymmetrical gene glow in a temperate walnut tree from East Asia. New Phytol. 2010;188:892–901.

    Article  PubMed  Google Scholar 

  55. Liu J, Möller M, Provan J, Gao LM, Poudel RC, Li DZ. Geological and ecological factors drive cryptic speciation of yews in a biodiversity hotspot. New Phytol. 2013;199:1093–108.

    Article  PubMed  Google Scholar 

  56. Wang IJ, Bradburd GS. Isolation by environment. Mol Ecol. 2014;23:5649–62.

    Article  PubMed  Google Scholar 

  57. Tian S, Kou YX, Zhang ZR, Yuan L, Li DZ, López-Pujol J, et al. Phylogeography of Eomecon chionantha in subtropical China: the dual roles of the Nanling Mountains as a glacial refugium and a dispersal corridor. BMC Evol Biol. 2018;18:20.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Fan DM, Sun ZX, Li B, Kou YX, Hodelb RGJ, Jin ZN, et al. Dispersal corridors for plant species in the Poyang Lake Basin of Southeast China identified by integration of phylogeographic and geospatial data. Ecol Evol. 2017;7:5140–8.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Koehler K, Center A, Cavender-Bares J. Evidence for a freezing tolerance-growth rate trade-off in the live oaks (Quercus series Virentes) across the tropical-temperate divide. New Phytol. 2011;38:962–81.

    Google Scholar 

  60. Cao YN, Wang IJ, Chen LY, Ding YQ, Liu LX, Qiu YX. Inferring spatial patterns and drivers of population divergence of Neolitsea sericea (Lauraceae), based on molecular phylogeography and landscape genomics. Mol Phylogenet Evol. 2018;126:162–72.

    Article  PubMed  Google Scholar 

  61. Nosil P, Funk DJ, Ortiz-Barrientos D. Divergent selection and heterogeneous genomic divergence. Mol Ecol. 2009;18:375–402.

    Article  PubMed  Google Scholar 

  62. Charlesworth B, Nordborg M, Charlesworth D. The effects of local selection, balanced polymorphism and background selection on equilibrium patterns of genetic diversity in subdivided populations. Genet Res. 1997;70:155–74.

    Article  CAS  PubMed  Google Scholar 

  63. Nielsen R. Molecular signatures of natural selection. Annu Rev Genet. 2005;39:197–218.

    Article  CAS  PubMed  Google Scholar 

  64. Via S, West J. The genetic mosaic suggests a new role for hitchhiking in ecological speciation. Mol Ecol. 2008;17:4334–45.

    Article  PubMed  Google Scholar 

  65. Saxe H, Cannell MGR, Johnsen B, Ryan MG, Vourlitis G. Tree and forest functioning in response to global warming. New Phytol. 2001;149:369–99.

    Article  CAS  PubMed  Google Scholar 

  66. Zhu H. Geographical elements of seed plants suggest the boundary of the tropical zone in China. Palaeogr Palaeoclimatol Palaeoecol. 2013;386:16–22.

    Article  Google Scholar 

  67. Lin SL, Zhao NX, Chen YZ, Yao JY, Jia XC. Distribution of figs (Ficus) in China and its significance in the issues for interspecific co-evolution. Acta Ecol Sin. 2007;27:4278–88.

    Google Scholar 

  68. Avise JC. Phylogeography: the history and formation of species. Cambridge: Harvard University Press; 2000.

    Book  Google Scholar 

  69. Nosil P, Harmon LJ, Seehausen O. Ecological explanations for (incomplete) speciation. Trends Ecol Evol. 2008;24:145–56.

    Article  Google Scholar 

  70. Wang Z, Fang J, Tang Z, Lin X. Patterns, determinants and models of woody plant diversity in China. Proc R Soc B. 2011;278:2122–32.

    Article  PubMed  Google Scholar 

  71. Huang JH, Chen B, Liu CR, Lai JS, Zhang JL, Ma KP. Identifying hotspots of endemic woody seed plant diversity in China. Divers Distrib. 2012;18:673–88.

    Article  Google Scholar 

  72. Doyle JJ, Doyle JL. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem Bull. 1987;19:11–5.

    Google Scholar 

  73. Kou YX, Zhang L, Fan DM, Cheng SM, Li DZ, Hodel RGJ, et al. Evolutionary history of a relict conifer, Pseudotaxus chienii (Taxaceae), in Southeast China during the late Neogene: old lineage, young populations. Ann Bot. 2020;125:105–17.

    Article  PubMed  Google Scholar 

  74. Jiang L, Bao Q, He W, Fan DM, Cheng SM, Lόpez-Pujol J, et al. Phylogeny and biogeography of Fagus L. (Fagaceae) based on 28 nuclear single/low-copy loci. J Syst Evol. 2020. https://doi.org/10.1111/jse.12695.

  75. Hall TA. BioEdit: a user-friendly biological sequence alignment program for windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41:95–8.

    CAS  Google Scholar 

  76. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28:2731–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.

    Article  CAS  PubMed  Google Scholar 

  78. Watterson GA. On the number of segregating sites in genetical models without recombination. Theor Popul Biol. 1975;7:256–76.

    Article  CAS  PubMed  Google Scholar 

  79. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. 1993;133:693–709.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Hudson RR, Kreitman M, Aguadé M. A test of neutral molecular evolution based on nucleotide data. Genetics. 1987;116:153–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Li HP. A new test for detecting recent positive selection that is free from the confounding impacts of demography. Mol Biol Evol. 2011;28:365–75.

    Article  PubMed  Google Scholar 

  83. Hubisz MJ, Falush D, Stephens M, Pritchard JK. Inferring weak population structure with the assistance of sample group information. Mol Ecol Resour. 2009;9:1322–32.

    Article  PubMed  PubMed Central  Google Scholar 

  84. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  CAS  PubMed  Google Scholar 

  86. Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8.

    Article  Google Scholar 

  87. Hey J. The divergence of chimpanzee species and subspecies as revealed in multipopulation isolation-with-migration analyses. Mol Biol Evol. 2010;27:921–33.

    Article  CAS  PubMed  Google Scholar 

  88. Hey J. Isolation with migration models for more than two populations. Mol Biol Evol. 2010;27:905–20.

    Article  CAS  PubMed  Google Scholar 

  89. Cornuet JM, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, Leblois R, et al. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics. 2014;30:1187–9.

    Article  CAS  PubMed  Google Scholar 

  90. Kuhner MK. LAMARC 2.0: maximum likelihood and Bayesian estimation of population parameters. Bioinformatics. 2006;22:768–70.

    Article  CAS  PubMed  Google Scholar 

  91. Heled J, Drummond AJ. Bayesian inference of species trees from multilocus data. Mol Biol Evol. 2010;27:570–80.

    Article  CAS  PubMed  Google Scholar 

  92. Myers EA, Xue AT, Gehara M, Cox CL, Rabosky ARD, Lemos-Espinal J, et al. Environmental heterogeneity and not vicariant biogeographic barriers generate community-wide population structure in desert-adapted snakes. Mol Ecol. 2019;28:4535–48.

    Article  PubMed  Google Scholar 

  93. Wright S. Isolation by distance. Genetics. 1943;28:114–38.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. McRae BH, Shah VB, Mohapatra TK. Circuitscape 4 User Guide. In: The Nature Conservancy; 2013. http://www.circuitscape.org.

    Google Scholar 

  95. Excoffier L, Laval G, Schneider S, Arlequin ver. 3.1: an integrated software package for population genetics data analysis. Evol Bioinformatics Online. 2005;1:47–50.

    CAS  Google Scholar 

  96. Wang IJ. Examining the full effects of landscape heterogeneity on spatial genetic variation: a multiple matrix regression approach for quantifying geographic and ecological isolation. Evolution. 2013;67:3403–11.

    Article  PubMed  Google Scholar 

  97. Conover W. Practical nonparametric statistics. New York: Wiley; 1999.

    Google Scholar 

  98. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190:231–59.

    Article  Google Scholar 

  99. Schoener TW. The Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology. 1968;49:704–26.

    Article  Google Scholar 

Download references

Acknowledgements

We are grateful to Dr. Bo Li and several previous postgraduates and undergraduates (Weidong Zeng, Tao Li, Bochao Li, Qinglin Meng, and Shupeng Dong) for their participating in the fieldwork.

Funding

This study was supported by grants from the National Natural Science Foundation of China (Nos: 31560064, 31660059, 31160043, 31160082) and Training Program for Academic and Technical Leaders (leading talents) in Major Disciplines of Jiangxi Province. The funding organizations provided financial support to the research projects, but were not involved in the design of the study, data collection, analysis of the data, or the writing of the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

D.M.F. produced data, analyzed data and wrote manuscript; H.L., Q.Y. and S.M.C. conducted experiments; S.Q.L. and Y.Y. collected samples and undertook the formal identification of the samples. Y.X.K. collected samples and analyzed data. Y.X.Q. wrote manuscript. Z.Y.Z. designed research, collected samples and wrote manuscript. All authors approved the final manuscript version.

Corresponding authors

Correspondence to Yixuan Kou or Zhiyong Zhang.

Ethics declarations

Ethics approval and consent to participate

The experiments did not utilize transgenic technology. All plant materials used in this study were obtained from wild plants and did not involve any endangered or protected species. Collection of plant material complied with institutional or national guidelines and was conducted in accordance with local legislation.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Fan, D., Lei, S., Liang, H. et al. More opportunities more species: Pleistocene differentiation and northward expansion of an evergreen broad-leaved tree species Machilus thunbergii (Lauraceae) in Southeast China. BMC Plant Biol 22, 35 (2022). https://doi.org/10.1186/s12870-021-03420-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-021-03420-9

Keywords

  • Population differentiation
  • Species diversity anomaly
  • Isolation-by-environment
  • Freezing temperature
  • Tropical niche conservatism
  • Machilus thunbergii
  • Southeast China