Genetic population structure of the alpine species Rhododendron pseudochrysanthum sensu lato (Ericaceae) inferred from chloroplast and nuclear DNA

Background A complex of incipient species with different degrees of morphological or ecological differentiation provides an ideal model for studying species divergence. We examined the phylogeography and the evolutionary history of the Rhododendron pseudochrysanthum s. l. Results Systematic inconsistency was detected between gene genealogies of the cpDNA and nrDNA. Rooted at R. hyperythrum and R. formosana, both trees lacked reciprocal monophyly for all members of the complex. For R. pseudochrysanthum s.l., the spatial distribution of the cpDNA had a noteworthy pattern showing high genetic differentiation (FST = 0.56-0.72) between populations in the Yushan Mountain Range and populations of the other mountain ranges. Conclusion Both incomplete lineage sorting and interspecific hybridization/introgression may have contributed to the lack of monophyly among R. hyperythrum, R. formosana and R. pseudochrysanthum s.l. Independent colonizations, plus low capabilities of seed dispersal in current environments, may have resulted in the genetic differentiation between populations of different mountain ranges. At the population level, the populations of Central, and Sheishan Mountains may have undergone postglacial demographic expansion, while populations of the Yushan Mountain Range are likely to have remained stable ever since the colonization. In contrast, the single population of the Alishan Mountain Range with a fixed cpDNA haplotype may have experienced bottleneck/founder's events.


