New insight into the phylogeographic pattern of Liriodendron chinense (Magnoliaceae) revealed by chloroplast DNA: east–west lineage split and genetic mixture within western subtropical China

Study species

L. chinense is a long-lived and monoecious tree species that inhabits subtropical forests and is the only Liriodendron species in Asia. Liriodendron trees are easily recognized by their distinctive leaves, with four lobes in most cases and a cross-cut notched or straight apex. Flowering occurs between April and May. L. chinense is pollinated by bees, flies, and beetles (Huang et al., 1999). Winged seeds from the aggregate fruit are dispersed by wind or gravity. The species is predominantly outcrossing, with a limited seed set and germination rate (Huang & Guo, 2002). Liriodendron trees prefer a temperate climate, sun or partial shade, and deep, fertile, well-drained, and slightly acidic soil.

Population sampling

Considering the conservation of the chloroplast genome, we collected samples from four to five individuals from 40 populations across the entire distribution area of L. chinense, which included 29 populations collected in a previous study (Yang et al., 2016) and 11 populations sampled in the present study (Fig. 1; Table 1). In addition, eight individuals from four populations (Missouri, Louisiana, South Carolina, and North Carolina in the USA) of L. tulipifera were sampled as the outgroup. Genomic DNA was extracted from young leaves using the CTAB method (Doyle & Doyle, 1987) and stored at −20 °C.

Genotyping and sequencing

We designed eight noncoding plastid DNA region primers as suggested by Shaw et al. (2007) based on the L. chinense chloroplast genome; amplification was performed using DNA from eight individuals from eight distant populations of L. chinense. Of the primers, psbJ-petA and rpl32-ndhF were found to be the most variable and stable. In addition, we also used three cpDNA primers (psbA-trnH, trnL intron-trnF, and trnK5’-matK) reported by Fetter (2014) for L. tulipifera. Only trnK5’-matK amplification produced polymorphic and stable results; in contrast, the psbA-trnH fragment contained too many SSRs, reducing the accuracy of the sequencing results. Additionally, trnL intron-trnF was too conservative. Hence, psbJ-petA, rpl32-ndhF, and trnK5’-matK were selected for a large-scale survey of haplotype variation in Liriodendron (Table 2). The cpDNA regions were amplified under the protocols reported by Fetter (2014). All sequences were generated using an ABI3730xl DNA sequencer (Applied Biosystems, Foster City, CA, USA), aligned using DNASTAR Lasergene V.7.1 and then proofed and corrected manually. All sequences were deposited in GenBank under accession numbers MH717109–MH717129 (psbJ-petA), MH717130–MH717142 (rpl32-ndhF), and MH717143–MH717150 (trnK5’-matK) (Table 2). The haplotype for each individual was identified by combining three cpDNA fragments.

Population genetic and phylogeographical analysis

The haplotype diversity (h; Nei, 1987) and nucleotide diversity (Pi; Tajima, 1983) were estimated for each population using DnaSP version 5.0 (Librado & Rozas, 2009). The species total genetic diversity (h_T and v_T) and population genetic diversity (h_S and v_S) were calculated by PermutCpSSR 2.0 (Pons & Petit, 1996). Genetic diversity between the scattered eastern mountain regions (E) and continuous western mountain habitats (W), including regions with genetic and habitat heterogeneity (northwest, NW; southwest, SW), were compared using a t-test or Mann–Whitney U test, depending on the distribution of each parameter, in the SPSS 13.0 program (SPSS Inc., Chicago, IL, USA). Nonparametric correlations (Spearman) between genetic diversity (h and π) in the western regions and latitude were also tested.

To quantify the variation in cpDNA sequences among populations and genetic clusters in the eastern and western regions, we performed hierarchical analysis of molecular variance (AMOVA) with 10,000 permutations in Arlequin version 3.1 (Excoffier, Laval & Schneider, 2005). Spatial analysis of molecular variance (SAMOVA) was applied to identify clusters of genetically similar populations. Geographically homogeneous populations were clustered into a user-defined number of groups (K) with a simulated annealing approach to maximize the proportion of total genetic variance observed between groups (F_CT). SAMOVA was conducted with SAMOVA 1.0 (Dupanloup, Schneider & Excoffier, 2002) by using 100 simulated annealing processes for each value of K from K = 2 to K = 20.

