If Dung Beetles (Scarabaeidae: Scarabaeinae) Arose in Association with Dinosaurs, Did They Also Suffer a Mass Co-Extinction at the K-Pg Boundary?

The evolutionary success of beetles and numerous other terrestrial insects is generally attributed to co-radiation with flowering plants but most studies have focused on herbivorous or pollinating insects. Non-herbivores represent a significant proportion of beetle diversity yet potential factors that influence their diversification have been largely unexamined. In the present study, we examine the factors driving diversification within the Scarabaeidae, a speciose beetle family with a range of both herbivorous and non-herbivorous ecologies. In particular, it has been long debated whether the key event in the evolution of dung beetles (Scarabaeidae: Scarabaeinae) was an adaptation to feeding on dinosaur or mammalian dung. Here we present molecular evidence to show that the origin of dung beetles occurred in the middle of the Cretaceous, likely in association with dinosaur dung, but more surprisingly the timing is consistent with the rise of the angiosperms. We hypothesize that the switch in dinosaur diet to incorporate more nutritious and less fibrous angiosperm foliage provided a palatable dung source that ultimately created a new niche for diversification. Given the well-accepted mass extinction of non-avian dinosaurs at the Cretaceous-Paleogene boundary, we examine a potential co-extinction of dung beetles due to the loss of an important evolutionary resource, i.e., dinosaur dung. The biogeography of dung beetles is also examined to explore the previously proposed “out of Africa” hypothesis. Given the inferred age of Scarabaeinae as originating in the Lower Cretaceous, the major radiation of dung feeders prior to the Cenomanian, and the early divergence of both African and Gondwanan lineages, we hypothesise that that faunal exchange between Africa and Gondwanaland occurred during the earliest evolution of the Scarabaeinae. Therefore we propose that both Gondwanan vicariance and dispersal of African lineages is responsible for present day distribution of scarabaeine dung beetles and provide examples.


Introduction
Recent explanations for the extraordinary diversity of insects, both at the species-and higher taxonomic levels, have largely centered on the role of co-diversification with flowering plants (angiosperms) [1].Significantly higher speciation rates within beetle groups associated with angiosperms compared to gymnosperm associated taxa [2] and correlated timing of crowngroup radiations of insect families with those of angiosperms [3][4][5][6] have been taken as strong evidence for the dependence of major insect radiations on angiosperm evolution.The majority of studies in this area, however, concern taxa with tight ecological associations with plants (i.e.herbivores or pollinators) and to date no model of the influence of angiosperm evolution on non-herbivorous insects has ever been proposed.Beetles are considered to be the most successful animal lineage in the world, occupying almost all terrestrial habitats, with the ~360-400,000 described species representing ~40% of known insect diversity [7].While numerous hyperdiverse lineages of herbivorous beetles exist, such as weevils (Curculionidae), leaf beetles (Chrysomelidae) and leaf chafers (Scarabaeidae: Melolonthinae), so do numerous hyperdiverse saprophagous and predatory lineages including rove beetles (Staphylinidae), ground beetles (Carabidae) and darkling beetles (Tenebrionidae).Even though non-herbivores represent ~40% of beetle genera [2], potential evolutionary drivers have largely gone unexamined (c.f.Hunt et al. [8]), but given that many saprophagous and predatory beetles do not have specialist diets (host or prey species) coevolution paradigms cannot easily be applied [9].
Scarab beetles (Scarabaeidae) are ideal for testing how lineages with different feeding biologies have radiated in response to major evolutionary events, as the family is composed of phytophagous and saprophagous lineages, with generalist and specialist adaptations to these ecologies.Boasting almost 30,000 species and 1,600 genera, Scarabaeidae is one of the largest beetle families but it is the subfamilies that form the clear division in predominate feeding ecology [10].The phytophages represent ~70% of the diversity and includes the richest subfamily Melolonthinae (11,000 spp.), that are predominantly leaf feeders, as well as the Rutelinae (4,000 spp.), Cetoniinae (3,300 spp.) and Dynastinae (1,500 spp.) which are more specialized feeders of fruit, flowers, pollen, sap, wood (including dead wood), and tubers, as well as leaves [10].The saprophagous group contains the Aphodiinae (3,300 spp.) that are more generalized detritivores feeding on dead or decaying matter including leaf litter, fallen logs and flowers, rotting fruits, mushrooms, dung and in some cases carcasses, and the Scarabaeinae (5,000 spp.) that are predominantly specialist dung feeders [10].Given the sheer diversity of species, it is unsurprising that exceptions to these general feeding biologies exist within lineages, however such specialist taxa represent a small proportion of diversity [10] with many of these feeding ecologies representing derived traits (e.g.pollen feeding Hopliini (Melolonthinae) [11] or necrophagy in Scarabaeinae [12]), thus making scarab beetles ideal for examining broad evolutionary trends.
The Scarabaeoidea first appear in the fossil record in the Middle Jurassic.Alloioscarabaeus cheni is a remarkably well preserved fossil from the Jiulongshan formation and is assigned to the superfamily due to the presence of antennae with lamellate clubs and dentate protibia with a terminal apical spur but due to unique wing development cannot be assigned to an extant family [13].The Jiulongshan formation, Daohugou Village, Inner Mongolia China is considered to represent late Middle Jurassic biota (~165 Ma) based on radiometric dating and biostratigraphic analyses [14][15][16][17].The age of Alloioscarabaeus provides a minimum age for the superfamily, however, early diversification of the Scarabaeoidea remains unclear due to gaps in the fossil record and the lack of other high quality fossils.Molecular studies have yet to yield consistent dates for the divergence of Scarabaeoidea, suggesting either a Lower Jurassic (~191.4Ma[8], ~174..9Ma[18]) or Lower Cretaceous (~141.11Ma [9]) origin.However, the latter estimate is younger than the age of Allioscarabaeus even when considering the 95% confidence interval (161.0-116.87Ma)[9] and may be underestimated.Further molecular evidence suggesting a Lower Jurassic origin of Scarabaeoidea can be taken from the phylogeny of stag beetles (Lucanidae), one of the earliest branching scarabaeoid families.The origin of the stag beetle crown group is estimated as ~160Ma, which is in line with a lower Jurassic origin of the superfamily [19].
To date, only one study has examined the evolution of scarab beetles (Scarabaeidae) using molecular divergence dating methods, and suggested that the timing of diversification tracks the sequential rise of angiosperms and mammals [18].The study found that the phytophagous lineage diverged ~108.9-128.1Ma(mean age of crown group tested using 6 alternative calibrations) while the major melolonthine tribes (leaf feeders) diverged in the Upper Cretaceous and the more specialized angiosperm feeding subfamilies diverged in the Paleogene.Due to these results, the diversification of plant feeding scarabs was attributed to the rise of angiosperms and a temporal lag [18].Within the saprophagous lineage, Ahrens et al. [18] estimated that the mean crown group ages of Aphodiinae and Scarabaeinae were 65.5-124.7 Ma and 86.6-100.2MArespectively, however no estimate of the origin of this saprophagous clade is explicitly given.Ahrens et al. [18] examine the origin of dung feeding in both lineages, and given that the timing was estimated at 72.6-85.5Mafor the Scarabaeinae and 43.4-78.1Mafor the Aphodiinae, they concluded that it was unlikely that dung feeding scarabs initially fed on dinosaur dung and was therefore associated with the rise of mammals.
The study of Ahrens et al. [18] provided the first insights in to the evolutionary history of scarab beetles, however, their 146 species, 4 gene data set could not recover a monophyletic Scarabaeidae.These results are in conflict with the results of a 2 gene, 282 taxa phylogeny of the Superorder Staphyliniformia which recovered a weakly supported sister relationship between the phytophagous and saprophagous scarab lineages [20].Despite conflicting evidence that the Scarabaeidae represents a single lineage, their monophyly is supported by numerous apomorphies including larval and adult morphological characters [10].As such we re-examine the phylogeny of scarab beetles to provide an independent test of the relationships and their divergence dates.We examine the effect of fossil calibration using two different fossil sets and the influence of maximum bound on estimated ages, allowing for a more rigorous test of the impact of potential co-radiations with angiosperms, dinosaurs and mammals on scarab beetle diversification.

Taxon Sampling
In total, 125 species from the superfamily Scarabaeoidea were collected from Australia and preserved in 96% undenatured ethanol.All necessary permits were obtained for the described study, which complied with all relevant regulations.Field permits were issued from the following State government institution: Department of Parks and Wildlife, Western Australia; Department of National Parks, Sports and Racing, and Department of Environment and Heritage Protection, Queensland; Department of Primary Industries; Office of Environment and Heritage, New South Wales; Department of Environment, Water and Natural Resources, South Australia, and Parks & Wildlife Service, Tasmania.No endangered or protected species were collected during this study.All voucher specimens are deposited at the Australian National Insect Collection, CSIRO, Canberra.Taxon sampling was relatively proportional to the diversity of the Australian scarabaeoid fauna.Additional scarabaeoid sequences representing the global fauna were downloaded from GenBank.Only 7 of 12 families recognized by Bouchard et al. [21] were included, however, the families absent from the study do not occur in Australia [22] so could not be collected during this study, and no sequences were available on GenBank at the time of analysis.The total data set comprised 450 taxa and included 45 nonscarabaeoid polyphagans and adephagans beetles as outgroups.A list of taxa used in analyses is provided in Table 1.