Background
The analysis of divergence patterns between species has been a major focus in systematics and molecular ecology [1][2][3][4][5]. Species descending from a common ancestor are expected to differentiate from each other and eventually attain reciprocal monophyly. Coalescence theory predicts that sister species are likely to be polyphyletic/paraphyletic in their genetic composition during the early stages after they split from their common ancestor [cf. [6]].
As sufficient time goes by, they become reciprocally monophyletic with the loss of complementary haplotypes. In contrast, lack of monophyly can be a result of several factors, such as incomplete lineage sorting, gene duplication, or hybridization [7][8][9][10][11][12]. Of them, hybridization leading to speciation has been known as one of the major forces for the diversification in plants [13]. As many hybridizing species remain morphologically discrete [14], hybrids usually share traits with parental species, often resulting in taxonomic difficulties. Interspecific gene flow usually leads to a mixed genetic composition in hybrids, but may have low impacts on parental species' genomes [cf. [15]]. Although molecular tools provide sufficient power in detecting genetic variation, empirical data on species delimitation are often difficult to interpret, since both gene introgression after speciation and shared ancestral polymorphisms can result in species paraphyly [15]. A complex of incipient species with different degrees of morphological or ecological differentiation provides an ideal model for studying species divergence.
Taiwan, an island of the island-arc system along the western edge of the Pacific Ocean, first emerged from the waters via collision between Eurasian and Pacific plates about 9 million years before present (Proto-Taiwan Stage), but attained its modern shape only 5-6 million years ago (MYA) [16,17]. Its rugged topography is characterized by hundreds of steep mountains, which in turn provide diverse habitats along distinct vegetation zones, including tropical, coastal evergreen forests, subalpine shrubs and alpine tundra along the Central Mountain Range [18]. Taiwan was connected to mainland SE Asia by a land bridge during lower sea levels of the last glacial maximum [19,20]; some species may have reached Taiwan by this route and subsequently become endemic [2,21].
Phylogeographical patterns and genetic characteristics of populations/species are shaped by the interaction between historical vicariance and recurrent genetic exchanges. These evolutionary events would leave evolutionary footprints in the genetic polymorphisms within and between populations across the distributional range. Of the geological events, regular historical glacial cycles in the Eurasia Continent had prevalent influences on survival and recolonization of populations/species [11]. Postglacial expansion of alpine species may follow the phalanx model, which describes the effects of slower expansions from refugia due to habitat constraints [2,22,23]. That is, during interglacial periods, alpine species moved upward to peaks and became fragmented. In contrast, down-slope movement of alpine species during glacial periods could connect populations that were previously isolated and allow for genetic exchange among populations [24]. On Taiwan island, only the peaks of high mountains over 3000 m in elevation were covered by ice sheets during the glacial maxima.
Rhododendron, one of the largest and most widespread woody plant genera, is distributed from the northern temperate zone, throughout tropical Southeast Asia, to northeastern Australia and consists of over 1000 species [25]. Hymenanthes, one of eight subgenera comprising Rhododendron, consists of 225 species, most of them distributed in the Himalaya-Southwest China region [25][26][27][28]. In Taiwan, 13 taxa occur within Rhododendron [29]. The R. pseudochrysanthum complex consists of R. pseudochrysanthum Hayata sensu stricto (s.s.), R. morii Hayata, and R. rubropunctatum Hayata, all in subgenus Hymenanthes subsection Maculifera, and all are phylogenetically related based on molecular analyses [30][31][32]. The remaining taxa of subsection Maculifera occur in west and central China, disjunctly from those in Taiwan [27]. R. pseudochrysanthum s.s. is an alpine species (3000-3900 m a.s.l.), while R. morii has a wider range from 1650-3600 m a.s.l., with some overlap in distribution. In contrast, R. rubropunctatum is restricted to elevations of 600-1200 m a.s.l. in northern Taiwan. The three species of R. pseudochrysanthum complex, all distributed along the central mountain range (CMR) and other mountain ranges west of the CMR, share morphological polymorphisms in leaf shape, petiole and flower shape, thus making taxonomy difficult, despite the differentiation of ecological distributions [33]. Possibly the closest relative of this species complex is R. hyperythrum Hayata, which is also endemic to Taiwan, occurring 1400-3700 m a.s.l. [25,28,29]. This species was previously placed in subsection Pontica [25,34], but phylogenetic work indicates that it is closely related to R. pseudochrysanthum [28]. These two might be sister species, but at present a phylogeny of subsection Maculifera is lacking, so this hypothesis cannot yet be tested. Rhododendron hyperythrum is also known to hybridise with members of the R. pseudochrysanthum complex [33]. Closely related to R. hyperythrum, the three species of the R. pseudochrysanthum complex with low morphological differentiations may be recognized as a single species (R. pseudochrysanthum sensu lato [s.l.]) [25,28,29,35].
Previous studies have investigated the phylogeography of the R. pseudochrysanthum in Taiwan, revealing close relationships among these taxa [30,32,51] and a lack of geographical subdivision. However, little attention has been paid to the demographic dynamics of R. pseudochrysanthum across the mountain ranges. In contrast to previous studies, sampling in this study covered vaster areas, with the aid of the bureau of the Yushan National Park. The present study analyzed genetic variation of the R. pseudochrysanthum s.l. at atpB-rbcL intergenic spacer of cpDNA and the internal transcribed spacer regions (ITS) of nrDNA. A phylogeographical approach was used to infer past evolutionary history and processes and to examine species delimitation and population structure of the species complex. The objectives of this study were to investigate the following:

