Cicada minimum age tree: Cryptic speciation and exponentially increasing base substitution rates in recent geologic time [version 1; peer review: awaiting peer review]

We developed a new time-calibrated tree incorporating primarily endemic along with some cryptic Ryukyu islands cicada data, following the recent publication of global cicada data by Marshall et al. (2018), Łukasik et al. (2018), Simon et al. (2019), Price et al. (2019), and Hill et al . (2021).  A total of 352 specimens were analyzed using BEAST v1. X software with a relaxed clock model. Fossil calibrations as old as Triassic were adopted largely following Johnson et al. (2018) and Moulds (2018), and a Quaternary geological event calibration was adopted following Osozawa et al. (2012, 2021b) and input into BEAST v1. X. Our timetree suggests that Tettigarctidae had a cicada basal lineage as old as 200.63 Ma, with Derotettiginae the next oldest lineage at 99.2 Ma. Tibicininae is a sister of the remaining subfamilies of Tettigomyiinae, Cicadettinae, and Cicadidae, and their species level differentiation and radiation began at 40.57 Ma. The Cicadinae clade consists of specific tribes with parapheletic relationship, and the vicariance and adaptive radiation generated many cryptic species in each tribe. We estimated base substitution rate as a function of age, and the result strongly indicates an exponential increase of base substitution rate in recent geologic time. The consequent increase in cicada biodiversity, including generation of cryptic species in the Ryukyu Islands and surroundings, may have been driven by the generation and spreading of C4 grasses and coeval Quaternary climate change.


Introduction
A phylogenetic tree of worldwide cicada was recently constructed by Marshall et al. (2018) and Simon et al. (2019) applying five concatenated sequences of mitochondrial COI and COII, and nuclear ARD1, EF-1a, and 18S rRNA, and by Łukasik et al. (2019) applying whole mitochondrial sequences for representative species in Marshall et al. (2018), and family level phylogenetic relation has been clarified. Although Tettigarctinae is an old diverged lineage and Derotettiginae may be next, their worldwide phylogenetic trees were not dated trees. Price et al. (2019;restricted to Platypleurini) and Hill et al. (2021; restricted to Asian Cicadinae) built partial (not worldwide) dated trees using BEAST v.2.5 (Bouckaert et al. 2014) applying COI and other sequences, but much of global cicada evolution has not been tied to absolute time.
The latest version of BEAST (Bayesian Evolutionary Analysis Sampling Trees v1. X; v1.10.4 2021; Suchard et al. 2018) released on 10 June 2018 has a clear and simple age calibration protocol and function, updated from BEAST v.1. 7 (v1. X ≒ v.1. 8). This calibration involves applying times of the most recent common ancestors (tMRCAs) of the ingroup species, i.e., applying the node age of a specific clade as a minimum age, in the associated software of BEAUti (Bayesian Evolutionary Analysis Utility; BEAST is the platform software). The maximum age constraint normally considered in MCMCtree (4.9e 2017; Yang 2007) was not clearly defined (Benton & Donoghue 2007;Marshall 2008;Hill et al. 2021), and simply handled by ignoring the maximum age in BEAST v.1. X calibration (Osozawa & Wakabayashi 2021;Osozawa et al. 2021a). We sought the oldest fossil of the corresponding node of specific clade with an assumption that the oldest fossil age was equivalent to the minimum age and equivalent to "tMRCA" in BEAST v1. X. Moulds (2018) reviewed the ages of cicada fossils. These redefined ages, ranging from 16.45 AE 0.45 Ma to 244.5 AE 2.5 Ma, were available for our fossil-based time calibrations in BEAST v1. X. Klopfstein (2021) suggested that recent node dating approaches including Misof et al. (2014) and Montagna et al. (2019) have a credibility problem: different studies using the same molecular data and even the same sets of fossils regularly arrive at drastically different age estimates. She showed that a major reason for these differences is well known: even welldated and firmly placed fossils can only provide a minimum age for a particular node. Therefore our fossil calibration applying solely minimum age (= tMRCA) was credible.
As shown by Osozawa et al. (2017a), Platypleura and some other endemic cicadas in Ryukyu Islands can be rigidly calibrated by a geological event calibration at 1.55 AE 0.15 Ma (Quaternary; Osozawa et al. 2012). As shown by Osozawa et al. (2021a), Meimuna opalifera and some other endemic cicadas on Hachijo-jima, Izu-Bonin islands, can be calibrated by a geological event calibration of emergent age at 0.24 Ma (Quaternary; Osozawa et al. 2021b).
Through these analyses, we corroborated the classification and some rearrangement of species into four subfamilies of Tibicininae, Tettigomyiinae, Cicadettinae, and Cicadinae included in a family Cicadidae by Marshall et al. (2018) and Łukasik et al. (2018), and then estimated the splitting dates of these subfamilies, tribes (especially Cicadinae tribes after Hill et al. 2021), and species (Figures 1-3). In the BEAST analyses, we included Derotettix, a relict species of new subfamily Derotettiginae with the oldest lineage in family Cicadidae (Simon et al. 2019), and attempted to estimate the crown age ( Figure 2). Comparison to the entire Hemipteroid insect timetree (Johnson et al., 2018) and entire insect timetree (Misof et al. 2014) could be conducted as an extension of this analysis, by adding other Hemiptera species as outgroup (Figures 1-3).
Our primary goal was to present the precise evolutionary history of all cicadas by constructing the BEAST timetree, and also taxonomic reconsiderations for Cicadinae tribes after Hill et al. (2021) and for Ryukyu endemic cicadas. Another BEAST v1. X function facilitates additional evaluation of the time variability of base substitution rates. Recent dating analyses employ a relaxed clock model, which allows each branch of a phylogenetic tree to have its own evolutionary rate (Drummond et al. 2012). Although the relaxed distribution can be set to lognormal in BEAUti, the rate of variability has not been documented prior to this study. The output figure of BEAST v1. X presents the base substitution rate and age at each node, and shows the acceleration of base substitution rates through the time.