The presence of phylogeographic structure was determined by whether the value of R_ST (considering the mutational distances between haplotypes) was significantly higher than that of G_ST (depending on the frequencies of haplotypes), with 1,000 random permutations in PermutCpSSR version 2.0 (Pons & Petit, 1996). Intraspecific genealogies of haplotypes were constructed using a median-joining network in Network 5.0 (Bandelt, Forster & Röhl, 1999), with all indels and single- or multiple-base substitutions treated as one-step mutations and weighted equally. As the SSR regions were too variable, SSRs were excluded from the network and the later maximum parsimony (MP) construction, though we included SSR regions in other analyses.

Two selective neutrality tests, namely, Tajima’s D (Tajima, 1989) and Fu’s F_S (Fu, 1997), were performed to infer potential population growth and expansion in Arlequin v3.1 (Excoffier, Laval & Schneider, 2005), and significantly negative values indicated population expansion (Fu, 1997). Pairwise mismatch distributions were analyzed to infer the historical demography of L. chinense and each (sub)clade. We calculated the expected frequency based on a population growth-decline model. The sum of squared deviations (SSDs) between the observed and expected mismatch distributions and the raggedness index (H_Rag; Harpending, 1994) were calculated with 1,000 parametric bootstrap replicates to test the null hypothesis of spatial expansion using Arlequin ver. 3.1. Furthermore, Bayesian Skyline Plots (BSPs) were used to test the hypothesis of demographic expansion by assessing the time variation in effective population size. The BSPs were constructed in BEAST ver. 2.1.3 (Bouckaert et al., 2014), and settings were referred to Amarilla et al. (2015). The substitutions rate was set as 1.0–3.0 × 10⁻⁹ per site per year (s/s/y). Two independent runs were carried out for a chain length of 1 × 10⁸ with a burn-in of 10%. The results of each run were visualized using Tracer v 1.7 (Rambaut et al., 2018) to ensure that stationarity and convergence had been reached (ESS > 200). All the demographic expansion tests were also performed for each of the three genetically heterogeneous regions.

Divergence between cpDNA lineages and ancestral area reconstructions

Phylogenetic relationships among the cpDNA haplotypes were reconstructed using the MP method with L. tulipifera as the outgroup. We conducted MP analyses in Mega 6 (Tamura et al., 2007) and used heuristic searches with the random addition of sequences (1,000 replicates) and tree bisection–reconnection branch swapping selected; gaps were treated as complete deletions. We calculated bootstrap values using 1,000 replicates of the original dataset, and five random additional sequences were added to each replicate.

A molecular clock model was applied to estimate the cpDNA lineage divergence time. The rate constancy of cpDNA haplotype evolution in L. chinense was first evaluated by Tajima’s relative rate tests (Tajima, 1993) in MEGA 6 (Tamura et al., 2007) using L. tulipifera as the outgroup. Then, the time since divergence (T) of the haplotype lineages was inferred from their net pairwise sequence divergence per base pair (d_A) using the Kimura 2-parameter model. T was calculated as T = d_A/2μ, where μ is the rate of nucleotide substitutions (Nei & Kumar, 2000); μ is derived from the pairwise genetic distance and an assumed T between L. chinense and L. tulipifera of 14.15 Mya (Nie et al., 2008).

