- Research article
- Open Access
Inferring speciation modes in a clade of Iberian chafers from rates of morphological evolution in different character systems
BMC Evolutionary Biology volume 9, Article number: 234 (2009)
Studies of speciation mode based on phylogenies usually test the predicted effect on diversification patterns or on geographical distribution of closely related species. Here we outline an approach to infer the prevalent speciation mode in Iberian Hymenoplia chafers through the comparison of the evolutionary rates of morphological character systems likely to be related to sexual or ecological selection. Assuming that mitochondrial evolution is neutral and not related to measured phenotypic differences among the species, we contrast hypothetic outcomes of three speciation modes: 1) geographic isolation with subsequent random morphological divergence, resulting in overall change proportional to the mtDNA rate; 2) sexual selection on size and shape of the male intromittent organs, resulting in an evolutionary rate decoupled to that of the mtDNA; and 3) ecological segregation, reflected in character systems presumably related to ecological or biological adaptations, with rates decoupled from that of the mtDNA.
The evolutionary rate of qualitative external body characters was significantly correlated to that of the mtDNA both for the overall root-to-tip patristic distances and the individual inter-node branches, as measured with standard statistics and the randomization of a global comparison metric (the z-score). The rate of the body morphospace was significantly correlated to that of the mtDNA only for the individual branches, but not for the patristic distances, while that of the paramere outline was significantly correlated with mtDNA rates only for the patristic distances but not for the individual branches.
Structural morphological characters, often used for species recognition, have evolved at a rate proportional to that of the mtDNA, with no evidence of directional or stabilising selection according to our measures. The change in body morphospace seems to have evolved randomly at short term, but the overall change is different from that expected under a pure random drift or randomly fluctuating selection, reflecting either directional or stabilising selection or developmental constraints. Short term changes in paramere shape possibly reflect sexual selection, but their overall amount of change was unconstrained, possibly reflecting their lack of functionality. Our approach may be useful to provide indirect insights into the prevalence of different speciation modes in entire lineages when direct evidence is lacking.
The study of speciation patterns has been traditionally dominated by the role of geography [1, 2], most likely due to the ready availability of species-level phylogenies and data on geographical ranges of species [3, 4]. Other speciation modes, such as sexual selection or ecological segregation, have usually been studied with phylogenies through the influence on species numbers or diversification rates [5–7], the evolution of ecomorphological characters in relation to habitat type and co-occurrence , the evolution of sexual characters [9–11], or the reconstruction of the ancestral niche . However, there are only few cases in which the prevalence of these three main speciation modes (geographical isolation, sexual selection, ecological segregation) has been assessed in a comparative framework, rather than testing their individual effect on a target group (but see e.g., [13, 14]).
In this work we use different morphological traits of a group of chafers of the genus Hymenoplia Eschscholtz to elucidate the prevalent speciation mode. The genus Hymenoplia is part of the tribe Sericini, a group of chafer beetles (Scarabaeidae), which are particularly diverse among the phytophagous scarabs (Pleurosticti) . Species of Hymenoplia are ecologically and morphologically very similar (see Methods), and their ranges have a large degree of overlap at larger geographical scales, although they generally do not co-occur at the local level. This pattern suggests either that there could be factors promoting speciation other than geographical isolation, or that the geographical signature of speciation has been lost by post-speciation range movements . Through the use of species-level phylogenies with multiple samples per species and different sets of morphological traits, we aim to detect the evidence for the influence of selection during trait evolution in order to discriminate between the dominant role of three general modes of speciation: geographical, ecological and sexual, taking into account infraspecific variation. Our null assumption is that if the rate of evolution of a character system across the phylogeny is proportional to time it could be considered "neutral", in the sense that it is not subjected to any strong directional or stabilizing selection (, see other examples in e.g.  for genital shape;  for ecological differences;  for general morphological traits).
To have a baseline for a neutral rate we use that of the mitochondrial DNA (mtDNA) [17, 20, 21]), and consider it proportional to time (despite possible exceptions, see Discussion), and not related to the phenotypic expression of any of the measured morphological characters. We compare this reference neutral evolutionary rate with that of three different morphological character systems: 1) shape of male genitalia as indicative of sexual selection; 2) general body size and shape and 3) morphological structural characters (including those used for species recognition), the last two indicative of a possible ecological partitioning. The character systems that show an evolutionary rate significantly correlated to that of the mtDNA could be said to be neutral and proportional to time, i.e. evolving under a Brownian motion model  or subjected to random fluctuating selection . Those differing significantly could be said to have evolutionary rates not proportional to time, likely to be due to stabilizing or directional selection . Specifically, we examine three hypotheses:
1) Geographical speciation: if the dominant mode of speciation in the group has been geographical isolation with subsequent drift for both molecular and morphological characters (e.g., ), all character sets could be expected to evolve in a predominantly neutral mode, and in consequence the rates of change be correlated among them as a result of their own correlation with time or to general pleiotropic effects.
2) Speciation driven by sexual selection (through divergence in shape and size of the genital structures): this should be reflected in the rates of evolution of the morphology of the intromittent portions of copulatory organs [9, 24]. If their shape and size are subjected to selection, their rates of evolution will not be proportional to time or its surrogate, the rates of mtDNA change. Given the very weak sexual dimorphism in species of Hymenoplia, only apparent by enlarged and lobiform anterior inner protarsal claws in the male, it may be expected that sexual selection did not have an impact on structural morphology or body shape.
3) Ecological segregation: although there is no information on the detailed ecology and biology of most species of Hymenoplia, it seems reasonable to assume that these potential differences could be reflected either in body shape and size or in structural external characters. Again, if these are subjected to selection, it is expected that their evolutionary rates will be de-coupled to that of the mitochondrial genome.
The decoupling of the molecular and morphological rates of evolution could be interpreted as indicative of the presence of directional or stabilizing selection. However, different types of selection are expected to result in contrasting patterns of morphological variation . Diversifying or directional selection should result in increasing divergence between species over evolutionary time, with a low ratio of intra- to interspecific variation . On the contrary, stabilizing selection may result in the random drift of the morphological change within some fixed boundaries [26–28]. In the later, rates of morphological evolution should be neutral and proportional to time within the allowed boundaries, but over longer evolutionary periods this relationship is not maintained, and rates should be decoupled from time (i.e. mtDNA rates) . Any sign of selection on the shape of the male copulatory organs is likely to be related to sexual selection and speciation, but the presence of selection (i.e. decoupling from neutral rates) in body morphospace or other morphological characters could be related to the same speciation process or to subsequet anagenetic evolution.
The comparisons of rates of change are performed at two levels [29, 30]: 1) Individual branches, i.e. internodal distances. We compare the change of each character system over all corresponding individual branches of the estimated topology. This gives an overall relationship of relative rates of change for multiple individual time periods. 2) Patristic distances between terminals. We computed a patristic distance matrix between all terminals for each character system, and then compared their association with the mtDNA distances through multiple Mantel tests. We also use pairwise plots of patristic distances to test for general differences between infraspecific and interspecific variation. Although our methods for hypothesis testing do not require the a-priori delimitation of species, which in groups with a complex taxonomy (like the genus Hymenoplia) is not a trivial issue, evolutionary processes may be expected to be different at infra and interspecific level under certain types of selection. The type of molecular data used (mitochondrial DNA) and the difficulties in the taxonomy of the group (see below) prevents the use of methods requiring the precise delimitation of species (e.g. Fontaneto et al. 2007), but the comparisons of the slopes of the regressions of intra- and among-species variation could provide insights into how variation is partitioned.
There is an extensive literature on the relationship between morphological and molecular evolutionary rates (e.g., [17, 19, 23, 25, 30–33]), although in most cases the rates were compared across different phylogenies, not between different morphological character systems in the same phylogenetic tree. By the use of parallel comparisons of morphological characters measured in the same specimens across the same phylogeny we avoid many of the problems associated with previous approaches (such as e.g. differences in taxon sampling, or the evolutionary background or the biology of the species; see  for a recent review). The existence of clear differences in the evolutionary rates among the studied character systems will provide insightful suggestions as to which has been the dominant speciation mode in this group of polyphagous scarab chafers, and thus contribute to understand which could have been the key factors for their diversification.
Species of Hymenoplia feed preferably on leaves and inflorescences of grasses as adults [34, 35] and presumably on humus and roots in the larval stages, like other members of the Pleurosticts . A recent phylogenetic analysis  showed that Hymenoplia, with ca. 46 taxa, is part of the least diverse (ca. 150 taxa) of the two sericine sister lineages of the Old World (with a combined diversity of ca. 3,200 species). Hymenoplia species occur exclusively in the Western-Mediterranean region, and a third of them in the Iberian peninsula . The traditional taxonomy of the genus Hymenoplia has been confusing, due to the lack of clear diagnostic characters, the variability among and within populations, and their largely overlapping geographical ranges [39, 40]. They do have similar ecologies, but, intriguingly, despite the general sympatry of their ranges (Figure 1), they very rarely co-occur syntopically in the same locality [35, 39]
Taxon sampling, DNA extraction and identification
Hymenoplia specimens were collected from grasses at 16 sites in the southern Iberian peninsula (Figure 1) in 2006. Among the hundreds of Hymenoplia specimens sampled by the first author from a vast range of study sites throughout all southern Spain (Figure 1), in only one site was more than one species found (D. Ahrens, unpublished data). Genomic DNA was extracted from thoracic muscle tissue with Charge Switch gDNA micro tissue kit (Invitrogen, Paisley, UK). Extracted specimens were dry mounted with minimal damage, allowing further morphological investigation. Vouchers are deposited in the D. Ahrens collection (NHM). Diagnostic characters to distinguish adult morphospecies were those traditionally used in taxonomic studies of the group, including body size and shape, coloration, surface sculpture and pilosity as well as male genital morphology [39–41]. All identifications were validated by comparison with the type specimens preserved in the collection of the MNCN. A representative subset of 38 male individuals of the eight named species for which we could obtain DNA samples (half of known Iberian Hymenoplia species ) plus one outgroup species (Additional file 1) were used for DNA sequencing and phylogenetic reconstruction.
Morphological characters and morphometric measurements
We studied three different morphological character systems:
1) General structural characters, including mainly features of cuticular integument such as body surface texture, pilosity and color, but also fine structures of legs and tarsomeres (Additional file 2), which are among those used for the species recognition. Nineteen discrete characters were scored for each adult male individual (Additional files 3, 4). The variation of structural morphological characters was visually checked based on a multiple scaling analysis performed in XLSTAT 2007 (Addinsoft, Paris) using the absolute morphological distances calculated from the original data in PAUP*4.0b10  (Additional file 5).
2) Morphometric characterization of the body shape and size. To characterize the morphospace defined by the studied species we measured total length, width of pronotum, length of elytra and width of elytra (Additional file 6). Due to the limited size range of the studied samples it was not necessary to transform the data, and we used the raw measurements for the analyses to obtain the matrix of Euclidean distances among specimens.
3) Sexual characters. We characterized the outline of the male left paramere (part of the intromittent genital organs) in lateral view. The shape analysis was performed with the software package SHAPEv.1.3  after the images of the lateral view of the paramere were converted into separate black-and-white bitmaps of the studied structures (paramere, phallobasis). SHAPE uses elliptic Fourier descriptors (EFDs) to analyze shape variation of two-dimensional outline data . Four coefficients and 20 harmonics were then extracted from shape outlines and treated as shape variables (Additional file 7). Chain coding, rotation and computation of harmonics were also carried out in SHAPE. Because the male genitalia of Sericini chafers have a complex three-dimensional structure (see e.g., , the shape analysis of a two-dimensional projection is likely to provide a more conservative, albeit more inaccurate, estimate of the total variation. The Fourier descriptors of the specimens were analyzed with Principal Components Analyses (PCA) in XLSTAT 2007 (Addinsoft, Paris), and converted into a distance matrix to be used to compute branch lengths in PAUP*4.0b10  (see below).
To determine the potential error associated with the use of the digital images for the morphometric study, we did five replicates of each image for each specimen; the images were taken by the same person, with a period of at least 2 to 24 h between replicates. The proportion of the total variation due to methodological error was quantified by dividing the trace of the pooled within-specimen covariance matrix by the trace of the total covariance matrix . In addition, we did a MANOVA test to assess whether interspecific variation was significantly higher than the measurement error. MANOVA tests demonstrated that the error was significantly smaller than interspecific differences (Wilk's Lambda P ≈ 0). The covariance structure of measurement error was not significantly correlated with the covariance structure of interspecific differences. All tests of correlated evolution are based on the mean shape values from the five replicates per structure.
The morphometric signal in body and paramere shape, i.e. the resulting principal component axes (Additional file 5) explaining most of the variation in the respective traits (see Results), was investigated independently for potential overlap and significance between infraspecific and interspecific variation using MANOVA and Canonical Variate Analysis in PAST .
Two mitochondrial gene regions were amplified and sequenced for the analyses: 1) cytochrome oxidase subunit 1 (cox1) and 2) 16S ribosomal RNA (rrnL), with the adjacent regions tRNA leucine (trRNA-Leu) and NAD dehydrogenase subunit 1 (nad1; mtDNA nomenclature follows ). The latter three gene fragments are hereafter referred as "rrnL-nad1". PCR and sequencing was performed using primers Pat (5'tccaatgcactaatctgccatatta) and Jerry (5'caacatttattttgattttttgg) for cox1 and 16SaR (5'cgcctgtttaacaaaaacat) and ND1A (5'ggtcccttacgaatttgaatatatcct) for rrnL-nad1 . Sequencing was performed on both strands using BigDye v.2.1 (Applied Biosystems, Carslbad, US) and an ABI3700 automated sequencer in the facilities of the CIB (CSIC, Madrid). Sequences were assembled and edited using Sequencher v. 4.5 (Genecodes Corp., Ann Arbor, USA).
There was no length variation in the protein-coding genes, and variation in the ribosomal genes was minimal (see Results), so progressive alignment procedures were conducted in ClustalX 1.83  under default gap opening penalty of 15 and extension penalty of 6.66. The sequences reported in this paper have been deposited in GenBank with accession numbers FJ847234-68/FJ956708-41 (Additional file 1).
To estimate the phylogenetic relationships among species of Hymenoplia we implemented Maximum Likelihood searches of the combined mitochondrial sequence in PhyMLv2.4.4  using a GTR+I+Γ model (as selected in Modeltest 3.06, ), with all parameters estimated from the data. To check the stability of the results we performed additional parsimony and Bayesian analyses. Parsimony tree searches were conducted using TNT 1.0 , with 10 ratchet iterations, 10 cycles of tree drifting and three rounds of tree fusing for each of 200 random addition sequences, coding gaps as 5th character and using five random addition sequences that included ten ratchet iterations, and three rounds of tree fusing using default settings. Bayesian analyses were conducted in MrBayes 3.12  using a GTR+I+Γ model, with four partitions (each codon for cox1 plus rrnL-nad1) and estimating all parameters independently in each partition. Partition homogeneity was tested in PAUP with 100 replicates. All trees were rooted with one species of the Mediterranean genus Paratriodonta Baraud, found to be closely related to (but clearly outside) Hymenoplia .
Node support was assessed by the node posterior probabilities in MrBayes, and by searching 100 pseudo-replicated data sets of nonparametric bootstrapping  both in PhyML and TNT (generated in the latter using the parsimony ratchet).
Comparison of the rates of evolution of the different character sets
To compare the rates of evolution of the different character systems we used as a reference tree the topology obtained with the ML search on the combined mtDNA dataset. We obtained the branch lengths of the four data sets to be compared (mitochondrial, structural characters, morphometry of the body shape, paramere outline) on this constrained topology using PAUP as follows (see e.g.,  for a similar approach):
1) The mtDNA branch lengths, used as a reference for "neutral" change, were obtained using a GTR+I+G model with parameters estimated from the data.
2) The length of the branches of the morphological structural character matrix (Additional file 4) were determined under two different approaches: a) with parsimony [Morphology (pars)] using accelerated transformation (ACCTRAN) character-state optimization, as the number of character state changes along the branch; and b) using a mean character difference distance matrix [Morphology (dist)] derived from the morphological structural character matrix (negative branches were set to be zero). We used a distance approach in addition to the parsimony branch lengths to make the three morphological data sets directly comparable (see below).
3) The length of the branches of the body morphospace and the outline of the male parameres were obtained using a Euclidean distance matrix, obtained from the raw measurements for the body morphospace and from the normalized EFDs of the shape analysis of the parameres (Additional files 6, 7). Negative branches were set to be zero.
We compared the morphological vs. the mtDNA branch lengths at two levels (see Introduction), individual branches and patristic distances between terminals.
1) Individual branches. We did multiple correlation tests (morphology, paramere shape, body shape vs. mtDNA branches) [29, 30] or ANOVA, with the amount of mtDNA change as explanatory variable, comparing it to different sets of branches (body shape, paramere shape, and morphology) grouped according to their amount of character state changes. We divided the mt branches in four (0, 0-0.001, 0.001-0.01, >0.01) or seven logarithmic categories (0, 0-0.0005, 0.0005-0.001, 0.001-0.005, 0.005-0.01, 0.01-0.05, >0.05), each having almost similar amounts of branch numbers. Since the distribution of the length of the branches of the morphological structural change was highly biased towards low values (Additional file 8) we also did the reverse ANOVA, i.e. with the length of the morphological branches as explanatory variable for the comparison with the DNA change. We divided the branches between those with no change (length = 0) and those with change (length > 0), and compared the length of the mitochondrial branches of the two groups with an ANOVA analysis. To have an overall comparison of the differences of the individual branches over the tree, we also used the K-score as measured by Ktreedist 1.0 . The software provides a single metric (the K-score) for the comparison between the branch lengths and the topological differences of the two trees, with no associated significance. In our cases, the topologies of all trees were identical, so K-scores were a measure of branch length differences only. Lower values of the K-score imply less differences from the mtDNA tree, with a higher degree of similarity between the branch lengths and, consequently, a predominantly neutral mode of trait evolution. The program first calculates the scale factor that minimizes differences between trees, and then computes the branch length distance (BLD, ) between the scaled comparison tree and the reference tree . The K-score is thus the minimum branch length distance between two trees, once one of them has been scaled. To estimate the probability that the observed K-scores were due to random processes, for each data set we obtained a null distribution of the morphological change in the same tree topology by randomizing the names of the terminals, obtaining the corresponding branch lengths in PAUP (which will now be random, with the constraint of the total amount of variation in the tree) and measuring the resulting K-score. This was computed using a Perl script (available on request), with 10,000 random replicas. To assess its significance, the observed K-score was compared to the null distribution of random K-scores .
2) Patristic distances. We computed the patristic distances between all terminals of the tree for the four character systems (mtDNA and three morphological data sets), and tested their association through multiple Mantel tests . Significance was assessed with 10,000 permutations of the observed matrix in the program zt v1.1 . Since diversifying and stabilizing selection can cause differences in the infraspecific and interspecific rates of change [e.g. ], we used pairwise plots of the patristic distances and least square regressions to test for possible differences in the slope of the relationship between infraspecific or interspecific morphological and mitochondrial distances (we did not use simple correlations with these data due to the non-independence of the individual measures, [17, 32].
The node density effect [62–64] could affect trees with either missing data or a punctuated form of evolution, introducing artifacts through the underestimation of the measure of root to tip distances in lineages with a lower number of intermediate nodes. We tested the presence of node density effect in the trees with branch lengths estimated for each character system (mitochondrial and morphological) using the "delta" test , implemented online in http://www.evolution.reading.ac.uk/pe/index.html. The strength of the effect is measured with the curvilinear relationship between the number of nodes and the root to tip distance, and it is considered to be significant if the strength of the relationship is significantly greater than 0, and the degree of curvature significantly greater than unity.
We obtained 826 bp of the 3' end of cox1 (without length variation) and a fragment of rrnL-tRNAleu-nad1 varying in length between 812-816 bp (ingroup). The combined molecular data matrix included 1,651 aligned positions, with 292 parsimony informative characters. The mean uncorrected p-distance between any two sequences was 0.104% and 0.048% for the cox1 and rrnL-nad1 partitions, respectively.
In all analyses Hymenoplia was split into three well supported clades (A, B and C) (Figure 2; see also Additional file 9). The sister group relationship between clade A and B had only low support. Within the clades A-C, all morphologically defined species were monophyletic and well supported, with the only exception of H. lineolata (Figure 2), which had low support in the ML analysis and was paraphyletic in the MrBayes and parsimony analyses due to the inclusion of H. fulvipennis (Additional file 9).
Except for structural morphology, the variation of body and paramere shape had no clear clustering according to the traditionally recognized species, with considerable overlap between them (Additional file 5). A box test (asymptotic approximation to χ2 and Fischer's F) with the measurement of the total body length was not significant (p = 0.3), revealing equal intra-class variance, i.e. the amount of variation within each species is similar. On the contrary, the Wilks' lambda test for the body length was significant (p: < 0.0001), detecting significant differences in the mean body size for at least one species. The Fisher distances for the species classes were significant for about half of the inter-species comparisons, and highly negatively correlated with their significance, as revealed by a mantel test; in other words, large differences between species had the tendency to be also significant. Most of the variation of body morphospace was represented by PC axis 1 (92.7%), with all measurements being strongly correlated with it (r>0.9).
For paramere shape, 95% of the variation was represented by PC axis 1-22. However, and despite the apparent overlap (Additional file 5), differences between the scores of the species in the PC axes corresponding to 95% of total trait variation were highly significant, as measured with MANOVA (paramere shape, Wilks' lambda95%: 4.37E-09; p: 2.028E-20; body shape, Wilks' lambda95%: 0.1994; p: 1.135E-05). The same was valid for the axes representing only 75% (i.e. PC axis 1-10) of shape variation in parameres (Wilks' lambda75%:0.0005526; p: 4.15E-16). The Canonical Variate Analysis applied to these PC axes shows a better species discrimination for paramere shape than for body shape (Additional file 10).
Comparison of molecular and morphological rates of evolution
For the comparison of the rates of evolution of the different character systems we used the topology of the tree obtained from ML search with the combined mtDNA (Figure 2). Branch lengths for the morphological rates of evolution were obtained in PAUP using this topology as a constraint (see Methods). None of the trees showed signs of the presence of node density effects, as measured with the "delta" test .
1) Comparison of individual branches
The length of the branches reflecting mtDNA change and both measures of morphospace (male parameres and general body shape; Figure 3) had approximate normal distributions (Additional file 8), and we used correlations (Pearson correlation) to estimate their relationship. All three correlations were very weak (r<0.1), and not significant (p > 0.05) (Table 1).
The branches with less mtDNA change were also significantly shorter for the structural morphological characters, as measured with both ANOVAs (using four or seven categories of mtDNA change, p4Cat = 0.015 and p7Cat = 0.0003 respectively). This was also the case for the reverse comparison, i.e. the branches with no morphological change were also significantly shorter for the mtDNA change (p < 0.0005, Table 2). On the contrary, there were no significant differences in the length of the morphospace or paramere outline branches between the four and seven categories of mtDNA change (Table 2). When the branches with morphological change were subdivided in three categories (length = 0, length = 1, length > 1) the results were similar (Table 2), with significantly less mtDNA change in the branches with less morphological character state changes, but no relationship with the other two character systems. The number of branches with more than two morphological changes was too low to allow for further subdivisions of the data (Table 2; Additional file 8).
We computed the null distribution of the values of the K-score between the tree with mitochondrial branch lengths and the trees with branch lengths reflecting the change in the three morphological character sets with 10,000 randomized replicas (Figure 4). Only in the comparison between the branches of the mtDNA and the paramere outline the observed value was within 95% of the null random distribution of values (Table 3; Figure 4), providing no evidence of a significant relationship between the evolutionary rates of the two character systems based on the K-score. For the body morphospace and structural morphology character systems (with both parsimony and distance branch lengths) the observed difference with the mtDNA branches was significantly lower than that expected at random, indicating a significant association between the evolutionary rates of these character systems (Table 3).
2) Comparison of patristic distances
The Mantel tests measuring the correlation between the distance matrices of the mtDNA and structural morphology (with both parsimony and distance branch lengths) and between mtDNA and paramere outline were highly significant. On the contrary, the correlation between the mtDNA and the body morphospace matrices was not significant (Table 1).
The pairwise plots of intraspecific patristic distances of the body morphospace and the paramere shape vs. mitochondrial distances show no significant relationship, with slopes of the regression of shape distance vs. sequence divergent not different from zero . The slope of the plot of the morphological vs. mtDNA distances was significantly negative, although with a very low correlation (R2 = 0.04; p < 0.05; slope -7.0, 95% confidence interval [-13.8, -0.1]; Figure 5). For interspecific distances, the slope of the plot of the paramere outline vs. mtDNA was not significantly different from zero, but the body morphospace had a negative, and the morphology a positive slope (body morphospace, R2 = 0.03; p < 0,0001; slope 0.8, 95% confidence interval [-1.2, -0.4]; morphology, R2 = 0.1; p < 0.0001; slope 11.0. 95% confidence interval [8.4, 13.51]; Figure 5).
In this work we exemplified the use of a novel approach to compare rates of evolution of different character systems with the study of speciation modes in a clade of Iberian Sericini chafers of the genus Hymenoplia. There are many factors known to influence character change and speciation, but the usual approach is to study their effect in isolation to increase the explanatory power of the statistical tests, avoiding being swamped by uncontrollable variation [10, 66–68]. The simultaneous consideration of alternative mechanisms is still uncommon [14, 23, 69], and there is still a lack of standard approaches to compare different sources of information in phylogenetic studies.
We obtained three sets of measures of different character systems potentially subjected to different selection pressures, and tested the basic assumption that evolutionary rates should be proportional to time (i.e. "neutral") when there is no such selective pressure, but decoupled from this neutral rate when there is directional or stabilising selection acting on them [1, 17, 18, 23]. The potential relationship between morphological and molecular evolutionary rates has been contentious, with successive re-examinations of the question in the light of new methodological approaches [30, 32, 33]. There are two main differences between ours and previous studies: 1) we assume that any possible correlation between the mtDNA and morphological rates is due to their independent correlation to time [17, 31]; and 2) we do not aim to establish a general correlation between "molecular" and "morphological" rates, but try to determine whether different morphological character systems evolve at a neutral rate proportional to time or, on the contrary, are decoupled from the rates of mtDNA change. Thus, other than indirect effects related to population size or general pleiotropic factors [32, 33], there should be no functional relationship between our measure of molecular evolution, based on mitochondrial genes, and any of the phenotypic traits considered. As we do not compare evolutionary rates among different phylogenies, we only need to establish a relative measure of the rates of different morphological character systems in the same phylogeny, minimizing problems associated with high type I errors or lack of statistical power, which are likely to be major drawbacks . The fact that for almost all type of tests we find both significant and not significant results (Table 4) suggests that our methods are powerful enough to detect a relationship, but with type I errors low enough to avoid generalized false positives.
Of the three tested character systems, the evolutionary rate of "structural morphology" was in all analyses proportional to that of the mtDNA (irrespective of the method used to reconstruct the branch lengths), and thus it could be assumed that structural morphological characters evolved at rates proportional to time. Structural characters are used regularly to separate morphospecies for classification purposes, and also for phylogenetic reconstruction. For most of the characters used here (Additional file 3) there is no evidence of their biological function (e.g. small differences in the density of the pubescence or the punctuation, surface integument, color), and thus how they could drive speciation through selection. The only putative "key innovation", the tarsal claws bearing a ventral membranous fringe (Additional file 3, 4), could be highly adaptive to feeding on grasses, but is not relevant here being a synapomorphy of Hymenoplia + Hymenochelus .
On the contrary, for both body shape and size and paramere outline some of the analyses showed a decoupling of the evolutionary rates with respect to that of the mtDNA. In the case of the body shape and size, the amount of change in the individual branches of the tree was more similar to that of the mtDNA than expected by chance, as shown with the null random distribution of the K-scores. However, when the accumulation of change over the tree was compared through the patristic distances, there was no correlation between the amount of change of the mtDNA and that of the body morphospace. This suggests that when measuring individual branches the deviation from the mtDNA rate is small enough not to be significant with the tests we used (ANOVA comparisons and the null distribution of the K-scores), but when measured over the whole tree with patristic distances the accumulation of changes in the body morphospace is decoupled from the amount of mtDNA change. The lack of correlation between body shape and size and mtDNA may be interpreted as the result of homoplasy, producing "saturation" in the signal due to the independent development of similar shapes and sizes in different lineages of the tree [17, 23]. This could be due to stabilizing selection favoring certain shapes and sizes (although weak enough not to be detectable when comparing individual branches, see above), or just to the existence of a limited range of possible body sizes and shapes for the whole lineage . In other words, there could be a random drift or fluctuating random selection at small temporal scales, detectable with the comparison of individual branches, but over longer periods measured with the patristic distances the limited morphospace (due to stabilizing selection, functional constraints or "saturation" of the morphospace, see ) would cause the individual trajectories to converge and lose their proportionality to mtDNA change . This is quite consistent with the dominance of the PC axis 1 for the morphospace (explaining 92% of variation) and that the discriminant analysis of the total body length revealed that general body size has some significant but not the only contribution to the variation of morphospace data.
The evolution of the shape of the paramere seems to be complementary to that of the body morphospace: when comparing individual branches the rates are not proportional to time, as would be expected under a neutral scenario [1, 70], suggesting strong changes driven by factors other than pure random drift. This directional change, together with high infraspecific variation (see below), are some of the marks of sexual selection . However, the patristic distances do seem to be to some extent proportional to genetic distance: even if at shorter term changes may be directional (and possibly driven by sexual selection), the accumulated change through several cladogenetic events has no directionality, becoming more like an unconstrained random drift and thus with a general correlation with genetic distance. The potential morphospace for the genital organs should be much larger than that of the body shape and size, given the wide range of different genital shapes known in Sericini in contrast with the limited variation in body shape and size (see e.g., ). The overall correlation of the patristic distances with the mtDNA rates suggests an unconstrained amount of change over longer periods of time, reflecting a lack of a functional role of the shape of the paramere. The unconstrained change of structural morphological characters and (to a lesser extend) paramere outline can be directly appreciated in the plots of the interspecific variation (Additional files 9, 10). On the contrary, body morphospace show a high degree of overlap, and poor discrimination among species.
Branch correlation and the assumption of neutral mtDNA evolution
Our conclusions are based on the assumption that the mtDNA evolves neutrally and proportionally to time. This is known not to be the general case, due to heterogeneity of rates among branches (e.g., [62, 72, 73]), the effect of incomplete sampling and node density [63, 64, 74], and in some cases active selection on the mitochondrial genome [75–79].
Our conclusions should be robust to random deviations from a clock-like behavior, as they should reduce, but not bias, the informative signal. It must be noted that our null hypothesis (i.e. neutral change) is based on a positive result, not on a negative one: we do find a significant correlation between the evolutionary rates of mtDNA and those of some character systems, a result highly relevant in itself [32, 33]. We did not find evidence of a node density effect, and the possible bias introduced by incomplete sampling could be expected to affect in a parallel way all character systems when plotted on the same tree topology, and thus not to introduce any strong bias in our results [19, 32, 80]. Although there are numerous examples of selection on mitochondrial genes (see  for a recent review), it is still unknown to what extend this is a widespread phenomenon that could undermine the general assumption of neutrality, especially among closely related species with a similar biology.
Methodological issues of evolutionary rate comparisons
We used three main statistical approaches to compare the evolutionary rate of the studied character systems (Table 4): 1 correlation and ANOVA; 2) comparison of the K-scores against a null distribution; and 3) Mantel test for the correlation of the matrices of patristic distances. The first two were measures of the comparison of individual branches of the tree, and the last one of the accumulated change over the pathway from the root to the tip of each terminal taxa (Omland 1997). Standard statistics can only be applied under certain restrictions, which were not met by most of our character systems (i.e. normality; uniformity of variance ; see Additional file 8). For those variables showing approximate normal distributions, correlations were all not significant (Table 1), possibly due to the heterogeneity in the distribution of the change across the branches (Additional file 8) or the limited number of branches involved . It seems reasonable to assume that the length of individual branches is subjected to multiple sources of variation, and that it would be difficult that simple correlations of a limited set of characters and specimens reflect any potential relationship. The comparison of the K-scores proved to be more powerful, discriminating between character systems. The absolute values of the K-score would suggest that the optimized tree based on paramere outline was closer to the mtDNA tree (K-score = 0.293) than that optimized on the body morphospace (K-score = 0.297) (Table 3). However, the randomization procedure revealed that the tree optimized on the paramere outline was not significantly different from a random tree, since the observed K-score was among the 95% of the null random distribution of values. Although the computation of the K-score includes a scaling factor, this depends on the distribution of the total amount of variation in the individual branches, so the final K-scores cannot be considered a non-dimensional number and thus comparable across character systems with different scale of change . It is thus necessary to associate the K-scores to a null distribution in order to obtain meaningful results in this comparative framework. For some trees their null distribution appeared clearly bimodal (Figure 4). This suggests the presence of two discrete "stable states" within the tree, determined by the presence of two long branches: when the random tree placed these two branches together, the overall K-score decreased in a marked step, creating a bimodal distribution. The effect of long branches in originating bimodal (or trimodal) null distributions can be confirmed by simulations of trees (J. Castresana and V. Soria-Carrasco, pers. com. 2009).
Trait change and hierarchical levels of the tree
A surprising result of our study was the amount of variation in the morphometric measures within what is usually recognized as a unique morphospecies, despite the reduced mtDNA distances (Figure 3; Additional files 5, 10). There was a perfect correspondence between morphologically recognized nominal species and mtDNA clusters of the sampled specimens (Figure 2), but some individuals within these clusters had relatively long branches for both the paramere outline and the body morphospace (Figure 3). As already noted, the difficult taxonomy of the group, and the limited amount of available data and specimens, makes difficult the rigorous assessment of the partition of the measured variation among and within species. We plotted the morphological vs. mitochondrial inter- and infraspecific distances as an exploratory tool ), and although the results have to be taken cautiously (due to the non-independence of the distances), they may provide useful insights into the data. The most striking difference was in the structural morphological characters, with a strong positive relationship with mtDNA distances at the interspecies level, but slightly negative (although significant) at the intraspecific level, possibly reflecting the choice of diagnostic characters of taxonomic use. It is also interesting to note the slightly (but significant) negative slope of the relationship between body morphospace and mtDNA distances, possibly reflecting the "saturation" of this character system (, see above) (Figure 5; Additional file 9, 10).
In any case, our basic approach does not require the delimitation of species, since it does not include species-counts or a different treatment of within- and between-species variation. It is thus not affected by potential taxonomic artefacts due to the difficult recognizability and complex taxonomy of the group [39, 40], and does not depend on the partition the morphological variation (e.g., ).
To compare the rates of change of the character systems over the tree we used statistics giving a global measure of difference, usually based on multiple pairwise comparisons (K-scores and Mantel test on patristic distance matrices). We did not study in detail the distribution of the change in different parts of the tree, such as e.g. between sister species, or between internal and terminal branches. The examination of the branch lengths in Figure 3 suggests some interesting patterns, such as the strong divergence between the syntopically co-occurring H. lineolata and H. fulvipennis, or the frequency of deep branches with length = 0 in the paramere outline. However, the number of taxa included was insufficient to draw any conclusion from these observations, which would require more comprehensive data to be tested .
In this work we outline an approach to compare evolutionary rates of different character systems in a clade of Hymenoplia chafers with the aim to study their speciation mechanisms. Although based on indirect evidence, we were able to obtain insightful conclusions on the prevalent mode of speciation, which could be the base for subsequent, more detailed ecological or biological studies. Thus, it seems likely that changes of the male genitalia are shaped by sexual selection, and that body morphospace is subject to either functional or developmental constraints with signs of stabilizing selection that cannot be unambiguously associated with speciation (Table 4). On the contrary, more structural morphological characters, including those used for species delimitation and recognition (other than male genitalia), seem to evolve at a rate proportional to time and are thus not related to speciation. These conclusions would be supported by the apparently widely occurring sympatry of the Hymenoplia species (Figure 1), although the taxa seem not to co-occur syntopically (see above). The possibility that sexual selection was an important factor in the diversification of the genus Hymenoplia is highly relevant, as much of recent biological diversity of scarabs in the Mediterranean region is attributed to vicariance only (e.g., [81, 82]). The approach presented here has the additional advantage of not being dependent on species definition, and could be of potential use for comparative studies not only related to speciation, but to any question requiring the comparison of the amount of change of different character systems over a phylogenetic tree.
Mayr E: Animal species and evolution. 1963, Cambridge, Mass: Harvard University Press
Coyne JA, Orr HA: Speciation. 2004, Sunderland MA, Sinauer Associates
Barraclough TG, Vogler AP: Detecting the geographical pattern of speciation from species-level phylogenies. American Naturalist. 2000, 155: 419-434. 10.1086/303332.
Phillimore AB, Orme CDL, Thomas GH, Blackburn TM, Bennett PM, Gaston KJ, Owens IPE: Sympatric speciation in birds is rare: Insights from range data and simulations. American Naturalist. 2008, 171: 646-657. 10.1086/587074.
Barraclough TG, Harvey PH, Nee S: Sexual selection and taxonomic diversity in passerine birds. Proceedings of the Royal Society London B-Biological Sciences. 1995, 259: 211-215. 10.1098/rspb.1995.0031.
Arnqvist G, Edvardsson M, Friberg U, Nilsson T: Sexual conflict promotes speciation in insects. Proceedings of the National Academy of Sciences of the United States of America. 2000, 97: 10460-10464. 10.1073/pnas.97.19.10460.
Morrow EH, Pitcher TE, Arnqvist G: No evidence that sexual selection is an 'engine of speciation' in birds. Ecology Letters. 2003, 6: 228-234. 10.1046/j.1461-0248.2003.00418.x.
Barraclough TG, Hogan JE, Vogler AP: Testing whether ecological factors promote cladogenesis in a group of tiger beetles (Coleoptera: Cicindelidae). Proceedings of the Royal Society London B-Biological Sciences. 1999, 266: 1061-1067. 10.1098/rspb.1999.0744.
Eberhard WG: Sexual Selection and Animal Genitalia. 1985, Harvard University Press, Cambridge, MA
Arnqvist G: Comparative evidence for the evolution of genitalia by sexual selection. Nature. 1998, 393: 784-786. 10.1038/31689.
McPeek MA, Shen L, Torrey JZ, Farid H: The tempo and mode of three-dimensional morphological evolution in male reproductive structures. Am Nat. 2008, 171 (5): E158-E178. 10.1086/587076.
Kozak KH, Wiens JJ: Does niche conservatism promote speciation? A case study in North American salamanders. Evolution. 2006, 60: 2604-2621.
Schneider CJ, Smith TB, Larison B, Moritz C: A test of alternative models of diversification in tropical rainforests: Ecological gradients vs. rainforest refugia. Proceedings of the National Academy of Sciences of the United States of America. 1999, 96: 13869-13873. 10.1073/pnas.96.24.13869.
Turgeon J, Stoks B, Thurm RA, Brown JM, McPeek MA: Simultaneeous Quaternary radiations of three damselfly clades across the Holarctic. American Naturalist. 2005, 165: 78-107. 10.1086/428682.
Ahrens D: The phylogeny of Sericini and their position within the Scarabaeidae based on morphological characters (Coleoptera: Scarabaeidae). Systematic Entomology. 2006, 31: 113-144. 10.1111/j.1365-3113.2005.00307.x.
Losos JB, Glor RE: Phylogenetic comparative methods and the geography of speciation. Trends in Ecology and Evolution. 2003, 18: 220-227. 10.1016/S0169-5347(03)00037-5.
Polly PD: On the simulation of the evolution of morphological shape: multivariate shape under selection and drift. Palaeontologia Electronica Polly PD: . 2004, 7 (7A): 28-[http://palaeo-electronica.org/2004_2/evo/issue2_04.htm]
Losos JB: Phylogenetic niche conservatism, phylogenetic signal and the relationship between phylogenetic relatedness and ecological similarity among species. Ecology Letters. 2008, 11: 995-1007. 10.1111/j.1461-0248.2008.01229.x.
Ricklefs RE: Time, species, and the generation of trait variance in clades. Systematic Biology. 2006, 55: 151-159. 10.1080/10635150500431205.
Kimura M: The neutral theory of molecular evolution. 1983, Cambridge Univ Press, Cambridge
Ballard JWO, Kreitman M: Is mitochondrial DNA a strictly neutral marker?. Trends in Ecology and Evolution. 1995, 10: 485-488. 10.1016/S0169-5347(00)89195-8.
Felsenstein J: Phylogenies and the comparative method. American Naturalist. 1985, 125: 1-15. 10.1086/284325.
Polly PD: On morphological clocks and paleophylogeography: towards a timescale for Sorex hybrid zones. Genetica. 2001, 112-113: 339-357. 10.1023/A:1013395907225.
Eberhard WG: Rapid divergent evolution of sexual morphology: comparative test of antagonistic coevolution and traditional female choice. Evolution. 2004, 58: 1947-1970.
Fontaneto D, Herniou EA, Boschetti C, Caprioli M, Melone G, Ricci C, Barraclough TG: Independently evolving species in asexual Bdelloid Rotifers. PLoS Biology. 2007, 5 (4): e87-10.1371/journal.pbio.0050087.
Lande R, Arnold SJ: The measurement of selection on correlated characters. Evolution. 1983, 37: 1210-1226. 10.2307/2408842.
Phillips PC, Arnold SJ: Visualizing multivariate selection. Evolution. 1989, 43: 1209-1222. 10.2307/2409357.
Wagner PJ: Exhaustion of morphologic character states among fossil taxa. Evolution. 2000, 54: 365-386.
Omland KE: Character congruence between a molecular and a morphological phylogeny for dabbling ducks (Anas). Systematic Biology. 1994, 43: 369-386. 10.1093/sysbio/43.3.369.
Omland KE: Correlated rates of molecular and morphological evolution. Evolution. 1997, 51: 1381-1393. 10.2307/2411190.
Wilson AC, Carlson SS, White TJ: Biochemical evolution. Annual Review of Biochemistry. 1977, 46: 573-639. 10.1146/annurev.bi.46.070177.003041.
Bromham L, Woolfit M, Lee MSY, Rambaut A: Testing the relationship between morphological and molecular rates of change along phylogenies. Evolution. 2002, 56: 1921-1930.
Davies TJ, Savolainen V: Neutral theory, phylogenies, and the relationship between phenotypic change and evolutionary rates. Evolution. 2006, 60: 476-483.
Galante E: Contribucíon al conocimiento de los escarabaeidos florícolas de la provincia de Salamanca (Col. Scarabaeioidea). Boletín de la Asociación Española de Entomologia. 1981, 5: 151-160.
Sanchez-Pinero F, Ruiz JL, Avila JM: Aportaciones a la distribucion y biologia de los Sericinae (Coleoptera: Scarabaeoidea, Melolonthidae) del sudeste de la Peninsula Iberica. Boletín de la Asociación Española de Entomologia. 1994, 18: 31-39.
Scholtz CH, Grebennikov VV: Scarabaeoidea Latreille, 1802. Coleoptera, Beetles. vol. 1: Morphology and Systematics (Archostemata, Adephaga, Myxophaga, Polyphaga partim). Edited by: Beutel RG, Leschen RAB, Berlin W. 2005, Handbook of Zoology. Arthropoda: Insecta part 38. New York: de Gruyter, IV: 367-426.
Ahrens D, Vogler AP: Towards the phylogeny of chafers (Sericini): analysis of alignment-variable sequences and the evolution of segment numbers in the antennal club. Molecular Phylogenetics and Evolution. 2008, 47: 783-798. 10.1016/j.ympev.2008.02.010.
Ahrens D: Taxonomic changes and an updated Catalogue of Palearctic Sericini (Coleoptera: Scarabaeidae: Melolonthinae). Zootaxa. 2007, 1504: 1-51.
Baguena Corella L: Hymenoplia Eschscholtz de la fauna Iberica. Eos. 1954, 30: 7-46.
Baraud J: Coléoptères Scarabaeoidea d'Europe (Faune de France, 78). Fédération française des Sociétés de Sciences naturelles et Société linnéenne de Lyon, Lyon. 1992
Baguena Corella L: Scarabaeoidea de la Fauna Ibero-Balear y Pirenaica. Madrid. 1967
Swofford DL: PAUP*: Phylogenetic Analysis using Parsimony. Version 4.0b. 2002, Sinauer Associates, Sunderland, MA
Iwata H, Ukai Y: SHAPE: a computer program package for quantitative evaluation of biological shapes based on elliptic Fourier descriptors. Journal of Heredity. 2002, 93: 385-385. 10.1093/jhered/93.5.384.
Kuhl FP, Giardina CR: Elliptic Fourier features of a closed contour. Computer Graphics and Image Processing. 1982, 18: 236-258. 10.1016/0146-664X(82)90034-X.
Ahrens D, Lago PK: Directed asymmetry reversal of male copulatory organs in chafer beetles (Coleoptera: Scarabaeidae) - implications on left-right polarity determination in insect terminalia. Journal of Zoological Systematics and Evolutionary Research. 2008, 46: 110-117. 10.1111/j.1439-0469.2007.00449.x.
Márquez EJ, Knowles LL: Correlated evolution of multivariate traits: detecting co-divergence across multiple dimensions. Journal of Evolutionary Biology. 2007, 20: 2334-2348. 10.1111/j.1420-9101.2007.01415.x.
Hammer Ø, Harper DAT, Ryan PD: PAST: Paleontological Statistics Software Package for Education and Data Analysis. Palaeontologia Electronica. 2001, 4 (1): 1-9.
Boore JL: Mitochondrial gene arrangement source guide, version 6.0. 2001, Department of the Environment, Joint Genome Institute, Walnut Creek, California, [http://www.uni-ulm.de/fileadmin/website_uni_ulm/iui.inst.190/Forschung/Projekte/GENESIS/MGA_Guide.pdf]
Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P: Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Annals of the Entomological Society of America. 1994, 87: 651-701.
Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG: The ClustalX windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Research. 1997, 25: 4876-4882. 10.1093/nar/25.24.4876.
Guindon S, Gascuel O: A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Systematic Biology. 2003, 52: 696-704. 10.1080/10635150390235520.
Posada D, Crandall KA: Modeltest: testing the model of DNA substitution. Bioinformatics. 1998, 14: 817-818. 10.1093/bioinformatics/14.9.817.
Goloboff P, Farris S, Nixon K: TNT (Tree analysis using New Technology). Cladistics. 2004, 20: 84-
Ronquist F, Huelsenbeck JP: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.
Felsenstein J: Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985, 39: 783-791. 10.2307/2408678.
Cooper N, Purvis A: What factors shape rates of phenotypic evolution? A comparative study of cranial morphology of four mammalian clades. Journal of Evolutionary Biology. 2009, 22: 1024-1035. 10.1111/j.1420-9101.2009.01714.x.
Soria-Carrasco V, Talavera G, Igea J, Castresana J: The K tree score: quantification of differences in the relative branch length and topology of phylogenetic trees. Bioinformatics. 2007, 23: 2954-2956. 10.1093/bioinformatics/btm466.
Kuhner MK, Felsenstein J: A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates. Molecular Biology and Evolution. 1994, 11: 459-468.
Gotelli NJ, Graves GR: Null models in ecology. 1996, Smithsonian Inst Press
Sokal RR, Rolf FJ: The principles and practice of statistics in biological research. 1995, New York: Freemann and Company
Bonnet E, Van de Peer Y: zt: a software tool for simple and partial Mantel tests. Journal of Statistical Software . 2002, 7 (10): 1-12. [http://www.jstatsoft.org/v07/i10]
Fitch WM, Margoliash E: A method for estimating the number of invariant amino acid positions in a gene using cytochrome c as a model case. Biochemical Genetics. 1967, 1: 65-71. 10.1007/BF00487738.
Sanderson MJ: Estimating rates of speciation and evolution: a bias due to homoplasy. Cladistics. 1990, 6: 387-391. 10.1111/j.1096-0031.1990.tb00554.x.
Venditti C, Meade A, Pagel M: Detecting the node density artifact in phylogeny reconstruction. Systematic Biology. 2006, 55: 637-644. 10.1080/10635150600865567.
Venditti C, Meade A, Pagel M: Phylogenetic mixture models can reduce node-density artifacts. Systematic Biology. 2008, 57: 286-293. 10.1080/10635150802044045.
Emerson SB: Male secondary sexual characteristics, sexual selection, and molecular divergence in fanged ranid frogs of Southeast Asia. Zoological Journal of Linnaean Society. 1998, 122: 537-553. 10.1111/j.1096-3642.1998.tb02162.x.
Monteiro LR, Gomes JL: Morphological divergence rates tests for natural selection: Uncertainty of parameter estimation and robustness of results. Genetics and Molecular Biology. 2005, 28: 345-355. 10.1590/S1415-47572005000200028.
Gosden TP, Svensson EI: Spatial and temporal dynamics in a sexual selection mosaic. Evolution. 2008, 62: 845-856. 10.1111/j.1558-5646.2008.00323.x.
Svensson EI, Eroukhmanoff F, Friberg M: Effects of natural and sexual selection on adaptive population divergence and premating isolation in a damselfly. Evolution. 2006, 60: 1242-1253.
Arnqvist G: The evolution of animal genitalia: distinguishing between hypotheses by single species studies. Biological Journal of Linnaean Society. 1997, 60: 365-379. 10.1111/j.1095-8312.1997.tb01501.x.
Ahrens D: The phylogeny of the genus Lasioserica inferred from adult morphology - implications on the evolution of montane fauna of the south Asian orogenic belt (Coleoptera: Scarabaeidae: Sericini). Journal of Zoological Systematics and Evolutionary Research. 2006, 44: 34-53. 10.1111/j.1439-0469.2005.00340.x.
Uzzell T, Corbin KW: Fitting discrete probability distributions to evolutionary events. Science. 1971, 172: 1089-1096. 10.1126/science.172.3988.1089.
Thorne JL, Kishino H, Painter IS: Estimating the rate of evolution of the rate of molecular evolution. Molecular Biology and Evolution. 1998, 15: 1647-1657.
Fitch WM, Bruschi M: The evolution of prokaryotic ferredoxins - with a general method correcting for unobserved substitutions in less branched lineages. Molecular Biology and Evolution. 1987, 4: 381-394.
MacRae AF, Anderson WW: Evidence for non-neutrality of mitochondrial DNA haplotypes in Drosophila pseudoobscura. Genetics. 1988, 120: 485-494.
Fos MA, Latorre A, Moya A: Mitochondrial DNA evolution in experimental populations of Drosophila subobscura. Proceedings of the National Academy of Sciences of the United States of America. 1990, 87: 4198-4201. 10.1073/pnas.87.11.4198.
Nigro L, Prout T: Is there selection on RFLP differences in mitochondrial DNA?. Genetics. 1990, 125: 551-555.
Rand DM: The Units of Selection of Mitochondrial DNA. Annual Review of Ecology and Systematics. 2001, 32: 415-448. 10.1146/annurev.ecolsys.32.081501.114109.
Dowling DK, Friberg U, Lindell J: Evolutionary implications of non-neutral mitochondrial genetic variation. Trends in Ecology and Evolution. 2008, 23: 546-554. 10.1016/j.tree.2008.05.011.
Webster AJ, Payne RJH, Pagel M: Molecular phylogenies link rates of evolution and speciation. Science. 2003, 301: 478-10.1126/science.1083202.
Martin-Piera F, Sanmartin I: Biogeografía de áreas y biogeografía de artrópodos Holárticos y Mediterráneos. Evolution and phylogeny of Arthropoda. Edited by: Melic A, De Haro JJ, Méndez M, Ribera I. 1999, Boletín de la Sociedad Entomológica Aragonesa, 26: 535-560.
Mico E, Sanmartín I, Galante E: Mediterranean diversification of the grass-feeding Anisopliina beetles (Scarabaeidae, Rutelinae, Anomalini) as inferred by bootstrap-averaged dispersal-vicariance analysis. Journal of Biogeography. 2009, 36: 546-560. 10.1111/j.1365-2699.2008.02010.x.
We are grateful to M. Alcobendas for assistance during laboratory work at the MNCN, and to I. Izquierdo and M. Paris for support during collection work in the MNCN. We thank S. Fabrizi for assistance during fieldwork; T. Branco and S. Montagud for collecting some of the samples; A. Cieslak, T. Barraclough, J. Castresana, V. Soria and two anonymous referees for helpful comments on the manuscript and discussions; and J. Castresana and V. Soria for providing a PERL script for the matrix randomization. Laboratory analysis and collection work at MNCN was supported by a SYNTHESYS grant (ES-TAF 2670) to DA, fieldwork in Spain and DA work was supported by the German Science Foundation (DFG/AH175/1), IR work was supported by grant CGL2007-61665.
DA collected the samples and performed the lab work, dissections, photography and morphometric analyses. DA and RI performed the statistical and phylogenetic comparative analyses of the data and contributed equally on the different drafts of the manuscript; figures and tables were done by DA. Both authors read and approved the final manuscript.