DNA amplification and gene sequencing
DNA was extracted from the head and thorax of specimens using a QIAGEN DNeasy tissue kit as per manufacturer protocols.Three mitochondrial genes, 16S rRNA, 12S rRNA and cytochrome c oxidase subunit I (COI), and the nuclear gene, 28S (LSU) rRNA, were amplified using primer pairs and amplification protocols as per Gunter et al. [23].PCR products were purified using EXOSAP-it (Affymetrix).Cycle sequencing reactions were performed using the BigDye Terminator v1.1 Cycle Sequencing Kit, the products of which were purified by alcohol precipitation and sent to the John Curtin School of Medical Research, Australian National University (ANU-JCSMR) for sequencing.Sequences were edited using Sequencher (v4.5;Gene Codes Corporation, Ann Arbor, MI, U.S.A.).Bidirectional sequences were aligned to form contigs and edited using Sequencher v. 4.5.

Multiple alignment and phylogenetic analysis
Once additional scarabaeoid sequences were downloaded from GenBank, sequences of each of the three fragments were aligned separately using default parameters of MAFFT [24] and Muscle [25](as implemented in Geneious ver.5.6 [26].Genes with limited or no length variability (28S, COI) were aligned with MAAFT as it produced a more reliable alignment (i.e. higher percentage pairwise identity and fewer inferred indels) than the MUSCLE alignments.The 16S-tRNA-Val-12S amplicon was aligned as three individual genes, 12S and tRNA-Val using MAFFT and 16S using MUSCLE as the high degree of length variability in this gene reduced reliability of MAFFT alignments.Each alignment was edited by eye before concatenation using Geneious into a final dataset spanning 4,584nt.Almost all taxa within the analysis contain data from at least three genes (96.8%) and 91.2% of taxa (n = 413) had complete gene coverage.
PartitionFinder [27] was used to determine the best partitioning strategy and nucleotide substitution models for phylogenetic inference.The optimal partitioning scheme divided the data into six partitions, separating each of the RNA genes except tRNA Val and 12S, and further subdividing COI by codon position.Bayesian Inference was conducted using MrBayes 3.2.1 [28][29].Each analysis consisted of 30 million generations with a random starting tree, and two simultaneous runs with four Markov chains sampled every 1000 generations were conducted with unlinked partitions.Stationarity in MCM chains was determined in Tracer [30] and burnin was set appropriately.A majority-rule consensus tree was obtained from the two combined runs to establish the posterior probabilities of clades.Maximum-likelihood analyses were performed using the RAxML [31] on the CIPRES portal [32] and the same partitioning strategies.

Fossil Selection (Table 2)
For age calibration, fossils and applied age constraints were selected based on recent reviews of the scarabaeoid fossil literature and relied heavily on the reviews of Krell [33][34] that critiqued           Table 1.Bolboceratinae is currently recognized as a subfamily of Geotrupidae, accordingly the fossils listed as "Bolboceratidae" and "Geotrupinae" are placed at these subfamily levels. doi:10.1371/journal.pone.0153570.t002 the reliability of taxonomic assignments made in initial fossil descriptions.The oldest scarabaeoid fossil, Alloioscarabaeus cheni, was not used as a fossil calibration point as it represents an extinct family, likely belonging to a primitive stem lineage [13], however, this fossil was used to cross-validate our results.The oldest unequivocal Scarabaeidae fossil, Juraclopus rohdendorfi cannot be unequivocally assigned to its proposed subfamily Aclopinae (as discussed in Bai et al. [35]), therefore both crown and stem positions for this fossil were tested.The phylogenetic position of Aclopinae is currently unknown and only Phaenognatha, which is currently classified in the Aclopinae (see Ocampo & Mondaca, [36]), is included in this analysis where it falls within the phytophagous clade.It has been suggested Phaenognathini may be a melolonthine tribe [22] but no classification changes have been made.Until Aclopus, the type genus, is sequenced the relationship between Aclopus and Phaenognatha will remain unclear and so we conservatively constrained the fossil within the phytophagous clade, instead of placing it on the Phaenognatha node.Accordingly, we tested two fossil calibration strategies for placement of the Juraclopus fossil: as either a crown or stem member of the phytophagous clade.For both fossil calibrated schemes (CS2: Juraclopus in the phytophagous crown-group, CS3: Juraclopus in the phytophagous stem-group), the root age of Polyphaga was constrained at 268-273 Ma on the basis age estimate of Polyphaga from Hunt et al. [8], which was deemed to be more consistent with the scarabaeoid fossil record than the estimate from McKenna et al. [9] (see above).These fossil calibration points were assessed against the fossil calibrations used by Ahrens et al. [18].The rationale behind testing alternative calibration schemes was to determine if the differences in age estimates between the present study and Ahrens et al. [18] were due to (i) the dataset and analysis methods, (ii) fossil calibration and priors, (iii) fossil selection, or (iv) influence of a molecular calibration point.The seven calibration schemes are summarized in Table 2 along with the details of the fossils used and the nodes constrained.

Calibration schemes
To explore the above divergence dating hypotheses, the following calibrations using the fossils selected by Ahrens et al. [18] were applied to our 445 taxon, 4 gene dataset: 1. hard bounds (minimum and maximum age) for all Ahrens fossil calibrations (excluding Aegialia which was not present in our taxon sampling) (CSi).Minimum and maximum ages reflect the upper and lower bounds of the dated fossil layer interval respectively (see Table 2).A direct comparison can be made to the results of Ahrens et al. [18] based on their 146 species, 4 gene data set providing a test of effects of dataset or analysis method.
Although hard bounds are stricter than the exponential priors used in Ahrens et al. [18], these comparison analyses serve to demonstrate the differences due to analysis calibration and not taxon or gene sampling 2. To test the effect of hard bounds on fossil calibrations, in particularly on Cenozoic fossils, we compare the results of CSi to the following two calibrations: CSii) using only Cretaceous fossils as calibrations with hard bounds (minimum and maximum age); CSiii) Cretaceous fossils only, hard bounds (minimum and maximum age) except for ingroup fossils which were only constrained with a minimum age.
3. To test if fossil selection is responsible for differences in age estimates, the results from the CS2 and CS3 can be compared to CSiii which all do not have Cretaceous or Cenozoic upper bounds applied within the ingroup.
4. The results from the program r8s are known to be influenced by the oldest calibration applied.All analyses of Ahrens et al. [18] only include fossil calibrations, the oldest dated from the Middle Jurassic.To determine if using the molecular age of the Polyphaga as a calibration point significantly influenced the estimated ages in CS2-CS3, we compare results of CSiii and CSiv which are identical except CSiv includes the molecular calibration point for the Polyphaga.Additionally, the tree was calibrated only using the root age of Polyphaga, 268-273 Ma as estimated in Hunt et al. [8], to provide a best guess estimate independent from fossil calibrations (CS1).