Results
Genetic variability of cpDNA and nrDNA and genetic differentiation in the Rhododendron pseudochrysanthum s.l The cpDNA atpB-rbcL intergenic spacer of R. formosanum, R. hyperythrum, and R. pseudochrysanthum s.l. ranged from 585 to 616 bp. The sequence of alignment of R. pseudochrysanthum s.l. had 585 bp, including 41 polymorphic sites, of which 15 were parsimony-informative. Most of these polymorphic sites were single-base and scattered along a sequence. The nrDNA ITS sequences of R. formosanum, R. hyperythrum, and R. pseudochrysanthum s.l. were aligned with a consensus length of 628 bp. No large indels or rearrangements were detected. The sequence alignment of R. pseudochrysanthum s.l. had 628 bp, including 133 polymorphic sites, of which 49 were parsimony informative.
With regard to cpDNA, 36 haplotypes (PH01-PH36; accession numbers: HQ850658 -HQ850693) were identified among 70 samples of the R. pseudochrysanthum s.l. (Figure 1, additional file 1: Table S1). Nucleotide diversity (π) detected in R. pseudochrysanthum s.l. was 0.0109 with high level of haplotype diversity (0.879) ( Table 1). High levels of nucleotide diversity (0.0079) and haplotype diversity (0.900) were detected in the populations of the Yushan Mountain Range, while the single population sampled from the Alishan Mountain Range was fixed at a single haplotype (Table 1). PH01 was the most common haplotype shared among the Alishan, Central, and Sheishan Mountain Ranges. In addition, no haplotypes in the Yushan Mountain Range (PH04-PH07, PH10-PH15, PH23-PH26) were shared by other mountain ranges. For nrDNA, 35 haplotypes (PR01-PR35; accession numbers: HQ850623 -HQ850657) were identified among 70 samples of the R. pseudochrysanthum s.l. (Figure 1, additional file 2: Table S2). Nucleotide diversity detected in R. pseudochrysanthum s.l. was 0.0114 with a high level of haplotype diversity (0.881). A high level of nucleotide diversity (0.0174) with a high level of haplotype diversity (0.906) were detected in the populations of the Yushan Mountain Range, while the lowest level of nucleotide diversity (0.0064) was detected in a single population of the Alishan mountains (Table 1).
To investigate the genetic structure of R. pseudochrysanthum s.l., a SAMOVA ( Table 2) was applied to define groups and to identify any locations of genetic uniqueness among the 14 populations. For cpDNA, an assumption with two groups (K = 2) displayed the greatest value of F CT and a maximal variance (F CT = 0.6167, P < 0.05). The groups at K = 2 exactly matched those found in the phylogenetic and network analyses, with those from the Yushan Mountain Range forming one group (YY, YP, YT) and those from the other three mountain ranges forming the other. Genetic structure was significant at three levels based on the spatial apportionment of genetic variation, i.e., among geographical groups (61.67%, P < 0.05), among populations within groups (4.93%, P < 0.05), and within populations (33.41%, P < 0.05).
The SAMOVA of nrDNA sequences also revealed the greatest value of F CT and a statistically maximal variance between two groups (F CT = 0.5889, P < 0.05) ( Table 2). Likewise, a large proportion of the genetic variability was partitioned among groups, although nonsignificantly (57.66%, P > 0.05). In contrast, when the three geographical groups of SL vs. SC vs. all others were compared, structuring with significant variabilities distributed among groups (54.78%, P < 0.05), among populations within groups (2.81%, P < 0.05), and within populations (42.41%, P < 0.05) was identified.
Significant genetic differentiation of cpDNA was detected between R. pseudochrysanthum s.l. Gene genealogies and phylogeographical analyses of the Rhododendron pseudochrysanthum s. l The ML tree of the cpDNA atpB-rbcL intergenic spacer haplotypes (Figure 2), rooted at R. formosanum, identified two major clusters I and II. A minimum-spanning network was also reconstructed based on the mutational steps between haplotypes of cpDNA ( Figure 3). Rhododendron pseudochrysanthum s.l. was identified with eight mutational steps from R. formosanum. Two major clusters I and II, corresponding to the ML tree, were identified with two mutational steps apart.
A ML tree and a minimum-spanning network were reconstructed of the nrDNA haplotypes, respectively (Figures 4,5). In total, clusters A-C were identified, none of which corresponded to any member of species.
We also utilized the *BEAST program to analyze the combined data of cp-and nrDNA sequences. A population tree of R. pseudochrysanthum s.l., rooted at R. formosanum (FH) and R. hyperythrum (HH), was reconstructed ( Figure 6). In contrast to the nuclear DNA tree, none of the outgroup populations were nested within R. pseudochrysanthum s.l.
Mismatch distribution analyses of cpDNA were perfomed to infer the long-term demographic history of populations. For R. pseudochrysanthum s.l., populations of the Central and Sheishan Mountain Ranges had a unimodal distribution (data not shown). Because the data used to produce mismatch distributions are not independent [52], Tajima's D was then used to detect departure from population equilibrium. Departures from the null model can be caused by many factors, such as changes in population size (e.g., expansion), which can lead to an excess of low frequency variants and negative D values. In contrast, processes such as population subdivision, or recent population bottlenecks can cause an excess of intermediate frequency variants leading to positive D values [53,54]. Most values of Tajima's D in R. pseudochrysanthum s.l. were negative, except in YT (0.325) and YY (1.057) of the Yushan Mountain Range (Table 1).
Using the program IM, we were able to refine six model parameters estimated in MCMC simulations and test for the isolation with migration model in R. pseudochrysanthum s.l. All of our parameter estimates had unimodal posterior distributions. The posterior distribution of θ other mtns (scaled effective population size of the Alishan, Central, and Sheishan Mountain Ranges to neutral mutation rate: 95.43) was larger than that of θ Yushan (scaled effective population size of the Yushan Mountain Range: 30.81), whereas the posterior distribution for θ A (scaled effective size of the ancestral population) peaked at 24.80 ( Figure 7A). Based on these θ values, the populations of the Alishan, Central and Sheishan Mountain Ranges were likely to grow substantially following divergence, whereas the populations of the Yushan Mountain Range remained stable or slightly grew relative to the ancestral population.
When the MLEs for the migration parameters were transformed into population migration rates, 2N e m = m* (θ/2) was estimated among the populations of different mountain ranges. Results from the IM revealed that the number of migrants between the Yushan Mountain Range and the other mountain ranges was effectively low. The posterior distribution of m (scaled effective migration rate) between the Yushan Mountain Range and the other mountain ranges peaked at 0.005, i.e., both m Yushan to other mtns = 0.005 (90% HPD = 0.005-1.035) and m other mtns to Yushan = 0.005 (90% HPD = 0.005-0.195), which were much lower than the m values among Alishan, Central, and Sheishan Mountain Ranges (m = 0.25-13.45) (Figure 7). All of the non-zero m values, nevertheless, implied the occurrence of historical gene flow between populations ( Figure 7B).
To elucidate the contribution of hybridization/introgression to the lack of species monophyly of R. formosanum, R. hyperythrum and R. pseudochrysanthum s.l., we  Table 1 for the detailed information of populations.

