A time scale for the origin and evolution of Carabus
The analyses of MIT, NUC and MIT-NUC combined datasets produce highly congruent topologies and high support for most nodes, including nodes T and K, representing the monophyly of Carabus and its sister relationship with Calosoma, respectively (Figure 3). In addition, all combined datasets show nearly complete overlap of the 95% HPD intervals on the estimated ages for these nodes; only the NUC dataset departs slightly from these values producing an older mean age estimate (Figure 2). The time scale obtained when all nuclear and mitochondrial genes are analyzed together (MIT-NUC dataset) situates the initial split between Carabus and Calosoma during the Oligocene, some 34 and 23 Ma, after the opening of the Atlantic Ocean and the split of the Nearctic and Palearctic regions . This timing is congruent with the observation that Carabus, essentially a flightless genus, is more diverse in the Palaearctic (more than 900 species) than in the Nearctic region (12 species), whereas Calosoma is slightly more diverse in the Nearctic (ca. 90 species) than in the Palearctic region (76 species). The evolutionary events that originated the main extant lineages of Carabus took place according to our data during the early Miocene, between 23 and 16 Ma (Figure 3).
These findings disagree with the previous hypothesis assigning an older Eocene origin to Carabus . The occurrence of North American endemic subgenera Tanaocarabus and Lichnocarabus was interpreted as evidence suggesting that the origin of the genus predated the opening of the Atlantic Ocean, considered the event responsible for the isolation of Nearctic from Western Palearctic lineages of Carabus. However, molecular data indicate that these Nearctic subgenera are instead related to Eastern Palearctic species [51, 52], which suggests an origin for the genus Carabus within the Paleartic region and a more recent dispersal event across Bering Strait land bridges, in agreement with the dates here provided.
Rates of molecular evolution in Carabus
The evolutionary rates for the genus Carabus as estimated here are, in general terms, higher than those reported in previous studies. The discrepancies are thought to be mostly related to the use of inappropriate calibration points in these studies, combined with simplistic corrections of genetic distances among species. Prüser and Mossakowski  in a study based on the nd1 gene in western Mediterranean species of the genus Carabus, calibrated the separation of six pairs of taxa with the opening of the Gibraltar Strait at the end of the Messinian (5.3 Ma). Five of these splits represented the separation of North African and European subspecies in C. (Macrothorax) morbillosus and C. (Macrothorax) rugosus, and the resulting rates ranged between 0.0020 to 0.0033 subs/s/Ma/l. However, these splits seemingly occurred well after the opening of the Gibraltar Strait (Andújar et al., unpublished data). The sixth node represented the split of subspecies of Carabus (Rhabdotocarabus) melancholicus from the Iberian Peninsula and North Africa, the only one with compelling evidence to represent a vicariant event resulting from the opening of Gibraltar Strait (Andújar et al., unpublished data). The rate estimated by Prüser and Mossakowski  was 0.0049 subs/s/Ma/l, much lower than our estimated rate for the slowest protein-coding gene--0.0113 (0.0081-0.0147) subs/s/Ma/l. Indeed, the divergence estimated by Prüser and Mossakowski  for these taxa, based on an uncorrected p-distance, was approximately half as low as distances corrected with appropriate evolutionary models (e.g., GTR + G model of evolution; ).
The degree of relaxation of evolutionary constraints acting on different genes or their portions ultimately reflects in differences on their inferred rate of molecular evolution [1, 53]. Our data showed a 16-fold difference between the slowest (rrnL) and the fastest (cob) evolving mitochondrial markers, in agreement with rate variation estimates between particular mitochondrial genes found by Pons et al.  for Coleoptera. This important disparity cautions against the extrapolation of rates of molecular evolution among different gene fragments or their combinations. On the other hand, the extrapolation of evolutionary rates for the same marker across taxa, at least when they are relatively closely related, appears as a safer assumption, although it is possible to find differences. Several studies using a similar approach as that described here but targeting different families (and suborders) of Coleoptera produced slightly different rates of evolution for the same gene fragments. Thus, for instance, Papadopoulou et al.  and Ribera et al.  estimated rates of 3.54% and 4.08% divergence/Ma in Tenebrionidae and Leiodidae, respectively, for a fragment homologous to the cox1-B fragment investigated here, which in our case yielded a slower rate of 2.90% divergence/Ma (0.0145 subs/s/Ma/l; 95% HPD:0.01-0.0198) in Carabus. Evolutionary rates obtained by Pons et al.  for the nd5 (0.0168; 95% HPD 0.0086-0.0279) and cob (0.0172; 95% HPD 0.0071-0.0311) genes are very similar to those we obtained for Carabus (Table 4), while they found an unexpected higher rate for the cox1 gene (0.0861), although their estimation appeared associated to a surprisingly wide HPD interval (95% HPD 0.0251-0.1760).
Heterogeneity in rates of molecular evolution
Invertebrate mitochondrial genomes have long been assumed to evolve at a standard molecular clock rate of 2.3% divergence/Ma . This rate was deduced from heterogeneous mitochondrial data, including restriction fragment length polymorphism, DNA-DNA hybridization and sequence data for several genes. Its utilization in phylogenetic studies has been frequent for datasets composed by individual and concatenated mitochondrial DNA gene fragments (e.g., only in the case of beetles and for concatenated data: [7, 54–61]). It is important to notice that in these studies similar rates for different taxa and different regions of the mitochondrial genome were assumed. This assumption has been lately refuted e.g., [9, 11, 12]. Moreover, in former studies where the standard rate was routinely applied the effects of methodological decisions on the calibration procedure had not been considered, a question that has been shown to be important [18, 23], as we also stress here.
There is fair variation in the rates of individual gene fragments in Carabus which are affected by methodological aspects of the calibration procedure. Major differences on both the rates and the TMRCA of Carabus are obtained in the case of the nuclear ribosomal genes, which reach a fourfold variation between analyses involving different clock models, as well as inclusion/exclusion of ambiguous characters and/or outgroups. Mitochondrial genes show comparatively less variation due to methodological decisions, and only the nd5 fragment showed a twofold difference depending on treatment (but see below). The most remarkable observation is that a calibration based on combined datasets results in the lowest effect of prior analytical decisions on both evolutionary rate and node age estimation, together with a reduction of 95% HPD intervals associated to these estimates. This effect hints at the importance of using multiple gene data on calibration exercises, both as a way to average across genes, diluting idiosyncratic behaviors of stand-alone markers, but also to help the analyses to converge into more precise estimates.
Clock calibration with protein coding genes
The response to changes in analytical conditions differs between protein coding mitochondrial markers and the HUWE1 gene fragment. Thus, for the latter, the selection of a suboptimal clock model produces marked differences in the estimation of both the rates of molecular evolution and the inferred TMRCA of Carabus, while this choice reveals no remarkable effects in the case of the mitochondrial genes. This behavior can be related to the underlying rate heterogeneity for each marker class: clock-like in the case of mtDNA data, but moderate in that of HUWE1. The selection of an unsuitable strict clock for the nuclear gene distorts branch length estimation and all parameters derived from these. In turn, the protein coding nuclear marker shows no apparent changes in the two parameters of interest through different partition schemes, while the mitochondrial genes prove more sensitive to complex evolutionary models, including the effect of adding or excluding an outgroup to the estimations. As with clock model selection, the partitioning scheme employed affects branch length estimation, and consequently the resulting inferred rate, a behavior already reported in other studies [23, 62]. Mitochondrial protein coding genes show a trend to slightly higher inferred evolutionary rates and TMRCA of Carabus with increasing partitioning; and the same is true for the MIT and MIT-NUC datasets. The influence of codon site-specific models is similar to that found by Papadopoulou et al. , in their study of darkling beetles from the Aegean islands based on two mitochondrial and two nuclear fragments. These authors found up to 11% discrepancy in inferred rates between NP and 3P partitioning, higher for the latter, both under ULN clock assumption and using relatively recent calibration nodes (9-12 Ma), as used in our study. Conversely, other studies found the opposite trend, with younger estimated ages when using complex codon partitioning models for mitochondrial data (e.g., Nearctic and eastern Palaearctic skinks ). This trend can be caused by the same underlying problem related to the specific way in which saturation effects are corrected, but combined in this case with the use of deep calibration points (96-148 Ma) to infer younger node ages.
There seems to be an effect of the relative position of calibrating nodes on the ages extrapolated along the tree. Saturation and incorrect branch length estimation in deeper parts of the tree despite complex model-based corrections could explain these differences. The use of relatively recent calibration nodes should produce reliable age estimates for other nodes in areas of the tree not affected by saturation, with error accumulating progressively for deeper nodes, possibly with underestimated branch lengths, and in this case providing with minimum age estimates. Depth constraints restrain branch stretching in most of the tree, and when deep nodes are affected by incorrect length estimation due to saturation problems, errors in resulting node ages will be extrapolated to more recent parts of the tree, resulting in older than true time estimates for such recent nodes .
Bayesian phylogenetic inference is known to produce incorrect long branch estimates at least in the case mitochondrial codon-partitioned datasets [17, 64]. These flawed estimations have been found to be stable across independent runs with fixed parameters, but also when Bayesian prior parameters are modified, an additional difficulty to detect them. It appears to be the case of the observed high rate increments when partitioning nd5 data by codon position (2P and 3P strategies) for both the ingroup and outgroup datasets. This type of misbehavior seems to be dependent on the characteristics of particular datasets, so that slight changes in taxon sampling can provide correct estimates . Taking advantage from the availability of nd5 sequences of Carabus in GenBank, we have conducted NP, 2P and 3P partitioning calibration analyses for the expanded nd5 dataset (58 sequences; Additional file 1: Table S2). These analyses produced moderate increments in the estimated evolutionary rate for nd5 as more partitions were considered, in agreement with the results from other mtDNA genes, supporting the reliability of these estimates.
Clock calibration with non-coding genes
The effect of methodological choices is markedly higher for non-coding gene fragments than it is for protein coding genes, on both mean rate and node age estimations, but also on their corresponding 95% HPD intervals. Indeed, calibration exercises based on non-coding genes frequently face alignment ambiguity and departures from the molecular clock, which can both jeopardize reliable branch length estimation. These problems increase with evolutionary distance between taxa, thus raising concerns about deep node calibration using slowly evolving ribosomal genes. Several alignment strategies have been proposed to deal with length variation in homologous DNA sequences (see an evaluation of methods in ), and discarding ambiguous DNA positions has been also proposed as a complementary method . We examined the effects of the latter procedure on global rate and node age estimation. Expectedly, character removal, especially when applying restrictive culling options (nogaps in Gblocks), has a marked effect which is negatively correlated with marker conservation. Highly variable gene fragments, with gappy alignments, can lose a significant amount of phylogenetic information when filtered out of ambiguous alignment positions, and consequently contribute with unreliable calibrations.
As with protein coding genes, there is a strong effect of clock model selection for nuclear ribosomal gene fragments as well. For these markers, SC was unfavored in most cases and its implementation resulted in lower mean values for inferred evolutionary rates and higher ages (up to fourfold) on the estimated TMRCA of Carabus, compared to the use of a favored ULN clock. While this behavior is not exclusive of non-coding genes, its prevalence for non-coding data cautions about their use for node age estimation without appropriate accounting for rate heterogeneity in calibration analyses.