Divergence time estimates
Divergence date estimates were calculated for all seven calibration schemes using a penalized likelihood (PL) method and we further tested four of these schemes (CS2, CS3, CSi and CSiii) under a Bayesian scenario using MCMCTree.PL divergence estimates were calculated with a truncated Newton algorithm using the software r8s, ver.1.7.1 [37].The consensus tree with branch lengths obtained from Bayesian analyses was used as a fixed input tree.The optimal level of rate smoothing was determined by cross-validation and a smoothing parameter was set accordingly.The use of PL in divergence dating analyses has been criticized as not as robust as newer, Bayesian methods.Bayesian divergence dating programs like BEAST [38], Multidivtime [39], and MCMCTree [40] have advantages such as being able to deal with uncertainties in topology and branch length and can allow for different types of prior distributions.Validation of our divergence estimates obtained from r8s was initially attempted in BEAST but given the size of the data set BEAST struggled to identify an appropriate seed tree and the analyses were aborted because the initial state had zero probability.This can occur if "the log likelihood of the tree is -Inf, [which] may be because the initial, random tree is so large that it has an extremely bad likelihood which is being rounded to zero" [38].Even when a fully bifurcating tree was used as input, BEAST continued to reject the log likelihood of the tree.Given the failure of the analysis to launch in BEAST, validation was next attempted in MCMCTree implemented in the PAML package [40].MCMCTree has been shown to handle analyses with much larger numbers of taxa than either BEAST or Multidivtime while still producing comparable time estimates [41][42].An estimation of substitution rates using a GTR model was calculated from the sequence alignment and input topology (the same output tree from the Bayesian analysis used in PL) with BASEML in the PAML package and used to adjust the rgene_gamma parameter accordingly in the control file.The Sigma2_gamma parameter was also set appropriately according to the estimated age between the tip and root of the tree (using the molecular divergence date estimated in Hunt et al. [8]).The gradient and hessian of branch lengths were calculated in MCMCTree and used in the final estimation of divergence times by Bayesian analysis using a modified fully bifurcating version of the original tree.The MCMC parameters were set so the first 50,000 iterations were discarded as burnin, and then sampled every 50 iterations until it had gathered 10,000 samples (burnin = 50,000 sampfreq = 50 nsample = 10,000), in total running for 550,000 iterations.All analyses were made with GTR as the model of substitution, with a flexible birth-death prior used (BDparas = 2 2 0.1).More uniform priors were also tested (BDparas = 1 1 0) but did not affect the posterior estimates significantly.Fine tuning parameters were set at 0.06, 0.5, 0.006, 0.12 and 0.4 for time, substitution rate, mixing, and substitution parameters, respectively.Acceptance proportions were checked after fine-tuning to ensure they ranged around ~30% and were always in the most reliable range of between 20-40% for MCMC analyses.

Diversification and speciation estimates
The chronograms resulting from analyses were used to generate lineage through time (LTT) plots using the packages gieger, ape, TreeSim and laser [43][44][45][46] in R. LTT plots were generated for specific clades through exclusion of the remaining taxa and pruned in Mesquite [47].For the purpose of the generalist versus specialist angiosperm feeding analysis the tribe Hopliini is included within the Melolonthinae (generalist, leaf feeding) LTT on the basis of current classification and while Hopliini feed on leaves and pollen [10], pollen feeding has been demonstrated to be a derived trait [11].LTTs were plotted in both linear and log scales as comparison of rates between groups with differing number of taxa is easier to visualize with the linear scales.
Clade-wide rates of diversification were estimated using the method of moments (MoM) detailed by Magallón & Sanderson [48] with three relative extinction rates (ε = 0/0.5/0.9).Species richness data for each clade is taken from Scholtz & Grebennikov [10] and the estimated age of each clade from PL analysis CS3.Relative diversification rates between clades did not differ depending on which of the calibration schemes were used.To examine speciation, extinction, net diversification and net extinction rates in a Bayesian framework, as well as potential mass extinctions, some of the resulting chronograms were assessed using the package TESS in R [49].Given that two thirds of genera were represented by a single species, we excluded the most recently diverging, species level sampling (45 or 50Ma depending on the analysis) to result in a topology with approximately 120-131 tips (representing ~2/3 tip sampling).A sampling fraction (rho) was then calculated based on the number of tips and nodes of each of the topologies at either 45 or 50Ma (depending on the analysis) based on the predicted diversity calculated using MoM detailed above (i.e. with no extinction Scarabaeinae richness was estimated as 1619 or 1457 at 45 and 50 Ma respectively).A summary of the analyses, excluded tips and sampling fractions are listed in S1 Table .A series of preliminary TESS / CoMET analyses was performed for each topology (with diversification set as 0.068635 for Scarabaeinae or 0.066875 for the Pleurosticti, as predicted by MoM), and held constant and extinction set at 0), to estimate the marginal posterior probability densities for the diversification-rate hyperpriors mu and sigma for both initial speciation and extinction rates for final TESS analyses.All TESS and CoMET analyses were run under the following conditions with unique combinations of prior settings for mu, sigma and rho: number of expected rate changes = log(2), expected survival probability of 25% (beta prior with shape parameters α = 0.25 and β = 0.75) consistent with prior estimated loss of 70-75% of contemporaneous terrestrial species [50][51], number of mass extinctions = 1, maximum number of MCMC iterations = 100000 with 4 chains, burnin = 50000, thinning each chain by sampling every 10 th state.MCMC simulations were terminated when the effective sample size reached 500 and we explored the maximum number of iterations and thinning to ensure ESS would reach 500 and the burnin and thinning setting were appropriate.
Within the Saprophagous clade, the aphodiines formed the earliest branching lineages representing Rhyparini (Rhyparus sp.only), Eupariini + Proctophanini + Psammodiini (PP = 0.84), and Aphodiini (PP = 1).Scarabaeinae was recovered as monophyletic (PP = 0.97) with Sarophorus + Coptorhina, and Odontoloma+ Dicranocara forming the earliest branching lineages followed by a comb containing most of the tribal diversity (PP = 0.97).Many monophyletic lineages representing the tribes were recovered while the tribes Coprini and Canthonini were recovered as paraphyletic (see Fig 1).Overall the topology is largely congruent with the traditional relationships based on morphology and recovers the highest molecular support for a monophyletic Scarabaeidae observed to date (PP = 0.78).

Divergence time estimates
We found a high degree of congruence between divergence dates estimated by penalized likelihood (PL) or Bayesian methods, and between similarly constrained analysis priors, independent of fossil choice or data set (Table 3).Fig 2 plots the accumulation of ingroup taxa (Scarabaeidae) using seven different calibration schemes analysed by PL and four calibration schemes analysed by Bayesian methods and highlights that the results of analyses CS2, CS3, CSiii and CSiv are highly congruent.
(i) Differences between Penalised Likelihood and Bayesian analyses.MCMCTree and r8s performed slightly differently due to the constrained topology i.e.MCMCTree requires a fully bifurcating tree as final input for age estimation.Fig 2B plots the accumulation of ingroup taxa (Scarabaeidae) of r8s and MCMCTree analyses and highlights the differences due to tree topology with Bayesian methods producing a smoother curve because of the absence of polytomies and diverging earlier yet accumulating more slowly than PL.Differences in the age estimates between r8s and MCMCTree are summarised in Table 3 and S2 and S3 Figs.In general age estimates from MCMCTree analyses were older than those estimated in r8s but not systematically so.The most marked difference was the age of the clade that corresponds to the large radiation of dung feeding tribes, which was unresolved in the underlying phylogenetic analysis.Interestingly the ages predicted for the other significant unresolved comb in the Melolonthine lineage was also significantly older in the MCMCTree analysis.This suggests that forced input of a fully bifurcating tree is not ideal, biasing results towards older age estimates.However, age estimates on resolved clades were in very close congruence between both methods, adding confidence to the overall performance of our PL analysis.
(ii) Differences between data set.To test if our dataset, alignment or dating methods significantly influenced the results, age estimates from CSi were directly compared to Ahrens et al. [18] as this calibration scheme was designed to replicate their constraints used in run 1 and 2. Results for clades ages listed in Ahrens et al. [18] are compared in Table 4. Overall results from CSi are similar to run 1 and 2 of Ahrens et al. [18], highlighting that our 450 taxa dataset analysed using r8s and MCMCTree produces comparable age estimates to Ahrens et al. [18] 146 taxa dataset analysed in BEAST, when the same fossil calibrations are used.
(iii) Differences between fossil selection.Results from CSiii can be compared to CS2 and CS3 to examine the influence of fossil selection.These three calibration schemes use only Mesozoic fossil calibration points and only outgroup taxa are constrained with maximum bounds, providing an opportunity to cross-validate fossil selection.Cross-validation of fossils is discussed in further detail below but in general, regardless of fossil set, when analyses are constrained using similar bounds the results are highly congruent for the majority of clades in both PL and Bayesian age estimates, ruling out major conflict between selected fossil sets (see Table 3 and Fig 2).
(iv) Use of a molecular age calibration point.As a test independent from fossil calibrations, CS1 was only constrained by the age of the Polyphaga estimated by Hunt et al. [8]; the resulting age estimates were older than those recovered when calibrated with fossils.Although the fossil record provides a direct timescale for the appearance of taxa, it is incomplete particularly for hexapods which are represented by a relatively poor fossil record [52].Given the poor fossil record of insects, it is likely that additional fossils from the earliest scarab lineages have yet to be discovered or remain unpreserved, therefore the estimates from CS1 are plausible if considering a gap in the fossil record (see Table 5).However, this analysis was conducted only as an upper most prediction from which we could cross-validate our fossil selection.
To test the effect of using the molecular age estimate of the Polyphaga from Hunt et al. [8] as a bound on PL analyses, the calibration scheme CSiv was designed to replicate CSiii with the addition of the Polyphaga being bound by the molecular age estimate above.The results demonstrate minimal differences in age estimates in the ingroup but age estimates were always older with the Polyphaga constrained (see Table 6).For example the estimated ages in analyses CSiii and CSiv for Scarabaeoidea were both in the Sinemurian (192.39Ma or 197.89Ma respectively), the Scarabaeidae in the Aalenian/Toarcian (170.36Ma or 174.78 Ma); the phytophagous clade in the Tithonian/Kimmeridgian (151.25Ma or 155.56 Ma) and the saprophagous clade in the Tithonian (147.25Ma or 151.52Ma).These results suggest that applying this bound at the root of the Polyphaga in CS2 and CS3 would have minimal effect on age estimations, allowing for comparison between calibration schemes.

