Cretaceous origin of giant rhinoceros beetles ( Dynastini ; Coleoptera ) and correlation of their evolution with the Pangean breakup

The giant rhinoceros beetles (Dynastini, Scarabaeidae, Coleoptera) are distributed in tropical and temperate regions in Asia, America and Africa. Recent molecular phylogenetic studies have revealed that the giant rhinoceros beetles can be divided into three clades representing Asia, America and Africa. Although a correlation between their evolution and the continental drift during the Pangean breakup was suggested, there is no accurate divergence time estimation among the three clades based on molecular data. Moreover, there is a long chronological gap between the timing of the Pangean breakup (Cretaceous: 110–148 Ma) and the emergence of the oldest fossil record (Oligocene: 33 Ma). In this study, we estimated their divergence times based on molecular data, using several combinations of fossil calibration sets, and obtained robust estimates. The inter-continental divergence events among the clades were estimated to have occurred about 99 Ma (Asian clade and others) and 78 Ma (American clade and African clade), both of which are after the Pangean breakup. These estimates suggest their inter-continental divergences occurred by overseas sweepstakes dispersal, rather than by vicariances of the population caused by the Pangean breakup.


INTRODUCTION
The tribe Dynastini, known as the giant rhinoceros beetles, includes some of the largest species within the Coleoptera, which is the most species-rich order among all living organisms.These species show remarkable sexual dimorphism, such as very long horns in the male's head and pronotum that are usually absent in females (Endrödi, 1985).Naturalists have been attracted for centuries by such marked morphological variety and excess male ornamentation (Rowland and Miller, 2012).It is important to clarify how this sexual ornamentation evolved, to gain a better understanding of natural selection and especially sexual selection.To elucidate this issue, reconstruction of a reliable phylogenetic tree is prerequisite.Nevertheless, no molecular studies of the Dynastini had been attempted to reconstruct their phylogenetic tree with extensive taxon sampling until recently (Rowland and Miller, 2012).Rowland and Miller (2012) inferred the phylogenetic tree on the basis of morphological characters as well as molecular data including mitochondrial and nuclear gene sequences.They indicated that this tribe can be divided into three subtribes, Dynastina (Dynastes, Golofa, Megasoma, Augosoma), Xylotrupina (Allomyrina, Xylotrupes, Trypoxylus, Xyloscaptes) and Chalcosomina (Haploscapanes, Chalcosoma, Beckius, Eupatorus, Pachyoryctes) (Supplementary Table S1).The first subtribe consists of the African clade (Augosoma) and the New World clade (Dynastes, Golofa, Megasoma).On the other hand, the latter two subtribes are distributed in Asia and Oceania, and constitute a large clade (hereafter called the Asian clade).From this categorization, Rowland and Miller (2012) suggested a correlation between Dynastini evolution and the Pangean breakup in the Early Cretaceous.
Indeed, the correlations between continental drift and phylogenetic branching patterns of terrestrial animals such as Eutherian mammals (Nishihara et al., 2009;dos Reis et al., 2012), Palaeognathae (ratite and tinamou: Cooper et al., 2001;Haddrath and Baker, 2012;Yonezawa et al., 2017), Pleurodira (side-necked turtle: Vargas-Ramírez et al., 2008), Anura (frogs: Bossuyt et al., 2006), Cichlidae (cichlid fish: Friedman et al., 2013) andPlecoptera (stoneflies: McCulloch et al., 2016) have been vigorously argued.Some of these studies (e.g., Haddrath and Baker, 2012;dos Reis et al., 2012;Yonezawa et al., 2017) suggested that biogeographical distribution patterns were not necessarily established merely by continental drift, but that the overseas 'sweepstakes dispersal route' (Simpson, 1940) should be considered.A sweepstakes dispersal route is "a possible route of faunal interchange which is unlikely to be used by most animals, but which will, by chance, be used by some" (Allaby, 1999).Since we focus on the establishment of intercontinental distribution patterns, hereafter we call them 'overseas sweepstakes dispersal' as an alternative to the hypothesis that vicariance was caused by the Pangean breakup.
In addition, no paleontological evidence supports the hypothesis that the breaking up of Pangea triggered the speciation of the subfamily Dynastinae because their fossil record can be traced back to after the Oligocene (Krell, 2000(Krell, , 2007)).Therefore, reliable divergence times within the Dynastini based on molecular data are important for a better understanding of the correlation between their evolution and geological events such as continental drift.Here, we present the first clear picture of the time-scale of giant rhinoceros beetle evolution and its correlation with geological events by detailed reanalyses of Rowland and Miller (2012)'s molecular and morphological data.

