Odonata Timetree; Exponentially Increased Base Substitution Rate Toward The Recent and Within the Carboniferous

Using BEAST v1.X, we constructed a credible timetree of 115 specimens of Odonata and ve species of Ephemeroptera (Paleoptera; Pterygota) and two species of Archaeognatha and three species of Zygentoma (Apterygota). 88 specimens we ourselves analyzed were collected from the Ryukyu islands, Taiwan, Japan, and China, and the resting sequence data were mostly from whole mitochondrial data found in GenBank / DDJB. The combined gene (not concatenated gene) analysis of the mitochondrial COI (795 bp), COII (548 bp), and 16S rRNA (517 bp), and the nuclear 28S rRNA (825 bp) were performed. Using the calibration function of BEAST v1.X, the timetree was constructed by applying a 1.55 Ma geological event (isolation of the Ryukyu islands from China), in addition to chronologically robust fossil dates ranging from 400 Ma for Archaeognatha, 300 Ma for Ephemeroptera, and 200 Ma for Odonata and to 1.76 Ma for Calopterygidae, for a total of 13 calibration points (event: 6, fossil: 12; Quaternary 7, pre Quaternary 11). The resultant timetree showed that molecular clock was not uniformly progressed, and the base substitution rate has exponentially increased from ca. 20 Ma to the Recent by over an order of magnitude. Our new and attractive nding indicates that the Quaternary severe climatic change including a start of glacial and interglacial cycle might have resulted in the extensive radiation and speciation of Odonata, and consequently increased the biodiversity. C4 pores generated in the Miocene effectively decreased atmospheric CO2, and triggered the Quaternary glaciation. Another peak of base substitution rate was found in the Carboniferous time around 320 Ma, and this may be analogous to the late Paleozoic icehouse. This glaciation has been triggered by the development of terrestrial plants to form thick coal layers, because this process also reduced the atmospheric CO2.