TESS analyses
No significant changes in speciation or extinction rates were detected in any of the TESS/ CoMET analyses conducted, including investigation of scarabaeine dung beetles or phytophagous Pleurosticti clade (see S4 Fig) .Within the Scarabaeinae, a potential mass extinction event in the Upper Cretaceous ~85-95Ma was detected in PL analyses of CS2, CS3 and CSiii but never in the MCMCTree analyses.We believe these differences are driven by topologies as the PL tree contained polytomies.

Systematics and Diversification
Divergence time estimates: Scarabaeoidea.The superfamily Scarabaeoidea was recovered as a strongly supported monophyletic lineage consistent with other studies [8-9, 18, 20].Given that Allioscarabaeus cheni represents an extinct family-level lineage within the superfamily, it was not used as a calibration point however is useful to cross-validate the age estimates from our analyses.We tested 7 calibration schemes, five of which (CS2-3, CSii-CSiv) predicted the mean age of origin to be in the Lower Jurassic (~176-203Ma), consistent with the age of Predicted ages of selected scarab lineages using two different fossil sets for calibration analyzed using penalized likelihood methods.See also Table 4 for additional node ages and Bayesian estimates of selected schemes and Fig 2 for differences in accumulation of ingroup taxa for all analyses.Generalized feeding biology is also listed.Exceptions to these generalized adult diets exist but only represent a small proportion of diversity and derived feeding ecologies doi:10.1371/journal.pone.0153570.t006 Allioscarabaeus cheni [13] and inline with the estimates of Hunt et al. [8] (~191.4Ma)and Ahrens et al. [18] (~174.3-190.9).Analysis CS1 was only calibrated with the estimated molecular age of the Polyphaga and predicted an Upper Triassic origin (~216Ma), while the most conservative analysis CSi with maximum bounds on Paleogene and Cretaceous fossils, predicted an Upper to Middle Jurassic origin (~160-170Ma).To date only the analysis of McKenna et al. [9] has predicted a Lower Cretaceous origin of Scarabaeoidea of 141.11Ma (161.0-116.87Ma),however this is in conflict with the Jurassic fossil record, which documents that at least three extant families (i.e.Lucanidae from Shara-Teg formation, Anda-Zhuduk and Daohugou Village Mongolia [53][54][55], Hyborsoridae from Karatau, Kazakhstan [56] and Alloioscarabaeidae from Jiulongshan formation, Mongolia [13]).
Given that calibration schemes CSi-iv use these lucanid and hyborsorid fossils as constraints we cannot cross-validate these node ages, however CS2-3 did not use these calibration points.The crown age of the Hyborsoridae in CS2-3 is estimated to have originated in Middle to Lower Jurassic (~172-177Ma) while the subfamily Hyborsorinae originated in the Upper Jurassic (~145-149Ma) (see Table 5) congruent with the oldest recorded Hyborsorinae fossils (Protohybosorus grandissimus from the Jurassic (155.7-164.7 Ma) [56] and Fortishybosorus ericeusicus from the Lower Cretaceous, Yixian formation [57]).Evidence for early hyborsorid diversification may be taken from the record of Mesoceratocanthus tuberculifrons (Hyborsoridae: Ceratocanthinae) from the Lower Cretaceous, Yixian formation [58].Our CS2-3 estimates predict the family Lucanidae originated in the Lower Jurassic (174-180Ma) and subfamilies began to diverge in the Middle to Lower Jurassic (~150-165 Ma) (see Table 5), these results are consistent with ages of known fossils of three lucanid subfamilies: Paralucanus mesozoicus (Paralucaninae) from Shara-Teg (145.5-150.8),Protolucanus jurassicus (Protolucaninae) from Anda-Zhuduk (145.5-150.8Ma)and Juraesalus atavus (Aesalinae) from Daohugou Village (159.8Ma) [54][55].Our results are also congruent with the estimates from a molecular phylogeny of the Lucanidae [19] which estimated the origin of the crown group Lucanidae sensu strictu at ~160Ma (154-171Ma) and the crown group Lucanidae sensu lato (includes Lucanidae s. s.Aesalinae, Diphyllostomatidae and the extinct lucanid subfamilies) at 167Ma (155-182Ma).Given that the fossil age of Protohybosorus grandissimus, Paralucanus mesozoicus, Protolucanus jurassicus and Juraesalus atavus are largely congruent with our CS2-CS3 estimates, this cross-validation highlights that our results do not appear to be grossly inflated.It also suggests the age of the Scarabaeoidea is in the Lower Jurassic as suggested by Hunt et al. [8], Ahrens et al. [18], and here in analyses CS2-3, and CSii-iv.
Divergence time estimates: Scarabaeidae.Ahrens et al. [18] did not recover a monophyletic Scarabaeidae so the timing of the origin of the family cannot be compared, however with the exception of CSi which was designed to replicate their analyses, their results consistently estimated much younger ages for major lineages within the Scarabaeidae than in the present study.The major differences between analyses lie in fossil selection and calibration method, i.e. selecting Cenozoic fossil and using exponential priors with "soft" maximum bounds.Exponential priors are best used when "there is a strong expectation that the oldest fossil lies very close to the divergence event being represented by the node, relative to a distant, soft maximum" [59].Given that the majority of fossils used in Ahrens et al. [18] were relatively young (Cenozoic) and were placed at nodes representing the most conservative interpretation of the taxonomic hierarchy (i.e.subfamily or tribe), which may have significantly underestimated the divergence dates when constrained using exponential priors with soft maximum bounds.To examine the effect of exponential priors with soft maximum bounds, we can compare analyses constrained with (CSi) and without (CSiii) maximum bounds on the ingroup taxa (Tables 3 and 6).The mean age of Scarabaeidae is estimated to have originated in the Aptian (113-125Ma) when maximum bounds are applied (CSi), or in the Aalenian (170.3-174.1Ma)when maximum bounds are excluded (CSiii), even considering the 95% confidence intervals of the Bayesian analyses there is no overlap between estimates (CSi: 99.2-156.7Ma;CSiii: 162.6-185.7Ma).This pattern is consistently recovered across ingroup clades, with lineage ages always estimated as significantly younger if maximum bounds are applied (CSi) than when only a minimum bound is applied (CSiii).Similarly the effect of removing calibrations with maximum bounds can be observed within Ahrens et al. [18], which excluded calibration points in some runs if they were identified as being in conflict with other settings in their BEAST analyses.For example the Aphodiinae calibration point "F" in runs 2-6 and the Aphodius calibration point "H" in run 4 were excluded, as a result the estimated age of Aphodiinae differs between run1 = 65.5Ma(57.3-73.3)and run2 = 69.2Ma(58.4-79.9)(includes calibration points F and H), and run3 (F is excluded) = 111.8Ma(94.8-129.7)or run4 (F and H are excluded) = 124.7Ma(108.0-142.7)[18].
Of the 10 ingroup fossils selected by Ahrens et al. [18] used to calibrate CSi, only one is Cretaceous, while nine are Cenozoic (1 Paleocene, 7 Eocene and 1 Miocene).Although these fossils are unequivocally assigned to the corresponding taxa, uncertainty exists as to whether they all represent the earliest members of their respective lineages, thus raising concerns regarding calibration using maximum bounds.Ahrens et al. [18] do not address the mismatch between their age estimates and other known reliable scarab fossils reviewed by Krell [33].When other known scarab fossils are considered, mismatches between node estimates inferred by Ahrens et al. [18] and in Csi are observed.For example, Ahrens et al. [18] recover an overall mean age of the tribe Sericini as ~93.4-111.3Ma(95% confidence intervals 81.3-124.7Ma)depending on various constraints including stem and crown placement of a Nearctic Sericini fossil from the Oligocene.However, the tribe Sericini is represented by multiple species from two genera in the Lower Cretaceous in the Zazunskaya and Leskovskaya Formations of Russia (~140-145Ma) and next appear 100my later in the fossil record, at Florissant, USA (33.9-37.2Ma)[33].Our estimate in CSi, also calibrated with this Oligocene fossil, recovers an even younger origin of the Sericini (88.39Ma (PL) or 74.42Ma(63.4-84.0)(Bayesian))inconsistent with the Cretaceous fossil record.Similarly, the age estimates for Melolonthinae recovered in CSi or Ahrens et al. [18] are not consistent with the fossil record.Given that Melolonthinae is recovered as a paraphyletic lineage within a monophyletic Pleurosticti clade, the crown age for the Pleurosticti can be used as a guide for the maximum age of stem group melolonthines.Ahrens et al. [18] estimate the crown age of the Pleurosticti clade in the Lower Cretaceous (mean age 108.9-128.1Ma;95% confidence intervals 95.9-142.1Ma)while our CSi estimates are again younger (101.4Ma(PL) or 90.86Ma (83.6-99.4Ma)(Bayesian)).Only when the maximum estimate from the 95% confidence intervals are considered, is one of the six runs of Ahrens et al. [18] congruent with the age of Cretomelolontha transbaikalica from the Baissa formation Russia (140.2-145.5Ma)[60].As is the case for tribe Sericini, a gap in the fossil record for the subfamily Melolonthinae exists, the next oldest melolonthine fossils are known from the Eocene (33.9-56.0Ma)[33].These extreme gaps in the fossil record highlight that the fossil record is poor for the Scarabaeidae and it is unlikely that the earliest branching members of major lineages have all been discovered.It is interesting to note that only trace fossils and no body fossils of scarabs are known from the Upper Cretaceous (66.0-89.8Ma),despite being represented in the middle Cretaceous by Aphodiinae (Cretaegialia spp.) and Scarabaeinae (Prionocephale deplanate) and in the Lower Cretaceous by Melolonthinae including, but not limited to, Sericini [33].Paleocene scarabaeoid fossils are also rare and represented only by Aphodiinae (Aphodius charauxi) and possibly Rutelinae (i.e.Anomalites fugitivus is reported from the Tertiary without specification [61]).Given the rarity of scarabaeoid fossils, including the absence of body fossils from the Upper Cretaceous, we believe that applying maximum bounds on fossils (particularly those from the Cenozoic) are inappropriate and we therefore consider the results of Ahrens et al. [18] and CSi to be underestimated.
When we consider all the other analyses inferred without maximum bounds applied to Cenozoic fossils (CS2, CS3, CSii, CSiii and CSiv) the age of the Scarabaeidae is consistently estimated to be in the Jurassic (see Table 6).Of these five analyses, only CSii has a maximum bound applied at a subfamily level (at the Scarabaeinae node: 83.5-92Ma); in all other analyses this node is only constrained by the minimum age bound.Consequently, the predicted ages of ingroup nodes are also always younger in CSii than for analyses without the maximum bound on the Scarabaeinae.Given the reported mismatch in estimated ages and the earliest known fossil, we consider the results of CSii are also likely to be an underestimate driven by an inappropriate analysis constraint.For the purpose of examining potential diversification hypotheses we consider the results of CS2, CS3, CSiii and CSiv to be plausible on the basis of crossvalidation of molecular estimates and the known fossil record, and because estimates are not driven by maximum bounds applied to ingroup taxa.Although the age of Scarabaeidae is constrained in CS3 (152-158Ma) the results are congruent with the results from other analyses and so we consider them plausible.
We tested two different fossil sets to calibrate the analyses, one based on the literature reviewed by Krell [33] and the second deemed to be unequivocally assignable by Ahrens et al. [18].Key morphological characters used to identify scarabaeoids are often poorly preserved in fossils so their exact systematic placement remains unclear, with various researchers taking more or less conservative views on their classification.One such example is the fossil of Juraclopus rohdendorfi Nikolajev, which was originally assigned to the tribe Aclopini (Scarabaeidae: Aclopinae) [62].Krell [33][34] accepts this placement and as such considers this the oldest known fossil from the family Scarabaeidae.Although Ahrens et al. [18] consider many of the fossils reviewed by Krell [33][34] to be unequivocally assignable, they provide no reasoning as to why Juraclopus is excluded from their fossil constraints.This fossil may have been excluded on the basis of uncertain placement within the phylogeny (Scarabaeidae is not monophyletic and no aclopines were included in their taxon sampling) or doubts over the identification of the fossil.Our analyses recover a monophyletic Scarabaeidae and include an aclopine that is recovered in the Pleurosticti clade, as such we test the placement of Juraclopus rohdendorfi at both the crown (CS2) and the stem (CS3) of the Pleurosticti clade.Interestingly, the estimated age for the origin of Scarabaeidae when not constrained with Juraclopus rohdendorfi (CSiii-CSiv) is highly congruent with the estimate obtained from CS2 suggesting that the fossil represents a member of the crown group Pleurosticti and also highlights that analyses are coming to a consensus for divergence estimates that predict the origin of the Scarabaeidae in the Aalenian-Toarcian (170.).
Another such fossil with debated identification is Prionocephale deplanate Lin, while Krell [33][34] and Ahrens et al. [18] list the fossil as assignable to Scarabaeinae, Tarasov & Génier [63] are less convinced.Tarasov & Génier [63] comment that the original description does not list any unambiguous characters of the subfamily, however, no specific explanation for or against the original identification is provided.It is clear further review of such Mesozoic scarabaeoid fossils is warranted; as such it is important to examine evolutionary hypotheses beyond the known and limited fossil record.Here we compare the timing of well-documented ecological events such as angiosperm and mammal evolution, continental drift and the Cretaceous-Paleogene mass extinction to explore our hypotheses on scarab beetle diversification.
Scarab diversification.The Scarabaeidae is divided into two strongly supported lineages representing major ecological groupings of phytophagous and saprophagous scarabs.Within the phytophagous lineage (the Pleurosticti), Melolonthinae is recovered as a paraphyletic grade with subclades representing major tribes, while Rutelinae, Cetoniinae and Dynastinae stem from a common ancestor (the CRD-clade) that diversified into the currently recognized subfamilies and tribes.Within the saprophagous clade, Aphodiinae form the most-early branching lineages and a monophyletic Scarabaeinae is supported.These results suggest both specialist phytophages and saprophages evolved from generalist ancestors, providing a compelling framework to examine feeding biology and the response to major biotic and abiotic evolutionary events.The divergence dates estimated in CS2-3 and CSiii-iv infer that the crown group Pleurosticti originated in the Upper Cretaceous.While the crown group of the saprophagous clade originated around the Jurassic-Cretaceous boundary, specialist adaptations to both of these feeding biologies originating in the middle of the Cretaceous.These age estimates significantly predate the mid-Cretaceous crown group origins of both clades predicted by Ahrens et al. [18], which as discussed above, may be an artifact of inappropriate calibration.Ahrens et al. [18] proposed that the evolution of scarab beetles tracks the sequential rise of angiosperms and mammals, however our results suggest an origin of the phytophagous lineage predating the ecological dominance of angiosperms and that specialist saprophages vastly predate mammal dominance.