Discussion
Systematic inconsistency between nuclear and chloroplast DNAs in Rhododendron pseudochrysanthum s.l Inconsistency in the cpDNA and nrDNA phylogenies (Figure 2 & 3) was detected in R. formosanum. The ML tree for the cpDNA agreed with the results of Chamberlain et al. [25] and Goetsch et al. [28] that R. formosanum is a monophyletic group relative to R. pseudochrysanthum s.l. In contrast, the ML tree of nrDNA revealed that one accession of R. formosanum fell within cluster B with the ingroup species. Lack of species monophyly may result from incomplete lineage sorting, which occurs when populations or species recently diverged, thus retaining ancestral genetic polymorphisms, or from genetic exchanges via hybridization and introgression that introduce foreign alleles into a species or population, resulting in a genetic admixture [10]. For given species, the time for a locus to acquire reciprocal monophyly is about 1.665 × θ/4U at cpDNA or 1.665 × θ/2U at nrDNA. Willyard et al. [55] suggest that ancestral polymorphisms would likely be shared between species if the coalescence time of the alleles exceeds the divergence time at one locus. In this study, the expected time to coalescence of R. pseudochrysanthum s.l.   [56,57] revealed that natural hybridization/introgression is not unusual in Rhododendron. Thus, interspecific hybridization/introgression would cause the lack of monophyly and prolong the duration of paraphyly, as frequently occurs in most plants, e.g., Ixeris [58], Begonia [59], and Ilex [60].
Genetic diversity and population structure of R. pseudochrysanthum s.l In this study, we investigated genetic diversity at the cpDNA and nrDNA loci of R. pseudochrysanthum s.l.
High levels of nucleotide diversity were detected at both cpDNA and nrDNA loci, agreeing with a pattern found in many East Asian plants, that is considered to have been shaped by a common glaciation history [cf. [11]]. During the Quaternary glaciations, Taiwan, as a refuge, provided shelters for species from a range of biomes and thus harbored much plant biodiversity. Examples include Cunninghamia konishii [21], Pinus luchuensis ssp. taiwanensis [2], Cycas taitungensis [41], and Cycas revoluta [1]. A potential refugium for R. pseudochrysanthum s.l. in Taiwan was inferred to harbor the genetic variations. Furthermore, highly diverse environments may also have created a wide spectrum of habitats to accommodate mutations [38].
The populations of the Yushan Mountain Range had the highest levels of genetic diversity (0.0079) in the cpDNA. A similar pattern was found in Taiwan fir (Abies kawakamii), a dominant species of subalpine forests at 2500-3700 m a.s.l. [61]. This common result is likely to be attributable to a stable population size in the Yushan Mountain Range. The habitat ranges of the Yushan Mountain Range, with the highest peak at 3952 m, are wider than those of the other mountain ranges (3886 m at most), providing diverse habitats for R. pseudochrysanthum s.l. Given a stable, effective population size under diverse habitats, the populations of the Yushan Mountain Range fostered a great number of genetic polymorphisms. In contrast, genetic diversity can be lost by chance in fragmented populations, as occurred in the population of the Alishan Mountain Range. PH01, the interior dominant haplotype, is shared among the Alishan, Central, and Sheishan mountain ranges. A single population of the Alishan Mountain Range was fixed for the dominant haplotype PH01, indicating a possible genetic bottleneck due to a founder effect and/or genetic drift during a period when population sizes were small [62].
Genetic structure of R. pseudochrysanthum s.l R. pseudochrysanthum s.l. inhabits mountaintops isolated by surrounding low-elevation habitats, forming small and fragmented patches. Genetic differentiation would be detected between these isolated populations of different mountain ranges. Nevertheless, genetic analysis of the cpDNA variation revealed very weak genetic differentiation among the populations within a mountain range. The short time of isolation since the last glacial retreat is seemingly not long enough for differentiating these populations within a mountain range and thus resulted in sharing ancestral genetic polymorphisms.
Specifically, spatial distribution of cpDNA revealed a noteworthy pattern of high differentiation between populations of the Yushan Mountain Range and the other mountain ranges (F ST = 0.56-0.72) ( Table 3). The divergence time for R. pseudochrysanthum s.l. is about 2.895-8.722 MYA, and given that within R. pseudochrysanthum s.l., the populations from different mountain ranges are more differentiated than the species, then it follows that this divergence time likely reflects an estimated divergence time between mountain populations. Molecular dating revealed that populations among different mountain ranges split some 2.895-8.722 MYA, a time range approximating or predating the formation of the Central Mountain Range of Taiwan [16]. Independent colonizations via individuals with two different subclusters, likely from adjacent mainland over Pleistocene, may result in high genetic differentiation within R. pseudochrysanthum s.l. Unexpected genetic structuring was detected with populations of the Yushan Mountain Range highly differentiated from other mountain ranges, despite this range lying between the Alishan and the Central Mountain Ranges. The single population of the Alishan Mountain Range was fixed at the dominant haplotype PH01, which was shared by the populations of the Central and Sheishan Mountain Ranges. Some other haplotypes may have been randomly lost from the Alishan population due to possible bottleneck/founder's events. Despite the fact that the Alishan Mountain Range is geographically distant to Central and Sheishan Mountain Ranges, sharing the dominant PH01 made these populations clustered together and genetically undifferentiated. Furthermore, subsequent glacial retreat and mountain uplifting triggered the geographical isolations and may have thereby enhanced the divergence of populations among different mountain ranges.
Today, over an interglacial period, Rhododendron populations of different mountain ranges represent   Table 1 for the detailed information of populations. some isolated patches. Given such geographical barriers, seeds are less likely to disperse across mountain ranges. Low capabilities of seed dispersal in current environments, plus possible independent colonization events followed by limited historical gene flow (Figure 7), may have together contributed to the high levels of genetic differentiation between populations of the Yushan Mountain Range and the other mountain ranges (F ST = 0.56-0.72). In contrast, low levels of genetic differentiation at cpDNA detected among the populations of Alishan, Central and Sheishan Mountain Ranges (F ST = 0.00-0.09) implied that there may exist historical gene flow during the glacial maxima, as the alpine populations migrated downwards and became connected to each other (Table 3).
Compared to the high genetic differentiation among populations between Yushan Mountain Range and the other mountain ranges based on cpDNA, nrDNA revealed much less differentiation among populations of different mountain ranges (F ST = 0.02-0.08) ( Tables 2,  3). Here, biparentally inherited nrDNA was used to estimate the extent of pollen flow. In contrast to the limited seed dispersal, pollen of Rhododendron species seemingly can move long distances [63,64]. Long-distance pollen dispersal of R. pseudochrysanthum s.l. across the discontinuity of mountain ranges reduced the genetic differentiation.
Ancestral population size, gene flow, and population dynamics of R. pseudochrysanthum s.l Given their different inheritances, nrDNA and cpDNA markers provide extensive aspects to the understanding of the evolutionary forces shaping fates of R. pseudochrysanthum s.l. Combining cpDNA and nrDNA data, IM analysis was used to evaluate the ancestral population size and gene flow. Populations of the Alishan, Central, and Sheishan Mountain Ranges were characterized by a recent population expansion as θ other mtns was approximately four times larger than ancestral θ A [65]. Star-like genealogies of cpDNA and negative Tajima's D of cpDNA and nrDNA (Table 1) also supported a scenario of recent demographic expansion in populations of the Central and Sheishan Mountain Ranges, a common pattern also occurring in Castanopsis carlesii [66], Cunninghamia konishii [67], and Trema dielsiana [68]. In contrast, the value of θ Yushan approximated ancestral θ A , indicating that the effective population size of the Yushan Mountain Range may have remained stable relative to the ancestral population.
Given low seed dispersal capabilities of Rhododendron in today's environments, non-zero migration rates (2Nem) detected among populations (m = 0.05-13.45) likely represented historical gene flow across mountain ranges. Here, compared to the gene flow among populations of Alishan, Central, and Sheishan Mountain

