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

Background Subtropical China is a global center of biodiversity and one of the most important refugia worldwide. Mountains play an important role in conserving the genetic resources of species. Liriodendron chinense is a Tertiary relict tree largely endemic to subtropical China. In this study, we aimed to achieve a better understanding of the phylogeographical pattern of L. chinense and to explore the role of mountains in the conservation of L. chinense genetic resources. Methods Three chloroplast regions (psbJ-petA, rpl32-ndhF, and trnK5’-matK) were sequenced in 40 populations of L. chinense for phylogeographical analyses. Relationships among chloroplast DNA (cpDNA) haplotypes were determined using median-joining networks, and genetic structure was examined by spatial analysis of molecular variance (SAMOVA). The ancestral area of the species was reconstructed using the Bayesian binary Markov Chain Monte Carlo (BBM) method according to its geographic distribution and a maximum parsimony (MP) tree based on Bayesian methods. Results Obvious phylogeographic structure was found in L. chinense. SAMOVA revealed seven groups matching the major landscape features of the L. chinense distribution area. The haplotype network showed three clades distributed in the eastern, southwestern, and northwestern regions. Separate northern and southern refugia were found in the Wu Mountains and Yungui Plateau, with genetic admixture in the Dalou Mountains and Wuling Mountains. BBM revealed a more ancient origin of L. chinense in the eastern region, with a west–east split most likely having occurred during the Mindel glacial stage. Discussion The clear geographical distributions of haplotypes suggested multiple mountainous refugia of L. chinense. The east–west lineage split was most likely a process of gradual genetic isolation and allopatric lineage divergence when the Nanling corridor was frequently occupied by evergreen or coniferous forest during Late Quaternary oscillations. Hotspots of haplotype diversity in the Dalou Mountains and Wuling Mountains likely benefited from gene flow from the Wu Mountains and Yungui Plateau. Collectively, these results indicate that mountain regions should be the main units for conserving and collecting genetic resources of L. chinense and other similar species in subtropical China.

INTRODUCTION for many plant species in subtropical China has accumulated (Qiu, Fu & Comes, 2011;Liu et al., 2012). Furthermore, there are no clear patterns of recolonization from the south but rather a pattern of range fragmentation (e.g., Wang et al., 2009, Wang et al., 2013, Feng et al., 2017 or only localized range expansion (e.g., Qi et al., 2012;Tian et al., 2015). Phylogeographic studies frequently suggest that the mountains in subtropical China functioned as refugia (i.e., Qiu, Fu & Comes, 2011;Liu et al., 2012;Qiu et al., 2017), dispersal corridors (e.g., Tian et al., 2015Tian et al., , 2018, or even genetic barriers to migration (Sun et al., 2014b;Zhang et al., 2016). As lineage differentiation and colonization mostly occurred far earlier than the last glacial maximum (LGM) (e.g., Cao et al., 2016;Zhang et al., 2016;Qiu et al., 2017), more ancient events should be considered in phylogeographic studies of subtropical organisms in China (Zhang et al., 2016). Moreover, the presence of northern refugia may have blurred a typical colonization pattern (Stewart et al., 2010), and individual biological characteristics and environmental differences may result in various phylogeographic patterns (Davis & Shaw, 2001;Petit, 2003). Hence, although the number of phylogeographic studies is rapidly increasing, given the complex topography and highly diverse flora in subtropical China, phylogeographic studies in this region remain preliminary.
The genus Liriodendron (Magnoliaceae) was widely distributed in the northern region of the Northern Hemisphere during the early to mid-Miocene but became extinct in Europe during the Pleistocene (Parks & Wendel, 1990), leaving the typical East Asian and eastern North American disjunct sister species Liriodendron chinense (Hemsl.) Sarg. and Liriodendron tulipifera Linn. Here, we focus on L. chinense, which naturally occurs in subtropical China and North Vietnam (Hao et al., 1995). In earlier research (Yang et al., 2016), we investigated the population genetics of central and peripheral populations (rear edge and leading edge) of L. chinense by using simple sequence repeat (SSR) markers from both the chloroplast and nuclear genomes to explore the roles of historical and present population genetics in determining the range limits of this species. However, the homology of chloroplast SSRs (cpSSRs) limited our ability to reveal the phylogeographic pattern of L. chinense, and a more complex phylogeographic pattern was found in the western region of L. chinense. Thus, in the present study, we applied three chloroplast DNA fragments (psbJ-petA, rpl32-ndhF, and trnK5'-matK) and increased sampling to 40 populations throughout the distribution area of L. chinense, especially surrounding the Sichuan Basin in the western regions, aiming to draw a clearer phylogeographic pattern of L. chinense and to examine the role of mountains in gene conservation in subtropical China. Here, we address the following questions: (1) Is there a clear relationship between the pattern of genetic divergence and the pattern of geographic distribution? (2) Is there an obvious lineage split or obvious genetic admixture between the continuous western mountain region and scattered areas in the eastern mountain region? (3) Are there any hotspots of haplotype diversity? If such regions do exist, what is the likely underlying process creating them? The results of this study should be valuable in furthering our understanding of the present genetic distribution of L. chinense and in helping to formulate appropriate strategies for the conservation of L. chinense and other plant species with similar distribution patterns in subtropical China.  Table 1 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.

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 −9 per site per year (s/s/y). Two independent runs were carried out for a chain length of 1 × 10 8 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 (p T = 1.756 × 10 −3 ), 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 −3 (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; p: 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  were divided into two groups (Fig. 2).
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 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.

DISCUSSION
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 westeast 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  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 ( 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., 2015Tian et al., , 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 Comes, 2011 andLiu 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.

CONCLUSION
A west-east lineage split across subtropical China was found for L. chinense. This pattern is different from the clear lineage divergence in the southwestern (southwestern Yugui Plateau)-central China split, which was mainly caused by the ancient uplift of the Himalayan-Tibetan Plateau and the formation of Asian monsoons (e.g., Wang et al., 2015;Cao et al., 2016;Zhang et al., 2016). The west-east split examined in the present study is more likely a process of gradual genetic isolation and allopatric lineage divergence, especially in some environmentally sensitive plant species, when the Nanling corridor was frequently inhabited by coniferous or evergreen forest during the Quaternary climate oscillation (Harrison et al., 2001;Qiu, Fu & Comes, 2011).
In the continuous western mountain region, genetic resources from the southern refugium and northern refugium combined in the Dalou Mountains and the Wuling Mountains. The contact regions became the regions of high genetic diversity and modern distribution centers of L. chinense. In the scattered eastern mountain regions, genetic differentiation likely occurred as the habitat gradually fragmented in each mountain region. This large-scale phylogeographical investigation demonstrates that L. chinense was structured into several distinct groups. Each group showed a marked mountain-associated pattern, falling almost entirely within the main mountain ranges, many of which were hypothesized to be areas of former glacial refugia (reviewed by Liu et al., 2012 andQiu et al., 2017) unique to the Wu Mountains, Yungui Plateau, Tianmu Mountains, Luoxiao Mountains, and Wuyi Mountains.
Collectively, these results indicate that mountain regions should be the main genetic resource conservation and collection units for L. chinense. Our results for L. chinense contribute to a growing body of evidence regarding the role of mountains in forming refuges for subtropical species and the different phylogeographic patterns in separate mountain regions. To obtain a better understanding of the phylogeographic pattern in subtropical China, we should emphasize studies on species with different biological characteristics in subtropical China.