Influence of Angiosperms on diversification
While estimates for the age of angiosperms vary, placing origins in the Jurassic or Triassic: 141-199Myr [64], 182-270Myr [65], 225-240Myr [66], or 221-275Myr [67], agreement exists that by the Albian-Turonian (90-110 Ma), angiosperms dominated the environment [68].Our estimates (CS2-3, CSiii-iv) predict the origin of the crown group Scarabaeidae in Middle to Lower Jurassic, with initial diversification of phytophagous scarabs at the Jurassic/Cretaceous boundary (see Fig 3).The oldest melolonthine lineages diverged ~130-150 Ma, suggesting that their initial evolution was either with the very earliest angiosperms or non-angiosperms as food plants.Their peak radiation was , lagging the origin of crown group angiosperms as has been reported in other leaf feeding beetles [5,69], including scarabs [18].In contrast, the CRD-clade (Rutelinae+ Cetoniinae + Dynastinae) that feed on angiospermspecific plant tissues originated 128-140Myr and split into the recognized subfamilies 92- 127Myr.The CRD-clade diversified at an accelerated rate from 90-105Myr (Fig 4B ), then accumulated at a similar rate to the leaf feeders until 70-80Myr, after which point they accumulated more rapidly than melolonthines.Ahrens et al. [18] estimate the mean age crown group origin of the Pleurosticti ~110-130Myr, with the specialist CRD clade originating ~90Ma (interpreted from Figure 1 [18] as not explicitly listed) and diverging into the recognized subfamilies ~46-72Ma.Although both our analyses and those of Ahrens et al. [18] suggest a co-radiation with angiosperms, our results predict a much smaller time lag, over 20Myr earlier than previously hypothesised with the timing of the first inferred peak of diversification of phytophagous scarabs corresponding with the ecological dominance of angiosperms ~90-110 Ma [68] (see Fig 3).
Interestingly, saprophagous scarabs underwent a similar diversification pattern to that observed in the phytophagous subfamilies, originating prior to or with the earliest angiosperms, and underwent a significant radiation 90-110Ma, corresponding to the major expansion of angiosperm diversity (Fig 3).This pattern is largely driven by the dung-specialist Scarabaeinae due to our limited aphodiine sampling.Ecologically, the main terrestrial vertebrate fauna of the Cretaceous were dinosaurs.Given that angiosperms only represented a small proportion of the floral diversity prior to the Albian (>113Ma), it is likely that angiosperms formed only a minor component of Lower Cretaceous dinosaur diet [70].The oldest direct evidence of angiosperm feeding is from the fossilized gut content of an Albian ankylosaur containing various angiosperm fruit [71], however coprolites and enterolites are rare [70].The incorporation of angiosperm tissue in to the diet of dinosaurs is supported by evidence from teeth and postulated jaw function of herbivorous dinosaurs, however it was likely that gymnosperms remained a significant component of their diet until their extinction [70,72].While faeces of herbivorous dinosaurs were certainly available before angiosperms originated, it may have been unsuitable for scarab feeding.Studies on scarabaeine dung beetle feeding show particle size is a critical determinant of dung selectivity [73].The incorporation of less fibrous angiosperms into the diet of herbivorous dinosaurs with grinding gastroliths [74] could have produced sufficiently fragmented dung suitable for scarab feeding.Non-angiosperm plants such as gymnosperms, ferns and allies are also of lower nutritional quality [75] so changes in floral composition were likely reflected in the available dung.We hypothesize that even if angiosperms were only a minor component of dinosaur diets, the less fibrous, more nutritious dung provided a new niche in which beetles first evolved.It also provides a model whereby the radiation of a non-phytophagous beetle group is driven by rise of angiosperms despite their ecological link being at one remove, i.e. through a dinosaur's digestive system.