To reconstruct the geographical ancestral area of L. chinense, a Bayesian binary Markov chain Monte Carlo (MCMC) (BBM) analysis was performed in RASP v3.0 (Yu et al., 2015; http://mnh.scu.edu.cn/soft/blog/RASP) using trees retained from Mega 6 (see above). L. tulipifera was used as an outgroup for our ancestral state reconstructions. Three geographic regions were defined according to the distribution pattern of L. chinense (Hao et al., 1995) and adjusted based on the presumed refugia, i.e., EA, East; SW, southwest; and NW, Northwest. The number of maximum areas at each node was set to four. One thousand trees produced by Mega 6 were used for BBM analyses. We set the root distribution to null, applied 10 MCMC chains with the JC + G model run for 106 generations, and sampled the posterior distribution every 100 generations.

Haplotype distribution and genetic diversity

For the cpDNA dataset, the psbJ-petA, rpl32-ndhF, and trnK5’-matK sequences, consisting of a consensus length of 2,080 bp with 23 nucleotide substitutions, were combined and aligned. In addition, one multiple-base substitution (12 bp in length), three indels (8–9 bp in length), and four SSRs were also identified (Supplementary data in Table S1). A total of 43 haplotypes were identified from these 32 polymorphisms (h1–h43, Table S1). Thirty-seven haplotypes were found in L. chinense, and six haplotypes were identified in L. tulipifera. No haplotype was shared by the two sister species. The geographical distribution of the haplotypes and their frequencies within populations are illustrated in Table 1.

CpDNA haplotypes were distributed relatively evenly over the entire range of L. chinense. No haplotype was obviously dominant, and 22 of the 37 haplotypes were unique to a single population (Table 1). In the western region, h38 was the most frequent haplotype in the northwest, and h30 was found only in the southwestern range. Each population in the Yungui Plateau contained only one haplotype, except for YJ, which was located in the northwestern Yungui Plateau and harbored three haplotypes. Population HF in the Wu Mountains and YY in the Wuling Mountains also harbored three haplotypes. In the eastern region, only h4 occurred in more than two populations. Most populations contained one to two haplotypes, except for LA, located in the Tianmu Mountains, which contained four haplotypes, largely due to the high degree of variation in the SSR regions (Table 1). At the species level, three cpDNA regions exhibited high haplotype diversity (v_T = 0.972, h_T = 0.970) and nucleotide diversity (π_T = 1.756 × 10⁻³), though the average population haplotype diversity was relatively low (h_S = 0.300, v_S = 0.195). The unbiased haplotype diversity (h) for each population was zero to one (LA), and the nucleotide diversity ranged from 0 to 0.128 × 10⁻³ (ZY and K). No significant deficiency in genetic diversity was detected in the scattered eastern mountain regions (P > 0.05, Table 1). In the continuous western mountain regions, genetic diversity was not significantly correlated with latitude (h: Spearman r = 0.340, P = 0.076; π: Spearman r = 0.318, P = 0.099), but hotspots of haplotype diversity occurred at mid-latitudes close to the current northern edge, around the Wu Mountains, Dalou Mountains, and Wuling Mountains. The number of haplotypes was clearly low in the southwestern range (Fig. S1).

Population genetic structure of L. chinense

AMOVA revealed strong population genetic structure at the species level (F_ST = 0.747, P < 0.001). Hierarchical AMOVA revealed that 17.74% of the variation between the two groups and 56.97% of the total variation in cpDNA of the species was distributed within the eastern and western regions (Table 3). By exploring the behavior of the indices F_CT (minimized) and F_SC (maximized) for different values of K, the best population grouping scheme for a set of sample populations can be identified in SAMOVA (Dupanloup, Schneider & Excoffier, 2002). In the present study, SAMOVA revealed a discontinuous increase in F_CT and a discontinuous decrease in F_SC when K increased from 2 to 20 (Fig. S2). However, F_CT reached its first peak and F_SC was relatively low when K = 7, and the trends in F_CT and F_SC reversed when K = 8. Although K = 11 was still acceptable, two groups were subdivided; thus, we used K = 7 as the best grouping scheme. The composition of groups corresponded well with the geographical distribution (Fig. 2). The first group had a geographically wide distribution and was composed of all the populations in the southwest. Specifically, the unique ‘island’ in the western range (Hao et al., 1995) of the Nanling Mountains was genetically consistent with populations in the continuous western region (the “belt”). Populations from the Daba Mountains composed group 2, and the three populations located in the Wu Mountains formed group 3. The populations in groups 2 and 3 exclusively formed the northwestern clade. Populations in the eastern region were obviously divided into four groups (4–7), with group 6 containing populations from the Tianmu Mountains in the northeastern region and group 7 including populations in the eastern Wuyi Mountains in the southeastern region. The populations located in the northern and central regions (Luoxiao Mountains) were divided into two groups (Fig. 2).

Table 3:

Results of AMOVA test of the population genetic structure in Liriodendron chinense.

Source of variation	d.f.	Percentage of variation	P	FST
Among groups	1	17.74	< 0.001
Among populations within groups	38	56.97	< 0.001
Within populations	124	25.29	< 0.05	0.747*

DOI: 10.7717/peerj.6355/table-3

Note:

* P < 0.001, 1,000 permutations.

Phylogeographic pattern of L. chinense

N_ST was significantly higher than G_ST (N_ST = 0.799, G_ST = 0.691; P < 0.001), indicating significant phylogeographic structure across the range of L. chinense. Four main clades (A–D) were revealed by the three cpDNA regions in the median-joining network (Fig. 3). Haplotypes from L. tulipifera formed clade D. Clade A was composed of all the haplotypes from the eastern and central populations in the eastern region. Haplotypes in clades B and C were specific to the western region. All populations in the southwestern region harbored haplotypes from clade B, whereas all populations from the northwestern region contained haplotypes in clade C; populations in the midwestern region usually contained haplotypes from both the B and C clades (Fig. 1; Table 1). Notably, clade A connected with clade D and showed a close relationship with L. tulipifera, and clade C appeared as a tip clade. Among these haplotypes, H2 in clade A was the fewest mutational steps from L. tulipifera, and this haplotype was widely distributed in the northern part of the eastern region, including the Tianmu Mountains, Dabie Mountains, and Lu Mountains.

The MP tree showed a pattern similar to that in the haplotype network (Fig. 4). Obvious eastern and western clades were found, and the western clade was further divided into two subclades: the southwestern and northwestern subclades. Relative rate constancy tests suggested equal rates among the cpDNA haplotypes of L. chinense and L. tulipifera (all P > 0.05), indicating that the clock-like evolution model (Tajima, 1993) was applicable for L. chinense. Based on a divergence time between the Liriodendron sister species of 14.15 Ma, a divergence time between the eastern and western regions of L. chinense of 0.433 Ma (95% CI = 0.0254–0.841 Ma) was suggested.

Demographic analysis and ancestral area reconstructions

None of the estimates of Tajima’s D were significant for any cpDNA (sub)clade of L. chinense. Although all Fu’s F_S values were positive, only the Fu’s F_S over the entire range (P < 0.01) was significant. The observed mismatch distributions of haplotypes for each (sub)clade failed to reject the spatial expansion model in most cases (SSD, H_Rag values, P > 0.05; Table 4), except for the SW lineage, with a significant SSD and H_Rag statistic (P < 0.05). The results of BSPs showing the estimated changes in effective population size over time for the whole and each (sub)clade were concordant with the results from the neutrality tests. They showed that L. chinense had a constant population size during the last 15,000 years (Figs. S3A–D).

Based on the topology of the interspecific chronogram, the BBM analysis of ancestral distribution areas (Fig. 4) supported a likely ancient vicariance event (node I, NA, and EA) of two sister species. A subsequent vicariance event (node II) between populations in eastern (EA) and southwestern (SW) China was found. This event was then followed by dispersal from the east (EA) to the southwest (SW). No obvious divergence event was detected between NW and SW; instead, several dispersal events were evident. Recent colonization events from the southwest (SW) to the northwest (NW) and the reverse (NW to SW) were inferred, based on the genetic admixture (H6 and H17) in the adjacent mountains (Figs. 1 and 4), with multiple expansions at different times.

East–west lineage split and the role of the Nanling Mountains

By investigating cpDNA sequence variation, an obvious west–east lineage split across subtropical China was found for L. chinense, and this L. chinense west–east genetic divergence was further shown by the phenotypic divergence between the eastern and western regions, such as divergence in petal color and lobed leaf number (Hao et al., 1995). Furthermore, according to the results of the molecular phylogenetic analysis based on RAD-seq, the eastern and western populations of L. chinense formed two clades (Zhong et al., 2019). These results, coupled with those from some other species (Shi et al., 2014; Sun et al., 2014a; Wang et al., 2015; Feng et al., 2017), emphasize the role of allopatric genetic divergence in shaping the phylogeographic patterns of plant species in subtropical China, providing new insight into the dual role of the Nanling Mountains in facilitating a west–east lineage split of temperate plant species in subtropical China.

Since the late Tertiary climate oscillation, repeated allopatry has caused frequent speciation (Hewitt, 2004). The breaking off of Beringia during the late Tertiary global cooling event induced species splits between Asia and North America (Milne & Abbott, 2002); L. chinense and L. tulipifera were among the species that diverged (Parks & Wendel, 1990). To the best of our knowledge, phylogeographic studies in subtropical China have reinforced the allopatric speciation hypothesis for species diversity, and the phylogeographic divergence in subtropical China is mostly attributed to the tectonic and climate changes during the Quaternary (Qiu et al., 2017). Specifically, the uplift of the Himalayan–Tibetan Plateau and the subsequent formation of Asian monsoons generated a split of the southwestern clade (southwestern Yungui Plateau) from the clades in other regions for many plant species, such as Tetrastigma hemsleyanum (Southwest China vs. Central and East China: 5.07 Ma; Wang et al., 2015), Taxus wallichiana (Sino-Himalayan vs. others: 6.58 Ma, Tonkin Bay region vs. Sino-Japanese forest subkingdom: 3.84 Ma; Gao et al., 2007), and Euptelea pleiosperma (West China vs. Central/East China: 1.66 Ma; Cao et al., 2016). Nonetheless, although L. chinense was also spread across different geological formations and two climate types in the Yungui Plateau, no lineage differences were found (Figs. 1 and 2). In contrast, the observed west–east divergence in L. chinense occurred in the eastern Yungui Plateau and Sino-Japanese forest subkingdom region (Wu & Wu, 1998), which have relatively consistent geological and climatic conditions. Compared with the long-term isolation and ancient allopatric divergence of Quercus spinosa (Southwest China vs. Central and East China: 25.67 Ma, Feng et al., 2017) in subtropical China, the probable divergence time of L. chinense between the east and west was estimated to have been approximately 0.433 Ma, which closely coincides with the early Mindel glaciation (approx. 0.30–0.45 Ma). Although recently occurring and not as pronounced, a west–east lineage split has also been detected in many other species, such as Fagus engleriana (Lei et al., 2012), Castanopsis eyrei (Shi et al., 2014), and Castanopsis fargesii (Sun et al., 2014a). We argue that the west–east lineage split was likely a process driven by gradual geographic isolation and allopatric genetic divergence. However, the process that underlies the divergence that occurs along a geographic feature, such as between the western (eastern Yungui Plateau) and eastern regions, remains unclear, though it is likely the result of the dual corridor and barrier roles played by the Nanling Mountains.

The Nanling Mountains, stretching from the west to the east for more than 1,000 km in subtropical China, have been recognized as an important corridor for the eastward migration of East Asiatic flora, referred to as the Nanling corridor (Wang, 1992). Recent phylogeographic studies have shown that the Nanling Mountains acted as a stepping stone and dispersal corridor for species, facilitating the exchange of genetic resources between the eastern Yungui Plateau and the eastern mountain regions in subtropical China during the Late Quaternary (e.g., Gao et al., 2007; Wang et al., 2009; Tian et al., 2015, 2018). In fact, according to the particular location of the Nanling Mountains and the westward dispersal event from the eastern regions inferred from the cpDNA-based BBM analysis, we assumed that the Nanling Mountains may have also contributed to the ancient dispersal of L. chinense to the eastern Yungui Plateau.

However, the Nanling corridor was frequently inhabited by coniferous or evergreen forest during the Quaternary glacial and interglacial periods (Harrison et al., 2001; Qiu, Fu & Comes, 2011). In the case of L. chinense, although there is a wide confidence range, the divergence time of the west–east lineage split (0.0254–0.841 Ma) occurred before the LGM (~0.018 Ma), and the Nanling Mountains were climatically unsuitable during the LGM (Yang et al., 2016). An inescapable consequence of the interruption of gene flow, the west–east lineage split may have developed due to gradual geographic isolation. In view of biological differences in the response to climate change among species, the Nanling Mountains may have functioned as a “filter” (such as the East China Sea land bridge during the Quaternary, Qiu et al., 2017) for different species. Hence, the west–east lineage split only partially occurred in some species, such as F. engleriana (Lei et al., 2012), C. eyrei (Shi et al., 2014), and C. fargesii (Sun et al., 2014a), but the split is not currently observed in species that are widely distributed in the Nanling Mountains and are able to produce sufficient gene flow between the eastern mountains and the Yungui Plateau (i.e., Taxus wallichiana, Gao et al., 2007; Sargentodoxa cuneata, Tian et al., 2015; and Tetrastigma hemsleyanum, Wang et al., 2015).

“Melting pot” due to admixture in the continuous western mountain region

Coexistence of the northern refugium in the eastern Daba/Wu Mountain regions and the southern refugium in the Yungui Plateau has also been found in many other subtropical species (i.e., Taxus wallichiana, Gao et al., 2007; Ligularia hodgsonii, Wang et al., 2013; Tetrastigma hemsleyanum, Wang et al., 2015; Euptelea, Cao et al., 2016; and Platycarya, Wan et al., 2017). Although not stressed by authors, genetic admixture of southern and northern refugia is frequently found, and the Wuling Mountains, which are located in the contact region, usually were the center of genetic admixture of many plant species, such as Taxus wallichiana (Gao et al., 2007), Tetrastigma hemsleyanum (Wang et al., 2015), and Emmenopterys henryi (Zhang et al., 2016).

For L. chinense, the Wuling Mountains and the Dalou Mountains are the present distribution center (Hao et al., 1995) and genetic diversity center (Fig. 1; Fig. S1), respectively. The high genetic diversity of this species may be explained by a center of ancient origin or genetic admixture from the northern and southern refugia. Coalescent theory predicts that older alleles will occupy the interior nodes of a haplotype network (Posada & Crandall, 2001). Nevertheless, none of the L. chinense populations there harbored ancient haplotypes from both clades. More precisely, the populations in the Dalou Mountains (NC, DZ, and TZ) contained the ancestral haplotype H17 from clade C (northwestern clade) and tip haplotypes (H3, H5, and H16) from clade B (southwestern clade). In contrast, populations in the Wuling Mountains (YY, K, and NS) all exclusively contained ancestral haplotype H13 from clade B and one tip haplotype, H6, from clade C (Figs. 1 and 3). Therefore, there were no areas where diversity was created. In fact, species usually shifted northward or southward to track warm interglacial or cold glacial conditions during the Quaternary climate oscillation (Hewitt, 2000). Increased diversity was more likely to occur in zones of admixture produced by redistribution of the genetic diversity created in refugia (Petit, 2003). For example, in E. henryi, genetic admixture of the northern lineage with the southern lineage is more likely explained by both introgressions following secondary contact and incomplete lineage sorting due to recent postglacial divergence during glacial periods (Zhang et al., 2016). Collectively, these findings suggest that the genetic diversity centers of L. chinense in the Wuling Mountains and Dalou Mountains were not the primary locations of ancient refugia but more likely a “melting pot” due to the admixture of divergent lineages from northern or southern refugia.

The next question, then, is what the underlying process is. In contrast to wave-like migration fronts, expansion from local refugia followed by subsequent population admixture is more likely to be characteristic of postglacial forest development (Hu, Hampe & Petit, 2009). In western subtropical China, the Dalou Mountains and the Wuling Mountains are located in the northern and northeastern portions of the Yungui Plateau, respectively, and they are geographically connected with the Yungui Plateau. We argue that these mountains can receive haplotypes from the eastern Yungui Plateau refugium by occasional capture because the expansion of L. chinense in the Yungui Plateau refugium appears to have been imperceptible. In view of the facts that clade C did not show a "star-like" pattern of population expansion and that the mismatch test and neutral test rejected an expansion pattern in the Yungui Plateau (Table 4), gene flow from the Yungui Plateau to adjacent mountains was weak and repeated several times, showing no obvious signals of sudden expansion. For example, despite the high cpDNA diversity in Populus adenopoda (Fan et al., 2018), the Yungui Plateau may have contributed little to its recolonization of northern and eastern areas. The most likely scenario for L. chinense was that the southern refugia resources had a pan-distribution pattern in the Yungui Plateau and the species reached its northern range limit in the Dalou Mountains (NC, DZ, and TZ) and Wuling Mountains (YY, K, and NS), where it encountered northwestern genetic resources (TZ, DZ, YY, K, and NS).

Glacial refugial populations at mid-to-high latitudes sometimes appear to play a greater role in the postglacial development of forest biomes than those inferred from fossil pollen records (Hu, Hampe & Petit, 2009). Similarly, in L. chinense, the southward introgression of northern refugial resources greatly contributed to the genetic admixture in the Dalou Mountains and Wuling Mountains. The northern refugium was centered in the Wu Mountains, where all the haplotypes in clade C were located (Table 1; Fig. 1). In contrast, in the Dalou Mountains (DZ and TZ), the haplotypes were almost exclusively H17, and the Wuling Mountains (YY, K, and NS) harbored only tip haplotypes from clade C. Although not significant, demographic analysis partially suggested weak expansion in the northwestern subclade (Table 4). Therefore, the Wu Mountains may have contributed to the southward colonization of L. chinense to the adjacent Dalou Mountains and Wuling Mountains (Figs. 1 and 2). This pattern was also found in other species, for example, in Tetrastigma hemsleyanum (Wang et al., 2015). Haplotypes in the Dalou Mountains (CC4 and CC5) and northern Wuling Mountains (CC6 and CC8) were also derived from the Wu Mountains (CC7, haplotype 9). Furthermore, genetic resources from the Wu Mountain region may have also contributed to genetic admixture in the Nanling Mountains (ZY, Fig. 1) or even farther east in the Luoxiao Mountains (Tetrastigma hemsleyanum, Wang et al., 2015). Given these results, we argue that genetic diversity in the Dalou Mountains and Wuling Mountains is largely due to the frequent contact of southern and northern refugium resources in response to past glacial and interglacial cycles.

Mountain refugia and implications for conservation

Evidence supporting the presence of multiple geographically isolated refugia of plant species in subtropical China without an obvious range shift to the south during the LGM has accumulated (reviewed by Qiu, Fu & Comes, 2011 and Liu et al., 2012), such as that found in E. cavaleriei (Wang et al., 2009), F. engleriana (Lei et al., 2012), Platycarya (Wan et al., 2017) and P. adenopoda (Fan et al., 2018). In the case of L. chinense, a “refugia within refugia” scenario (c.f. Gómez & Lunt, 2007) of multiple localized glacial refugia across the range has been suggested (Yang et al., 2016). Similarly, in the present study, the results of three cpDNA sequence variations indicated the presence of five mountain-based refugia across the range of L. chinense. These refugia were located in the Wu Mountains (northwest), Yungui Plateau (southwest), Luoxiao Mountains (center), Tianmu Mountains (northeast) and Wuyi Mountains (southeast), falling almost entirely within the main mountain ranges (Figs. 1 and 2).

These mountain regions reportedly act as important refugia for many plant species. For instance, the Wu Mountains were a large endemism center (López-Pujol et al., 2011) and functioned as an important refugium for relict tree species (e.g., E. cavaleriei, Wang et al., 2009; Tetracentron sinense, Sun et al., 2014b; and E. henryi, Zhang et al., 2016). The Tianmu Mountains formed a northern refugium counterpart to the Wu Mountains, especially for many relict species (e.g., Ginkgo biloba, Gong et al., 2008; and F. engleriana, Lei et al., 2012). The Wuyi Mountains were previously reported to serve as an important refugium (e.g., Wang et al., 2009; Sun et al., 2014a; Wan et al., 2017). We also found divergent haplotypes in the Wuyi Mountains regions, and special attention should be paid to the haplotype divergence of L. chinense between the eastern and western ranges of the Wuyi Mountains (Figs. 1 and 3), which has also been found in E. cavaleriei (Wang et al., 2009) and E. henryi (Zhang et al., 2016).

From the perspective of conservation biology, regions possessing different haplotypes should be targeted to maintain the maximum amount of genetic diversity. Hence, the main mountain refugia located in the scattered eastern and continuous western mountain regions should be included. Furthermore, in view of their high degree of harbored genetic diversity, the Daba and Wuling Mountains need special protection. Although L. chinense is listed as an endangered species in China (Fu & Jin, 1992), fortunately, most populations are currently protected in nature reserves or areas used as fengshui forests by local villagers (Tang et al., 2013). Due to the high genetic heterogeneity of L. chinense, in situ conservation priorities should target populations from the five main mountain refugial regions. For ex situ conservation, reintroduction, seed banking, and collection of germplasm resources should be implemented to reserve the unique germplasm resources of L. chinense.