Cicada minimum age tree: Cryptic speciation and exponentially increasing base substitution rates in recent geologic 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.

Taxon sampling

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 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).

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 ± 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, 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, and Hill 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 ± 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.

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?

Regarding BEAST2 (= *BEAST2, StarBEAST2), our approach diverged from previous studies such as Osozawa et al. (2016), Price et al. (2019), and Hill et al. (2021), who utilized BEAST v2.5 (Bouckaert et al., 2014). Instead, we opted for BEAST v1.8 and subsequently v1. X. While the calibration function in BEAUti of BEAST v2.5 bears similarities to BEAST v1. X, there are notable differences. In BEAST v2.5, the “Partition” tab only permits the input of individual sequence data. Consequently, if the sequence data are not concatenated, separate BEAST runs must be conducted for each set of applied sequence data (e.g., mitochondrial COI and nuclear 28S rRNA), as demonstrated by Osozawa et al. (2016). The resulting tree files from these runs must then be combined into a single file using LogCombiner. However, when merging these tree files, the branches in the resultant tree become folded, reflecting the incongruent topology arising from different sequence data sources, such as mitochondrial COI and nuclear 28S rRNA. To mitigate this issue, Osozawa et al. (2016) employed DensiTree to obscure the foldings. Consequently, we discourage the usage of BEAST v2.5 due to the inconvenience and potential confusion caused by folded branches in the combined tree.

In the case of BEAST v2.6, which was released in May 2019, and BEAST v2.7, released in 2023 (Bouckaert et al., 2019), significant changes were made to the protocols. A tutorial for these versions can be found in https://taming-the-beast.org/tutorials/starbeast2-tutorial. Notably, the inclusion of cladistic data alongside molecular data became possible with the implementation of total-evidence dating (Zhang et al., 2016). For extinct species, tip dating is set at their youngest fossil age, while for extant species, it is set at zero age. However, the fossil age is often poorly constrained, with minimum and maximum age ranges typically used. In BEAST v2.6 and v2.7, the calibration and node dating function that was implemented in BEAST v2.5 was abandoned, and node dating for extant species is solely based on applying and assuming the base substitution rate.

In the context of BEAST v2.6 and v2.7, it is important to clarify that the term “tip” does not refer to terminal nodes for extant species. Instead, it refers to the tip node representing extinct fossil species from ancient times (c.f., https://beast.community/first_tutorial). The tip date for fossil species is inferred from the fossil age, and it is worth noting that the age assigned is not necessarily the minimum age for the oldest fossil, but rather the youngest fossil, which is often poorly constrained. Additionally, it is crucial to ensure that these fossil species are indeed extinct, and determining their relative placement in relation to the lineage of extant species can be problematic, as it involves the concept of ghost lineage. It is important to understand that tip dating does not contribute to the quality of node dating.

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 18S rRNA gene sequences appeared in the Partition box (For Figure 3; COI file only for Figures 1 and 2). Note that COI and 18S 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.

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, calibration dates were set in Priors bellow.

Tips and Traits: Default.

Sites: Substitution Model: HKY (Hasegawa, Kishino and Yano) model, Base frequencies: Empirical, Site Heterogenety Model: Gamma, Number of Gamma Categories: 4, Partition into codon positions: Off. The GTR model generates similar topology.

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).

Trees: Tree Prior: Speciation: Yule Process.

Priors: tMRCA (time of MRCA) was input from the calibration point date as Prior Distribution: Normal, and as the Mean and Standard deviation. See bellow for Priors as detailed setting of age calibration.

Operators: Default.

MCMC: Length of chain: 10,000,000.

Running BEAST was done by incorporating xml input file made by BEAUti. The consequent tree was drawn by FigTree v1.4.2, for that, the tree files were input into TreeAnnotator. The 95% highest posterior density for confidence intervals of ages can be output in FigTree, but not shown in 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).