Origin of dung feeding
The co-evolution of scarabaeine dung beetles with dinosaurs is controversial.Trace fossils within herbivorous dinosaur coprolites attributed to tunneling dung beetles appear in the Upper Cretaceous, the oldest from the Campanian [72,76], younger than the oldest Scarabaeinae fossil (Prionocephale deplanate Lin) [33].Conversely, it has been argued that the mixing of uric acid with feces common to all non-mammals would have made dinosaur dung unsuitable [77], a hypothesis supported by the rarity of extant dung beetles utilizing reptile or avian dung.Molecular-clock based estimates for the crown age of marsupials and placental mammals are as old as 82 and 101Ma, with peak mammal diversification between 75-85Ma [78].Fossil evidence however is inconsistent with the molecular-clocks, suggesting a mammal crown group age of 65Ma and peak diversification over the next 10Ma [79].Even if fossil evidence underestimates the age of mammals, they were not a significant component of the Cretaceous fauna and our evidence shows that both the origin and crown group radiation of most dung beetle tribes (~119-130Ma and 85-98Ma respectively) predates the major mammal radiation.Furthermore, it is unlikely that coprophagy by scarabaeines evolved in association with the small insectivorous stem-group mammals present in the Mesozoic [78] as very few extant scarabaeines feed on small dry dung pellets with only specialists using insectivore dung as a resource [73].
A strong mismatch between the inferred divergence times of both mammals and dung beetles provides the first molecular evidence that the initial evolution of coprophagy in scarabaeines was associated with dinosaurs and not mammals.Ahrens et al. [18] estimated the mean age of dung feeding Scarabaeinae to be ~72-85Ma, thus attributing their origin to co-evolution with mammals alone.However our reanalysis suggests that inappropriate fossil calibrations with maximum bounds are responsible for these younger ages.Historical biogeographical scenarios also support our significantly older age of scarabaeine dung beetles than the past molecular estimates.Two of the oldest dung feeding tribes, Dichotomini and Canthonini, have strong Gondwanan distributions [63,[80][81][82], resulting from either vicariance or dispersal.Our results suggest that with the origin of dung feeding Scarabaeinae is congruent with the continental breakup of Gondwanaland at ~95-100Ma with the Gondwanan tribes diverging ~85-95Ma (Table 3).The Gondwanan distribution of Dichotomini and Canthonini is improbable with the origin of dung feeding Scarabaeinae at ~76Ma.

Biogeography of Scarabaeinae
The geographic origin of dung beetles remains controversial with many reviews of biogeographic evidence suggesting a Mesozoic origin of scarabaeines from Gondwanan ancestors [80,[82][83].Alternatively, Sole & Scholtz [84] proposed an African origin of dung beetles based on a 3-gene phylogeny of 25 African genera from the tribes Dichotomini and Canthonini.Using molecular clock methods, the age of the origin of the two tribes was estimated at around 56Ma.Sole & Scholtz [84] concluded that given this inferred age and as representatives of these early branching scarabaeines are found in Africa, the fauna therefore originated and radiated on Africa and that dispersal is responsible for their current distribution.A similar study of African Scarabaeinae that sampled 4 genes from 55 of 105 African genera suggested the earliest split in the subfamily in Africa occurred between 42 and 27 Ma based on molecular clock estimates of COI data only [85] and tested two differing substitution rates used in previous analyses of scarabaeine divergence [84,86].Mlambo et al. [85] question the validity of their older estimates due to poor resolution of a number of parameters, and favour the higher arthropod rate of 0.012 mutations per million years to suggest an origin of the subfamily at 33.9Ma (as listed in Fig 6 of [85]).To date the most extensively sampled molecular phylogeny is that of Monaghan et al. [87] which sampled 3 genes from 214 species representing all continents (excluding Antarctica) to explore systematic relationships of the Scarabaeinae, evolution of nesting strategies, and biogeographic scenarios.The Dispersal-Vicariance (DIVA) analysis supported an out of Africa hypothesis, with long distance dispersal to all other continents with little back migration [87], however no timeline was given for the origin or dispersal of scarabaeine lineages.Phillips et al. [81] suggested a Mesozoic/Cretaceous origin of scarabaeine dung beetles with subsequent diversification in the Tertiary, and hypothesized that a combination of vicariance and dispersal (predominately from Africa) accounted for the present distribution of dung beetle tribes.While long distance dispersal is certainly likely for the distribution of some tribes (e.g.Onthophagini, Oniticellini and Sisyphini), long distance dispersal is unlikely for clades containing both Old World and New World taxa due to their limited flight abilities [63].The biogeographic distribution of scarabaeine tribes was discussed in detail by Tarasov & Génier [63], under two evolutionary scenarios, Cenozoic vs. Mesozoic origins, who suggested that the distribution of dung beetles better fit a vicariance hypothesis.Tarasov & Génier [63] hypothesized that if the origin of the Scarabaeinae coincided with the separation of Gondwana from Africa (~160Ma) and clades containing both Old World and New World taxa originated while Africa and South America remained connected (~110-95Ma) then the phylogenetic pattern recovered within their morphological phylogeny would support a vicariant pattern.
Previous studies of African dung beetles [84][85][86] rely on molecular clock methods due to the limited fossil record.The use of molecular clocks remains controversial as 'standard' substitution rates are often calculated from small samples of closely related species and given the influence of calibration choices on calculated rates [88] 'standard' substitution rates may not be realistic due to rate variation rates across time scales [89] or among invertebrate species [90].For example, Wirta et al [86] and Mlambo et al. [85] both use a standard range of rates of 0.0075 and 0.012 substitutions/site/Ma, yet the results from these two studies are incongruent with each other.Wirta et al. [86] hypothesises that Malagasy Oniticellini/'Helictopleurini' diverged from the African Oniticellini either 44(29/64) or 28(18-39) Ma for the rates 0.0075 and 0.012 substitutions/site/Ma, but Mlambo et al. [85] in replicating Wirta et al.'s [86] methods, estimated the earliest split in the Scarabaeinae at 42(32/53) or 27(20/35) Ma and origin of Africa Oniticellini, a derived tribe, at ~13.2Ma.Such inconsistencies highlight the limitations of molecular clock methods.To date, only the divergence dating analyses of Ahrens et al. [18] and those presented here are calibrated using fossil data and samples from beyond the Scarabaeinae and Aphodiinae to allow the incorporation of more distant calibration points.The origin of Scarabaeinae inferred by these later studies [present study, 18], significantly predate the molecular clock based estimates [84][85], which may be attributed to fossil calibration, or the broader taxonomic and geographic sampling.Our various calibration schemes recovered the origin of scarabaeine dung beetles in the Lower Cretaceous (mean age Aptian-Berriasian, minimum to maximum 95% confidence 116.5-152.4Ma),roughly consistent with the split of Africa from Gondwana.The opening of the South Atlantic Ocean and the split from Africa and Gondwana is poorly understood, however the general consensus based evidence from plate tectonics is that the intercontinental rift was initiated ~140Ma and that complete separation of the two continents occurred just prior to the Aptian, ~126Ma [91][92][93].This timing is consistent with hypotheses derived on the similarities in Vertebrate fauna, that South America, Africa, Madagascar, India and Australia remained partially connected until the Aptian [94][95][96][97][98].It has also been proposed that faunal exchange between Africa and South America may have persisted until the Albian-Cenomanian through a connecting land bridge [99], a chain of islands or through Antarctica and Australia [96][97].Given the inferred age of Scarabaeinae in the Lower the major radiation of dung feeders prior to the Cenomanian, and early divergence of both African and Gondwanan lineages, we hypothesis that the faunal exchange between Africa and Gondwanaland occurred during the earliest evolution of the Scarabaeinae.Therefore we propose that both Gondwanan vicariance and dispersal of African lineages is responsible for present day distribution of the Scarabaeinae.
Support for a vicariance pattern can be observed in the Canthonini clade 2 representing Australian, Malagasy and New Caledonian fauna, estimated to have originated ~85-105Ma with the new Caledonian genera diverging ~50-55Ma and splitting from a clade containing the Australian genera Demarziella and Diorygopyx ~76-82Ma (Table 3).These ages are broadly consistent with the break-up of the Gondwana landmass with New Caledonia (+ New Zealand) splitting from Gondwana ~80Ma.Unfortunately no New Zealand dung beetles were included in our analysis to test the subsequent splitting of New Caledonia and New Zealand ~30-40 Ma.Further investigation is needed to confirm our hypothesis that the origin of Scarabaeinae occurred while Africa and South America were still connected but our results indicate that biogeographic distribution and the estimated origin of Canthonini fits a vicariance scenario.
In terms of endemic scarabaeine radiations, the origin of the Malagasy Oniticellini (often classified as Helictopleurini) is interesting to consider in a biogeographic context.The lineage was paraphyletic in our original Bayesian analysis but the 5 Helictopleuris spp.were recovered as monophyletic in the fully-bifurcating tree used in MCMCTree analyses.The mean age of crown group "Helictopleurini" in MCMCTree analyses was estimated at 61.16-66.17Mawhile PL estimated a minimum age of 56.18-60.92Ma in PL for the clade containing 3 of 5 Helictopleuris species sampled.Given that Madagascar split from India ~80Ma [100], long distance dispersal best explains the biogeographic pattern of Malagasy Oniticellini.This dispersal hypothesis is in agreement with the findings of a comprehensive molecular phylogeny of "Helictopleurini" that dated the origin of the lineage ~23-37 Ma using molecular clock methods [86].Wirta et al. [86] explored a scenario that radiation of "Helictopleurini" was triggered by the arrival and radiation of lemurs and other Malagasy mammals, stating that their estimates broadly corresponded to this hypothesis.Interestingly our estimates are remarkably congruent with the arrival of lemurs in Madagascar ~55-65Ma [101][102][103].While our sampling is too limited to examine Malagasy radiations in detail, we can compare mechanisms for successful long distance dispersal to Madagascar from Africa in the Paleocene.The detailed study of spatial and temporal arrival patterns of Madagascar's vertebrate fauna, tectonic history and oceanic currents demonstrated that Early Cenozoic surface currents were periodically conducive to rafting from Africa [104].Successful transoceanic dispersal by rafters was highest in the Paleogene, with success decreasing over time, reaching its lowest levels in the mid-Miocene [104].Although we do not rule out flight as the method of colonization, if Oniticelli arrived in Madagascar in the Paleocene, then transoceanic dispersal via rafting was possible.As such, we propose an alternative method of long distance dispersal to Madagascar given the limited flight ability of scarabaeine dung beetles and rarity of medium distance dispersal events [63], although rafting is less probable with an Oligocene origin as suggested by Wirta et al. [86].

The impact of the K-Pg event and radiation in the Paleogene
The mass-extinction of non-avian dinosaurs in the Cretaceous-Paleogene (K-Pg) event, 65-66Ma would likely have had an impact on dung beetle diversity given the loss of an important resource, dinosaur dung.A rate shift is observed in the Scarabaeinae LTT plots around this period and provides evidence for co-extinction (Fig 4C).In contrast to the LTT analyses, TESS and CoMET analyses could neither identify a mass extinction, nor any significant rates of diversification or extinction both within the Scarabaeinae and phytophagous clades.This failure to identify an increase in diversification within the lineage coinciding or with a sequential lag with the origin of angiosperms is surprising.Unfortunately TESS and CoMET methods are reliant on sampling fraction and with the diversity of these scarab lineages (~5000 Scarabaeinae and ~19800 Pleurosticti), our sampling fractions were always between 2-8%, raising questions as to whether these methods are suitable for such limited taxon sampling.Accordingly, the following discussion will focus on the LTT results as they are more robust to the sampling used in this study.
The composition of extant scarabaeine dung beetle guilds provides insight into their survival through the K-Pg extinction.Most extant dung beetles display clear preferences for dung based on host diet, digestive processes and dropping size [73], although a few extant species are extreme generalists exist which may feed on herbivore, carnivore and omnivore dung, bird and reptile droppings, insect frass, carcasses and dead insects (e.g.[105]).It is unknown if Cretaceous dung beetles were generalists or specialists, however if particle size is a critical determinant to dung feeding, it is unlikely Cretaceous scarabaeines were able to feed on both large droppings of mega-herbivore dinosaurs and small, dry pellets of the primitive insectivorous mammals.However many small (<100kg) herbivorous and omnivorous dinosaurs [74] existed during the Cretaceous and it is possible that dung from these animals could have been utilized by small, generalist dung beetles that could also feed on dung from primitive mammals and so survive the loss of dinosaur hosts.
The shape of our LTT plots (Fig 4C) is consistent with a partial extinction of Scarabaeinae around the K-Pg event [106], likely the loss of those species that were dependent on dung of large dinosaurs such as those recorded in coprolite trace fossils.The extant dung beetle fauna is likely to be descended from species either already adapted to Cretaceous mammal dung or generalists capable of utilizing the dung of both small dinosaurs and early mammals, or alternative resources.For ~20 Ma after the plateau in LTT plots, dung beetle diversity accumulates at its most rapid rate (Fig 4C ).The radiations during the Paleogene could be attributed to the diversification of mammals, or indirectly to the rise of grasses ~70-60Ma [107].The evolution of grasslands in the mid-Eocene and in-turn grazing mammals [79], may have provided a second important ecological opportunity for scarabaeine radiation, but at present, species-level sampling is not sufficient to disentangle the causes of the most recent radiations.However, the role that mammals play on species-level diversification and specialization of feeding in scarabaeines can be explored through richness and abundance of native mammals compared to dung beetle food resources.The islands of Mauritius, New Caledonia and New Zealand lack indigenous mammals (with the exception of native bats) and the species richness of their dung beetle fauna is also low (Mauritius n = 5 spp.; New Caledonia = 13 spp.;New Zealand n = 14 spp.[80]).These faunas are all hypothesized to have east Gondwanan origins, being closely related to Malagasy, or Australian and New Guinean faunas respectively [80], which in comparison exhibit more diverse mammalian and dung beetle faunas (Madagascar: 88 mammals excluding 29 bats [108], ~200 scarabaeines [80]; Australia + New Guinea: ~315 mammals excluding ~120 bats [109], ~430 scarabaeines [80]).A detailed study of food resources of the native dung beetles of New Zealand demonstrated that they have evolved a generalist diet of dung (native reptile, bat, bird and insect, and non-native mammals) and carrion [105].Similarly, derived feeding biologies of Scarabaeinae have been connected to the abundance of large mammals in an ecosystem.Halffter & Matthews [12], noted that copro-necrophagous and necrophagous dung beetles were most common in the Neotropics, proposing a connection between the limited abundance of large mammals compared to open areas with large mammals (e.g. the Afrotropics) where there is almost a complete absence of necrophagous species.Necrophagy in Scarabaeinae is a derived trait, with the Neotropical tribes being either exclusively coprophagous, or predominately coprophagous and occasionally necrophagous [12].Furthermore, almost all genera that have necrophagous are predominately coprophagous or copro-necrophagous with the exception of Deltochilum which is fundamentally necrophagous but occasionally copro-necrophagous, coprophagous [12], and even predatory [110].Unfortunately, species-level sampling is not sufficient to explore evolution of derived feeding ecologies in Scarabaeinae, however it is evident that the local and regional ecological factors potentially driving such specialization warrant further investigation.
Comparison of diversification rates in scarabaeine dung beetle tribes that originate during the mid-Cretaceous versus end-Cretaceous/early-Paleogene provides further evidence for our hypothesis of co-extinction with dinosaurs (Fig 5).The Onthophagini represent ~40% of scarabaeine richness, yet originated approximately 20Ma after the oldest tribes Coprini, Dichotomiini and Canthonini (which represent 8%, 15% and 20% of richness respectively) (Table 7).The high diversification rate of Onthophagini (0.1102 lineages My -1 ) provides strong evidence of a rapid radiation in the Paleogene.Additionally, two of the three tribes that originated in the Paleogene display higher diversification rates than the average rate of the subfamily (Sysiphini = 0.083382, Phanaenini = 0.10314, vs. Scarabaeinae = 0.068635).Extant dung beetles from all tribes are intimately tied to feeding on mammalian dung and there is no reason to suggest that only these younger tribes underwent a rapid radiation with mammals in the Paleogene.Our LTT plots indicate that Scarabaeinae accumulating quickly throughout the Upper Cretaceous, as such, the lower species richness and inferred diversification rates for the early-diverging tribes is best explained by extinction at the K-Pg boundary.
The impact of the K-Pg event on insect diversity is considered to be relatively minor with no family-level extinctions, rather species-level turnover whose intensity varied significantly between different regions [111].Recently, K-Pg mass extinctions have been inferred for butterflies [112][113] and xylocopine bees [114] using similar methods to those used here.In each case the mass extinction of herbivorous insects is attributed to a loss of angiosperm diversity at the K-Pg.However, angiosperms only suffered a modest impact from the K-Pg event: negligible loss of family-level diversity globally and only moderate (18-30%) loss of families/genera at regional scales [50].As such the cause of these mass extinctions of butterflies and bees is harder to directly infer given the minimal loss of angiosperm diversity.Within scarab beetles the significance of the K-Pg event may not be confined only to dung beetles (Scarabaeinae).The accumulation of Melolonthinae appears to halt from 55-65Ma before continued diversification at the same rate to those prior to the interruption (Fig 4A).Interestingly this pattern is not observed in the "CRD clade" of scarabs that are specialised to feed on angiosperm-specific tissues, thus highlighting that perhaps more complex interactions may be the cause of the loss of diversity of herbivores and pollinators at the K-Pg rather than direct loss of angiosperm diversity.

Conclusions
In the absence of a clear and rich fossil record, divergence dating analyses of invertebrates will remain a complex task.We are no closer to reaching a consensus on evolutionary history of scarab beetles in particularly scarabaeine dung beetles.Past studies have estimated the mean age of Scarabaeinae at 33.9Ma [86], 56Ma [85], 86.6-100.2Ma[18] and 118.8-131.6Ma[PL methods presented here].The differences in estimates can be attributed to divergence dating methods (including calibration points, priors and programs), taxon or geographic sampling.We present the first evidence of the mid-Cretaceous origin, and Upper Cretaceous radiation of dung beetles, the timing of which is consistent with major global events including rise of the angiosperms and continental breakup of Gondwana.
The complex associations between plants and animals make it difficult to determine if single or multiple factors drove the evolution of contemporary biotas.This study of phytophagous and saprophagous scarabs revealed that angiosperms played a major role in the diversification of scarabs as both a direct and indirect food source.We hypothesise that a change in host diet, to incorporate more nutritious and less fibrous angiosperm foliage, provided a palatable dung source for dung beetle feeding thus creating a new niche for diversification.This indirect influence of angiosperms on dung beetle diversification opens up our understanding of how nonherbivorous insects (approx.50% of species) evolved following the Cretaceous terrestrial revolution.The confirmation that dung-feeding scarabaeines evolved in the mid-cretaceous and rapidly diversified prior to the availability of a suitable mammalian dung source, provides the first inference of initial origin of any insect group to novel ecological interactions with dinosaurs and radiation due to specialization on a dinosaur resource, i.e. their dung.Our results reveal the complex evolutionary forces shaping these lineages, particularly scarabaeine dung beetles, which have undergone co-radiations with angiosperms and mammals plus potential co-extinction with dinosaurs.The reconstructed history of the Scarabaeidae highlights the complexity of evolutionary diversification and ecological adaptation and provides evidence for an extinction event not captured in the fossil record.As the loss of dung beetle diversity at the K-Pg can be readily associated with the extinction of dinosaurs, this is the first insect mass extinction for which we can confidently infer the cause.within the interval, with dashed line representing significance thresholds of 2lnBF = 2, 6, or 10 (BF = Bayes Factor).Time scale excludes most 45 or 50 Ma depending on analysis that corresponds to species level sampling (see S1 Table for parameters).(PDF) S1 Table .TESS analysis settings.Penalized likelihood (PL) output from r8s, Bayesian output from MCMCTree.Sampling fraction (rho) is calculated from the number tips at the corresponding cut time and compared to the predicted diversity calculated using Method of Moments [48].The most recent 45 or 50 Ma were excluded to minimize effect of limited species level sampling.(DOCX) Figulus sp.
TAS Bold represent newly generated sequences, non-bold represent sequences from GenBank with TaxonomyID numbers.Collection locality is provided for samples collected for this study.doi:10.1371/journal.pone.0153570.t001 Fig 1. Bayesian Phylogeny of scarab beetles (Scarabaeidae).Taxa color-coded by scarab subfamily, with outgroups in grey (superfamily) and black (other beetles).Grey circles indicate polyphagous (P) and saprophagous (S) lifestyles.White circles represent the node priors A-S as per Table 2.* represents nodes for which divergence dates are inferred.See also S1 Fig for posterior probability and terminal names.doi:10.1371/journal.pone.0153570.g001

Fig 2 .
Fig 2. Comparison of accumulation of Scarabaeidae in different divergence dating analyses.(A) Penalized Likelihood (PL) method estimated in r8s using seven different calibration schemes (B) Bayesian methods estimated in MCMCTree using four different calibration schemes (bold) and shadowed by their corresponding PL estimate (pale).doi:10.1371/journal.pone.0153570.g002

Fig 3 .
Fig 3. Lineage through time (LTT) plots of scarab lineages.LTTs are compared to the proportion of extant crown group angiosperm lineages as per Schneider et al. [115].Shaded areas represent congruent estimates from analyses CS2, CS3, CSiii and CSiv estimated by Penalized Likelihood methods with divergence maxima plotted from CS2 and minima from CSiii.doi:10.1371/journal.pone.0153570.g003

Fig 5 .
Fig 5. Scarab beetle diversification rates.Diversification rates calculated using the method of moments [48] and three extinction rates (ε = 0, 0.5, 0.9) of A) dung beetle tribes * The tribes Coprini and Canthonini were not recovered as monophyletic so the node age of Scarabaeinae B was used to estimate diversification rate; B) scarab clades and subfamilies.* The subfamily Rutelinae was not recovered as monophyletic so the node age of Rutelinae+ Dynastinae was used to estimate diversification rate.F = Family clade, P = Phytophagous clade, S = Saprophagous clade.doi:10.1371/journal.pone.0153570.g005