Conclusions
In this study, inconsistent gene genealogies of cpDNA and nrDNA were detected. In the cpDNA tree, agreeing with Chamberlain et al. [25], R. pseudochrysanthum s.l. is phylogenetically more closely related to R. hyperythrum than to R. formosanum. In contrast, nrDNA failed to identify reciprocal monophyly of ingroup vs. outgroup, indicating possible incomplete lineage sorting and interspecific hybridization/introgression. For R. pseudochrysanthum s.l., the spatial distribution of cpDNA revealed a noteworthy pattern of high differentiation between the populations of the Yushan Mountain Range and the populations of the other mountain ranges. Molecular dating revealed that the split of the populations at different mountain ranges approximated or predated the formation of the Central Mountain Range of Taiwan island, likely suggesting independent colonizations via individuals from previously differentiated populations on the mainland. At the population level, the populations of Central, and Sheishan Mountain Ranges may have undergone postglacial demographic expansion, while populations of the Yushan Mountain Range are likely to have remained stable ever since the colonization. In contrast, the single population of the Alishan Mountain Range with a fixed cpDNA haplotype may have experienced bottleneck/founder's events.

Methods
Sample collection, DNA extraction, PCR, and DNA sequencing We assessed the levels of genetic variation of cpDNA (atpB-rbcL intergenic spacer) and nrDNA ITS sequences among populations of Rhododendron pseudochrysanthum s.l. (Table 1). In addition, R. formosanum Hemsl. (subgenus Hymenanthes subsection Argyrophylla [25]), and R. hyperythrum (subgenus Pontica subsection Maculifera) were selected to serve as outgroups. In the field, R. hyperythrum is sympatric with R. pseudochrysanthum s.l., while R. formosanum is allopatric from two other species. Morphologically, R. hyperythrum and R. formosanum are distinguishable from R. pseudochrysanthum s.l. According to Goetsch et al. [28], R. pseudochrysanthum s.l. is phylogenetically more related to R. hyperythrum than to R. formosanum. A total of 70 R. pseudochrysanthum s.l. individuals were collected from 14 populations throughout their natural range, mostly from four major mountain ranges in Taiwan, i.e., Alishan, Central, Sheishan and Yushan Mountain Ranges (Figure 1 and Table 1). One to 12 individuals about 100 m apart were sampled based on the natural size of each population. Young, healthy leaves were collected in the field, rinsed with tap water, and dried in silica gel. All samples were stored at -70°C until processing.
Leaf tissue was ground to powder in liquid nitrogen and stored at -70°C. Total genomic DNA was extracted from the powdered tissue using a cetyltrimethyl ammonium bromide (CTAB) procedure [69]. PCR amplification was carried out in a 100-μL reaction. The reaction was optimized and programmed on a MJ Thermal Cycler (PTC 100) (MJ, Alameda, California, USA) as one cycle of denaturation at 95°C for 4 min; 30 cycles of 45 s denaturation at 92°C, 1 min 15 s annealing at 52°C, and 1 min 30 s extension at 72°C; followed by 10 min extension at 72°C. Template DNA was denatured with reaction buffer, MgCl 2 , NP-40 and ddH2O for 4 min (first cycle), and cooled on ice immediately. Primers, for the ITS (ITS5: 5'-GGAAGTAAAAGTCGTAACAAGG-3' and ITS4: 5'-TCCTCCGCTATATGATATGC-3') [70] or cpDNA atpB-rbcL intergenic spacer (atpB-1: 5'-ACATCKAR-TACKGGACC AATAA-3' and rbcL-1: 5'-AACACCA GCTTTRAATCCAA-3') [71], dNTP and Taq polymerase (Promega, Madison, Wisconsin, USA) were added to the above ice-cold mix. The reaction was restarted at the first annealing at 52°C. PCR products were purified by electrophoresis in 1.0% agarose gel using 1X TAE buffer. The gel was stained with ethidium bromide, and the desired DNA band was cut and eluted using agarose gel purification (Qiagen, Valencia, California, USA). Both direct and subcloning sequencing was conducted to check PCR error. Chloroplast DNA sequences were determined by direct sequencing, while nrDNA was sequenced via subcloning. For subcloning, the desired DNA fragments were cloned into the pGEM-T vector system (Promega, Madison, Wisconsin, USA). Seven clones of the plasmid DNA were selected randomly and purified using a plasmid mini kit (Qiagen). DNA sequencing in both directions was conducted with an Applied Biosystems Model 377A automated sequencer (Applied Biosystems, Foster, California, USA). No within-individual variation was detected for either DNA marker.