Figure 4. Number of base changes of transition and tansversion vs corrected pairwide distance diagram for whole mitochondrial gene.

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 specific calibration points are as follows: tMRCA of Mogannia minuta (Q1), M. hebes (Q2), Cryptotympana facialis (Q3), dark winged Platypleura (Q4) and right winged Platypleura (Q5; Osozawa et al. 2017a), Graptopsaltria nigrofuscata + G. bimaculata (Q6), Meimuna kuroiwae (Q8), M. oshimaensis + M. iwasakii + M. goshizana (Q9), Tanna japonensis + T. japonensis ishigakiana + T. sozanensis + T. sp. (Q10) Euterpnosia chibensis + E. chibensis daitoensis + E. chibensis okinawana (Q11), E. iwasakii + E. viridifrons + E. olivacea + E. gina + E. sp. (Q12): The date of the geological event, which records the isolation of the Ryukyu Islands from the Chinese mainland by the opening of the Okinawa trough that began (i.e., islands had separated from mainland and each other by this time) at 1.55 ± 0.15 Ma (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 ± 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 ± 2.5 Ma).

A2: The oldest fossil Belostomatidae was reported from the Zagaje Formation, Poland (Popov 1996) of Hettangian age (200.3 ± 1.0 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 ± 1.4 Ma (He et al. 2006).

A4: Fossil Gerridae were recovered from French amber (Perrichot et al. 2005) of Albian age (107 ± 6 Ma).

A5: Fossil Delphacidae were found in the Green River Formation, USA (Grande 1980). Ar-Ar dating on silicic tuff within the formation yields ages of 53.5 – 48.5 Ma (weighted average age of 51.25 ± 0.31 Ma; Smith et al. 2003).

A6: Fossil Flatidae were found in the Maíz Gordo Formation, northwest Argentina (Petrulevičius 2011) of Paleocene age (61 ± 5 Ma).

Calibration point B: Stem Tettigarctinae: Oldest fossil of Tettigarctinae were found in strata Dorset, England (Whalley 1985) of Hettangian age (203.1± 1.0 Ma).

Calibration point C: Stem Derotettiginae: The preferred food of Derotettix mendosensis is Amaranthaceae in Argentina (Simon et al. 2019), and this worldwide C4 plant was phylogenetically studied by Piirainen et al. (2017). This plant fossil was reported by Zucol et al. (2018), and the fossil-bearing horizon was dated by the Ar-Ar method at 49.512 ± 0.019 Ma (Eocene; Woodburne et al. 2014). However, fossil Burmacicada protera were found from Burmese amber (Poinar & Kritsky 2011). Detrital zircons from the amber bearing matrix yielded a maximum depositional age U-Pb age of 98.79 ± 0.62 Ma, that was interpreted to closely approximate the actual depositional age on the basis of geologic relationships and associated fossils (Shi et al. 2012). We applied this older date of Burmese amber for stem Derotettiginae or crown Cicadidae.

Calibration point D: Stem Platypedia putnami (= crown Tibicininae): Fossil Platypedia primigenia were found in the Florissant Formation, Colorado, USA, and the associated strata was dated by the Ar-Ar method at 35.15 ± 1.65 Ma (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 ± 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 E: Crown Cicadettinae: Paracicadetta oligocenica (Boulard & Nel 1990) were recovered from deposits of Céreste, France, and this famous fossil locality was considered to be of Ruperian age (31 ± 2.9 Ma; Ducreux et al. 1985).

Calibration point F: Stem Lyristes plebejus: Fossil Lyristes sp. were reported from Seifhennersdorf, Germany (Tietz et al. 1998), and associated strata was dated by the K-Ar method as 30.44 ± 1.52 Ma (Walther & Kvacek 2007).

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 (16.45 ± 0.45 Ma; Roček et al. 2011).

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 ± 0.4 Ma (Fujiwara et al. 2008).

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 (COI) and Figure 3 (COI + 16S rRNA), the following description follows Figure 2 with 352 specimens. Our analyses was concordant to the subfamily classification of Marshall et al. (2018), Łukasik et al. (2019), and Simon et al. (2019). Figures 1 and 2 also include data in Price et al. (2019) and Hill et al. (2021).

Hemiptera, including Cicadoidea, has a single common ancestor of 242.96 Ma, as calibrated by the 244.5 ± 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 ± 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 ± 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).

Simon et al. (2019) proposed a new subfamily Derotettiginae consisting of a single species of Derotettix mendosensis, which is a sister of the remaining Cicadidae species and the oldest lineage species in Cicadidae dated at 99.2 Ma, as calibrated by point C at 98.79 ± 0.62 Ma. Łukasik et al. (2018) showed such a basal lineage of D. mendosensis in Cicadidae.

Our timetree showed that Tibicininae is a sister of Tettigomyiinae + Cicadettinae + Cicadinae and differentiated at 66.15 Ma, and Tibicininae started differentiation at 57.31 Ma, as calibrated by point D at 57.6 ± 1.6 Ma. Tettigomyiinae is a sister of Cicadettinae and differentiated at 35.46 Ma, Tettigomyiinae + Cicadettinae is a sister of Cicadinae differentiated at 40.57 Ma. Cicadettinae started differentiation at 30.85 Ma, as calibrated by point E at 31 ± 2.9 Ma. Cicadinae started differentiation at 38.25 Ma. Differentiation of Tettigomyiinae + Cicadettinae took place simultaneously after 35.46 Ma.

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 ± 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.

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). The date of 84 Ma can be considered to be a starting date of continent level vicariance, which may have triggered the Tibicininae differentiations relative to especially Cicadinae shown in Figures 1 and 2.

