Phylogeny and divergence times inferred from rps16 sequence data analyses for Tricyrtis (Liliaceae), an endemic genus of north-east Asia.

BACKGROUND AND AIMS: Tricyrtis is a genus of monocots with attractive and sophisticated flower shapes and colours, endemic to north-east Asia. There are 18 known species. The highly restricted geographical distribution of the genus is of great interest in terms of both abiotic (continental drift) and biotic (long-distance dispersal) impacts on monocot plant speciation events and their timing, and of evolutionary patterns of diversification leading to the extant taxa. The aims of this study were to (i) predict the time of speciation (divergence) events at infraspecific levels of Tricyrtis, (ii) estimate the rate of evolution of the genus and (iii) provide information on an excellent plant model system in terms of studying loss of biodiversity or extinction of organisms in the dynamic earth environment.
METHODOLOGY: To investigate the divergence time and evolution rate of Tricyrtis, Bayesian Markov chain Monte Carlo (MCMC) analyses were performed by calculating the mean branch lengths of evolutionary paths based on base substitution variations between rps16 intron nucleotide sequences from the 18 known species. PRINCIPAL
RESULTS: Based upon the relaxed molecular clock model test data, a Bayesian phylogenetic inference tree is presented, and the divergence times and rate of evolution of Tricyrtis were estimated. The analyses also suggest that evolution is occurring at the infraspecific level of the genus in a manner that is not strictly clock bound.
CONCLUSIONS: Continental drift may have been the main speciation process giving rise to the current distribution of the taxa of Tricyrtis. The single-locus gene sequence data presented here are a significant step towards an improved future understanding of the molecular evolution of Tricyrtis via multi-locus evaluation.

Introduction

Investigations of evolutionary diversification, combining current phyto-geographical information with geological data (Särkinen ) and biotic and/or abiotic elements (Tiffney and Mazer 1995), or predictions of biodiversity loss (Cox 2001; Takahashi ) can be of great interest to plant systematists. In the study of evolution, estimates of divergence time can provide crucial information on speciation events (Barraclough and Vogler 2000; Barraclough and Nee 2001; Nee 2001; Nichols 2001; Espert and Burghardt 2010). The high endemism of the monophyletic genus Tricyrtis in north-east (NE) Asia is of particular interest to both molecular evolutionists and plant systematists (Hong 1999; Hong ; Maki ; Fay ; Rudall ; Takahashi and Maki 2007; Takahashi ) in the context of speciation mode and rate, adaptive radiation, and diversification routes and direction. There are 18 known species in the genus, most being found in wet and shaded areas at relatively high altitude in Nepal, China, South Korea, Japan and Taiwan (Takahashi 1980, 1987; Peng and Tiang 2007).

The geographical mode of speciation (e.g. allopatric versus sympatric) of both ancestral and extant species can be deduced by examining the current distribution patterns of sister taxa using interspecific phylogeny approaches (Losos and Glor 2003). Depending upon whether the distribution pattern of sister taxa overlaps geographically or not, either sympatric or allopatric speciation modes can be assumed (Dieckmann and Doebeli 1999; Losos and Glor 2003). It has also been suggested that the driving force underlying the geographical mode of speciation (adaptive evolution, vicariance) can be inferred from phylogeny (Avise and Walker 1998; Losos and Glor 2003).

Ecological and evolutionary factors including both biotic and abiotic elements (Popp ) have been explored by researchers to explain the relationships between areas of endemism and evolutionary history (Ronquist 1997; Särkinen ). A co-evolutionary role of pollinators (a biotic element) in flower evolution is relatively well known (Pérez ; Kirchoff ). In particular, bee or hummingbird pollinators in Schizanthus (Solanaceae) have been described as important factors in adaptive evolutionary processes of this genus in the Chile region (Pérez ).