Introduction
A reliable timetree (dated phylogeny) for the entirety of Odonata (dragon y and damsel y) have not yet been reported to date, except for the very recent report in bioRxiv (Kohli et al., 2020). We have presented the timetree of the east Asian Coeliccia damsel y (details of the vicariance in Osozawa et al., 2017c) calibrated solely by the Quaternary date of geological event (1.55 ± 0.15 Ma; isolation of the Ryukyu-Japan-Taiwan islands from China; triggered the vicariant speciation onto the most recent common ancestor of the present each island and Chinese endemic species; no geological evidence of land bridge for dispersal; Osozawa et al. 2012). A problem with that analysis was that calculated basal node age of several million years was younger than what would be considered reasonable. This chronological problem also existed in analyses of other East Asian Odonata timetrees (Anotogaster dragon y, Osozawa This means that the BEAST v1.X calibration is an "active" one. If the timetree has unreasonable shape and very distinct input and output age, we reset the calibration and repeatedly run BEAST until to get reasonable one. In our previous analyses such as for Coeliccia damsel y in Osozawa et al. (2017c), we calibrated by geological event of 1.55 Ma and applied strict clock model in BEAST v1.8.2 analyses with the same calibration function to v1.X, and represented the base substitution rates of rapid mitochondrial and slower nuclear genes respectively. The rate was comparable to the previously reported rates such as by Brower (1994), although the rate was represented as relatively high (solution of the above problem is addressed below). Our present BEAST v1.X analyses applied relaxed clock model alternative to the strict clock model in Osozawa et al. (2017c), and the rate median at each node is not constant but valuable.
We made rate vs age diagram (Figs. 1 and 2 inset) to verify the variability, and as a main result we indicate an order of magnitude increase of base substitution rate in Odonata evolutional history found in the Quaternary and also in the Carboniferous. We propose here that the consequent increased biodiversity of Odonata was driven by severe environmental change that began in the Quaternary, and by global glaciation during the Carboniferous. We also discussed the ultimate factor of these glacier events, which was connected to the plant evolution.

Timetree (Figs. 1 and 2)
Insecta basal node was estimated at 413 Ma, Archaeognatha was represented as the oldest lineage in insecta, and Zygentoma is second. Ephemeroptera is a sister of Odonata differentiated at 347. 56 Ma, and this Palaeoptera is a sister of Apterygota (Archaeognatha + Zygentoma).
Anisoptera is a sister of Zygoptera, and each forms distinct major clade. Although Epiophlebia superstes is included in the Anisoptera clade, it can be considered to represent a distinct clade as an intermediate Anisozygoptera, concordant to Hasegawa and Kasuya (2006)  Note that ancient divergence times were successfully estimated by our mitochondrial gene application, with no sign of saturation of mutations, and application of nuclear gene was not affected on the topology and divergence time. Actually the base substitution rate of combined gene (Figs. 1 and 2 inset) is comparable to that of solely mitochondrial gene (Osozawa et al., 2017c).
In introduction, we noted a problem that the basal node age was estimated much younger than expected age when calibrated solely by the Quaternary geological event at 1.55 Ma. This problem was solved by contemporaneously applying old ages calibrated by fossil calibration points from A (400.45 Ma; Devonian) to I (1.76 Ma; Quaternary) in addition to geological event at 1.55 Ma (Quaternary) from Q1 to Q6, and the consequent dated tree (Figs. 1 and 2) is reliable.
Epiophlebia superstes can be considered to be included in the Anisoptera major clade, and is a sister of the other Anisoptera, common to Kohli et al. (2020). The stem age was estimated to be 189.5 Ma, but younger than 200.3 ± 1.0 Ma of calibration point C. In the Anisoptera clade, E. superstes persists as an independent lineage, so classi cation to Anisozygoptera was possible, and the species age is as old as 180.17 m.y. Accordingly, E. superstes can be considered a living fossil. Genus Epiophlebia is known from isolated areas of the Japan islands, China (the same species to China was reported also from North Korea; Gunther et al., 2013), and the Himalaya, but the vicariance was indicated to be very mild, suggesting that the isolation in these areas occurred recently and is unrelated to its very old lineage (Busse et al., 2012).
Petaluridae including Tanypteryx pryeri and Petalura gigantea was dated at 149.7 Ma (early Cretaceous). Ware et al. (2014) proposed that Petaluridae is an old species dating back to the Jurassic, re ecting the breakup of the supercontinent Pangaea. As we note in methodology, the Pangaea breakup and initiation of the Atlantic Ocean was recorded in the Santana Formation, and may at 125 -100 Ma (Martill, 2007), which is younger than the stem node age of Petaluridae. However, Petaluridae can be considered to be primitive and extant species (Ware et al., 2014).
Family level differentiations for Anisoptera were occurred during early Cretaceous time around 130 Ma, and those for Zygoptera were occurred during Paleogene to early Neogene between 72.64 to 22.19 Ma. So family age of Zygoptera is younger than that of Anisoptera. Subspecies level differentiations of Gomphidae were during Paleogene. Species level differentiation was after Paleogene, but mostly Neogene, and endemic species calibrated by 1.55 Ma as points Q1 to Q6 were differentiated during the Quaternary.
Yu and Bu (2011) cladistically analyzed damsel ies, and showed that Mesopodagrion and Rhipidolestes possess in distinct clade, and Megapodagrionidae was divided into two distinct families. Two distinct clades of Megapodagrionidae in Figs. 1 and 2 re ect this subdivision.
Vicariance increased biodiversity Q1: Stylogomphus, Q2: Asiagomphus, Q3: Anotogaster, Q4: Chlorogomphus (Osozawa and Wakabayashi, 2015), Q5: Rhipidolestes, Q6: Coeliccia (Osozawa et al., 2017c) were calibrated by the 1.55 ± 0.15 Ma geologic event, and each of the multi-furcation or polytomy indicates the generation of multi-endemic species (Fig. 2). Isolation of islands by the opening of the Okinawa trough is a physical process, but such process clearly triggered increasing biodiversity. These endemic species might have experienced a severe bottleneck by the isolation of habitat, but the bottleneck might have effect on species the low diversity and the low nucleotide substitution rate (Zhai et al., 2017), which is con ict to the present case.
Calibration point I: Calopteryx + Matrona was calibrated by fossil date of 1.76 ± 0.22 Ma, and splitting between Calopteryx japonica (Japan) and Matrona basilaris (Amami-Okinawa; failed to amplify Taiwan specimen) appears to be a consequence (or a precursor) of 1.55 Ma vicariance (Fig. 2).
In general, back arc spreading to form continental islands may trigger vicariance that increases extent of biodiversity. There are many back arc basins in the western Paci c Ocean, and some are still active, including the Okinawa trough. The high diversity of terrestrial organisms observed on continental islands separated from continental landmasses by sea-oor spreading of this sort may be a consequence of the physical isolation of these islands by the rifting process and the consequent vicariance.
In contrast, the oceanic islands are also generated back arc spreading by building of volcanic edi ces in the ocean, but these had no terrestrial life prior to their formation of islands because they emerged above sea level as a result of progressive volcanism. Such islands initially acquired their terrestrial species by dispersal from other land masses, and the species diversity associated with such islands may be expected to increase by vicariance after initial colonization by dispersal (cf., Osozawa, K. et al., 2016).
We showed that lotic damsel ies tended to speciate vicariantly in contrast to the lentic damsel y (Osozawa et al., 2017c), and the similar report was by Letsch et al. (2016). However, as for Octogomphinae, lotic Davidius and lentic Trigomphus show no such tendency (Fig. 2). Libellulidae is mostly lentic, and vicariance may be negligible (= same sequence between island populations), and resulting in low species diversity, but we failed to amplify enough island specimens of Libellulidae to fully evaluate this premise. This is a case also for mostly lentic Coenagrionidae damsel y by lacking corresponding sequence data.  Osozawa and Wakabayashi, 2015). However, ages of older nodes were underestimated, and branch lengths were shortened, when using this Quaternary geologic calibration alone.