In Cicadettinae, Prasiini is a sister of Cicadettini. Muda kuroiwae in Prasiini (Hayashi & Saisho 2011) is endemic and restricted to Okinawa-jima and Kume-jima, and represents as a sister of the similar species of Katoa taibaiensis on the Chinese mainlamd.

In the Moganniini clade, Nipponosemia terminalis (Matsumura 1913) (synonym: Vagitanus terminalis) is a sister of Mogannia spp. N. terminalis has been documented from the Yaeyama islands and Miyako-jima (endangered and protected), Ryukyu, and Taiwan (Figure 2 from the Taiwan specimen; another species of N. virecens is known from the Kaoshun peninsula, southern most Taiwan; Lee & Hayashi 2004), but Yang & Wei (2013) reported N. terminalis and other three Nipponosemia from China. A detailed phylogenetic study for these Nipponosemia species would be useful. The genitalia and morphological character are similar to Mogannia (Hayashi & Saisho 2011), concordant to the sister relationship with Mogannia. See Osozawa et al. (2021a) for the Mogannia minuta vicariant speciation on the Ryukyu Islands and the accidental typhoon dispersals in recent and also ancient times. Mogannia hebes in northern Taiwan and in southern China has sister relationship reflecting vicariance by the Taiwan strait (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). Unipomponia decem (Psithyristriini?) was renamed Pomponia decem (Lee & Sanborn 2010).

Megapomponia (Lee & Sanborn 2010) was associated with the genera from Dundubiina (Hill et al. 2021), and included in the Dundubiini clade. Oceanian Cosmopsaltriini is a sister of Asian Dundubiini, reflecting large scale vicariance driven by bio-geographic barrier of the Wallacea line, as well as endemism within the islands by Oceanian arc fragmentations (Boer & Duffels 1996). In the Dundubiini clade, see Osozawa et al. (2021a) for the vicariant and cryptic speciation of Meimuna kuroiwae on the Ryukyu islands and the accidental typhoon dispersals in recent and ancient times (including Meimuna boniensis on the oceanic Bonin islands). See Osozawa et al. (2021a) for the vicariant and cryptic speciation of Meimuna opalifera on the oceanic Hachijo-jima island (Osozawa et al. 2021b) by the accidental typhoon dispersal from the Japan continental islands. Meimuna mongolica in Korea and China is the basal lineage relative to the sympatric M. opalifera. Meimuna oshimensis endemic on the Amami and Okinawa islands (cryptic species), Meimuna iwasakii endemic on the Yaeyama islands and Taiwan (no specimen collected from Taiwan specimen; cryptic species), and Meimuna goshizana and M. gakokizana, other endemic species on Taiwan, were vicariantly speciated or adaptively radiated in Taiwan.