MATERIALS AND METHODS
Phylogenetic inference and divergence time estimations Detailed procedures for tree inference and divergence time estimation were as follows.The nucleotide sequence data published by Rowland and Miller (2012) (cytochrome c oxidase II, 16S rRNA, histone III and arginine kinase) were retrieved from NCBI.Rowland and Miller (2012)'s data cover all genera among the tribe Dynastini.In addition, the corresponding nucleotide sequences of Tribolium castaneum (family Tenebrionidae) were also included as an outgroup.The accession numbers are presented in Supplementary Table S1.Retrieved sequences were automatically aligned by the MAFFT program (Katoh and Standley, 2013), and carefully checked by eye.The impact of the partition model on both tree inference and divergence time estimation was discussed in detail in our previous studies (Omoto et al., 2009;Wu et al., 2014).The unlinked partition model (Pupko et al., 2002) was applied to take account of the different tempo and mode of nucleo-tide substitutions among genes and codon positions.The best partition was found by the PartitionFinder program (Lanfear et al., 2012).The best partition is shown in Supplementary Table S2, and it was applied for both the phylogenetic inference and the divergence time estimations under the GTR + Γ model.The phylogenetic tree was inferred by the maximum likelihood (ML) method using the RAxML program ver.8.0.20 (Stamatakis et al., 2008).The rapid bootstrap method (Stamatakis et al., 2008) was applied to evaluate the confidence of the internal branches with 1000 replications.Divergence times were estimated by the Bayesian relaxed clock method (Thorne et al., 1998) with the MCMCTREE program implemented in PAML ver.4.7 (Yang, 2007;Inoue et al., 2010).To reduce the computational burden, the normal approximation method was applied.The autocorrelated rates model (Thorne et al., 1998) and the independent model (Rannala and Yang, 2007) were used for the drift model of the evolutionary rates along the branches.The prior distributions of the root rates were set to 4 and 2.26 for the α and β parameters of Γ distribution, respectively, and the parameters of σ 2 were set to 1 and 0.2.Other parameters such as the Birth and Death rate were set to default.The condition of the MCMC was as follows: the first 200,000 generations were discarded as burnin; then, 20,000 trees were sampled per 200 generations.
The history of morphological character traits was inferred using MESQUITE ver.3.3 (http://mesquiteproject. org) with the maximum parsimony method on the basis of the 17 characters reported by Rowland and Miller (2012).
Fossil calibrations Reliable fossil records of the family Scarabaeidae (Krell, 2000(Krell, , 2007) ) were used as calibration points.The age of the oldest known record (OKR) of a fossil in a particular crown group does not necessarily correspond to the exact timing of its emergence because of the incompleteness of the fossil record.Therefore, the emergence timing of the crown group should be assumed between the OKR of crown taxa and the OKR of the stem taxa (Yonezawa et al., 2007).Here we suggest the calibrations of divergence times (maximum and minimum boundaries of the nodes) in the Dynastini based on this idea.The first calibration is the split between the subfamilies Dynastinae and Melolonthinae (node 2 in Fig. 1).The subfamily Dynastinae (Dynastini species, Orizabus and Cyclocephala) is closer to the subfamily Cetoniinae (Cotinus mutabilis) than to the subfamily Melolonthinae (Polyphylla decemlineata).This relationship is supported by morphological (Krell, 2000) and molecular data (Rowland and Miller, 2012, this study).The reliable OKR of Melolonthinae (Cretomelolontha transbaikalica: Fossil 1 in Fig. 1) is from Baysa, Zazinskaya Formation, Transbaykalia, Russia (Nikolajev, 1998;Krell, 2000Krell, , 2007)).This formation belongs to the Valanginian stage (132.9-139.8Mega annum (Ma; million years ago)) of the Lower Cretaceous.On the other hand, the OKRs of Dynastinae and Cetoniinae can be traced back to only after the Oligocene (Fossil 2 in Fig. 1).Therefore, the minimum boundary of the Dynastinae-Melolonthinae split (node 2) should be 132.9Ma.Generally, the maximum boundary of the divergence time is difficult to set, because the OKRs of the stem taxa do not guarantee the emergence of the crown taxa, especially when a long ghost lineage exists or a stem taxon survives for a very long time period as a relic (Yonezawa et al., 2007;Benton et al., 2009).The OKR of the whole family Scarabaeidae is Antemnacrassa spp. or Juraclopus rohdendorfi (Nikolajev, 2005; see also Krell, 2007; Fossil 3 in Fig. 1) from Karatau, Kazakhstan, and these fossils are thought to be of Kimmeridgian age (152.1 to 157.3 Ma) in the Late Jurassic.Therefore, the tentative maximum boundary of this split (node 2) was assumed to be 157.3Ma.The second calibration is the split of the family Scarabaeidae and the family Tenebrionidae (Tribolium castaneum) (node 1 in Fig. 1).Because the family Elateridae is closer to the family Tenebrionidae than to the family Scarabaeidae (Hunt et al., 2007), the OKR of Elateridae (Elaterophanes: Fossil 4 in Fig. 1) from Dorset, U.K., identified as Hettangian age (199.3-201.3Ma) in the Early Jurassic (Whalley, 1985;Hunt et al.,  1 and 2, Supplementary Fig. S1 and Supplementary Table S3.H. JIN et al. 2007), can be used as the minimum boundary.On the other hand, the OKR (Fossil 5 in Fig. 1) of the Coleoptera (Cupedidae) is known from Madygen, Fergana, Kyrgyzstan and identified as Ladinian age (237-242 Ma) in the Middle Triassic (Arnol'di et al., 1992;Hunt et al., 2007).This age can be tentatively used as the maximum boundary of this node.These calibrations are summarized in Table 1.
Evaluation of the prior distributions by Bayes factor Although the prior distributions of the divergence times (e.g., drift model for changing of the evolutionary rates along the branches, fossil calibrations) are important, normally these are evaluated by circumstantial evidence, especially consistency and inconsistency with the fossil record (e.g., Zhong et al., 2009;Wu et al., 2015).Lepage et al. (2007) suggested the application of the Bayes factor to evaluate drift models (autocorrelated model and independent model).Here, we expand this idea to evaluate not only the drift models, but also fossil calibrations by the Bayes factor with Tracer ver.1.6 (http://tree.bio.ed.ac.uk/software/tracer/).

RESULTS AND DISCUSSION
Phylogenetic tree The ML tree of the Dynastini is shown in Supplementary Fig. S1.The branching patterns are essentially the same as those of Rowland and Miller (2012)'s tree except for the basal position of Chalcosoma within the subtribe Chalcosomina in this study, and the nodal support values are generally higher.The evolutionary history of 17 morphological characters described by Rowland and Miller (2012) was traced on the basis of this ML tree topology (Supplementary Figs.S1 and  S2).Among the 17 characters, dorsal integument color (character 17 in Rowland and Miller, 2012) evolved most frequently (5 steps).Except for this character, male ornamentation (2.44 steps on average) evolved significantly faster than other characters (1.43 steps on average) (t-test: P = 0.037).It is not surprising that sexual selection should be the strongest driving force of the morphological variety in this tribe.At the same time, the reconstructed ancestral states suggest that convergent evolution frequently occurs in sexual ornamentation.
There was no significant correlation between the morphological evolutionary rates and the molecular evolutionary rates along the branches (data not shown).
Divergence times The result of the divergence time estimation is shown in Fig. 1 and Table 2.The autocorrelated model was selected for the drift model of the evolutionary rate (log Bayes factor (BF) ≈ 40: strong support), but selection of these two models gave a very limited effect on the estimates (Table 2).The choice of fossil calibrations did not have any considerable effect on the estimates.Exclusion of the calibration on the Scarabaeidae-Tenebrionidae split (node 1 in Fig. 1) made the marginal likelihood score higher (ln BF = 1.05 -1.39), but ln BF values of < 3 are considered negligible (Kass and Raftery, 1995).It suggests that the assumption of the maximum boundary on their split is not appropriate, and that the actual divergence time of this node (Scarabaeidae-Tenebrionidae split) is older than this maximum boundary.
To test the hypothesis that the subfamily Dynastinae emerged in the Oligocene, we also estimated the divergence times under each of the following two assumptions.The first assumption was that the divergence between the subfamily Dynastinae and its sister subfamily Cetoniinae occurred during the Oligocene, and the second assumption  was that the crown group of the subfamily Dynastinae emerged during the Oligocene (the split of Orizabus and other members of the Dynastini was assumed to have occurred in the Oligocene (23.0 Ma to 33.9 Ma)).These assumptions made the marginal likelihood scores lower (ln BF = 3.87 -4.80), indicating that the exclusion of such assumptions results in positive effects (20 > ln BF > 3: Kass and Raftery, 1995).These results suggest that our evaluation of the fossil calibration is objective.
The effect of taxon sampling on divergence time estimation was also investigated.Sparser taxon sampling typically had no considerable effect (except for the estimate of Dynastes and Megasoma in data set 2) (Supplementary Table S3).However, in one case, when the taxon samplings were dense in the outgroup and sparser in the ingroup, the estimates became grossly younger (please refer to Supplementary Table S3 and Supplementary Table S4 for details of the members of the ingroup and outgroup).Since the nucleotide sequence data of Dynastinae species other than Dynastini are insufficiently accumulated in public databases, the effect of taxon sampling in the outgroup should be examined in a future study.In every case of prior distributions in terms of rate drift models as well as fossil calibrations (except for the 'bad' assumption that the subfamily Dynastinae emerged in the Oligocene), the estimated divergence time (node 6 in Fig. 1,98.5 Ma, between Dynastina (African-New World clade) and Xylotrupina + Chalcosomina (Asian clade), which is also the time of the most recent common ancestor (tMRCA) of the tribe Dynastini, falls within the late Early Cretaceous to the early Late Cretaceous (94.1-114.0Ma).The divergence time between the African clade and the New World clade (node 18 in Fig. 1) was estimated to be 56.9-99.8Ma (the Late Cretaceous).Even considering the instability caused by taxon sampling, these dates are within the late Early to the Late Cretaceous, and the Late Cretaceous to the Paleocene, respectively.Since the breakup of Pangea  (Krell, 2000(Krell, , 2007)), as mentioned above, there are very long ghost lineages.It is notable that the phylogeographical branching pattern and divergence times are compatible with those of the Eutherian mammals (Fig. 1).Overseas sweepstakes dispersal just after the Pangean breakup can be considered to be a major driving force for the establishment of a distribution pattern.However, considering the large standard errors in our estimates for the separation between Asian and African + New World lineages (node 6, 76.6-120.2Ma) and that between African and New World (node 18, 56.9-99.8Ma) lineages, the breaking up of Pangea itself also should be taken into consideration.
According to Allaby (1999), sweepstakes dispersal "requires a major barrier that is occasionally crossed.Which groups cross and when they cross are determined virtually at random."Since the distances among Africa, South America and Eurasia in the Late Cretaceous were shorter than they now are (Nishihara et al., 2009), it is plausible that distribution areas were established by overseas sweepstakes dispersal across the short intercontinental distances during the Cretaceous.As several species of the Asian clade (Xylotrupes gideon, Allomyrina pfeifferi, or the genus Haploscapanes) are distributed in Australasia beyond the Wallace Line, if the distances between the landmasses were sufficiently short, these species could expand their distribution areas.However, at the same time, if such overseas dispersal occurred frequently, the phylogeographic pattern would be disturbed.Considering the very clear relationships between the phylogenetic branching patterns and geographic distribution patterns, such overseas dispersal events are likely to have been very rare.
Our results suggest that all major Dynastini lineages emerged in the Paleogene, especially during the Paleocene and Eocene (Fig. 1), despite their poor fossil record.The Paleocene and Eocene can be characterized as the warm period in the Cenozoic era (Zachos et al., 2001).Because most of the Dynastini species are adapted to tropical regions, it is plausible that they adapted to various empty ecological niches after the K-Pg extinction, and diversified at the morphological level thereafter.Interestingly, although morphological evolutionary rates did not show significant correlation with estimated time intervals (Supplementary Fig. S3), the distribution pattern of the changes in the time-scaled tree (Fig. 1) suggests that the rates accelerated episodically in the Paleocene to Eocene.These accelerated rates in morphological evolution are harmonious with the episodic adaptive radiation after the mass extinction.Explosive emergence of higher taxonomic groups after the mass extinction, at both molecular and morphological levels, is known in several clades such as the Eutherian mammals (molecular level: dos Reis et al., 2012;morphological level: O'Leary et al., 2013) and the Aves (molecular level: Jarvis et al., 2014;morphological level: Clark et al., 2005).These findings will contribute to a better understanding of faunal successions after mass extinctions.
Besides the tribe Dynastini, the subfamily Dynastinae includes the seven tribes Agaocephalini, Oryctini, Phileurini, Pentodontini, Oryctoderini, Cyclocephalini and Hexodontini.Some tribes (e.g., Oryctini, Pentodontini) are globally distributed, and some (e.g., Oryctoderini, Hexodontini) are regionally distributed.Phylogenetic relationships among these tribes are poorly understood (but see the morphological analysis by Endrödi, 1985).In addition, how their sexual dimorphism and ornamentation evolved, and how they established their distribution area, are interesting issues.Molecular phylogenetic analyses based on more extensive taxon sampling of this subfamily will shed more light on these issues in the future.

Fig. 1 .
Fig. 1.The result of divergence time estimation.A time-calibrated phylogenetic tree of the Dynastini as inferred from mitochondrial genes (COII and 16S) and nuclear genes (H3 and ArgKin).The time-scaled tree and divergence times were estimated with the autocorrelated model.The horizontal axis indicates geological time in Ma.Nodal bars indicate the 95% CIs of the divergence time estimations.Numbers in circles near the nodes are the individual nodal numbers of the Dynastini that are used in Tables1 and 2, Supplementary Fig.S1and Supplementary TableS3.The species names of Asian clade, African clade and American clade members are indicated in the pale green box, the pale blue box and the pale orange box, respectively.Nodal bars colored in purple indicate the divergence events among Asian, African and American clades, while those colored in green and in orange indicate the emergence times of the Asian and American clade, respectively.Details of Fossils 1-5 as well as the OKR of the Cetoniinae are described in the main text, and also correspond with Table1.Gray circles on branches indicate changes in morphological character states (see Supplementary Fig.S1for details).The outline of the time-calibrated phylogenetic tree of Eutherian mammals (dos Reis et al., 2012) is also shown for comparison, and the colors of the nodal bars are concordant with those of the Dynastini.Fossil 4 (Elaterophanes) belongs to the family Elateridae, and not the family Tenebrionidae; however, since Elateridae are phylogenetically closer to Tenebrionidae than to Scarabaeidae, Fossil 4 represents the Elateridae + Tenebrionidae clade, and is therefore indicated on the branch of Tenebrionidae.The image of Tribolium castaneum (top) is from Wikipedia.The copyright of the Dynastini images belongs to Mushi-sha Ltd., and the images are used with their permission.The photos of the Eutherian mammals (Panthera tigris, Loxodonta africana and Myrmecophaga tridactyla) were taken by M. Hasegawa.
Fig. 1.The result of divergence time estimation.A time-calibrated phylogenetic tree of the Dynastini as inferred from mitochondrial genes (COII and 16S) and nuclear genes (H3 and ArgKin).The time-scaled tree and divergence times were estimated with the autocorrelated model.The horizontal axis indicates geological time in Ma.Nodal bars indicate the 95% CIs of the divergence time estimations.Numbers in circles near the nodes are the individual nodal numbers of the Dynastini that are used in Tables1 and 2, Supplementary Fig.S1and Supplementary TableS3.The species names of Asian clade, African clade and American clade members are indicated in the pale green box, the pale blue box and the pale orange box, respectively.Nodal bars colored in purple indicate the divergence events among Asian, African and American clades, while those colored in green and in orange indicate the emergence times of the Asian and American clade, respectively.Details of Fossils 1-5 as well as the OKR of the Cetoniinae are described in the main text, and also correspond with Table1.Gray circles on branches indicate changes in morphological character states (see Supplementary Fig.S1for details).The outline of the time-calibrated phylogenetic tree of Eutherian mammals (dos Reis et al., 2012) is also shown for comparison, and the colors of the nodal bars are concordant with those of the Dynastini.Fossil 4 (Elaterophanes) belongs to the family Elateridae, and not the family Tenebrionidae; however, since Elateridae are phylogenetically closer to Tenebrionidae than to Scarabaeidae, Fossil 4 represents the Elateridae + Tenebrionidae clade, and is therefore indicated on the branch of Tenebrionidae.The image of Tribolium castaneum (top) is from Wikipedia.The copyright of the Dynastini images belongs to Mushi-sha Ltd., and the images are used with their permission.The photos of the Eutherian mammals (Panthera tigris, Loxodonta africana and Myrmecophaga tridactyla) were taken by M. Hasegawa.

1
Nodal numbers correspond to those of Fig. 1 and Supplementary Fig. S1. 2 Numbers of the fossils are concordant with Fig. 1. 3 ln Marginal likelihood scores were estimated by the harmonic mean.

Table 2 .
(Nishihara et al., 2009)es thought to have occurred around 110 to 148 Ma (Early Cretaceous), and more specifically 120 Ma(Nishihara et al., 2009), the initial diversification of the tribe Dynastini seems to have occurred during or just after the Pangean breakup.The tMRCA of the New World clade (node 19) is 46.6-87.6Ma, and that of the Asian clade (node 7) is 49.1-88.6Ma.This means that the tribe Dynastini had already diversified in the Late Cretaceous to the Eocene.Since fossils of the subfamily Dynastinae first emerged in the Oligocene is