Winding up the molecular clock in the genus Carabus(Coleoptera: Carabidae): assessment of methodological decisions on rate and node age estimation
© Andújar et al; licensee BioMed Central Ltd. 2012
Received: 1 December 2011
Accepted: 28 March 2012
Published: 28 March 2012
Rates of molecular evolution are known to vary across taxa and among genes, and this requires rate calibration for each specific dataset based on external information. Calibration is sensitive to evolutionary model parameters, partitioning schemes and clock model. However, the way in which these and other analytical aspects affect both the rates and the resulting clade ages from calibrated phylogenies are not yet well understood. To investigate these aspects we have conducted calibration analyses for the genus Carabus (Coleoptera, Carabidae) on five mitochondrial and four nuclear DNA fragments with 7888 nt total length, testing different clock models and partitioning schemes to select the most suitable using Bayes Factors comparisons.
We used these data to investigate the effect of ambiguous character and outgroup inclusion on both the rates of molecular evolution and the TMRCA of Carabus. We found considerable variation in rates of molecular evolution depending on the fragment studied (ranging from 5.02% in cob to 0.26% divergence/My in LSU-A), but also on analytical conditions. Alternative choices of clock model, partitioning scheme, treatment of ambiguous characters, and outgroup inclusion resulted in rate increments ranging from 28% (HUWE1) to 1000% (LSU-B and ITS2) and increments in the TMRCA of Carabus ranging from 8.4% (cox1-A) to 540% (ITS2). Results support an origin of the genus Carabus during the Oligocene in the Eurasian continent followed by a Miocene differentiation that originated all main extant lineages.
The combination of several genes is proposed as the best strategy to minimise both the idiosyncratic behaviors of individual markers and the effect of analytical aspects in rate and age estimations. Our results highlight the importance of estimating rates of molecular evolution for each specific dataset, selecting for optimal clock and partitioning models as well as other methodological issues potentially affecting rate estimation.
KeywordsMolecular clock Rates of molecular evolution Deep node ages Partitioning model Clock model Outgroup selection Gblocks Mitochondrial genes Nuclear genes Coleoptera Carabus
Time calibration of phylogenetic trees is a key factor to reconstruct the evolutionary history of taxa . This is usually accomplished by extrapolating the known age of a node (e.g. based on fossil data) to the remainder of the tree, and assuming a molecular clock. Some animal groups, such as mammals or birds, are especially suited for time estimation based on a rather complete fossil record, but other organisms, typically species-rich groups of invertebrates with poor or inexistent fossil record, may represent more of a challenge for a similar exercise. In the case of insects, for instance, reliable paleontological evidence is frequently lacking, which leads to applying a proposed standard rate of 2.3% divergence/My for the insect mitochondrial genome . However, the development of independent calibration analyses, mainly based on the age of geologic phenomena that underlie the origin of particular cladogenetic events, have found rates either slower e.g., [3–5] or faster e.g., [6–9] than this standard. These studies illustrate how often the routine application of a standard rate may lead to incorrect inference of evolutionary histories.
The discrepancies in the rates of molecular evolution can be attributed to lineage specific effects or to the molecular marker employed [1, 10–12]. However, they could also reflect biases in the calibration procedure . Decisions relative to the analytical procedure with potential effects on estimated rates include the suitability of selected lineage splits and the strategy used to enforce ages to nodes [14–16], methodological aspects such as the method for branch length estimation (i.e., maximum likelihood vs. Bayesian methods; [17, 18]), model of among-branch rate variation (i.e., strict vs. relaxed clock models; [19–22]), selection of evolutionary model e.g., , partitioning of data e.g., [12, 23], taxon sampling e.g.,  or inclusion of ambiguously aligned regions . But the effects of the methodology and specifically the combined effect of choices relative to these methodological aspects have not been fully explored. The increasing number of studies that rely on calibration analyses highlight the importance of investigating how these factors influence evolutionary rate and node age estimation in real datasets.
Within the Coleoptera, the family Carabidae has been the focus of several molecular clock calibration attempts. For instance, Contreras-Diaz et al.  and Ruiz et al.  estimated cox1-cox2 rates of 3.04% and 0.92% divergence/My for Canarian species of Trechus and Sphodrini ground beetles, respectively. Prüser and Mossakowski  used a strict global clock method and the opening of the Gibraltar Strait at the end of the Miocene to calibrate hypothetic vicariance events between Iberian and North African populations of Carabus species. They found rates between 0.39 and 0.98% divergence/My for nd1 data. Su et al.  and Tominaga et al.  investigated nd5 rates for two endemic subgenera of Japanese Carabus using a similar approach and the isolation of Japan from the continent at 15 Mya. These authors found a very low evolutionary rate for this gene, 0.28% divergence/My.
The genus Carabus is a Holarctic taxon that includes about 950 species, currently divided in eight main divisions . It is highly diversified in the Palaearctic, where it is distributed throughout continental Eurasia, but is also present in Japan, Iceland, the Canary Islands, North Africa and several Mediterranean islands. Indeed, the genus Carabus represents a good research subject to conduct comparative calibration analyses, since there are a fair number of available DNA sequences, a relatively solid systematic knowledge of its limits and major splits , and the evolutionary history of these apterous beetles can be linked to geologic events that provide multiple potential calibration hypotheses about the origin of clades.
Here we perform calibration analyses on nine gene fragments that include protein coding and ribosomal genes belonging to both mitochondrial and nuclear genomes. Calibration analyses are based on the ages of the most recent common ancestor (TMRCA) for three cladogenetic events, including the origin of subgenera Mesocarabus, Macrothorax and the sister pair Eurycarabus and Nesaeocarabus, as estimated on an nd5 phylogeny using eight calibration points based on dates for fossil and past geologic events. The procedure takes into account the uncertainty intervals of these calibration points, to avoid a false perception of precision on age estimation in subsequent calibration analyses. Overall, we have conducted a total of 152 independent calibration analyses on individual and concatenated DNA matrices, using different outgroups, clock models, partition schemes and alternative treatments of ambiguous characters. We used these data to address several specific aims: (i) to obtain a reliable time scale for the origin and evolution of the genus Carabus and discuss the obtained rates of molecular evolution with those reported in other studies, and most critically (ii) to evaluate the effect of methodological decisions relative to calibration analyses in the resulting calibrated phylogenies.
Taxon and gene sampling
Species, locality data, voucher reference and accession numbers for each specimen and sequence.
Pirin Mts., Bulgaria
C. (Archicarabus) nemoralis
C. (Platycarabus) irregularis
C. (Rhabdotocarabus) melancholicus
C. (Rhabdotocarabus) melancholicus
C. (Rhabdotocarabus) melancholicus
C. (Carabus) deyrollei
C. (Eurycarabus) famini
C. (Eurycarabus) famini
El Alia, Tunisia
C. (Nesaeocarabus) abbreviatus
C. (Nesaeocarabus) abbreviatus
C. (Morphocarabus) monilis
C. (Mesocarabus) dufourii
C. (Mesocarabus) lusitanicus
Ciudad Real, Spain
C. (Mesocarabus) macrocephalus
La Coruña, Spain
C. (Mesocarabus) riffensis
C. (Chrysocarabus) auronitens
C. (Chrysocarabus) rutilans
C. (Lamprostus) coriaceus
C. (Macrothorax) morbillosus
C. (Macrothorax) morbillosus
C. (Macrothorax) morbillosus
C. (Macrothorax) morbillosus
El Alia, Tunisia
C. (Macrothorax) rugosus
C. (Macrothorax) rugosus
C. (Megodontus) violaceus
Popovi Livadi, Bulgaria
C. (Limnocarabaus) clatratus
C. (Tachypus) cancellatus
C. (Tachypus) cancellatus
Each specimen was characterized for nine DNA fragments corresponding to seven different ribosomal and protein coding genes from mitochondrial (cox1-A, cox1-B, nd5, cytb, rrnL) and nuclear (LSU-A, LSU-B, ITS2, HUWE1) genomes (Additional file 1: Table S1). The sequences of HUWE1 are homologous to the Anonymous gene described by Sota and Vogler  for the genus Carabus, and to the predicted HUWE1 gene as identified by BLAST searches against the Tribolium castaneum genome. PCR reactions were made using PuReTaq Ready-To-Go PCR beads (GE Healthcare, UK) or Qiagen Taq Polymerase with 39 cycles at 50-54°C for primer annealing. Purification of PCR products and sequencing in both directions with the same primers used for PCR was performed by Macrogen Inc. (Seoul, Korea). Sequence accession numbers are given in Table 1.
Mitochondrial protein coding genes were unambiguously aligned and checked for their correct translation to amino acids using Mega 4 . The 5'-end of the HUWE1 fragment was also unambiguously aligned and correctly translated to amino acids, while the 3'-end showed a length-variable intron sequence. All ribosomal markers (rrnL, LSU-A, LSU-B and ITS2) also showed length variation and required objective alignment prior to phylogenetic analysis.
Variable-length DNA fragments for the dataset including outgroups were aligned under each combination of five iterative refinement methods (FFT-NS-i, E-INS-i, G-INS-I, L-INS-I and Q-INS-I) and three scoring matrices (1-PAM, 20-PAM and 200-PAM) in Mafft 6.240 [33, 34]. Each individual alignment was assessed for congruence with respect to a combined matrix including unambiguously aligned regions of every gene. To get this combined matrix, every fragment was independently aligned in MAFFT with FFT-NS-i parameters, and local ambiguities were removed with Gblocks  with the No-gaps option and other default parameters; the resulting fragments were concatenated. Congruence was measured using both the incongruence length difference index (ILD; ) and the rescaled ILD . ILD values were estimated from parsimony-based tree lengths in every case using PAUP* 4.0 . The alignment conditions maximizing character congruence for every length-variable marker, i.e. producing the lowest rescaled ILD value, were objectively selected as those generating the best homology hypothesis for these data and employed in subsequent analyses.
Favored alignments were used to produce three concatenated matrices, including all mitochondrial fragments (MIT), all nuclear fragments (NUC) and both datasets (MIT-NUC). In the case of ribosomal genes, selected alignments were previously processed with Gblocks  using the All-gaps option and default parameters; only selected positions were included in the concatenated matrices. These alignments are available at Treebase.org (submission number 12410: http://purl.org/phylo/treebase/phylows/study/TB2:S12410)
All phylogenetic and calibration analyses described below were run independently for both the complete dataset (outgroup dataset) and a subset of data representing only the 28 ingroup Carabus sequences (ingroup dataset).
Bayesian phylogenetic inference for each individual and concatenated datasets, without specifying partitions and without clock assumptions, were run with MrBayes 3.1 [39, 40] under a GTR + G + I model to assess the reliability of nodes to be used in subsequent calibration tests. The GTR model of nucleotide substitution was identified by jModeltest  as the best fitting for nd5, cox1-B, LSU-A, LSU-B, ITS2 and all concatenated datasets. For the remaining genes (cox1-A, cob, rrnL and HUWE1) a simpler model was favored, but the GTR was second best. The parameter gamma (G) was favored for all data sets but not invariant characters (I) in the case of nuclear ribosomal genes, although this parameter was included in the second best fitting model. To homogenize this part of the methodology, a complex GTR + G + I model was used in every analysis as their parameters are co-estimated with the tree and they can match eventually any simpler model at expense of increasing perhaps the variance of estimates. Analyses enforced two independent runs, each with three hot and one cold chain, for 20,000,000 generations, whereby trees were sampled every 1,000 generations. Convergence of independent runs was checked in Tracer 1.5  and their half compatible consensus tree was calculated excluding 10% of initial trees, after the plateau in tree likelihood values had been reached. Trees were visualized using FigTree 1.1.2  and node posterior probabilities were interpreted as support values.
Calibration hypotheses employed to time-calibrate the initial phylogeny of the genus Carabus and related taxa based on the nd5 gene
Evolutionary event: node
Event age (Ma)
Prior 95% age interval
Split between two Canarian endemic species: Carabus (Nesaeocarabus) coarctatus
and C. (N.) abbreviates: C
of Gran Canaria
(a = 0, b = 14.5)
Carabus (Autocarabus) cancellatus fossil: F
of Cantal (France)
Lognormal (μ = 25,σ = 1.5, offset = 5)
Radiation of Damaster: J1
Final disconnection of Japan from mainland
Truncated Normal (μ = 3.5,σ = 1, a = 0.1, b = 1000)
Radiation of Leptocarabus: J2
Final disconnection of Japan from mainland
Truncated Normal (μ = 3.5,σ = 1, a = 0.1, b = 1000)
Radiation of Ohomopterus: J3
Final disconnection of Japan from mainland
Truncated Normal (μ = 3.5,σ = 1, a = 0.1, b = 1000)
Split between subgenus Isiocarabus and Ohomopterus: J4
Initial disconnection of Japan from mainland
Normal (μ = 15,σ = 1)
Split between Carabus (Eurycarabus) genei from Corsica and North African Eurycarabus: M2
Opening Gibraltar strait
Exponential (μ = 0.5, offset = 5.3)
Split between two Carabus (Rhabdotocarabus) melancholicus subspecies: M3
Opening Gibraltar strait
Exponential (μ = 0.5, offset = 5.3)
Calibration points employed to time-calibrate molecular phylogenies of single and combined datasets in Carabus
AGE and 95% HPD interval (Ma)
Split between Carabus (Macrothorax) rugosus and C. (Macrothorax) morbillosus
Shape: 63.787; Scale: 0.118
Split of Carabus (Mesocarabus) riffensis from European Mesocarabu
Shape: 68.490; Scale: 0.1604
Split between the subgenera Eurycarabus and Nesaeocarabus
Shape: 66.361; Scale: 0.144
Calibration analyses of each marker and their combination were conducted in BEAST 1.5.4 based upon four independent runs of 50 million generations each, sampling every 2,000th generation, and using a Yule tree prior and a GTR + G + I evolutionary model, with ten categories for the gamma distribution. Samples from these independent runs were compared, checked for convergence and combined after conservatively removing 10% of initial trees in Logcombiner 1.5.4 , drawing one sample every 8,000th generation. Mean, standard error, highest posterior density intervals (95%HPD) and effective sample size of likelihood, evolutionary rates and the TMRCA of Carabus were inspected using Tracer 1.5. Consensus trees were obtained in TreeAnnotator 1.5.4  using the mean age option.
Calibration analyses on individual protein coding genes. Protein coding genes, including the 5'-end of the HUWE1 fragment, were analyzed under different clock assumptions including SC and ULN clocks, as well as different codon partition schemes (1P, 2P and 3P), for both the outgroup and ingroup datasets. Additionally, the complete HUWE1 gene fragment (i.e. including the non-coding 3'-end) was analyzed without codon partitioning. Results were analyzed to simultaneously select the best clock and partition model using BF as above.
Calibration analyses on individual ribosomal genes. Calibration analyses for trees based on ribosomal genes were conducted considering (i) all positions (complete dataset), (ii) only positions selected by Gblocks with the No-gaps option and other default parameters (nogaps dataset), and (iii) only positions selected with the All-gaps option (allgaps dataset). Each individual analysis was done under different clock assumptions (SC and ULN) and for the outgroup and ingroup matrices. The best clock model was assessed in every case using BF comparisons as above.
Calibration analyses on concatenated datasets. Outgroup and ingroup datasets for the MIT, NUC and MIT-NUC concatenated matrices were analyzed in BEAST under different clock assumptions (SC and ULN) following different partition schemes, which included no data partitioning (NP), partitioning by gene but not by codon in the case of protein coding genes (G-1P), partitioning by gene, and each protein coding gene with two codon partitions where first and second positions are considered together (G-2P), and partitioning by gene and by codon position (G-3P). Among-branch rate variation was always treated as linked among partitions, and BF comparisons were used again to select for the best clock and partitioning models.
The effect of different treatments of data was explored using Wilcoxon signed rank test on the values of two relevant parameters, including the estimated rate of molecular evolution for the marker or markers investigated and the estimated age of the ingroup node (TMRCA of Carabus).
Phylogenetic framework for calibration tests
The alignment of choice based on congruence optimization with unambiguously aligned data was that obtained implementing the Q-INS-i algorithm for all nuclear genes, except for the LSU-A fragment in the case of the ingroup dataset. The latter, together with the mitochondrial ribosomal gene rrnl data, were optimal using the E-INS-i method (Additional file 1: Table S3).
Calibration analyses on the expanded nd5 gene dataset used a strict clock and 2P codon partitioning as selected by BF, resulting in a rate of molecular evolution of 0.0154 (95% HPD 0.0112-0.0198) substitutions per site per million years per lineage (subs/s/Ma/l). The time calibrated phylogeny obtained is shown in Figure 1. Nodes A, B and C were recovered with high support (posterior probability of 1.0) with median ages ranging from 7.5 Ma (node A) to 11.9 Ma (node B; Table 3).
Selection of optimal analytical conditions
For each individual and combined dataset and alternative partition and clock model schemes, the four independent BEAST runs resulted in similar global rates, TMRCA of Carabus and likelihood values, reaching stationary equilibrium and ESS values always higher than 500 for the likelihood parameter and higher than 200 for all the other parameters, with very few exceptions. The results from these independent runs and for each dataset were thus pooled together to generate samples of 180 million generations, with ESS values always higher than 200.
BF comparisons resulted in the selection of the 2P codon partition strategy for mitochondrial coding genes, the combination of all mtDNA markers, and the combination of all markers, while no data partitioning was selected for the nuclear HUWE1 fragment alone and for all nuclear markers combined (Additional file 1: Table S4). The strict clock was favored for all protein coding genes (including the entire HUWE1 gene fragment), the combined mtDNA including outgroup taxa, and also for some ribosomal genes using the ingroup dataset and after gap exclusion. Otherwise, the relaxed (ULN) clock was preferred for all other individual and combined markers (Additional file 1: Table S4). Relative to the effect of ambiguously aligned characters or outgroup inclusion/exclusion, no direct BF comparisons were possible, and a pragmatic decision was taken in every case. Nogaps datasets were discarded due to their major effect on rates and age estimations (see below). Allgaps and complete matrices showed similar rates and ages, and the first option was selected. Finally, we discarded outgroups for estimation of molecular rates and TMRCA of Carabus since this genus was not found monophyletic in most individual marker analyses, producing a net overestimation of the age for this node, taken as the one including all ingroup taxa but not only. When monophyly was recovered, such as in concatenated datasets, the differences in both rates and ages between outgroup and ingroup datasets was low.
Evolutionary rates and ingroup ages
Rate of molecular evolution and TMRCA of Carabus for each individual fragment and combined datasets under optimal analytical conditions
The mean estimated age of the ingroup oscillated between 13.4 (8.35-24.85) Ma as estimated for LSU-A data, and 31.2 (16.8-52.47) Ma, in the case of ITS2, whereas for most genes this value was between 20 and 30 Ma with widely overlapping 95% HPD intervals (Figure 2, Table 4). The combined analyses of all genes resulted in a TMRCA of Carabus of 25.2 (18.41-33.04) Ma. The dispersion of values around the mean was higher for genes analyzed under a relaxed clock. The split between Carabus and Calosoma (node K) was estimated to have occurred at around 30.6 Ma (95% HPD 24.27-37.65) with combined genes and the dataset including outgroup taxa (Figure 3).
Effect of partitioning scheme
Effect of ambiguously aligned characters
The simultaneous analysis of non-coding sequence information with exon information for the nuclear HUWE1 gene had little effect on the estimation of evolutionary rates, although the estimation of the ingroup age decreased or increased when assuming strict or relaxed clocks, respectively (Figure 6).
Effect of clock model
The choice of strict versus relaxed clock had effects on the estimation of parameters of interest, generally associated with the actual clock model best fitting the data. Individual and combined genes in which the strict clock was preferred showed null to low effect of clock model on the estimation of rates and ingroup ages, with both parameters showing at most a trend to slightly higher values when a relaxed clock was enforced (Wilcoxon signed rank test on mean rate, P < 0.01; Wilcoxon signed rank test on median TMRCA of Carabus, P < 0.01). In all other cases, except for the total evidence dataset, the use of the suboptimal strict clock resulted in lower rate estimates and higher ingroup ages (Figures 4, 6, 5).
Effect of outgroup
The dataset including outgroups generally resulted in higher rates of molecular evolution compared with the analyses using ingroup data only (Wilcoxon signed rank test, P < 0.01) (Figures 4, 6, 5). Exceptions to this pattern affected the fast evolving ribosomal markers LSU-B and ITS2, as well as the combination of nuclear genes and the total evidence dataset when investigated under a strict clock model. In turn, for the TMRCA of Carabus no general trend could be identified (Wilcoxon signed rank test, P = 0.1143), despite it was generally retrieved as older for most treatments when including outgroups, for instance with individual mitochondrial genes (Wilcoxon signed rank test, P < 0.01), a fact associated to most individual gene fragments failing to recover the monophyly of Carabus. Exceptionally, this age was younger for all nuclear markers independently (except LSU-B) and for their combination (NUC dataset) when the unfavored strict clock was enforced.
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.
Our study centred on the genus Carabus illustrates some of the common problems but also alternative solutions and suitable strategies to molecular clock calibration analyses, exploring the possibility of synergic effects of wrong methodological decision. Mitochondrial genes generally fit the strict clock model of evolution, providing an adequate and simple frame to conduct calibration analyses. However, the combination of several genes including nuclear and mitochondrial fragments results in the best strategy to minimise the effect of both the idiosyncratic behavior of individual markers and analytical aspects on the estimation of evolutionary rates and node ages. But even for mitochondrial genes or combined datasets, rate differences between markers, together with the effect of specific analytical decisions on their estimation, advise against extrapolating rates between different studies. As a summary of the findings reported here, it should be always desirable to check for the appropriate data partition scheme, and for the underlying evolutionary model for each partition; it is also necessary to investigate whether the data really fit a molecular clock previous to any calibration study, applying suitable corrections if needed. The analysis of time on trees should avoid as much as possible regions with dubious homology due to alignment ambiguity, or keep them to a minimum, since gap exclusion has a dramatic effect on branch length estimation and concomitant rates and node ages. Finally, the analyses benefit from enquiry circumscribed to ingroup taxa (and using relatively recent calibration nodes), minimizing the biases introduced by highly divergent lineages and saturation.
Natural logarithm of the Bayes factors
Highest posterior density
Time to the most recent common ancestor
Two codon partitions: with 1st and 2nd position together
Three codon partitions
Partitioning by gene with no codon partition
Partitioning by gene: with two codon partition of coding genes with 1st and 2nd position together
Partitioning by gene: with three codon partition of coding genes
Uncorrelated log normal clock
Dataset excluding ambiguous character with the Nogaps option and default parameters in Gblocks
Idem: applying Allgaps option
Concatenation of mitochondrial genes
Concatenation of nuclear genes
Concatenation of all genes
This research was funded by projects CGL2006/06706, CGL2009-10906 (JS) and CGL2008-00007 (JGZ) of the Spanish Ministry of Science and Innovation (the latter with support of the European Regional Development Fund, "Una manera de hacer Europa"), as well as project 08724PI08 of the Fundación Séneca (Murcia) (JS). CA received the support of FPU predoctoral studentship of the Spanish Ministry of Education. Thanks are due to Obdulia Sánchez, Ana Asensio (University of Murcia) and Gwenaelle Genson (CBGP Montpellier) for their technical assistance. Carlos Ruiz (University of Murcia), Ignacio Ribera (IBE, Barcelona) and especially Paula Arribas (University of Murcia) helped with their comments and advise. Jean-Yves Rasplus (CBGP Montpellier) facilitated C. (Tachypus) cancellatus and Achille Casale (University of Sassari) helped with the identification of some taxa. MrBayes was run using the resources of the Computational Biology Service Unit from Cornell University, which is partially funded by Microsoft Corporation; most final computing was carried out using the facilities of the Ben Arabi supercomputer of the Fundación Parque Científico de Murcia.
- Bromham L, Penny D: The modern molecular clock. Nat Rev Genet. 2003, 4: 216-224. 10.1038/nrg1020.PubMedView ArticleGoogle Scholar
- Brower AVZ: Rapid morphological radiation and convergence among races of the butterfly Heliconius erat inferred from patterns of mitochondrial-DNA evolution. Proc Nat Acad Sci USA. 1994, 91: 6491-6495. 10.1073/pnas.91.14.6491.PubMedPubMed CentralView ArticleGoogle Scholar
- Pruser F, Mossakowski D: Low substitution rates in mitochondrial DNA in Mediterranean carabid beetles. Insect Mol Biol. 1998, 7: 121-128. 10.1046/j.1365-2583.1998.72056.x.PubMedView ArticleGoogle Scholar
- Gómez-Zurita J, Juan C, Petitpierre E: The evolutionary history of the genus Timarch (Coleoptera, Chrysomelidae) inferred from mitochondrial COII Gene and partial 16S rDNA sequences. Mol Phylogenet Evol. 2000, 14: 304-317. 10.1006/mpev.1999.0712.PubMedView ArticleGoogle Scholar
- Ruiz C, Jordal B, Serrano J: Molecular phylogeny of the tribe Sphodrini (Coleoptera: Carabidae) based on mitochondrial and nuclear markers. Mol Phylogenet Evol. 2009, 50: 44-58. 10.1016/j.ympev.2008.09.023.PubMedView ArticleGoogle Scholar
- Luchetti A, Marini M, Mantovani B: Mitochondrial evolutionary rate and speciation in termites: data on European Reticulitermes taxa (Isoptera, Rhinotermitidae). Insectes Soc. 2005, 52: 218-221. 10.1007/s00040-005-0806-5.View ArticleGoogle Scholar
- Contreras-Diaz HG, Moya O, Oromi P, Juan C: Evolution and diversification of the forest and hypogean ground-beetle genus Trechus in the Canary Islands. Mol Phylogenet Evol. 2007, 42: 687-699. 10.1016/j.ympev.2006.10.007.PubMedView ArticleGoogle Scholar
- Ribera I, Fresneda J, Bucur R, Izquierdo A, Vogler A, Salgado J, Cieslak A: Ancient origin of a Western Mediterranean radiation of subterranean beetles. BMC Evol Biol. 2010, 10: 29-10.1186/1471-2148-10-29.PubMedPubMed CentralView ArticleGoogle Scholar
- Papadopoulou A, Anastasiou I, Vogler A: Revisiting the insect mitochondrial molecular clock: the mid-Aegean trench calibration. Mol Biol Evol. 2010, 27: 1659-1672. 10.1093/molbev/msq051.PubMedView ArticleGoogle Scholar
- Kumar S: Molecular clocks: four decades of evolution. Nat Rev Genet. 2005, 6: 654-662. 10.1038/nrg1659.PubMedView ArticleGoogle Scholar
- Mueller RL: Evolutionary rates, divergence dates, and the performance of mitochondrial genes in bayesian phylogenetic analysis. Syst Biol. 2006, 55: 289-300. 10.1080/10635150500541672.PubMedView ArticleGoogle Scholar
- Pons J, Ribera I, Bertranpetit J, Balke M: Nucleotide substitution rates for the full set of mitochondrial protein-coding genes in Coleoptera. Mol Phylogenet Evol. 2010, 56: 796-807. 10.1016/j.ympev.2010.02.007.PubMedView ArticleGoogle Scholar
- Schwartz RS, Mueller RL: Variation in DNA substitution rates among lineages erroneously inferred from simulated clock-like data. PLoS One. 2010, 5: e9649-10.1371/journal.pone.0009649.PubMedPubMed CentralView ArticleGoogle Scholar
- Graur D, Martin W: Reading the entrails of chickens: molecular timescales of evolution and the illusion of precision. Trends Genet. 2004, 20: 80-86. 10.1016/j.tig.2003.12.003.PubMedView ArticleGoogle Scholar
- Heads M: Dating nodes on molecular phylogenies: a critique of molecular biogeography. Cladistics. 2005, 21: 62-78. 10.1111/j.1096-0031.2005.00052.x.View ArticleGoogle Scholar
- Inoue J, Donoghue PCJ, Yang Z: The impact of the representation of fossil calibrations on Bayesian estimation of species divergence times. Syst Biol. 2010, 59: 74-89. 10.1093/sysbio/syp078.PubMedView ArticleGoogle Scholar
- Brown JM, Hedtke SM, Lemmon AR, Lemmon EM: When trees grow too long: investigating the causes of highly inaccurate Bayesian branch-length estimates. Syst Biol. 2010, 59: 145-161. 10.1093/sysbio/syp081.PubMedView ArticleGoogle Scholar
- Schwartz RS, Mueller RL: Branch length estimation and divergence dating: estimates of error in Bayesian and maximum likelihood frameworks. BMC Evol Biol. 2010, 10:Google Scholar
- Aris-Brosou S, Yang Z: Effects of models of rate evolution on estimation of divergence dates with special reference to the metazoan 18S ribosomal RNA phylogeny. Syst Biol. 2002, 51: 703-714. 10.1080/10635150290102375.PubMedView ArticleGoogle Scholar
- Drummond AJ, Ho SYW, Phillips MJ, Rambaut A: Relaxed phylogenetics and dating with confidence. Plos Biol. 2006, 4: 699-710.View ArticleGoogle Scholar
- Battistuzzi FU, Filipski A, Hedges SB, Kumar S: Performance of relaxed-clock methods in estimating evolutionary divergence times and their credibility intervals. Mol Biol Evol. 2010, 27: 1289-1300. 10.1093/molbev/msq014.PubMedPubMed CentralView ArticleGoogle Scholar
- Battistuzzi FU, Billing-Ross P, Paliwal A, Kumar S: Fast and slow implementations of relaxed clock methods show similar patterns of accuracy in estimating divergence times. Mol Biol Evol. 2011, 28: 2439-2442. 10.1093/molbev/msr100.PubMedPubMed CentralView ArticleGoogle Scholar
- Brandley MC, Wang Y, Guo X, Montes de Oca AN, Fería Ortiz M, Hikida T, Ota H: Accommodating heterogenous rates of evolution in molecular divergence dating methods: an example using intercontinental dispersal of Plestiodon (Eumeces) Lizards. Syst Biol. 2011, 6: 3-15.View ArticleGoogle Scholar
- Linder HP, Hardy CR, Rutschmann F: Taxon sampling effects in molecular clock dating: An example from the African Restionaceae. Mol Phylogenet Evol. 2005, 35: 569-582. 10.1016/j.ympev.2004.12.006.PubMedView ArticleGoogle Scholar
- Lemmon AR, Brown JM, Stanger-Hall K, Lemmon EM: The effect of ambiguous data on phylogenetic estimates obtained by maximum likelihood and Bayesian inference. Syst Biol. 2009, 58: 130-145. 10.1093/sysbio/syp017.PubMedView ArticleGoogle Scholar
- Su ZH, Tominaga O, Okamoto M, Osawa S: Origin and diversification of hindwingless Damaster ground beetles within the Japanese islands as deduced from mitochondrial ND5 gene sequences (Coleoptera, Carabidae). Mol Biol Evol. 1998, 15: 1026-1039.PubMedView ArticleGoogle Scholar
- Tominaga O, Su ZH, Kim CG, Okamoto M, Imura Y, Osawa S: Formation of the Japanese carabina fauna inferred from a phylogenetic tree of mitochondrial ND5 gene sequences (Coleoptera, Carabidae). J Mol Evol. 2000, 50: 541-549.PubMedGoogle Scholar
- Deuve T: Illustrated Catalogue of the Genus Carabus of the World (Coleoptera: Carabidae. 2004, Sofia-Moscow: PENSOFTGoogle Scholar
- Sota T, Ishikawa R: Phylogeny and life-history evolution in Carabus (subtribe Carabina: Coleoptera, Carabidae) based on sequences of two nuclear genes. Biol J Linn Soc. 2004, 81: 135-149. 10.1111/j.1095-8312.2004.00277.x.View ArticleGoogle Scholar
- Arndt E, Brücker M, Marciniak K, Mossakowski D, Prüser F: Phylogeny. The genus Carabus L. in Europe, a synthesis. Edited by: Turin H, Penev L. 2003, Sofía-Moscow-Leiden: PENSOFT & European Invertebrate Survey, 307-325.Google Scholar
- Sota T, Vogler AP: Incongruence of mitochondrial and nuclear gene trees in the Carabid beetles Ohomopterus. Syst Biol. 2001, 50: 39-59.PubMedView ArticleGoogle Scholar
- Tamura K, Dudley J, Nei M, Kumar S: MEGA4: Molecular evolutionary genetics analysis (MEGA). Mol Biol Evol. 2007, 24: 1596-1599. 10.1093/molbev/msm092.PubMedView ArticleGoogle Scholar
- Katoh K, Misawa K, Kuma K, Miyata T: MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002, 30: 3059-3066. 10.1093/nar/gkf436.PubMedPubMed CentralView ArticleGoogle Scholar
- Katoh K, Asimenos G, Toh H: Multiple alignment of DNA sequences with MAFFT. Bioinf DNA Seq Anal. 2009, 537: 39-64. 10.1007/978-1-59745-251-9_3.View ArticleGoogle Scholar
- Castresana J: Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 2000, 17: 540-552.PubMedView ArticleGoogle Scholar
- Farris JS, Kallersjo M, Kluge AG, Bult C: Testing significance of incongruence. Cladistics. 1994, 10: 315-319. 10.1111/j.1096-0031.1994.tb00181.x.View ArticleGoogle Scholar
- Wheeler WC, Hayashi CY: The phylogeny of the extant chelicerate orders. Cladistics. 1998, 14: 173-192. 10.1111/j.1096-0031.1998.tb00331.x.View ArticleGoogle Scholar
- Swofford DL: PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). 2000, Sunderland, Massachusetts: Sinauer AssociatesGoogle Scholar
- Huelsenbeck JP, Ronquist F: MrBayes: Bayesian inference of phylogenetic trees. Bioinformatics. 2001, 17: 754-755. 10.1093/bioinformatics/17.8.754.PubMedView ArticleGoogle Scholar
- Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.PubMedView ArticleGoogle Scholar
- Posada D: jModelTest: Phylogenetic model averaging. Mol Biol Evol. 2008, 25: 1253-1256. 10.1093/molbev/msn083.PubMedView ArticleGoogle Scholar
- Rambaut A, Drummond AJ: Tracer v1.5 2007. [http://beast.bio.ed.ac.uk/Tracer]
- Rambaut A: FigTree v.1.1.2 2008. [http://tree.bio.ed.ac.uk/software/figtree/]
- Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7:Google Scholar
- Pagel M, Meade A: A phylogenetic mixture model for detecting pattern-heterogeneity in gene sequence or character-state data. Syst Biol. 2004, 53: 571-581. 10.1080/10635150490468675.PubMedView ArticleGoogle Scholar
- Miller KB, Bergsten J, Whiting MF: Phylogeny and classification of the tribe Hydaticini (Coleoptera: Dytiscidae): partition choice for Bayesian analysis with multiple nuclear and mitochondrial protein-coding genes. Zool Scr. 2009, 38: 591-615. 10.1111/j.1463-6409.2009.00393.x.View ArticleGoogle Scholar
- Ho SYW, Lanfear R, Bromham L, Phillips MJ, Soubrier J, Rodrigo AG, Cooper A: Time-dependent rates of molecular evolution. Mol Ecol. 2011, 20: 3087-3101. 10.1111/j.1365-294X.2011.05178.x.PubMedView ArticleGoogle Scholar
- Rambaut A, Drummond AJ: TreeStat v1.6.1: tree statistic calculation tool. 2010Google Scholar
- Ziegler PA: Evolution of the Arctic-North Atlantic and the Western Tethy. 1988, Tulsa, OK: American Association of Petroleum Geologists MemoirGoogle Scholar
- Turin H, Penev L, Casale A: The Genus Carabus in Europe. A synthesis. 2003, Sofía-Moscow-Leiden: PENSOFT & European Invertebrate SurveyGoogle Scholar
- Su ZH, Imura Y, Okamoto M, Kim CG, Zhou HZ, Paik JC, Osawa S: Phylogeny and evolution of Digitulati ground beetles (Coleoptera, Carabidae) inferred from mitochondrial ND5 gene sequences. Mol Phylogenet Evol. 2004, 30: 152-166. 10.1016/S1055-7903(03)00163-5.PubMedView ArticleGoogle Scholar
- Su ZH, Imura Y, Zhou HZ, Okamoto M, Osawa S: Mode of morphological differentiation in the Latitarsi-ground beetles (Coleoptera, Carabidae) of the world inferred from a phylogenetic tree of mitochondrial ND5 gene sequences. Genes Genet Syst. 2003, 78: 53-70. 10.1266/ggs.78.53.PubMedView ArticleGoogle Scholar
- Ohta T: Slightly deleterious mutant substitutions in evolution. Nature. 1973, 246: 96-98. 10.1038/246096a0.PubMedView ArticleGoogle Scholar
- Ribera I, Hernando C, Aguilera P: Agabus alexandrae sp n. from Morocco, with a molecular phylogeny of the Western Mediterranean species of the A. guttatus group (Coleoptera: Dytiscidae). Insect Syst Evol. 2001, 32: 253-262. 10.1163/187631201X00191.View ArticleGoogle Scholar
- Barraclough TG, Vogler AP: Recent diversification rates in North American tiger beetles estimated from a dated mtDNA phylogenetic tree. Mol Biol Evol. 2002, 19: 1706-1716. 10.1093/oxfordjournals.molbev.a003993.PubMedView ArticleGoogle Scholar
- Ribera I, Vogler AP: Speciation of Iberian diving beetles in Pleistocene refugia (Coleoptera, Dytiscidae). Mol Ecol. 2004, 13: 179-193. 10.1046/j.1365-294X.2003.02035.x.PubMedView ArticleGoogle Scholar
- Gómez-Zurita J, Garneria I, Petitpierre E: Molecular phylogenetics and evolutionary analysis of body shape in the genus Cyrtonus (Coleoptera, Chrysomelidae). J Zool Syst Evol Res. 2007, 45: 317-328. 10.1111/j.1439-0469.2006.00393.x.View ArticleGoogle Scholar
- Balke M, Gómez-Zurita J, Ribera I, Viloria A, Zillikens A, Steiner J, Garcia M, Hendrich L, Vogler AP: Ancient associations of aquatic beetles and tank bromeliads in the Neotropical forest canopy. Proc Nat Acad Sci USA. 2008, 105: 6356-6361. 10.1073/pnas.0710368105.PubMedPubMed CentralView ArticleGoogle Scholar
- Gómez-Zurita J, Funk DJ, Vogler AP: The evolution of unisexuality in Calligrapha leaf beetles: Molecular and ecological insights on multiple origins via interspecific hybridization. Evolution. 2006, 60: 328-347.PubMedView ArticleGoogle Scholar
- Hidalgo-Galiana A, Ribera I: Late Miocene diversification of the genus Hydrochus (Coleoptera, Hydrochidae) in the west Mediterranean area. Mol Phylogenet Evol. 2011, 59: 377-385. 10.1016/j.ympev.2011.01.018.PubMedView ArticleGoogle Scholar
- Faille A, Casale A, Ribera I: Phylogenetic relationships of Western Mediterranean subterranean Trechini groundbeetles (Coleoptera: Carabidae). Zool Scr. 2011, 40: 282-295. 10.1111/j.1463-6409.2010.00467.x.View ArticleGoogle Scholar
- Brown JM, Lemmon AR: The importance of data partitioning and the utility of bayes factors in Bayesian phylogenetics. Syst Biol. 2007, 56: 643-655. 10.1080/10635150701546249.PubMedView ArticleGoogle Scholar
- Zheng Y, Peng R, Kuro-o M, Zeng X: Exploring patterns and extent of bias in estimating divergence time from mitochondrial DNA sequence data in a particular lineage: A case study of salamanders (Order Caudata). Mol Biol Evol. 2011, 28: 2521-2535. 10.1093/molbev/msr072.PubMedView ArticleGoogle Scholar
- Marshall DC: Cryptic failure of partitioned Bayesian phylogenetic analyses: Lost in the land of long trees. Syst Biol. 2010, 59: 108-117. 10.1093/sysbio/syp080.PubMedView ArticleGoogle Scholar
- Golubchik T, Wise MJ, Easteal S, Jermiin LS: Mind the gaps: Evidence of bias in estimates of multiple sequence alignments. Mol Biol Evol. 2007, 24: 2433-2442. 10.1093/molbev/msm176.PubMedView ArticleGoogle Scholar