Ethical approval
The present study did not concerned invertebrate experiments and did not involve endangered or protected species. We obtained permission of collection in the Taroko National Park, Taiwan, from the director (No. 0990012881;August 1~11, 2010), with a help of Bor-ming Jahn, and permission of collection in the Tokara islands, from the Toshima village headman, from August 29~September 8, 2010. Collection in the Ryukyu islands was before the designation of National Park since 2016. No specific permission was required outside the national parks and private areas. Marshall et al. (2018), Simon et al. (2019), Łukasik et al. (2019) included comparatively few Asian cicada species in their analyses. We have previously published 70 isolate data from Platypleura primarily from the Japan, Ryukyu, and Taiwan islands (Osozawa et al. 2017a; our aim was the vicarince acted on each island population started at 1.55 Ma and the cryptic speciation), and 21 of these data were used in the present analyses by excluding duplicated sequence data. We also collected and analyzed cicada specimens, adding isolate data from 92 specimens. Accordingly, our own data total Figure 1. Simplified cicada timetree built by BEAST v1.X, applying a 1,534 bp in maximum COI sequence. Inserted figure: Base substitution rate (= rate median shown at each node; substitutions per site per million year; s/s/ myr) vs age (= posterior age shown at each node) diagram. Red approximate curve with its formula was drawn by an Excel function, with the intersection for the curve = 0.0128 s/s/myr, the rate median shown on Tracer. 21 + 92 = 113 specimens (Table 1). Note that we collected all the 35 species from Japan including the Ryukyu Islands, but excepting severely protected Platypleura albivannata (see Osozawa et al. 2017a; may be extinct without DNA sequence data) and Meimuna boninensis (see Osozawa et al. 2021a).