The vascular plants occurring in the EAS (Eastern Asia)–ENA (Eastern North America) disjunct floras (Wen and Zimmer 1996; Qian and Ricklefs 2000; Xiang ) provide intriguing potential insights into monocot phylogenetic history. Qian and Ricklefs (2000) showed that ecology and age are major elements in the generation of species-rich ‘diversity anomalies’ in these regions, which include the Japanese Islands and the Korean Peninsula. The divergence time of the EAS–ENA floras, based upon rbcL sequence data using Cornus and other sister taxa as a model phylogenetic system (Xiang ), has been provided as supportive evidence for paleontological estimates (Miocene–Pliocene).

Reconstructions of plant phylogenetic trees are carried out by estimating clade ages (divergence time) considering non-neutral change rates of DNA sequence data (‘relaxed clock’) as priors (Kumar and Gadagkar 2001) as well as biogeographic patterns and fossil data as posteriors, using Bayesian inference-derived software programs (Renner 2005; Pirie ; Antonelli ). These include BEAST (Bayesian Evolutionary Analysis by Sampling Trees), MrBayes, Mesquite and Bali-Phy (http://evolution.genetics.washington.edu/phylip/software.html).

Douzery measured the evolutionary time of the major clades using a large dataset of 129 proteins from 36 eukaryotes. They analysed these datasets taking a Bayesian relaxed molecular clock approach combined with fossil calibrations and global molecular clock constraint conditions. Their estimates using the relaxed clock model suggested more recent divergence events when compared with the global molecular clock models. Clearly, divergence time calculations can be clock model dependent.

In the current study, all 18 known taxa of Tricyrtis have been included in DNA sequence data analysis using the rps16 intron (non-coding plastid DNA) region. This comprehensive approach encompassing the whole genus (even though initially only a single locus has been studied) has not been reported previously in estimates of divergence time and evolutionary rates using Bayesian Markov chain Monte Carlo (MCMC) analyses. Therefore, in this study, divergence time and rate data were calculated and used to infer the geographical mode of speciation, and to interpret the roles of biotic factors (interactions with endemic insects) and abiotic factors (break-up of tectonic plates) in the formation of current areas of endemism of Tricyrtis.

Materials and methods

Taxon sampling

Eighteen DNA sequences were obtained from the rps16 intron region (∼729 base pairs) using fresh and dried specimens of Tricyrtis and its outgroups (North American genera, Calochortus and Scoliopus) shown in Table 1. The world distribution map of Tricyrtis is shown in Fig. 1. Most specimens in the study were checked directly with herbariumType specimens from the Royal Botanic Gardens, Kew, UK, and the Natural History Museum, London, UK. The primers for flanking the rps16 intron region of Tricyrtis were rpsF (GTGGTAGAAAGCAACGTGCGACTT) and rpsR2 (TCGGGATCGAACATCAATTGCAAC). For the polymerase chain reactions, the conditions were 97 °C (50 s), 50 °C (50 s) and 72 °C (1 min + 50 s) at each step. Thirty cycles were carried out and the delay time was 7 min at 72 °C. The total genomic DNAs of Scoliopus (Trilliaceae) and Calochortus (Calochortaceae) were obtained directly from the Jodrell Laboratory, Royal Botanic Gardens, Kew, UK.

Table 1

Plant material source information.

Taxa	Localities	Collection informationa	Accessionb (NCBI/GenBank)
Sect. I. Brachycyrtis
T. macranthopsis	Japan, Honshu, Wakayama Prefecture	BMA (fresh)	JN561092
T. ishiiana var. surugensis	Japan, Honshu, Shizuoka Prefecture	BSWJ (fresh)	JN561093
Sect. II. Flavae
T. flava	Japan, Kyushu, Miyazaki Prefecture	BMA (fresh)	JN561095
T. nana	Japan, Honshu, Mie Prefecture	BSWJ (fresh)	JN561094
T. perfoliata	Japan, Hyugu Province, Miyazaki Prefecture, Mt Osuzu, alt. 500 m	HK: 1989–3971	JN561097
T. ohsumiensis	Japan, Osumi Province, Kagoshima Prefecture, Mt Hoyoshidake	BMA (fresh)	JN561096
Sect. III. Hirtae
T. formosana var. amethystina	Taiwan, Ta Ta Chia, alt. 2650 m	BSWJ 237	JN561099
T. hirta var. hirta	Japan, Honshu, Idzu Province, 250 m Mt Amagi	BSWJ (fresh)	JN561100
T. suzukii	Taiwan (wild)	JCW (fresh)	JN561101
Sect. IV. Tricyrtis
T. affinis	Japan, Honshu, Yamashiro Province, Mt Hieizan	HK: 1989–3903	JN561105
T. latifolia	China, Hubei, alt. ∼1300 m Japan, Honshu, Arakawamachi in Echigo	HK: 1981–4281	JN561102
T. macropoda	S. Korea, Wando Is., alt. 40 m	HK: 1989–3904	JN561104
T. maculate	Nepal, Arun Valley	HK: 1992–202	JN561103
T. setouchiensis	Japan, Tokushima Prefecture, Mt Kotsusan	BMA (fresh)	JN561098
Outgroups
Scoliopus bigelovii		Chase 441K (1976–1909)	JN561106
Calochortus minimus		Ness 606, PUA	JN561107

aThe abbreviations refer to the following: BMA, Brian Mathew (Kew); BSWJ, Bleddyn & Sue Wynn-Jones (Wales); HK, Living collections, Royal Botanic Gardens, Kew; JCW, Jenn-Che Wang (Taiwan).

bRefers to GenBank (EMBL/NCBI) accession number for rps16 intron (chloroplast DNA) sequences.

Fig. 1

The current world distribution map of Regions A and B refer to the localities of sect. Brachycyrtis (blue boundary) and sect. Flavae (red), respectively. Those areas are also magnified in the box. The letter C indicates the localities of sect. Hirtae (green), and D shows the region where sect. Tricyrtis plants are known to grow (purple).

Bayesian analyses

BEAST v1.6.1 (Drummond , ) was used for the Bayesian MCMC inferred analyses of the nucleotide sequence data (Drummond and Rambaut 2007). BEAUti (Bayesian Evolutionary Analysis Utility version) v1.6.1 (Drummond , ) was utilized to generate initial xml files for BEAST. A Yule process of speciation (‘a pure birth’ process) was used as a tree prior for all the tree model analyses. The Yule tree prior is widely recognized as giving the best fit model for trees describing the relationships between different species. The parameter can be regarded as explaining the net speciation rate (Nee 2006).

For the MCMC posterior analyses, the length of chain was 10 000 000. After 100 tree burn-in processing, 10 000 trees were used for the analyses (Table 2). The BEAUti xml file (Drummond , ) was run in the BEAST v1.6.1 program (Drummond , ) and the maximum clade credibility (MCC) chain generations were repeated five times for each molecular clock model with independent runs to ensure suitable convergence and adequate mixing. The MCC tree in Fig. 2 was generated under the relaxed clock model (HKY substitution).

Table 2

Overall estimated evolutionary rate using molecular clock models.

Taxa	Substitution modela	Site heterogeneity model	Clock model	Mean rateb (SSMY)
Tricyrtis	HKY	G + I	Strict clock	0.000418
Outgroup	HKY	G + I	Relaxed clock (uncorrelated lognormal)	0.000380
Tricyrtis	GTR	G + I	Strict clock	0.000418
Outgroup	GTR	G + I	Relaxed clock (uncorrelated lognormal)	0.000381

aThe abbreviations refer to the following: HKY, Hasegawa, Kishino, Yano (1985); G + I, gamma + invariant sites; GTR, general time reversible.

bSSMY refers to the unit of mean rate of evolution as base substitutions per site per million years.

Fig. 2

An MCC tree generated via Bayesian analyses. The divergence times are shown in scales (Myr). The grey bars at the nodes refer to the 95% HPD range for estimating divergence times showing node age errors. If the posterior probability is >0.5, only those nodes were labelled with the posterior probability to show dating with confidence. The branch colour coding refers to the following rate variations: (i) a posterior rate > average rate, shown in red; (ii) a posterior rate < average rate, shown in blue; (iii) a posterior rate is close to the average, shown in black; (iv) a posterior rate is slightly lower than the average, shown in brown. Underneath the scale at the bottom of the tree, the number 4.0 means 4 Myr in length.

LogCombiner v1.6.1 (Rambaut and Drummond 2010, ) was used to combine the log files from the independent BEAST runs. Tracer v1.5 software (Rambaut and Drummond 2009) was used for the output of the model parameters to examine the sampling and convergence results obtained from BEAST. TreeAnnotator v1.6.1 software (Rambaut and Drummond 2010, ) was used to annotate the phylogenetic results generated by BEAST as a form of single ‘target’ tree. On the target trees are shown summary statistics of posterior probabilities of the nodes: the 95% highest posterior density (HPD) limits of the node heights, rates and the posterior estimates. For the annotated BEAST MCC tree output analyses, the FigTree v1.3.1 (Rambaut 2009) program was also used (Fig. 2). The posterior probability limit (Drummond ) was set to 0.5 (Fig. 2). This is equivalent to the bootstrapping value in PAUP (Phylogenetic Analysis Using Parsimony analysis) analyses.

The topologies of the phylogenetic trees, and 95% HPD clade support statistic values, were examined using various Bayesian sampling tree analyses. The phylogenetic reconstruction and the topology of the MCC trees were compared with the previous studies by the current authors using maximum parsimony approaches (unpublished data).

Estimating the evolutionary rate

The ages of internal nodes in this study were calibrated based on the information shown in rbcL coding cpDNA sequence data analyses by Vinnersten and Bremer (2001). According to these authors, the age of the genus Tricyrtis (Eurasian origin, Tricyrtis affinis as a representative taxon) was estimated at ∼22 million years (Myr), and they also reported the ages of Calochortus and Scoliopus as 40 and 5 Myr, respectively, based on their fossil and molecular data (Vinnersten and Bremer, 2001). In the current analyses, both clock.rate and meanRate parameters were used in molecular clock models to calculate the rate of evolution (Table 2, Fig. 2). The units for rate of evolution in the strict molecular clock model using the clock.rate parameter are substitutions per site per million years (Drummond ). The equivalent meanRate parameter is used in the relaxed molecular clock model, and this was calculated using Equation (1) (Drummond ): For r= rate on the ith branch and t = length of time (time units for the ith branch), then b is equivalent to r(the branch length in substitutions per site). The meanRate of this example can be derived from Equation (1).

Evolutionary models

As the substitution models used sequence data from a total of 18 species, the HKY85 (Hasegawa, Kishino and Yano 1985) and GTR (general time-reversible) substitution models were utilized and compared (Table 2). The site heterogeneity model was G (gamma) + I (invariant sites) in both strict and relaxed clock models, allowing rate variations among the sites in the alignment (Table 2; Drummond ).

Divergence time estimation

In the BEAUti v1.6.1 program, the ‘relaxed phylogenetics’ models were tested with various evolutionary rates for divergence time estimation (Drummond ). A strict clock model was also tested for comparison with the relaxed clock model (Table 2). For the strict clock model, the clock rate was set to unity. The assumption underneath the strict clock model is that a global clock rate has no variation among lineages (Drummond ). On the other hand, the relaxed clock model (BEAST software, Drummond , ) requires no a priori rate correlation between its ancestor and the lineages. In the current analyses, the uncorrelated lognormal distribution model was taken under the relaxed clock model in BEAST v1.6.1 (Drummond , ). This model was used to test whether the current data are consistent with a clock-like fashion.

Results

Estimated mean rate of evolution under the molecular clock models

The overall estimated mean rate of evolution is presented in Table 2. When two molecular clock rate variation models (strict versus relaxed clock) were tested, the mean rate of the individual clock model did not show any variation regardless of whether different substitution models were used, including HKY and GTR (Table 2). For instance, the mean rate of evolution was 0.000418 substitutions per site per million years (SSMY, 95% HPD: 0.000296–0.000553) under the strict clock model, and 0.000380 SSMY (95% HPD: 0.000267–0.000505) when the relaxed clock model was tested, as shown in Table 2. The 95% HPD is regarded as a Bayesian representation of confidence interval (Drummond and Rambaut 2007). How calibration and estimation were done for the mean rate of evolution is explained in Materials and methods.

The mean coefficients of variation (σγ) under the relaxed clock model were 0.99 (HKY) and 1.00 (GTR) (Table 2), indicating that Tricyrtis DNA sequences are evolving in a non-clockwise fashion. This also shows that a significant level of rate heterogeneity between lineages is present (Drummond ).

Branch rate distribution

The mean rate was in the range of 0.0003–0.0015 SSMY at the infraspecific level of Tricyrtis (Table 3, Fig. 3). The fastest rate was found in the clade branch in sect. Flavae, where Tricyrtis ohsumiensis shows robust support at the internal node having a posterior probability value of 1.0 (Tables 2 and 3; Figs 2 and 3). The rate of evolution in that branch was 0.0015 SSMY (95% HPD: 0.0003–0.0034). Considering the mean rate estimates of between 0.0038 (95% HPD: 0.0027–0.0051) and 0.0042 (95% HPD: 0.0003–0.0055) SSMY (Table 2), the rate of evolution above is 3.6–4.0 times faster than the average. The next fastest rate of evolution was found at the external node, where the clade of sect. Hirtae diverged from the rest of the clades of Tricyrtis (Table 3, Figs 2 and 3). At that node, the posterior probability was 1.0 and the rate of evolution was 0.0012 SSMY (95% HPD: 0.0003–0.0025). This rate is 2.9–3.2 times faster than the average.

Table 3

The node age and rate of evolution with the posterior probability densities are shown.

Node numbera	Mean age	Mean rateb	95% HPDc
	(node, Myr)	(SSMY)	(SSMY)
1	3.2767	0.0015	0.0003–0.0034
2	19.9648	0.0012	0.0003–0.0025
3	9.1026	0.0004	0–0010
4	3.7729	0.0004	0–0011
5	6.6387	0.0004	0–0011
6	1.3031	0.0003	0–0008
7	33.1201	0.0005	0–0014

aOnly those nodes with posterior probability > 0.5 are numbered in Fig. 3.

bSSMY refers to the unit of mean rate of evolution as base substitutions per site per million years.

cPresents lower–upper 95% HPD intervals, respectively.

Fig. 3

A Bayesian inference tree with estimated divergence times. Where the posterior probability at a node is above the threshold value of 0.5, this is shown in the figure. The grey bars at the nodes refer to 95% confidence intervals for estimating divergence times. The branch colour coding refers to the following rate variations: (i) a posterior rate > average rate, shown in red; (ii) a posterior rate < average rate, shown in blue; (iii) a posterior rate is close to the average, shown in black; (iv) a posterior rate is slightly lower than the average, shown in brown. Underneath the scale at the bottom of the tree, the number 4.0 means 4 Myr in length.

The external node of sect. Hirtae, where it branches from the clade of sect. Tricyrtis and sect. Brachycyrtis, shows the rate as 0.0004 SSMY (95% HPD: 0–0.0010) with robust support (Table 3, Figs 2 and 3). Within the sect. Hirtae, the rate of Tricyrtis amethystina and T. formosana at the internal node was also estimated as 0.0004 SSMY (95% HPD: 0–0.0011), indicating close to the average rate (Table 3, Figs 2 and 3). The sister clade of Tricyrtis maculata and T. latifolia shows the average rate 0.0004 SSMY (95% HPD: 0–0.0011; Table 3, Figs 2 and 3). Another sister clade of Tricyrtis perfoliata and T. flava gives a relatively slower rate than that of the average as 0.0003 SSMY (95% HPD: 0–0.0008; Table 3, Figs 2 and 3). At this node, the posterior probability was only at the threshold level, 0.5 (Fig. 2).

Divergence time and age estimation

In Table 3 and Fig. 3, the estimated mean age is presented with 95% HPD intervals in million years at those nodes with posterior probability above the threshold level. The posterior distribution of the evolutionary root age was similar to the prior distribution pattern, which is an indication of an appropriate choice of parameters.

When the strict clock model (HKY) was tested, the posterior distribution pattern of the mean age of the root of the tree was 19.91 Myr (95% HPD: 18.93–20.88) in Tricyrtis, and the root age of the outgroup (Calochortus) was 40.01 Myr (95% HPD: 39.05–41.01; Table 4). In the relaxed clock model test (HKY), the posterior distribution pattern of the mean age of the root of the tree was 19.96 Myr (95% HPD: 19.02–20.95) in Tricyrtis, and the root age of the outgroup was 40.0 Myr (95% HPD: 39.03–40.97; Table 4). These posterior distribution values are very similar to the specified prior distributions in the analyses. All the effective sample sizes were more than 500, indicating a good mixing of tree sampling.

Table 4

Estimated divergence time and age using molecular clock models.

Taxa	Substitution modela	Clock model	Mean ageb	95% HPDc
Tricyrtis	HKY	Strict clock	19.91 (0.0053d)	18.93–20.88
Outgroup	HKY	Strict clock	40.01 (0.0053)	39.05–41.01
Tricyrtis	HKY	Relaxed clock (uncorrelated lognormal)	19.96 (0.0057)	19.02–20.95
Outgroup	HKY	Relaxed clock (uncorrelated lognormal)	39.99 (0.0053)	39.03–40.97
Tricyrtis	GTR	Strict clock	19.91 (0.0054)	18.96–20.87
Outgroup	GTR	Strict clock	40.01 (0.0052)	39.00–40.97
Tricyrtis	GTR	Relaxed clock (uncorrelated lognormal)	19.97 (0.0055)	18.96–20.90
Outgroup	GTR	Relaxed clock (uncorrelated lognormal)	39.98 (0.0053)	38.99–40.92

aThe abbreviations refer to the following: HKY, Hasegawa, Kishino, Yano (1985); GTR, general time reversible.

bDenotes that the unit of mean age is in million years.

cRefers to lower and upper 95% HPD intervals, and the units are in million years.

dIndicates standard error of the mean (SEM).

When the strict clock model using the GTR substitution model was tested, the posterior distribution pattern of the mean age of the root of the tree was 19.91 Myr (95% HPD: 18.96–20.87) in Tricyrtis, and the root age of the outgroup (Calochortus) was 40.01 Myr (95% HPD: 39.0–41.0; Table 4). In the relaxed clock model test, the posterior distribution pattern of the mean age of the root of the tree was 19.96 Myr (95% HPD: 18.96–20.90) in Tricyrtis, and the root age of the outgroup was 39.98 Myr (95% HPD: 38.99–40.92; Table 4).

Discussion

The estimated divergence times and rates of evolution of Tricyrtis point to the primary cause of vicariance of the genus being an abiotic factor, continental drift, rather than a biotic one, long-distance dispersal. Furthermore, the centre of adaptive radiation of Tricyrtis may have been the Japanese Islands, during the Miocene–Pleistocene period, and the direction of diversification can be seen (Figs 1 and 3) to be moving towards the Himalayan region in the continuing macroevolutionary process.

The estimated divergence time of most clades of Tricyrtis at 11 Myr ago (Miocene) agrees approximately with the separation times from their parent tectonic plates of the Japanese Islands and Taiwan, at 3 and 20 Myr ago, respectively (Taira 2001; Huang ). The speciation events leading to the split of NE Asian Tricyrtis from Calochortus and Scoliopus may have occurred during the Eocene period (more than 40 Myr ago; Fig. 3). It is unclear how the EAS–ENA disjunct distribution pattern of ancestral clades of Tricyrtis and Calochortus was initially formed at the generic level.

Paleogeologically, the Japanese Islands originated from the North American/Pacific/Philippine Sea tectonic plates, and the separation of Japan from the Eurasian plate (including China and the South Korean Peninsula) has been estimated as occurring about 30 Myr ago (Taira 2001). This suggests that the split between Tricyrtis and Calochortus occurred before the break-up of these tectonic plates. However, the possibility of molecular data giving an erroneously early date for this split cannot be ruled out and in some instances the fossil record may indicate a more recent origin than suggested by the molecular data (Douzery ).

The Tricyrtis and Calochortus genera show both similarities and dissimilarities (Patterson and Givnish 2002; Givnish ) in various character states, and they are assumed to share a common ancestor in the deep node (Fig. 3). The geographical distribution patterns of ancestral groups of the two taxa, Calochortus and Scoliopus, were investigated by reconstructing trees (Patterson and Givnish 2002). According to previous reports (Patterson and Givnish 2002; Givnish ), both Tricyrtis and Calochortus share character states such as having passively dispersed fruits and large, showy flowers. However, they also show differences in having rhizome (Tricyrtis) or bulb (Calochortus)-type storage organs, and narrow and parallel versus broad and reticulate leaves (Patterson and Givnish 2002). Even though Tricyrtis produces septicidal capsule-type fruits throughout the whole genus, long-distance dispersal may not be the major factor explaining the geographical mode of speciation at infraspecific level. However, at the generic level (between Tricyrtis and Calochortus), seed dispersal type may have played a significant role at the time of speciation before the tectonic break-up.

Long-distance dispersal by wind may not be a primary factor explaining the mode of geographical distribution of Tricyrtis, considering the preferred habitats of most Tricyrtis species in wet areas and that they have no clear-cut means of wind dispersal for crossing abiotic barriers. However, based upon previous reports, it is very tempting to relate the speciation events to co-evolution with insect pollinators. The pollination mechanisms of Tricyrtis have been reported (Takahashi 1987, 1994; Maki ). The pollen- and nectar-foraging visitors of Tricyrtis flowers (Takahashi 1987, 1994) include (i) Ceratina japonica and Tetralonia sp. (Anthophoridae), (ii) Bombus diversus diversus (Apidae), (iii) Megachile tsurugensis (Megachilidae), (iv) Lasioglossum occidens, L. mutilum, L. duplex (Halictidae), etc. Interestingly, the insects listed above show highly restricted distribution patterns (Global Biodiversity Resources Discovery System) even though insects can in principle overcome abiotic barriers either by flying around them or being carried by other objects. For instance, L. mutilum is distributed only in southern parts of the Japan Islands (Kyushu, Shikoku, Honshu), and M. tsurugensis also has a highly restricted occurrence in areas in southern parts of Japan as well as in the Himalayas (Takahashi 1987). In addition, L. occidens shows a disjunct distribution pattern between North America and South East Asia. Therefore, it is reasonable to assume that limited and specific visitors of Tricyrtis flowers may have played significant roles in the speciation and diversification of the genus.

A revisionary taxonomic study of Tricyrtis was reported by Takahashi (1987), and the monophyletic genus was divided into four sections: Brachycyrtis; Flavae; Hirtae and Tricyrtis (Takahashi 1980, 1987, ). Classical taxonomic and other approaches (population genetics, genetic diversity, pollination mechanisms and geographical distribution patterns) by Takahashi since the 1970s have provided a solid basis for using Tricyrtis as an excellent model system for the study of speciation rates. However, he and others have not carried out any molecular evolutionary studies using a species-level phylogeny of the genus from DNA sequence data. Therefore, the current molecular estimation of divergence time and rate of Tricyrtis provides additional information for measuring the speciation rate of the genus and for evaluating its macroevolution.

The current divergence time data of Tricyrtis species fit relatively well with the accepted tectonic ‘break-up’ time of the Japanese Islands and Taiwan areas during the Miocene–Pleistocene periods (Tables 3 and 4). Based upon the current data in this study, clades of Tricyrtis at sectional levels may have gone through a speciation process from their common ancestor, resulting in two different groups (sect. Flavae and other remaining sections; Figs 2 and 3) after the genus split out from other close relatives of monocot genera. Considering the divergence time (∼13–14 Myr ago) of sect. Flavae, which appeared later than that of other remaining sections (Fig. 3), species in this section (endemic to Japan) may have undergone an allopatric speciation process caused by the break-up of tectonic plates. The mid-Miocene period is when the Southwest Japan Arc is assumed to have split from the Eurasian Plate (Taira 2001). This Southwest Japan Arc area includes South-West Kyushu, Shikoku and Honshu where the extant sect. Flavae species (Tricyrtis nana, T. flava, T. perfoliata, T. ohsumiensis) occur endemically. The common ancestor of other remaining taxa (sects. Hirtae, Brachycyrtis and Tricyrtis) may have gone through a sympatric speciation process considering its earlier time of appearance during the Miocene, close to the common ancestor of all the Tricyrtis species (Fig. 3). Tricyrtis nana has been suspected of being the ancestor of the other species in sect. Flavae due to its low genetic diversity and relatively broad distribution patterns compared with the three other species (Maki ; Takahashi ). The current data give supportive evidence that T. nana may have originated ∼3Myr before those other taxa (Fig. 3).

Using the current data, the taxa endemic in Taiwan, such as T. amethystina and T. formosana, are found to have diverged from their most common ancestor around 3.8 Myr ago (Pleistocene; Table 4, Fig. 3), and this fits well with the predicted time of tectonic separation of Taiwan from the Eurasian plate (including China). In the case of Tricyrtis hirta, which occurs in both the Japanese Islands and Taiwan, this must have diverged from their most common ancestor ∼9 Myr ago (Table 4, Fig. 3). The speciation mechanism of T. hirta is not clear based upon the current data, and should be explored further in the near future.

Tricyrtis maculata is the only known species that occurs in the Himalayan region, and its divergence time is shared with T. latifolia at ∼6.6 Myr (Table 3, Fig. 3). This is intriguing evidence, which can be related to the current spatial distribution pattern of those two taxa. Tricyrtis latifolia grows in China and Taiwan and seems likely to be older than other Taiwanese taxa (Table 4, Fig. 3). This means that the diversification of the extant Taiwanese taxa may have occurred after the tectonic separation of Taiwan from the Eurasian plate. Furthermore, T. maculata may have gone through its speciation process much later than the time when the most common ancestor of other Tricyrtis taxa had split, indicating that the diversification direction is from Japan towards the Himalayas. Tricyrtis macropoda diverged from the common ancestor slightly earlier than either T. maculata or T. latifolia (Table 4, Fig. 3) and they all share a Eurasian geographical origin. The extant T. macropoda is known to occur currently in China, Korea and Japan (Takahashi 1987). This too implies that the diversification direction of the genus may be moving from Japan towards the Himalayas.

Most Tricyrtis species occur in a limited region of the islands of Japan—Kyushu, Shikoku and Honshu (Takahashi 1987)—and there is no evidence that the plants grow in Hokkaido. The area hit by the tsunami in March 2011 is in Honshu, and that region includes habitats for the rare endemic T. nana (sect. Flavae) as well as T. affinis and T. macropoda (sect. Tricyrtis; Takahashi 1987, ). This type of natural disaster may bring about the extinction of rare species, including Tricyrtis taxa, reinforcing warnings on speciation events and possible sudden loss of biodiversity in the very dynamic and active Earth environment.

Conclusions and forward look

Bayesian MCMC analytical methods are increasingly being used to estimate posterior distributions, thereby leading to a better informed choice of evolutionary model for the study of plant systematics, molecular evolution, medicine and bioinformatics. The current study of MCC trees using endemic Tricyrtis provides a good example of how these methods can advance the understanding of molecular evolution in plants. Bayesian analyses are shown to be preferable to other clock-model-based dating methods for analysing branch lengths, topology and temporal rate estimations. Their use in the current study allowed estimations of the divergence time and molecular evolutionary rate of Tricyrtis.

The current geographical distribution patterns of Tricyrtis, together with the divergence time and rate estimation data in this study, imply allopatric speciation caused by an abiotic factor, continental drift. However, sympatric speciation may also be operating in local regions with a faster recent speciation rate. Insect pollinators of Tricyrtis may also have played important roles in co-evolution, leading to the diversification of the groups.

Different genes can evolve at different rates under different selection pressures. Because of this, other cpDNA regions will be tested in addition to the rps16 intron cpDNA used here. This will refine our estimates of evolutionary rates, divergence times and molecular dating.

Sources of funding

This work was funded by Korea Forest Service, grant 2010325 (2010–2011).

Contributions by the authors

All the authors contributed to this paper to a similar extent.

Conflict of interest statement

None declared.