Post-Triassic Spermatophyta Timetree Adding the Quaternary Radiated Asarum Wild Gingers

base Abstract Background Angiospermae radiation was known as the mid-Cretaceous event, but adaptive radiation of Asarum is also expected in the Quaternary. In order to know such the Angiospermae evolutionary history through the time, we constructed a whole Spermatophyta timetree employing BEAST v1. X associated with robust fossil calibration function.


Background
A goal of botany may be the correlation of botanical evolutionary events with the timeline of Earth history (Wilf & Escapa, 2015). For this purpose, for the rst trial, we prepared a whole Spermatophyta Bayesian inference (BI) tree constructed using the latest and the most advanced version of BEAST (v1. X;Suchard et al., 2018). This BI tree is a timetree simultaneously dated by multipoint calibration function associated with BEAST V1. X, and precisely represented the Angiosperm history since the Jurassic to the Quaternary (duration: 182.99 m.y.; Fig. 1). We included Ginkgo biloba (Ginkgoidae) and representative species of Pinales (conifers) and Cycadales for Gymnospermae, and also included Amborella trichopoda for the oldest species of paleodicots of Angiospermae.
Angiosperm was diversi ed in the Cretaceous time, and the adaptive radiation might have increased the base substitution rate (Magallon et al., 2015). They also predicted that while substantial amounts of angiosperm morphological and functional diversity had deep evolutionary roots, the extant species richness was probably acquired later (= Quaternary). Ho et al. (2005), developers of BEAST v.1.3 (Drummond et al., 2006), noted higher mutation rate after 1 to 2 Ma (Quaternary) compared to the prior lower rates. To revisit the issue of accelerated base substitution rates requires calibration with younger ages such as those of Quaternary age as well as the older aged calibration.
Genus Asarum, section Heterotropa is a suitable target for the Quaternary calibration, and we included in our BEAST analyses. Heterotropa is diversi ed in the Japan, Ryukyu, and Taiwan islands (Figs. 2 and 3; Saruma henryi is a sister of of the Asarum species as shown by Kelly (1998) and Sinn et al. (2015ab). Section Asiasarum is a sister of section Euasarum + section Heterotropa, and section Euasarum is a sister of section Heterotropa. In section Asiasarum, A. sieboldii is not differentiated among collection sites in Japan, but is differentiated relative to the Korean specimen. In section Euasarum, inter-species and intra-species differentiation is recognized, although with stolons for plant propagation (no stolon for Heterotropa).
The Chinese Heterotropa excluding A. forbesii (= section Longistylis; Sinn et al., 2015b) is a sister of A. forbesii + the Japan-Ryukyuan-Taiwan Heterotropa, and A. forbesii is a sister of the Japan-Ryukyu-Taiwan Heterotropa. Four Chinese species with rbcL data in the DDBJ / GenBank are included in our analyses, but an additional 14 Chinese species were also included in section Longistylis clade based on solely ITS and matK data.
Endemic species of section Heterotropa were severely differentiated on the Japan-Ryukyu-Taiwan islands and the divergence time was considered to be the geological event of 1.55 Ma. The differentiation tends to occur within islands, and therein constitutes distinct clades. The species of each clade tends to have common characteristics of calyx and leaf ( Fig. 2; Table 1). However, sympatric occurrence of genetically distinct species within a single island without a sister relationship, differentiation into a distinct clade on the same island, and inter-island sister relationships were found. Viola (Malpighiales; rodids) has crown node of 15.49 Ma. Endemic Viola okinawensis is a sister of widely distributed V. grypoceras. Viola amamiana, Amami Oshima, is a sister of V. tashiroi, Iriomote-jima. Viola yedoensis is differentiated in each island of Ryukyu, Tanwan, and HongKong. Viola betonicifolia is a similar species with and a sister of sympatric V. yedoensis. Sister relation was found for additional calibrated species by 1.55 ± 0.15 Ma for nodes of A5 of Lilium (monocots) and A6 of Pinus (conifers).

Exponential increase in base substitution rate
The base substitution rate has exponentially increased since 10 Ma (Miocene; Neogene), and the heavy red trendline and equation was shown on inset in Fig. 1. All the Heterotropa rates were less than 0.01 as shown in the lower rates on right hand gure of Heterotropa (Fig. 1), and the 1.55 Ma calibration for Heterotropa did not affect the increasing rata in the Quaternary time.
We also drew a "trendline" by smoothly connecting plots with maximum rate through the time (light blue curve on inset in Fig. 1). This curve shows another peak of rate around 130 Ma (mid-Cretaceous).

Discussion
Higher level phylogeny, re ecting mid-Cretaceous Angiospermae radiation Spermatophyta (Gymnospermae + Angiospermae) phylogeny comparable to ours was presented by Rydin et al. (2002). We estimated the crown age of Ginkgo biloba (Ginkgoales) at 166.63 Ma, suggesting actually a living fossil.
They showed order-family stem ages were between 140 and 120 Ma applying their calibrations comparable to ours. According to their statements: "Eudicots are morphologically characterized by tricolpate pollen (or derived from this condition; Walker & Walker, 1984;Donoghue & Doyle, 1989), which probably evolved on the stem lineage of this clade. Tricolpate pollen grains are morphologically distinctive, can easily become preserved as fossils, and unequivocally indicate membership to a single clade, thus providing an exceptionally good calibration. The eudicots were estimated as being nearly contemporaneous with paleodicots and monocots; therefore, the major components of angiosperm extant diversity began to diversify in the middle Cretaceous." Note that Ceratophyllales expressed the basal eudicots ( Fig. 1) is characterized by tricolpate.
In our timetree of Fig. 1 as noted above, crown age of Angiospermae was estimated at 132. 44 Ma,paleodicots;131.53 Ma,monocots at 116.68 Ma,and eudicots: 120.31 Ma, and after the almost contemporaneous differentiation, the order level diversi cation was followed until Paleogene (red diamonds in Fig. 1), supporting Magallon et al. (2015) comments.
Family-species level phylogeny, re ecting the relatively recent Angiospermae radiation Ran et al. (2018) showed that Pinales (Gymnospermae) consists of two clades of Cupressaceae and Pinidae. Nagalingum et al. (2011) made timetree of Cycadales consisting of the Cycadaceae and Zamiaceae clades, and pointed out that Cycadales is a living fossil but each speciation was a relatively young event in the Neogene. We estimated genus level differentiation of Cycadaceae after 39.56 Ma and that of Zamiaceae after 34.35 Ma (Paleogene). In a case of Pinales, genus level differentiation of Cupressaceae was after 26.32 Ma, and that of Pinaceae was after 21.79 Ma.
Smith & Brown (2018) made Spermatophyta large dated tree following divergence-time estimate by Magallon et al. (2015). They showed that number of the lineage was increased since ca. 130 Ma, and it also exponentially increased toward the Recent since ca. 20 Ma, in their Fig. 2, inset b.
Yellow diamond in Fig. 1 shows that family level differentiation was since the Paleogene but continued to the Neogene, followed by genus-species level differentiation including extensive Asarum radiation mentioned below.
Allopatric speciation is expected to have acted on allopatric island populations of Heterotropa, but these endemic species are not always allopatric but partly sympatric within an island, such as on the Amami Oshima islands (Maeda, 2013). Amami Oshima yields many Heterotropa species, and we analyzed six endemic species. All these differentiated Amami species splitting into three distinct clades ( Fig. 1 right hand enlarged gure; red colored; Amami) may not have been generated solely by the allopatric speciation (cf., Osozawa et al., 2014Osozawa et al., , 2015ab, 2016Osozawa & Wakabayashi, 2015).
The geological event that formed isolated islands of Ryukyu at 1.55 Ma is also tied to climatic and other environmental changes. The isolation of Ryukyu islands from continental Asia also prevented detrital sediment from the Yangtze and Yellow rivers from reaching these locations, so coral reefs began to form around each island in the absence of this detrital input. The separation of the Ryukyu islands from the mainland also resulted in the in ow of the Kuroshio current through the trough (Fig. 3), and this current, thus rerouted, owed into and warmed the present East China Sea (Osozawa et al., 2012). The rifting event resulted in a sea or seaways where the mainland once was, and this affected storm patterns. Typhoon tracks once passed over land and lost energy, but after the rifting remained over the sea and maintained their strength further along their course (Osozawa & Wakabayashi, 2015). In the case of the Honshu of the Japan islands, the Tsushima warm current, a branch of the Kuroshio current, began to ow into the Japan Sea as a result of the rifting event (Fig. 3). This current warmed the seawater during the winter, and the consequent wet, northwestern seasonal wind caused heavy snowfall along the Japan Sea coast ( Fig. 3; Osozawa et al., 2012).
As a consequence, the Heterotropa endemic speciation and severe diversity in the Japan-Ryukyu-Taiwan islands (= the simultaneous genetic and morphological expansion; Figs. 1 ~ 3) might have rather resulted from the mechanism of adaptive radiation from a single ancestral Asarum inhabiting the eastern margin of original Asia continent (Fig. 3 left), similar to that observed for the Hawaiian silversword alliance (Baldwin and Sanderson, 1998;Blonder et al., 2015; cf., Alpine plants in Boucher et al., 2015).
On Amami Oshima, adaptive radiation may be in progress (Matsuda et al., 2017). A. fudsinoi there has many morphotypes of calyx and leaf as if distinct species, and A. fudsinoi is genetically similar to other A. gusk, A. celsum, and A. pellucidum with distinct morphology (Fig. 2 right center).
Similar to Fagus crenata ( Fig. 1; although the MRCA age is estimated at 20.39 Ma), A. megacalyx might have adapted to a heavy snowy environment and radiated in northern Honshu (Fig. 3). Hiura (1978) originally proposed the idea of adaptation to a snowy climate, and this included the premise that the increased size of calyx during the several month growth period indicated adaptation to the heavy snow accumulation.
Asarum forbesii is a sister of the island Heterotropa, and it differentiated at 1.91 Ma; these have a sister relation to the other Chinese Heterotropa (section Longistylis) differentiated at 3.72 Ma ( Fig. 1; Longistylis species were also radiated since 2.63 Ma). At 3.72 Ma, a single common ancestor of A. forbesii inhabited along the Yangze river and the Heterotropa colonized the eastern area of the Asian continent partly overlapped in habitat (Fig. 3), and started to differentiate at 1.91 Ma (Fig. 1). This date may be related to start of glacier and inter-glacier climatic change of the Quaternary time since 2.58 Ma.

Triggers for the Recent and mid-Cretaceous increases in base substitution rate
We noted that Spermatophyta differentiation was recognized within the middle Cretaceous and the Quaternary (Fig. 1) as pointed out by Magallon et al. (2015) and suggested by Smith & Brown (2018). We discuss here such the differentiation that we also noted above was tied to the increasing base substitution rate ( Fig. 1 inset), and the basic factors.
A strict clock model assumes that every branch evolves according to the same evolutionary rate (Drummond et al., 2006), but the rate and branch length is clearly not proportional, and thus we should apply relaxed clock model . We can now discuss the variable rate, through building a credible Spermatophyta timetree by fully applying multiple point calibrations and partition functions of BEAST v1. X.
The molecular clock hypothesis states that DNA sequences evolve at a relatively constant rate over time even in a case of relaxed molecular clock, but the absolute clock needs to be calibrated (Ho, 2008). We showed that the base substitution rate of Spermatophyta has increased exponentially toward the Recent as shown by the heavy red trendline (Fig. 1 inset). This phenomenon had been shown for taxa such as primates (Ho et al., 2005(Ho et al., , 2011 and insects (Papadopoulo et al., 2010), but generated little attention. BEAST v1. X is possible to run simultaneously applying multiple calibration points, but BEAST v.1.3 (Ho et al., 2005) and BEAST v.1.4 (Papadopoulo et al., 2010) needed to run repeatedly by applying a date at every calibration point, and a rate at every run is assumed to be constant at each node (like as strict clock model, although relaxed clock model was applied). They found that a Quaternary date calibration produced a rapid rate. Note that we instantaneously calibrated by multi dates including older dates such as Jurassic as well as younger dates of Quaternary of 1.55 Ma, not solely calibrated by the Quaternary date, and the reason of increasing rate is unrelated only to the Quaternary calibration (Asarum rates were estimated less than 0.01; Fig. 1 right hand enlarged gure). Also note that combined gene analyses were not possible in the old versions of BEAST above, and such the users needed to discuss the rate for every each gene or concatenated gene. Therefore, although their trendlines and the equations are similar to ours, their timetrees are actually incorrect by not re ecting drastic changing of base substitution rates through the time, but by assuming a constant rate as if strict clock.
We justi ed that the recently increasing base substitution rate is an actual event. A possible trigger of increasing rate may have been the start of glacial and interglacial cycle and the severe environment change in the Quaternary time. As noted above, this might be a factor in Heterotropa radiation, although physical isolation of islands was also acted on. Feedback from biologic development might have also played a role in triggering the onset on glaciations. The Quaternary glaciations may have been triggered by expansion of land grasses prevailing Poales, and the late Paleozoic glaciations (Montañez et al., 2013) might have been triggered by the development of terrestrial tree ferns to form thick coal layers, because these processes increased carbon xation that consequently decreased atmospheric CO2 concentration (Taira, 2007). C4 plants are e cient for such CO2 xation (Saga, 2004), and C4 Poales plants appeared and started diversi cation from 15.58 Ma ( Fig. 1; C4 Amaranthaceae, eudicots, from 25.93 Ma). Carbon isotope ratios from mammalian fossil tooth enamel showed that diet change into C4 plants started at 9.9 Ma in eastern Africa (Uno et al., 2011). Also isotopic analysis of the mammalian fossil tooth from the sub-Himalayan Siwalik Group, Pakistan, showed that C4 savannah replaced C3 forest and woodland between 8.5 and 6.0 Ma (Badgley et al., 2008). The expansion of C4 grasses was a global phenomena including North America and South America beginning in the late Miocene and persisting to the present day, concerning the present glacial-interglacial period (Cerling et al., 1997). Note that increasing rate started at 10 Ma, not at 2.58 Ma of the Quaternary (Fig. 1 inset). Nagalingum et al. (2011) showed that cycads underwent a near synchronous global diversi cation beginning in the Miocene. Our result revised the estimate for the initiation of this genus level diversi cation to 39.56 Ma (Oligocene; Fig. 1). Such a date coincided with a large and rapid drop in atmospheric CO2 and the onset of cooling also at the Eocene to Oligocene transition (33.9 Ma; c.f., Pound & Salzmann, 2017). Note that conifer radiation was also estimated at 37.99 Ma (Eocene; Fig. 1).
Fossil record indicates a dramatic increase in phylogenetic diversity and ecological abundance of angiosperms around the mid-Cretaceous (Friis et al., 2006), and this event including noted above might be tied to another rate peak around 130 Ma as shown by the light blue trendline (Fig. 1 inset), which might express a start of order level radiation of Angiospermae (Fig. 1). Thereafter, both the paleodicots and monocots were extensively differentiated in order level, and then family level, but for eudicots until the Neogene (Fig. 1). Insect beetles have been clari ed to have a rate peak also around middle Cretaceous time (Hunt et al., 2007;McKenna et al., 2015;Gunter et al., 2016;Toussaint et al., 2016;Zhang et al., 2018), and the increasing rate should be tied to the co-radiation of Angiospermae as the food plant. Also broad-leaved angiosperm plant species, especially dicots, initiated a global ecological transformation towards the Cretaceous biodiversity (Jan de Boer et al., 2012).

Materials And Methods
Taxon sampling (Table 1) Amborella trichopoda, New Caledonia origin, was offered by the Koishikawa Botanical Garden, the University of Tokyo. Some Angiospermae specimens of which sequence data were unavailable in DDBJ / GenBank were collected from Experimental Station for Medical Plant Studies, Graduate School of Pharmaceutical Sciences, Tohoku University (Table 1). Endemic Viola and others were collected (Table 1), collaborated the Asarum sampling below.
The genus Asarum consists of three sections of Heterotropa (evergreen, main section), Euasarum (deciduous, stolon; 16 species in Kelly, 1998;previous Asarum;renamed by Sinn et al., 2015b), and Asiasarum (deciduous; not known from the Ryukyu and Taiwan islands; two species). The phylogenetic relationship of these sections, focused on section Euasarum, was studied by Kelly (1998), Sugawara et al. (2005), and Sinn et al. (2015ab). The phylogeny of section Asiasarum, distributed in East Asia, including the Japan islands, was studied and subdivided into species by Yamaji et al. (2007). Section Heterotropa has not been analyzed over the entirety of the Japan, Ryukyu, and Taiwan islands (Fig. 3) except for a recently published paper of Takahashi & Setoguchi (2017). We submitted our Heterotropa sequence data to DDBJ / GenBank in 2015, and they were released in 2019 (Table 1).
On the basis of ower structure, the pollinators of Heterotropa are probably tiny ies and millipedes (Hiura, 1978;Sinn et al., 2015b), and seeds may be carried/transported by ants (Maeda, 2013;Sinn et al., 2015b). The dispersal ability of Heterotropa may be small (Sugawara, 2006), and this may be a factor in generating endemic Heterotropa species (Takahashi & Setoguchi, 2017). Note, however, that some islands, such as Amami Oshima island, yield plural endemic species (Fig. 2; Table 1), and these are mostly sympatrically distributed within such an island including Amami Oshima (Maeda, 2013). Therefore, allopatric speciation, proposed by Takahashi & Setoguchi (2017), cannot be simply applied to explain the extraordinary diversi cation of Heterotropa in the Japan-Ryukyu-Taiwan islands.
Collected leaf samples were stored in a desiccator with silica gel lling the space under the platform. Voucher Asarum specimens were also prepared. Specimen numbers (= isolate), localities (= country), species, and collectors are registered in the DDBJ / GenBank. These accession numbers, and additional information on the character of calyx, leaf, and stem are shown in Table 1. Some representative calyx photos are shown in Fig. 2. Many endemic Asarum species are also known from China, and the Chinese sequence data have recently become available in the DDBJ / GenBank. DNA extraction, polymerase chain reaction ampli cation, and sequence alignment A problem applying molecular phylogenetic method to plants was the lack of suitable primers to amplify recognizable variation (Kanno et al., 2004;Fazekas et al., 2008). However, the chloroplastic matK and rbcL genes were recommended for use as markers by the CBOL Plant Working Group (2009). In addition, the nuclear ITS gene is known to be an additional suitable marker (Asarum: Kelly, 1996;Sugawara et al., 2005;Yamaji et al., 2007;China Plant BOL Group, 2011;Sinn et al., 2015ab;Takahashi & Setoguchi, 2017), although the APG (APG, 1998;APG II, 2003;APG III, 2009; APG IV, 2016) did not include the ITS data in their analyses. We analyzed by the combining these genes (not by the concatenated data set) in BEAST v1. X also as the rst trial, and found that the resolution was much improved. DNA extraction was done using the GenElute™ Mammalian Genomic DNA Miniprep Kit by Sigma-Aldrich. Before this operation, crushed young leaf specimens were cleaned with 1% Tris and EDTA solution in 1.5 mL tube.
Ampli cation was done by GoTaq G2 Green Master Mix, Promega, and the temperature of incubation was 94˚C for 60 seconds, with denaturation at 94˚C for 30 seconds, annealing at 53˚C for ITS (56˚C for matK; 54˚C for rbcL) for 60 seconds, extension at 72˚C for 60 seconds, cycled 35 times, and nal extension was at 72˚C for 5 minutes.
The PCR product was puri ed using Wizard SV Gel and PCR Clean-Up System, Promega. Sequencing was done by Macrogen Japan.
Sequence alignment was done using ClustalW incorporated in MEGA X (Kumar et al., 2016). No gap was observed in aligned sequences for Asarum, although gapes were found in ITS region for the resting sequences even for Aristolochia species relative to the Asarum sequences. The ITS gene (694 bp), the chloroplastic matK gene (833 bp), and the chloroplastic rbcL gene (544 bp) were readable for Asarum and the number of base pair were similar to the resting. Codon translation for the matK and rbcL sequences was checked using the ExPASy-Translate tool, Bioinformatics Resource Tool. These sequences were also checked by BLAST (Basic Local Alignment Search Tool) offered by DDBJ, and such data are re ected in our registered data in the DDBJ / GenBank.

Phylogenetic analyses associated with fossil + geological event calibration by BEAST
A Bayesian inference (BI) tree ( Fig. 1; Heterotropa tree extends to the right in this diagram) was constructed using the software BEAST v1. X released on 10th June 2018 (Bayesian Evolutionary Analysis Sampling Trees; Suchard et al., 2018), running BEAUti (Bayesian Evolutionary Analysis Utility), BEAST, TreeAnnotator, and FigTree, in ascending order. Before operating the BEAST software, the BEAGLE Library must be downloaded.
For graphic explanation of the operation of this software, see the "BEAST operating manual" at: http://kawaosombgi.livedoor.blog/ We surveyed literatures including fossil calibration reviews such as found in Smith et al. (2010), and employed chronologically reliable calibration dates bellow. Some are based on radio-isotopic dating the fossil-bearing strata, whereas others are based on biostratigraphy assigned to an age/stage on the geologic time scale, for which absolute age ranges are generally based on radio-isotopic dates of associated strata in key global localities (Wilf & Escapa, 2015). This time scale has been standardized by International Commission on Stratigraphy (ICS) (www.stratigraphy.org) and the most recent version of the time scale is available at http://www.stratigraphy.org; as of this writing the most recent version available is v v2020/03, and the explanatory paper related to the generation of the time scale is Cohen et al. (2013).
Calibration point A: Gonkgo biloba is often described as a living fossil and an extant genus, and Gonkgo fossil wood was reported from the Liaoning province, northern China (Jiang et al, 2016). Chang et al. (2009Chang et al. ( , 2014 reported Ar-Ar ages of 160.7 ± 0.4 Ma and 166.7 ± 1.0 Ma, respectively, and we adopted the latter as the prior input and stem age (Fig. 1).
Calibration point C: Modern but also fossil Ceratozamia was reported from the European Oligocene basinal strata (Kvacek, 2014;Kovar-Eder, 2016), and corresponds to calcareous nannoplankton zone NP23 (31.8 ± 2.2 Ma) and planktonic foraminifera zone P18 (32.9 ± 0.9 Ma; Fig. 1). Fossil localities of other modern cycads such as Australia are not well constrained by geochronology. Note that extinct cycad species such as Nilssonia cannot be adopted for calibration dates.
Calibration point D1: Archaefructus from the Jehol Biota, northeast China, is an extinct but earliest known genus of Angiospermae. The lower Archaefructus fossil horizon within the Yixian Formation (Sun et al., 2011) was dated by the Ar-Ar method applied to intercalated silicic tuff, and the date was 130.7 ± 1.4 Ma (He et al. 2006). We applied this date for the crown age of paleodicots (monophyletic), rather than the whole Angiospermae. Calibration point D2: No fossil Amborella has not been found, but this species is the oldest lineage of paleodicots in the APG system, and the stem age was applied also at 130.7 ± 1.4 Ma.
Calibration point E: Leefructus mirus also from from the Jehol Biota is included in a basal eudicot family of Ranunculaceae (Sun et al., 2011). The fossil horizon within the middle Yixian Formation (Sun et al., 2011) was dated by the Ar-Ar method applied to sanidine included in intercalated silicic tuff, and the date was 124.60 ± 0.25 Ma ).
Calibration point F: Palynological data of the Early Cretaceous continental sequences of western Portugal revealed that monocot radiation was in Aptian time (119 ± 6 Ma) (Hochuli et al., 2006). They stated that it preceded the radiation of dicots by at least 10 m.y., but it con icts to the above calibration of D and E and the APG system, and we do not follow this statement.
In BEAUti, the following software settings were used. Partitions: fasta les were converted into nexus les by ClustalW offered by DDBJ, and the loading was by using the Import Data or plus button. Partitions de ned by the ITS, matK, and rbcL gene sequences appeared in the Partition box.
Taxa: Loading of taxa as ingroup was by using the plus button. The left screen: Taxon Set (monophyletic boxes were checked for all, and stem box were checked in case by case), and the right screen: Included Taxa. As output in Fig. 3, calibration dates were set in Priors. Clocks: Clock Type: Uncorrected relaxed clock, Relaxed Distribution: Lognormal. Uncorrelated relaxed clocks allow each branch of a phylogenetic tree to have its own evolutionary rate under log-normal distribution, and the node rate is the rate median of three branches (Drummond et al., 2006).
Priors: tmrca (time of MRCA) was input from the calibration point date noted above as Prior Distribution: Normal, and the Mean and Standard deviation.
Running BEAST was done by incorporating xml input le made by BEAUti. The consequent tree was drawn by FigTree v1.4.2, for that, the tree les were input into TreeAnnotator. The 95% highest posterior density for con dence intervals of ages can be output in FigTree, but not shown in Fig. 1 to avoid confusion. In FigTree, posterior probability ("posterior"), posterior age ("Node ages"), and "rate median" (not constant) can be output, and these are shown at each node in Fig. 1. This function was not fully used in any previous paper, and we found in this paper the inconsistent rates through the time as suspected by the relaxed clock model of BEAST (Drummond et al., 2012). Consequently, we made base substitution rate ("rate median" shown at each node in FigTree) vs age ("Node age" shown at each node in FigTree) diagram (Fig. 1 inset)  Monte Carlo methods; ITS:internal transcribed spacer; matK:chloroplastic maturase K; rbcL:ribulose-1, 5-bisphosphate carboxylase/ oxygenase large subunit; DDBJ:DNA Data Bank of Japan; APG system:Angiosperm Phylogeny Group system.

Declarations Competing interests
The authors declare that they have no competing interests.
Authors' contributions SO collected samples and coordinated the research, carried out the DNA analyses, and wrote this paper, CN helped to collect specimen from Ryukyu (nomenclator of Viola okinawensis), and JW (native Englishspeaking American geologist) edited the writing.
All authors have read and approved the manuscript.

Availability of data and materials
The sequence data of ITS, matK, and rbcL are available in GenBank / DDBJ, and the accession numbers are in Table 1. Blank: Uploaded and released after acceptance.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.