Increase base substitution rate and biodiversity
Base substitution rates estimated by old fossil dates of A to I for calibration are < 0.01 s/s/myr, and these resulted in older nodes and extension of branches. Both the fossil calibrations including the Quaternary of point I and the Quaternary geologic event calibrations of point Q1 to Q6 appear to be robust, so the discrepancy might be the result of changing base substitution rates that may have been relatively constant from 200 Ma at slower rate, but exponentially increased since at ca. 20 (Papadopoulo et al., 2010) needed to run repeatedly by applying a date at every calibration point, and a rate at every run was assumed to be constant at each node (= strict clock model). They found that a Quaternary date calibration produced a rapid rate. Note that we instantaneously calibrated by multi dates including older dates up to the Devonian as well as younger dates of the Quaternary, rather than solely calibrating by the Quaternary date, and the increasing rate is not solely an artifact of the Quaternary calibration. It should also be noted that combined gene analyses were not possible in these older versions of BEAST, and such the users needed to discuss the rate for every each gene such as COI vs 28S rRNA. 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 governed by the strict clock.
Increased base substitution rate may be expected to increase biodiversity, as is the case for the 1.55 Ma vicariance-induced biodiversity mentioned above. Vicariance is partly explained as radiation for Rhipidolestes amamiensis as forming distinct clades (Fig. 2). An alternative is that greater apparent biodiversity in geologically recent time may be an artifact of the greater availability (preservation) and consequent greater number of studies of recent geological sections (Rohde and Muller, 2005 A possible forcing factor of exponentially increasing evolution rate and probable biodiversity toward the Recent may be the start of glacial and interglacial cycle and the environment change in the Quaternary time since 2.58 Ma. The Quaternary glaciations may have been triggered by expansion of C4 land grasses and evolution of sea diatoms, because this process increased the carbon xation that consequently decreased atmospheric CO2 concentration (Taira, 2007). The expansion of C4 grasses was a global phenomena also including North America and South America beginning in the late Miocene and persisting to the present day, concerning the present glacial-interglacial period with a time lag (Cerling et al., 1997).
We found another peak of base substitution rate in the Carboniferous time (300 -350 Ma; 0.063 s/s/myr at the node between Ephemeroptera and Odonata; Figs. 1 and 2 inset), and this may be analogous to the late Paleozoic icehouse (= the Karoo ice age in South Africa; Montañez et al., 2013). Whereas glacial episodes in Earth history appear to be in uenced by landmass con guration also in this period, feedback from biologic development may have also played a role in triggering the onset on this glaciation. The late Paleozoic glaciation had been triggered by the development of terrestrial plants to form thick coal layers in the Carboniferous time (Franks et al., 2014), because this process also increased carbon xation that effectively reduced atmospheric CO2 (Taira, 2007).

Conclusion
Constant evolution rates that are the basis of the strict molecular clock model do not apply in the cases of Odonata and may not apply in a broader sense. The present biodiversity including Odonata may be a consequence of a geologically recent exponential increase in base substitution rate. This may re ect the onset of Quaternary glacial and inter glacial cycles associated with severe climatic change, triggered adaptive radiation of Odonata. Another peak was found in the Carboniferous time, also tied to the late Paleozoic icehouse.

Materials
East Asia Odonata sampling Sequence data including whole mitochondrial data for several species of Archaeognatha, Archaeognatha, Ephemeroptera, and Odonata that we did not analyzed or failed to get sequence data were found in GenBank / DDBJ and included in our BEAST v1.X analysis (Fig. 2).
Location map is in Fig. 3, and representative collected specimens are shown on Fig. 4. DNA sequence data are newly registered in the DDBJ / GenBank, excluding previously uploaded data of Chlorogomphus by Osozawa and Wakabayashi (2015) and Coeliccia by Osozawa et al. (2017c). Species isolate, country, accession number, collection date, and collector data newly obtained in this paper are shown in Table 1. We collected lotic Davidius nanus, D. moiwanus, and Lanthus fujiacus (Octogomphinae; Gomphidae; Japan; Kiyoshi and Sota, 2006), and lentic Octogomphinae species of Trigomphus citimus tabei and T. interruptus (southwest Japan), and T. melampus (northeast Japan). We will compare the extent of vicariance of these species with contrasting habitat, as analyzed in Osozawa et al. (2017c).
Sieboldius albardae (Lindeniinae; Gomphidae; Japan and Korea; lotic) is characterized by small head relative to large body.
Tanypteryx pryeri (Petaluridae; endemic in Japan) was considered to be morphologically and phylogenetically an extant, relict, primitive species with the other Petaluridae originated from the Pangean super continent, and the larvae live in fens and boggy seeps taking several years to mature (Ware et al. 2014). We analyzed other three Petaluridae species including P. gigantea, the world's largest dragon y found in Australia, using corresponding sequence data in Ware et al. (2014), although their data were not complete.
We collected many Aeshnidae species, but we only obtained the COI-COII sequence data from Polycanthagyna melanictera (Okinawa; lentic), due to mis-or no ampli cation. P. melanictera is mostly known from Japan, as well as such as Tobi-shima off shore of Sakata, northern Honshu, and also Hachijo-jima of the Izu-Bonin oceanic islands, indicating its strong dispersal ability by ying. Anaxparthenope is a very common species in the coastal area near Sendai (northeast Honshu), persisting after the 2011 Tsunami, but only A. imperator and A. junius data were available in GenBank / DDJB. We also failed to obtain reliable sequence data of Oligoaeschna pryeri (Japan), considered to be a primitive species in Aeshnidae (Ozano et al., 2012).
Although we collected Corduliidae and some are endemic species in the Ryukyu islands, we failed to amplify the sequences. However, corresponding sequence data for three Corduliidae species were available in GenBank / DDJB.
Owing to mis-ampli cation of COI-COII regions, our analyses of Libellulidae were limited to sg38 Orthetrum japonicum ( Fig. 2; Table 1). To facilitate calibration using the age of fossil Libellulidae, we also used GenBank / DDJB sequence data of whole mitochondrial. Libellulidae species are mostly lentic and commonly found on even small islands of Ryukyu, but an exception is lotic Sympetrum pedemontanum elatum (Japan; no COI-COII sequence data). However, we included S. eroticum (lentic) using GenBank / DDBJ data.
We collected Epiophlebia superstes (Anisoptera; lotic) from the Sendai hill (newly found habitat by the senior author) of northern Honshu. The larva spends 5 to 8 years in small streams, and leaves the water just before the emergence (Ozono et al., 2012).

DNA Extraction and Polymerase Chain Reaction Ampli cation, Sequence Alignment
We used the sequence data of the mitochondrial COI gene (3P; 795 bp), COII gene (5P; 548 bp), 16S rRNA (517 bp), and the nuclear 28S rRNA gene (825 bp). See our newly obtained data with accession numbers shown in Table 1. See Osozawa et al. (2017c) for analytical details, including the primer set used.
Note that these COI and COII gene data are intersecting the nonprotein cording tRNA-Leu region (Osozawa et al., 2017c), and the corresponding region data are seldom found in the DDBJ / GenBank, but we included sequence data in our analyses from DDBJ / GenBank as noted above (no isolate and only species name in Fig. 2).
Quaternary geologic event: Ryukyu continental islands formed at 1.55 Ma for BEAST v1.X analyses Osozawa et al. (2012) showed that the Ryukyu islands formed by back-arc spreading of the Okinawa trough that separated them from the Chinese mainland. This sea-oor spreading resulted in separation of the islands by 1.55 ± 0.15 Ma and this separation has increased progressively since then. Subsidence to form each island progressed rapidly, and these islands were simultaneously separated from the China mainland and isolated from one another by the development of the Okinawa trough, as well as other major seaways between the islands, including the Tsushima and Taiwan straits (Fig. 3). These straits, including minor ones in some cases, are expected to have acted as migration barriers to trigger vicariance, so 1.55 Ma is a geologically robust calibration date applied to our phylogenetic analyses in BEAST v1.X for six Odonata groups of Stylogomphus (Gomphidae; calibration point Q1), Asiagomphus They stated that their calibration date strictly follows the protocol proposed by Parham et al. (2012), that is geologically robust in the opinion of the present authors. We agree with their statements that "Review the published geological age and the stratigraphic range of the fossil and synchronize the data with the most updated standard global geochronology". "Review the quality of the stratigraphic information for each fossil. In particular, amber deposits such as those from the Baltic or Dominican Republic could not be adequately correlated with the geological time scale". Therefore, they have largely excluded amber fossils from their analysis, with an exception of Burmese amber, which was dated at 99 Ma using U-Pb method (Shi et al., 2012). "Biostratigraphic information has been updated and adapted to the current geological time scale" (the most recent International Chronostratigraphic Chart by the International Commission on Stratigraphy). Misof et al. (2014) used the calibration date for fossil egg insertion scars into host plant by damsel y (Moisan et al., 2012) and considered the basal node of Odonata to be late Triassic (235-221 Ma). However, the Madygen Formation in Kyrgyzstan (Voigt et al. 2006) that yielded this fossil has associated no constraining radioisotopic age data and no reliable biostratigraphic information (no index fossils that can be related to robust absolute age determination), so we consider this late Triassic date to be poorly supported, and we used alternatively reliable early Jurassic date for damsel y (see below). However, there is some uncertainty of the fossil ages (NOT fossil identi cations), so we reevaluated the geologic constraints associated with the fossil locations. We adopted several pre-Quaternary calibration dates, after evaluating their geological reliability, including, but not limited to, whether the fossil ages were constrained by modern and more accurate radioisotopic dating or by stratigraphical correlation. We applied 11 fossil calibration points from A to I (Figs. 1 and 2), and resulting in robust time calibration as old as 400 Ma. Note that the Devonian calibration point A was for Archaeognatha, and the Carboniferous calibration point B was for Ephemeroptera (our own biblographic survey and reevaluation; Figs. 1 and 2).
The labels preceding each text section below correspond to that given in the following section as calibration points on actual BEAST v1.X analysis (Figs. 1 and 2). These fossil dates are evaluated in order of decreasing age from A to I. As will be seen in the subsequent section, the ages with larger uncertainties were not used in BEAST analyses for generating Figs C: Wing pattern of fossil Liassophlebia resembles to that of Epiophlebia superstes included in the former Anisozygoptera (Nel et al., 1993), and the basal node represents crown Epiprocta (Anisozygoptera + Anisoptera) (Kohli et al., 2016). Liassophlebia fossil was found from the Bayreuth Formation between the Keuper (one of the three Triassic strata in Germany) and the Gryphaeensandstein Formation (Kohli et al., 2016  Libellulidae is also found from the famous fossil locality of Cereste, southern France (Nel and Paicheler, 1993), including Lethe butter y (Nel et al., 1993;Pfretzschner, 1998) , 1985), which represents a large relative uncertainty, so this age is not adopted for calibration, also considering the younger age than the Green River Formation.
I: Wing fossils of Calopteryx japonica (=Matrona basilaris japonica; analyzed and shown in Fig. 2) and Oligoaeschna pryeri (no reliable COI-COII sequence data were obtained) were found from the uvial to deltatic strata of the Beppu-Shimabara graben (Esaki and Asahina, 1957), northeastern extension and branch of the Okinawa trough ( Fig. 1; cf. Phylogenetic analyses by BEAST v1.X; Age calibration Bayesian inference (BI) trees were constructed using the software BEAST v1.X. The BEAGLE Library must be downloaded beforehand to operate the platform software of BEAST. The details are described in In BEAUti, "Taxa" and "Priors" (date setting is in "Prior"; date should be set in box of "tMRCA") were set as follows for ingroup species relative to four partitions of COI, COII, 16S rRNA, and 28S rRNA. These partitions automatically appear in "Partitions" if these sequence data are uploaded by using the plus button or by "Import Data". The protocol is that a MRCA at each calibration of point of A to I and Q1 to Q6 initiated to differentiate into the ingroup species. In FigTree, posterior probability ("posterior"), posterior age ("Node ages"), and "rate median" (not constant because we applied relaxed clock model) can be output, and these are shown on each node in Fig. 3. 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.3 left bottom inset). Figure 1 Simpli ed BI tree for combined genes of Odanata species, calibrated by both the 1.55 Ma event and older fossil dates, using BEAST v1.X. Calibration point H, for example, represents a crown node of a common ancestor of Libellulidae ingroup species, but the stem node can be speci ed if it is needed. In BEAUti, Page 23/26 monophyletic can be speci ed for Libellulidae ingroup species. Note that the ingroup of calibration point Q6, for example, contains Chinese Platycnemididae species as well as Ryukyu-Taiwan species. Inserted Figure: Base substitution rate vs age diagram. Red thick curve with equation: Exponential trendline drawn by Excel function, Red thin curve: Free handed trendline drawn by connecting plots with maximum rate through the time around 320 Ma. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 2
Page 24/26 BI tree for combined genes of Odanata species, calibrated by both the 1.55 Ma event and older fossil dates, using BEAST v1.X. Inserted Figure: Base substitution rate vs age diagram. Red thick curve with equation: Exponential trendline drawn by Excel function, Red thin curve: Free handed trendline drawn by connecting plots with maximum rate through the time around 320 Ma.

Figure 3
Geographical map of the Ryukyu and Japan Islands, Taiwan, Korea, and China. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 4
Anisoptera species analyzed in this paper.