Data analysis
Nucleotide sequences were aligned with the program MAFFT version 6 [72] and later adjusted visually. After the alignment, maximum-likelihood (ML) analyses (HKY model) were performed with PhyML v3.0 [73], and bootstrap consensus values were calculated with 1000 bootstrap pseudoreplicates. Nucleotide substitution models were determined with the Akaike Information Criterion by ModelTest 3.7 analysis [74,75]. The HKY + I + G model with 6 substitution categories was determined to be the most suitable model by Modeltest and was used for all subsequent nucleotide analyses. The number of mutations between DNA sequence haplotypes in pairwise comparisons was used to construct a minimum-spanning network with the aid of the MINSPNET [76] in a hierarchical manner [77].
To reconstruct the phylogeny of populations, the combined data of both cp-and nrDNA sequences were analyzed by the program *BEAST v1.5.4 (Bayesian Evolutionary Analysis Sampling Trees) [78]. This method uses coalescent theory to provide joint inferences of a species tree topology, divergence times, population sizes and gene trees from multiple genes sampled from multiple individuals across a set of closely related species.
Levels of genetic diversity within populations and species were quantified with measures of nucleotide diversity (π) [79] using DnaSP Version 5.1 [80]. Patterns of geographical subdivision and levels of genetic differentiation among populations were estimated hierarchically with the aid of DnaSP. Mismatch distribution analyses and Tajima's D were also conducted [81].
SAMOVA (spatial analysis of molecular variance) was applied to identify groups of geographically adjacent populations of R. pseudochrysanthum s.l. that were maximally differentiated based on sequence data [82]. We performed the analyses based on 100 simulated annealing steps and examined maximum indicators of differentiation (F CT values) when the program was assigned to identify K = 2-6 partitions of populations.