Table 1 .
Taxon and GenBank accession numbers within the phylogenetic analysis.

Table 2 .
Calibration points used for the estimation of the divergence times.
Type of calibration refers to the age estimated from either a molecular estimate (Mol.) or from a fossil (Fos.) and listed as Ma.Analyses using calibration points are highlighted, "minmax" represent clades where both minimum and maximum age constraints are applied, "min" represents only minimum age constraints, empty cells represent a constraint is not applied.Calibration scheme as follows CS1: molecular constraint only; CS2-3 molecular constraint plus selected fossil constraints with Juraclopus rohdendorfi placed at either the crown (CS2) or stem (CS3) of the Pleurosticti clade; CSi-CSiv: fossil selections based of previous study CSi: fossil constraints only, hard bounds applied to all fossils; CSii: fossil constraints only, but limited to only Cretaceous fossils, hard bounds applied; CSiii: fossil constraints only, limited to only Cretaceous fossils, hard bounds only applied to outgroup taxa while ingroup taxa constrained only by minimum age; Civ: replicate of CSiii with molecular age constraint also applied.#Baissa formation is noted as being "probably pre-Barremian" and most paleoentomologists date the Zaza formation as Valanginian-Hauterivian 1

Table 3 .
Predicted ages of scarab lineages using different fossil calibrations and divergence dating programs.

Table 5
summarises the congruence between the estimated ages of clades from our analyses and fossils used as calibration points.Cross-validation analyses highlight that divergence estimates are strongly influenced by how fossil calibration points are treated.Mean age estimates from CSi (maximum bounds applied to 9 ingroup taxa) suggest that the Scarabaeidae diverged

Table 5 .
Cross validation of divergence dates estimated in r8s.Analyses CS1-3 and CSi-iv represent different fossil sets for calibration.Bold text represents nodes constrained with minimum bounds only, bold and underlined represent constrained nodes with minimum and maximum bounds.Oldest known fossils only from the Mesozoic are listed.

Table 6 .
Crown group divergence times.

Table 7 .
[48]e diversification rates estimated using the methods of moments[48].Clade age is the midpoint taken from molecular age estimates when calibrated using the minimum and maximum fossil age of Juraclopus rohdendorfi placed at the crown of the polyphagous scarab clade (E in Fig1).* represents taxon that were not recovered as a monophylum so the age of the node separating the oldest split was used to estimate diversification rates. doi:10.1371/journal.pone.0153570.t007