Taxon sampling
We incorporated representative sequence data from the GenBank/DDJB. This is because Tettigarctinae, Derotettiginae, Tibicininae, and Tettigomyiinae are not known from East Asia, and Cicadettinae has only two species of Kosemia in the Japan main islands. Thus to extend our analyses beyond East Asia, the Marshall et al. (2018) and Łukasik et al. (2019) data were essential for us. We combined our data from 113 East Asian specimens with data from 75 specimens from the studies of Marshall et al. (2018; including 20 East Asian specimens), and Łukasik et al. (2018;including 27 East Asia specimens). In addition, we incorporated data of 15 Platypleurini (other than Platypleura) from Price et al. (2019), and data of 149 Asian Cicadinae from recently published Hill et al. (2021). Accordingly we analyzed sequence data from 113 + 75 + 15 + 149 = 352 specimens.
Platypleura cicada (Osozawa et al. 2017a) experienced vicariance triggered by the 1.55 AE 0.15 Ma isolation of the Ryukyu, Japan, and Taiwan islands from Chinese continent (Osozawa et al. 2012), and we collected specimens from each island population for each Platypleura species. Similarly, we collected cicada specimens for the present analyses from each island population of Mogannia (Cicadettinae), and Cryptotympana, Graptopsaltria, Hyalessa, Pomponia, Figure 2. Cicada timetree built by BEAST v1.X, applying 1,534 bp COI sequence. OUTs with isolate number: our own analyzed specimens shown in Table 1, and others: from GenBank/DDJB. In outgroup Hemiptera; #: analyzed family by Johnson et al. (2018); % analyzed family by Misof et al. (2014). Inserted figure: Base substitution rate (= rate median shown at each node; substitutions per site per million year; s/s/myr) vs age (= posterior age shown at each node) diagram. Red approximate curve with its formula was drawn by Excel function, with the intersection for the curve = 0.0128 s/s/myr, the rate median shown on Tracer.
Meimuna, Tanna, and Euterpnosia (Cicadinae). Hyalessa maculaticollis was known to be affected by vicariance within China and Japan (Liu et al. 2018). Our 113 East Asian specimens consist primarily of these endemic and cryptic species inhabiting Japan, Taiwan, and the Ryukyu islands.
DNA sequence COI and 18S rRNA sequence data from our collected 113 isolates, including Platypleura in Osozawa et al. (2017a), are shown in Table 1. Primers used, amplifications, and sequencing are given in Osozawa et al. (2017a). These sequence data were aligned by ClustalW in MEGA 5 (Tamura et al. 2011). The COI sequence data comprise 1,534 bp, and the 18S rRNA sequence 874 bp, with high enough resolution to construct a phylogenetic tree, as we showed previous analyses of Platypleura (Osozawa et al. 2017a). We did not analyze calmodulin in Osozawa et al. (2017a), because the resolution was insufficient. The COI data in Marshall et al. (2018) comprised 1,485 pb, comparable to ours. The COI data in Price et al. (2019) comprised 940 bp, andHill et al. (2021) comprised 648 bp, comparable to ours, so we incorporated these data into our present analyses. Nuclear 18S rRNA shows less variation with much slower base substitution rate compared to mitochondrial COI (Osozawa et al. 2017a; COI: 0.0270 substitutions/site/myr, 18S rRNA: 0.000492 s/s/myr; strict clock model; solely calibrated by 1.55 AE 0.15 Ma following Osozawa et al. 2012). The tree topology was unaffected by 18S rRNA (Osozawa et al. 2017a), but 18S rRNA was included in the analyses in this paper.  Table 1, and others: from GenBank/DDJB. In outgroup Hemiptera; #: analyzed family by Johnson et al. (2018); % analyzed family by Misof et al. (2014). Inserted figure: Base substitution rate (= rate median shown at each node; substitutions per site per million year; s/s/myr) vs age (= posterior age shown at each node) diagram. Red approximate curve with its formula was drawn by Excel function, with the intersection for the curve = 0.0114 s/s/myr, the rate median shown on Tracer. Note that this rate is a little slower than that solely of COI in Figures 1 and 2, reflecting slower rate of 18S rRNA than COI (see Osozawa et al. 2017a). We used COI and 18S rRNA sequence data, from 352 total specimens (239 from GenBank/DDBJ + 113 of our own) for the COI timetree in Figures 1 and 2, and 155 total specimens (42 from GenBank/DDBJ + 113 of our own) for the COI +18S rRNA timetree in Figure 3. The COI and 18S rRNA data in table 1 in Marshall et al. (2018) contain missing and incomparable data, so some of their GenBank/DDBJ data were not applicable for our analyses. Whole mitochondrial sequence data by Łukasik et al. (2019) are included in our analyses as corresponding COI regions. Within COI sequence data in Marshall et al. (2018), 21 data for Cicadettinae and 14 data for Cicadinae were incorporated into our analyses. 18S rRNA sequence in Marshall et al. (2018) was used for only Nablistes heterochroma (Tettigomyiinae) and Platypedia putnami (Tibicininae). Only the COI sequence data in Price et al. (2019) for Platypleurini and in Hill et al. (2021) for Asian Cicadinae were applied to our study.
North American Cryptotympanini were analyzed by Hill et al. (2015), applying 1,467 bp of COI and 783 bp of nuclear EF-1a with sufficient resolution. Cicadettini, primarily from Australia, was analyzed by Marshall et al. (2016), applying 1,492 bp of COI and 1,047 bp of nuclear EF-1a also with sufficient resolution. Some of these COI sequence data were included in our analyses.
For our initial analysis, we constructed a minimum age tree solely applying COI sequence data (Figures 1 and 2; 352 specimens) that covers Tettigomyiinae and Tibicininae species. Following this analysis, we constructed a minimum age tree by applying both COI and 18S rRNA sequences ( Figure 3; 155 specimens, i.e., 352 À 155 = 197 specimens lack 18S rRNA sequences). These analyses showed that topology and ages associated with the analyses were not impacted by inclusion or exclusion of 18S rRNA sequence data.
Why was BEAST2 not used ? Price et al. (2019) and Hill et al. (2021) employed BEAST v.2.5 (Bouckaert et al., 2014). We employed BEAST v.2.5 in Osozawa et al. (2016), but thereafter we changed to BEAST1. The calibration function in BEAUti in BEAST2 is similar to BEAST1 (see next section bellow), but note that "Taxa" in BEAST v1. X was combined with "Priors" in BEAST v.2.5. However, In BEAST2, applied sequence data should be concatenated as done in Price et al. (2019) and Hill et al. (2015Hill et al. ( , 2021. Otherwise, as in Osozawa et al. (2016), the calculation in BEAST2 for each applied sequence data (e.g., COI and nuclear 16S rRNA) is done separately (in this example, two BEAST runs are needed). The consequent two tree files must be combined into one by LogCombiner, but in the combined tree, branches are folded reflecting the discordant topology made by e.g., mitochondrial COI and nuclear 16S rRNA sequence data. BEAST v1. X is recommended to do combined gene analyses, and in general, concatenated gene analyses are not recommended.
The most recent BEAST 2.6 released in May 2019 did not change the protocols. A misunderstanding about tip dating is as follows (https://beast.community/first_tutorial). The "Tip Dates" panel in BEAST2 or the "Tips" panel in BEAST1 allows giving the taxa dates. However, taxon dates (or 'tip dates') are only important in certain cases, e.g., when sampled from rapidly evolving viruses or ancient fossil DNA material. In the case of the cicadas, we are analyzing a tree that represents millions of years of evolution, so the dates of the tips can be assumed to be zero. This should be the default, i.e., the taxa all have a date of zero and the "Use tip dates box" should not be selected.
Phylogenetic analyses by BEAST v1. X A Bayesian inference (BI) tree (Figures 1-3) was constructed using the software BEAST v1. X, running BEAUti, BEAST, TreeAnnotator, and FigTree, in ascending order. Before operating the BEAST software, the BEAGLE Library must be downloaded. Tracer v.1.6 was applied for checking the calculation status and estimating the median base substitution rate.
For graphic explanation of the operation of this software, see Osozawa (2021a; BEAST v1.X tutorial, in a case of four cicada genera) at: dx.doi.org/10.17504/protocols.io.bq6mmzc6.
In BEAUti, the following software settings were used (Appling Appendix BEAUti file, readers may run the platform software BEAST and check the protocol and reliability).
Partitions: Loading fasta files was by using the Import Data or plus button. Partitions defined by the COI and 16S rRNA gene sequences appeared in the Partition box (For Figure 3; COI file only for Figures 1 and 2). Note that COI and 16S rRNA partitions automatically appear in Partitions without employing PartitionFinder, and the partitioning is performed simply by applying each COI and 18S rRNA sequence, instead of the concatenating of genes by SeaView (Gouy et al. 2010) as done by Price et al. (2019) and Hill et al. (2021). Additional partitioning by PartitionFinder 2 (Lanfear et al. 2016) in MCMCTree and BEAST2 analyses is not required in the present BEAST1 analyses.
Roček et al. 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; see Table 2), and the right screen: Included monophyletic Taxa (= specific clade) and the resting Excluded Taxa in the central screen. As input in Figures 1-3 Figures 1-3 to avoid visual complexity. 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 Figures 2 and 3. This rate related function was not used in any previous paper, and we found in this paper variable base substitution rates the time as suggested 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 (Figures 1-3 inset) using a function of Excel.
The inset in Figures 1-3 shows that the base substitution rate was relatively slow until the Quaternary higher rate. To evaluate whether the slow rate reflected saturation, we examined the relation between pairwise distance and number of transition or transversion for each gene, using the MEGA5 function (Tamura et al. 2011; Figure 4).
Fossil and geological event calibrations by BEAST v1. X Calibrations points are shown on minimum age trees in Figures 1-3, and these dates were input in "Priors" in BEAUti as noted above; they are summarized below ( Table 2). As noted above, corresponding ingroup species were included in ingroup taxa (= leaf node taxa in a specific clade) by Taxon Set on the Taxa screen in BEAUti. Fossil calibrations are after Johnson et al. (2018) and Moulds (2018) ( Table 2; Figures 1-3). For these fossil calibrations, some are based on radio-isotopic dating of 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. This time scale has been standardized by the International Commission on Stratigraphy (ICS) (www.stratigraphy.org) and the most recent version of the time scale is available at http://www. stratigraphy.org/index.php/ics-chart-timescale, and the explanatory paper related to the generation of the time scale is Cohen et al. (2013).
Calibration points Q1 to Q6 and Q8 to Q12 are after our geological event calibration that adopts a 1.55 Ma date (Osozawa et al. 2012). This geologic event calibration was used in previous studies of Platypleura cicadas (Osozawa et al. 2017a) and four cicada groups (Osozawa et al. 2021a).
The  (Osozawa et al. 2012). The age assignment is from multiple biostratigraphic and radio-isotopic ages from the oldest marine strata on the landward side of the islands as well as the sides facing other islands, so that the age of such strata constrains the physical separation of the islands from the mainland and each other. There is no geologic evidence for land bridges that could have aided dispersal in the Ryukyu Islands.
Calibration point Q7 (Meimuna opalifera) is distinct from the above 1.55 AE 0.15 Ma event calibration. Hachijo oceanic island is a part of the Izu volcanic arc, and we recently estimated the emergence time of Hachijo-jima as an island at 0.24 Ma (Osozawa et al. 2021b). This date is applicable for crown Meimuna opalifera on the Hachijo-jima + the Japan-Tokara islands (= Stem Meimuna opalifera on Hachijo-jima).
With the assumption that the oldest fossil age is equivalent tMRCA (= minimum age), the specific fossil calibration points and associated dates are as follows: Calibration point A1: Crown Hemiptera: Fossils of Aphidoidea were reported from the French Bundsandstein (Szwedo & Nel 2011;Bashkuev et al. 2012) of Anisian age (244.5 AE 2.5 Ma).
A3: Fossil Ledridae (Zhang 1997) and fossil Cercopidae (Hong 1982) were recovered from the Jehol Biota of northern China. The Jehol Biota horizon has been dated by the Ar-Ar method on associated silicic tuff at 130.7 AE 1.4 Ma (He et al. 2006 (Mcintosh et al. 1992). However, we used an older crown date for crown Tibicininae based on fossil Davispia bearcreekensis that were found in the Fort Union Formation, Montana, USA (Cooper 1941). The age of the enclosing strata has been considered Thanetian in age (57.6 AE 1.6 Ma) (Flores & Bader 1999). Crown Cryptotympanini: Fossil Hadoa grandiose were also found in the Florissant Formation, Colorado, USA, but this calibration generated an unreasonable tree and was not adopted.
Calibration point G: Crown Cryptotympana: Fossil Cryptotympana incasa and C. miocenica (G1), and also Hyalessa lapidescens (G2) were found in Shanwang, Shandong, China (Zhang 1989;Zhang et al. 1994), and these strata are considered to be time correlative to the European MN5 mammalian unit ( Calibration point H: Crown Meimuna spp.: Fossil Meimuna protopalifera were found in the Itamuro Formation, Tochigi, Japan (Fujiyama 1969;Yoshikawa 2005), and the zircon fission track age of correlative terrestrial strata of the Nashino Formation of the Sendai area is 6.4 AE 0.4 Ma (Fujiwara et al. 2008).

Results
Hemiptera minimum age tree (Figure 1) Our timetree spans a range as old as ca. 250 Ma, and there is no evidence of saturation of mutations (Figure 4), suggesting our minimum age tree is robust and reliable.
Because the topology is concordant between Figures 1 and 2  Hemiptera, including Cicadoidea, has a single common ancestor of 242.96 Ma, as calibrated by the 244.5 AE 2.5 Ma age reviewed above as A1. The dated tree of the outgroup Hemiptera calibrated by A1 to A6 was concordant to Johnson et al. (2018) and Misof et al. (2014).
In the Cicadoidea ingroup, Tettigarctidae was an old lineage that differentiated from Cicadidae at 200.63 Ma, as calibrated by 200.3 AE 1 Ma (calibration point B), so Tettigarctidae is essentially a living fossil that has persisted since 200.63 Ma. We estimated a date of the common ancestor of two extant species of Tettigarcta tomentosa (Tasmania) and T. crinita (southeast Australia) at 13.96 Ma, and the youngest fossil of Tettigarctinae was reported from the Aquitanian (21.735 AE 1.295 Ma), southern New Zealand (Kaulfuss & Moulds 2015). However, Tettigarctidae includes 19 extinct genera according to Kaulfuss and Moulds (2015) and with many more genera according to Moulds (2018). A single common ancestor of Cicadidae except Derotettiginae started differentiation and speciation into Tibicininae, Tettigomyiinae, Cicadettinae, and Cicadinae at 66.15 Ma. Although the pre-Miocene fossil Cicadidae collectively include ten extinct genera, comprising Davispia and Lithocicada for Tibicininae, Paracicadetta, Paleopsalta, Minyscapheus, and Miocenoprasia for Cicadettinae, and Burmacicada, Camuracicada, Tymocicada, Dominicicada for Cicadinae, the remaining 23 genera post-Oligocene fossil cicadas are extant (Moulds 2018). Cicadidae, consisted of only one species but coexisted with a Tettigarctidae species between 200.63 and 66.15 Ma, and cicada biodiversity was extremely low during this period except for extinct species and D. mendosensis.
In the Cicadettinae major clade, each tribe constitutes a distinct clade. In the Cicadinae major clade, apart from older five tribe clades containing only one specimen, six tribe clades of Platypleurini, Cryptotympanini, Psithyristriini, Dundubiini + Cosmopsaltriini, Polyneurini + Sonatini, and Leptopsaltriini + Gaeanini are recognized. Discrepancies are addressed by reconsideration of taxonomy in the discussion.
The geologic calibration points Q1 to Q12 at 1.55 AE 0.15 Ma (and 0.24 Ma) apply to multi furcations that were recognized for Mogannia minuta and other cicadas endemic to in the Ryukyu Islands and Taiwan (and in Hachijo-jima) as noted above. Each island or island group population was mostly genetically distinct, endemic, and cryptic, as shown for Platypleura in Osozawa et al. (2017a). This also applies to Meimuna opalifera on Hachijo oceanic island (Osozawa et al. 2021a,b). However, note that some cicadas were accidentally dispersed by super typhoons up to 1,000 km in modern and ancient times including Meimuna boninensis (Osozawa et al. 2021a).

Inconsistent cicada base substitution rate (Figures 1-3 insets)
Comparing base substitution rate vs age shows that the rate has not been constant; the rate appears to have exponentially increased into the Holocene. The data points, approximate curve, and associated equation are shown on the insets of Figures 1-3. The curves and associated rates are similar for analyses based on COI alone (Figures 1 and 2 insets), and combined COI + 18S rRNA (Figure 3 inset). Figure 4 shows that even mitochondrial COI gene with rapid base substitution rate (Osozawa et al. 2017a) is never saturated toward the ancient time up to ca. 250 Ma.

Discussion
Taxonomic implications from the dated tree Tibicininae is solely from North and South America with an exceptional occurrence from the Mediterranean region, but absent from Asia and Africa (+ Australia). The stem age is estimated at 66.15 Ma (Figures 1 and 2), and if we assume that Tibicininae was generated by vicariance its differentiation may have been influenced by the formation of the Atlantic Ocean. Marine magnetic anomalies on the Atlantic Ocean floor can be used to ascertain spreading history and separation of continents that resulted from this spreading. The configuration at Chron34 (84Ma) after the Cretaceous magnetic quiet zone (long normal polarity epoch; superchron K-T at 118-84 Ma) was shown by Moulin et al. (2010), and the south Atlantic Ocean spread over 500 km (minimum distance between Africa and South America) at Chron 34 (84 Ma (Osozawa et al. 2011), and this species in southern Taiwan was differentiated relative to the northern Taiwan species reflecting vicariance triggered by the physical barriers of the Yilan basin and Lanyang valley (Osozawa et al. 2017b; http://kawaosombgi.livedoor.blog/?p=26 and others).
We combined East Asian Platypleura data after Osozawa et al. (2017a) with mostly African Platypleurini data excluding Platypleura after Price et al. (2019), and the terminal node of the East Asian Platypleura in the Platypleurini clade suggests the possibility of a Gondwanan origin and dispersal to Far East of Japan and Ryukyu, and Taiwan (Price et al. 2019). See Osozawa et al. (2017a) for the Platypleura vicariant speciation (see below for the cryptic speciation) on the Ryukyu Islands.
Tacua speciosa (Tacuini) represents the basal lineage of Cryptotympanini concordant with Marshall et al. (2018). In the Cryptotympanini clade, Auritibicen in Japan is the basal lineage, and Lyristes plebejus (synonym: Tibicen plebejus) in Croatia is the next. Asian Cryptotympana species is a sister of North American Noetibicen species (Hill et al. 2015).
Intercontinental dispersal by way of a Bering land bridge during Oligocene to Miocene climatic optima (Wu et al. 2015) was proposed for Papilionoidea butterflies feeding on Magnoliidae.
Zammara smaragdina, Costa Rica, represents the basal lineage of the remaining major clades that are paraphyletic each other. Distantalna splendida, renamed from Tosena splendida by Boulard (2009;referred in Hill et al. 2015), is represented by Tosenini, as a next basal lineage, distinct from another Tosenini of Tosena melanopteryx in the Psithyristriini composite clade. Tosena (Tosenini) and Pomponia (Psithyristriini) has a sister relationship, and these are similar tribes (species level transfer may be needed; Duffels & Hayashi, 2006). Pomponia backanensis, northern Vietnam, was described by Pham et al. (2015). Pomponia linearis on the Yaeyama islands and Taiwan, mildly differentiated each other as cryptic species, was renamed P. yayeyamana based on Kato (1932). The original P. linearis was reported from primarily Indochina, and has been treated as the P. linearis complex, including cryptic Chinese and Indian populations (Hayashi & Saisho 2011 (Chen 2011), and may have adaptively radiated within the island.

Recently increased cicada biodiversity
Hemipteroid insects of Psocodea, Thysanoptera, and the subject of this study, Hemiptera, include 120,000 described species which comprise over 10% of known insect diversity; they date back to 400 Ma (Hemiptera: 300 Ma; Johnson et al. 2018). Johnson et al. (2018) estimated that differentiation into species took place primarily in the Cretaceous, including Cercopoidea, Gerriidea, Flatidae, and Cicadoidea, which are common to our analyses. However, they analyzed only two to nine taxa, in contrast to the 344 taxa of Cicadoidea and 8 taxa for other Hemiptera of our analyses. Misof et al. (2014) estimated mostly pre-Paleogene dates of differentiation into species including Cercopoidea, Aphididae, and Delphacidae, concordant with our analyses, with less than 13 taxa analyzed. Their higher-level phylogeny suggested long branches and an old lineage of each super family species concordant with ours, but did not suggest the geologically recent increase in insect diversity apparent from our analyses of 352 Hemiptera taxa ( Figure 2).
In Figure 2,  Cryptic species on each island of Ryukyu chain are typical examples of increased biodiversity. For example, Platypleura kaempferi in the Amami and Okinawa islands has light colored wings, contrasting with dark colored wings in Japan-Korea-China and Taiwan, and the clades are distinct from each other (Osozawa et al. 2017a). P. kaempferi is not a single species but includes at least two cryptic species of light or dark winged Platypleura. Cicadas calibrated by other Quaternary calibration points include cryptic species, which also contributed to increasing biodiversity. The Okinawa trough is currently spreading (widening) and the Ryukyu islands are separating from the Chinese mainland. Accordingly vicariant speciation and radiation is in progress, which is also contributing increasing biodiversity. On the Chinese mainland, Hyalessa maculaticollis and Platypleura hilpa extensively radiated to form cryptic species (Liu et al. 2018(Liu et al. , 2020. Exponentially increased base substitution rate as a factor of Hemiptera diversity, and their possible causes Figures 1-3 insets show a large range of base substitution rates for different time periods, at variance with the constant molecular clock hypothesis (relatively constant rate over time; Ho 2008). The trend in base substitution rates shows an exponential increase into the Holocene.
Such an increase in base substitution rate was first shown for taxa such as primates by Ho et al. (2005) who showed that a Quaternary calibration date resulted in a more a rapid base substitution rate than that associated with an older calibration date. They employed an older version of BEAST (v1. 3; Drummond & Rambaut 2003) that required repeated runs, applying a date at each calibration point. In contrast, BEAST v1. X, used in our analyses, can simultaneously apply multiple calibration points, as we have done using dates ranging from the Triassic to the Quaternary. As a result, the calculated increasing rate of base substitution in our analyses is not an artifact of a Quaternary calibration, but is constrained by multiple age calibrations across a wide range of geologic time. Therefore, although the base substitution rate trendlines and associated equations of Ho et al. (2005) are similar to ours, their timetrees do not reflect the changing of base substitution rates through time, but rather reflect a constant base substitution rate as if constrained by a strict molecular clock. A similar analysis was done for beetles in the Aegean region by Papadopoulo et al. (2010).
The increasing base substitution rate is apparently associated with the recently increasing cicada diversity, expansion, and radiation (in Figure 2 timetree) that started at 40.57 Ma. The timing of the most rapid diversification coincides with Quaternary environmental change, marked by the start of glacial-inter glacial cycles. The initiation of Quaternary glaciations may have been triggered by rapid expansion of land grasses (Poales), that led to increased carbon fixation that decreased atmospheric CO 2 concentrations, because of the high efficiency of CO 2 fixation of such C4 plants (Sage 2004;Taira 2007). C4 Poales appeared and began diversification during the Oligocene (23 -33.9 Ma) based on molecular clock approach, and after 14.5 Ma based on fossil evidence (Sage 2004). We estimated 20.35 Ma by our Angiospermae timetree, that also employed BEAST v1. X with robust plant fossil calibrations (Osozawa et al. 2021c).
Food plants of D. mendosensis are, however, C4 dicots of Amaranthaceae (see figure 9 in Sage 2004) and Chenopodiaceae in degraded salt-plain habitats in arid regions of central Argentina (Simon et al. 2019). These dicot fossils and C4 monocot fossils of Poales grass (Chloridoidae) were reported from the Eocene in Patagonia by Zucol et al. (2018), and the fossil horizon was dated by the Ar-Ar method at 49.512 AE 0.019 Ma (Woodburne et al. 2014). The C4 photosynthetic pathway began at ca. 50 Ma in South America, earlier than elsewhere. For Chloridoideae, however, transition from C3 to C4 photosynthesis occurred in the Oligocene (23~33.9 Ma) as reported by Christin et al. (2008), consistent with our estimate for C4 dicots at 31.92 Ma (Osozawa et al. 2021c), whereas Sage (2004) suggested a fossil date at 14.5 Ma as noted above.
The trigger of increasing biodiversity may have been the generation and radiation of C4 plants and development of grass lands on the Earth since the Oligocene or perhaps more definitively since middle Miocene, by decreasing atmospheric CO2 concentrations. This may have led to the start of Quaternary ice ages and resultant adaptive radiation and increasing base substitution (≒ mutation) rates. Thus, biologic activity, including spreading C4 grasses may have significantly impacted Earth's environment.

Data availability
Sequence data in Table 1 are found in GenBank/DDBJ by incorporating the accession number.