Isolation with migration and estimating of ancestral population size
We used the simulation program IM [83] to investigate the scenario of isolation with migration between populations. By applying coalescence simulations and Bayesian computation procedures, IM yielded six model parameters, including the population-split time (t = tμ), scaled migration rates (m = m/μ), and scaled effective population size (θ = 4Neμ) for the ancestral and two current populations. The posterior probability densities of these parameters were generated by Markov chain Monte Carlo (MCMC) simulations, and three of these simulations were run with individual simulations being updated 50 million times. Within each simulation, we used a procedure to swap among 10 heated chains (Metropolis coupling). Each simulation yielded a marginal density histogram for the six parameters of interest, and the peaks of the resulting distributions were considered as the maximum likelihood estimates (MLE) of the parameter with credibility intervals equaling the 90% highest posterior density (HPD) intervals. Here we compared populations of the Yushan Mountain Range with other populations of three different mountain ranges based on the results of F ST analyses. Parameters in the IM model are scaled by the mutation rates per generation for the loci, i.e., μ = UG, where U is the geometric mean of mutation rates of all the loci and G is the generation time.

Estimation of divergence times of cpDNA and nrDNA lineages
In the present study, we estimated the divergence time of the sister lineages of R. pseudochrysanthum s.l. split for both cp-and nrDNAs using a Bayesian relaxed clock method with BEAST. We excluded sequences of 18S, 5.8S and 26S RNA regions of nrDNA from the calculation of the divergent time. To estimate the divergence between lineages a well-supported mutation rate is required. In this study, applying evolutionary rates of 1.0-3.0 × 10 -9 [84] and 5.0-7.8 × 10 -9 [85] substitutions per site per year for synonymous sites of cpDNA and nrDNA, respectively, we estimated the divergence times of cpDNA and nrDNA lineages for R. hyperythrum and R. pseudochrysanthum s.l. Additionally, Bayesian estimates of the divergence time of the MRCA of the R. pseudochrysanthum s.l. clusters were obtained using BEAST v. 1.3; available from http://beast.bio.ed.ac.uk/ Main_Page [86]. Following the result of ModelTest analysis [75], the HKY model of nucleotide substitution with estimated base frequencies, gamma shape distribution (with six categories), proportion of invariant sites, and a relaxed clock with uncorrelated log normal distribution of branch lengths were chosen. Posterior estimates of the mutation rate and divergence time were obtained by Markov chain Monte Carlo analysis, with samples drawn every 500 steps over a total of 1 million steps. All operators were automatically optimized. Convergence of parameters and mixing of chains were followed by visual inspection of parameter trend lines and checking effective sampling size (ESS) values with three pre-runs. The ESS parameter was found to exceed 100, which suggests acceptable mixing and sufficient sampling. Adequate sampling and convergence to a stationary distribution were checked using TRACER v. 1.3 [87].
There are two primary processes that may cause a locus to show non-monophyly for a group of closely related species, i.e., incomplete lineage sorting and hybridization/introgression. Distinguishing between these two processes is difficult, nevertheless, estimates of divergence time can provide a means of assessing the probability that incomplete lineage sorting alone could account for the observed non-monophyly. The expected time for a locus to become monophyletic following divergence from a sister lineage is 1.665 × NeG and 1.665 × 2NeG for cp-and nuclear DNAs, respectively [55,88]. According to IM [83], the effective population size can be estimated as Ne = θ/(4UG). Therefore, the former is then equal to 1.665 × θ/4U, while the latter is 1.665 × θ/2U. Here, θ values were obtained from IM analyses, and U is the geometric mean of mutation rates of cp-and nrDNAs.