Hyalessa is in Sonatini was renamed from Oncotympana in the distinct Oncotympanini. Hyalessa maculaticollis in Japan and China (Liu et al. 2018) is deeply differentiated. Sonatini is a sister of Polyneurini (once included in Tosenini; noted in Hayashi & Saisho 2011), and these constitute the Polyneurini + Sonatini clade. Graptopsaltria nigrofuscata in Japan and G. bimaculata on the Amami and Okinawa islands were vicariantly speciated.

Terpnosia cf. graecina was synonymised with Leptopsaltriina (discussed in Hill et al. 2021) constitutes the basal lineage with Thailand Leptopsaltriini in the Leptopsaltriini-Gaeanini major clade. Kalabita operculata, Malaysia, was a member of above-mentioned Platypleurini known to have diversified in Africa (Price et al. 2019), but this Asian species is included in the Leptopsaltriini + Gaeanini clade. Trilar (2006) presented an adult photo of K. operculata showing that it lacks the pronotum that characterizes Platypleurini. Furthermore, the spectrogram - oscillogram of its song is similar to those of Euterpnosia spp., Leptopsaltriini, shown in Hayashi and Saisho (2011). Platypleurini both in African and Asia are monophyletic as noted above, and an exception is unreasonable. Hill et al. (2021) transferred Kalabita Moulton, 1923 from Platypleurini to Leptopsaltriini. Tosena paviei is a member of Tosenini, but is a sister of Callogaeana guanxiensis in Gaeanini, and renamed as Vittagaeana paviei in Gaeanini by Hill et al. (2021). Both the species are from Vietnam, and constitute a Gaeanini clade with another Vietnam species of Balinta cf. tenebricosa of Gaeanini. An interesting result is that Gaeanini including T. paviei is not monophyletic, and paraphyletic in the Leptopsaltriini major clade. Hill et al. (2021) showed wing phenotypes of Gaeanini- Tosenini, which are distinct from phenotypes of Leptopsaltriini. Note that Gaeana (Gaeanini) is also a sister of Tanna (Leptopsaltriini) with different wing and chest phenotypes. Tanna japonensis is differentiated between Japan and the isolated population on the Amami-Oshima island, and further differentiated from the isolated Tanna japonensis ishigakiana population on Ishigaki-jima island. Taiwan yields Tanna sozanensis (sister of T. japonensis ishigakiana) and the other seven Tanna species (Chen 2011), might have adaptively radiated within the island. Tanna in China is a sister of T. sozanensis. According to Hill et al. (2021), Cicadmallus micheli is characterised by an unusual ‘hammer-head’ morphology but otherwise bears morphological relationships to Leptopsaltriini, and represents the basal lineage of Indochina Terpnosia (Thai & Yang 2009) and East Asian Euterpnosia of the Leptopsaltriini clade. See Osozawa et al. (2021a) for the vicariant and cryptic speciation for Euterpnosia for the northern population on Japan-Amami-Okinawa and the southern population on Yaeyama-Taiwan, as well as accidental typhoon dispersal in ancient times (Euterpnosia chibensis daitoensis on the oceanic Daito islands from Tokuno-shima continental island). Taiwan yields Euterpnosia gina, E. olivacea, E. viridifrons, and other 12 Euterpnosia species (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, ingroup Cicadidae, excluding Derotettiginae, underwent extensive differentiation into 341 taxa after 66.15 Ma, mostly after 40.57 Ma, leading to increasing biodiversity of Cicadidae, although Price et al. (2019) suggested that number of lineages saturated in the Pleistocene. Cicadidae consisted of only two species including D. mendosensis between 99.2 and 66.15 Ma, although Cicadidae contains many extinct species that remain to be identified as fossils (Moulds 2018).

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, 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 rapid diversification of Delphacini is linked to shifts in host plants, particularly from C3 to C4 grasses within the Poaceae family (Urban et al. 2010). 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₂ concentrations, because of the high efficiency of CO₂ 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 ± 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.