During the last decade, molecular phylogenetics based on the comparative analysis of markers from mitochondrial and nuclear genomes allowed for the revision of mammalian systematics [1,2]. Sequencing efforts initially focused on mitochondrial DNA (mtDNA) due to the availability of conserved primers, the presence of rapidly evolving sites, and a reduced effective population size inducing the rapid fixation of variants between subsequent speciation events. Some protein-coding genes such as cytochrome b [CYB] have emerged as central tools in investigating intraspecific [e.g. [3-6]] to ordinal-level [e.g. [7-9,6]] evolutionary relationships. To date (January 2006), ca. 5,000 complete CYB entries are available for mammals in public data bases, and furthermore, another mitochondrial gene – the cytochrome c oxidase I – has been proposed as the marker for animal DNA barcoding [10]. However, the strong nucleotide saturation encountered at third codon positions as well as a high sensitivity to taxon sampling soon caused the questioning of the utility of mitochondrial protein-coding genes in general, and CYB in particular for resolving deep phylogenies [11,12]. Owing to their different evolutionary patterns – a less biased base composition and lower saturation than mitochondrial genes [13,14] – nuclear genes have represented a reasonable alternative to mtDNA for reconstructing deep-level mammalian phylogenies [[15-21]; but see [22]]. By providing complementary information, mitochondrial and nuclear genes are thus generally employed at different levels, low-level phylogeny and taxonomy for the former, and deeper-level for the latter.

Several mtDNA-based studies of mammalian systematics have shed light on the difficulties experienced in using CYB to resolve phylogenies at lower taxonomic levels. Several works at the family level recovered multifurcations among genera or even species, which led to the conclusion of rapid, near-simultaneous divergences of multiple lineages (= star-phylogeny), without any time left for synapomorphies to accumulate in mtDNA. This interpretation has been put forward for some speciose clades of rodents, such as Ctenomyidae [23], Echimyidae [24], or Sigmodontinae [25]. The star-phylogeny hypothesis was also supported by the lack of resolution for relationships among genera, contrasting with well-defined nodes above and below multifurcations [23,24,26]. However, an alternative hypothesis suggests that CYB sequences have undergone substitution saturation throughout the course of speciation events, leading to a loss of the original phylogenetic signal, and producing a soft-polytomy [27]. In this case, the use of slower-evolving nuclear DNA (nuDNA) might represent a better choice for resolving such recent radiations. Supporting this hypothesis, sequences from such markers, either exons or introns, have recently improved among-families [28-30], or even among-genera [31,32] phylogenetics of spiny rats, squirrels, whales, weasels, and spiral-horn antelopes.

In the present study, we focused on the evolutionary history of voles, lemmings, and muskrats. These animals belong to the Arvicolinae, one of the six Cricetidae subfamilies of the highly-diverse Muroidea, which encompasses one third of all rodent species. The Arvicolinae represent themselves one of the most impressive placental radiations in the Northern Hemisphere, consisting of 151 species and 28 genera [33]. As stated by [33], the explosiveness and recency of arvicoline evolution can be dramatically highlighted by the more than 60 species of Microtus and the inconsistency of their systematic treatment. Phylogenetic reconstructions based on mitochondrial sequences (CYB and NADH-dehydrogenase 4 [ND4]) provided unresolved topologies at two different levels: (1) among arvicoline genera, and (2) among species of Microtus [34-36]. These results reinforced observations previously made by paleontologists [37,38] who concluded that two successive pulses of speciation occurred during the evolutionary course of arvicoline rodents over the last 10 million years (Myr).

Following the encouraging results from recent studies [28-32], we thus compared the relative contribution of a nuclear marker – exon 10 of the growth hormone receptor (GHR) gene – to one of the mitochondrial CYB for inferring phylogenetic relationships among major Arvicolinae lineages. We performed analyses in a probability framework because (i) maximum likelihood (ML) and Bayesian methods are indeed based upon explicit models of sequence evolution, (ii) they allow the definition of independent models to reflect the contrasting substitution patterns among the three codon positions of GHR and CYB, and (iii) they have proven to be robust to a number of systematic biases during phylogenetic reconstruction [39-41]. In addition to such inter-sites comparisons, we evaluated the amount of saturation at each codon position. Using the variable-length bootstrap method, we also estimated the amount of sites required by mitochondrial, CYB-like, and nuclear, GHR-like data to resolve recent and deeper nodes of the evolution of voles and lemmings.

Figure 1 Maximum posterior probability trees reconstructed from the mitochondrial CYB (left, A) and nuclear GHR (right, B) sequences. Two reliability indices are given on nodes: the Bayesian posterior probabilities/the maximum likelihood bootstrap percentages. (more ...)

We also scanned the GHR alignment for potential indel signatures. A 3-bp deletion at position 1,925 of the GHR sequence of Rattus norvegicus (Accession number NM_017094) is diagnostic for the monophyly of Arvicolinae. A 3-bp in-frame insertion – a direct AGC repeat – at position 1,778 of Rattus norvegicus (Accession number NM_017094) is shared by Neodon + Phaiomys, and independently by Blanfordimys + Lasiopodomys.

Independent paleontological and molecular studies respectively estimated that divergence times among some arvicoline genera occurred 3–5 Myr ago (Mya) [37] to 5–9 Mya [30,34,42]. The adaptive radiation of Microtus has been dated to approximately 2 Myr by most paleontologists [37,43] but there is molecular evidence for splits between 2.6 to 4.4 Mya [34]. The GHR topology is here more resolved than the CYB one for intergeneric speciation events, corresponding to mean uncorrected pairwise divergences of 3.8 % ± 0.9 (GHR) versus 14.2 % ± 1.6 (CYB). For the 2–5 Myr old Microtus radiation, corresponding to divergences of 2.2% ± 0.6 (GHR) and 11.7 % ± 1.2 (CYB), the GHR topology is at least as resolved as the CYB one. These results emphasize previous conclusions [44] that highlighted the superiority of nuclear genes over mitochondrial genes even for divergences spanning the last 5–10 Myr.

Figure 2 Maximum posterior probability tree reconstructed from the combination of the mitochondrial CYB and nuclear GHR sequences. Three reliability indices are given on nodes: the Bayesian posterior probabilities/the bootstrapped Bayesian posterior probabilities/the (more ...)

Clades moderately to strongly supported in the GHR tree data benefit from similar support values after the combination of GHR and CYB sequences. Actually, a substantial increase of support is found for nodes defining the basal position of Prometheomys, and the monophyly of Microtus s.l. is reinforced. On the opposite, the sister-clade relationship recovered between Myodini and Arvicolini is less supported in the combined tree than in the GHR one. Expectedly, the M. arvalis + M. guentheri clade, recovered by CYB data, is similarly supported in the combination.

To summarize, the combination of CYB and GHR sequences provides a phylogeny where two-third of the nodes benefit from a medium to high support (Figure 2). The addition of nuclear sequences also allows the division of a single, multi-taxa polytomy into two noticeably smaller multifurcations. The first includes Arvicolini, Ondatra, Lemmus, Dicrostonyx, and Phenacomys, although the latter two genera are consistently found as sister-group across the different reconstructions performed. The second multifurcation includes phylogenetic relationships among Microtus s.l. subclades with evidence for a deeper divergence of the Neodon + Phaiomys lineage.

In addition, a first molecular evidence is here provided for the phylogenetic position of enigmatic genera suspected to have differentiated early within the arvicoline radiation. The long-clawed mole-vole, Prometheomys, currently localised in the Caucasus Mountains, was considered as a relic of an archaic lineage formerly widespread throughout most of Eurasia [33,48-51]. This view is corroborated by its first emergence in the GHR and combined phylogenetic trees (Figures 1B and 2). The collared lemming, Dicrostonyx as well as the heather vole, Phenacomys, exhibit plesiomorphic morphological traits which have compromised their systematic affiliation to other arvicolines. Dicrostonyx was initially grouped with true lemmings whereas Phenacomys was included in either Arvicolini or Myodini [33]. By obtaining a moderately supported sister-taxa relationship between Dicrostonyx and Phenacomys, our study is in accordance with previous result from the analysis of highly repetitive DNA (LINE-1) elements [52]. Although unexpected, this clade might find a biogeographical explanation as these two taxa are primarily known – from a palaeontological viewpoint – in the Arctic region of Eurasia and North America [37,53,54].

To summarize our results from a taxonomic standpoint, the phylogenetic relationships discussed above agree with the tribal recognition of Lemmini (Lemmus), Ondatrini (Ondatra), Dicrostonychini (Dicrostonyx), Phenacomyini (Phenacomys), Prometheomyini (Prometheomys), Myodini (Myodes and Eothenomys) and Arvicolini (Arvicola, Chionomys, Microtus, Neodon, Phaiomys, Blanfordimys, and Lasiopodomys) [33].

Actually, three lineages are clearly identified within Microtus. Firstly, strong evidence is provided for a Neodon + Phaiomys subgenera clade, contradicting former studies [33] which, using molar and other external and cranial contrasts, rejected any close phylogenetic affinity between them. Secondly, Microtus kikuchii, M. oeconomus, and M. middendorffi form a well-supported association. Independent source of phylogenetic data – here nuclear GHR sequences – thus confirm the "Asian lineage" recently identified on the basis of mitochondrial phylogenies and chromosomal data [35,36,55], as well as the monophyly of the subgenus Alexandromys, as redefined in [33]. The five nearctic species included in our analysis and representatives of various subgenera (Aulacomys, Mynomes, Pedomys, as well as Microtus longicaudus and M. chrotorrhinus whose systematic status is unclear) are associated with some support in the combined analysis. The recognition of an American clade within Microtus strengthens previous results from biochemical analyses [37] and mitochondrial phylogenies [35,36]. Third, some support for a Terricola + Microtus clade is recovered, suggesting close phylogenetic affinities among Western Palaearctic taxa, and reinforcing the conclusions of [36].

Table 1 Bayesian estimates of DNA substitution model parameters for CYB, GHR, and their codon positions.

Table 2 Akaike Information Criterion (AIC) sensitivity analyses about the effect of model parameters on log-likelihood gains provided by codon partitions.

Table 3 Maximum likelihood bootstrap support for a selection of nodes, computed according to codon positions 1, 2, or 3 of CYB and GHR.

Figure 3 Saturation plots of the number of observed differences as a function of the numbers of inferred substitutions for each pair of sequences at each codon position of CYB and GHR genes. The Y = X straight line corresponds to the situation where there is no (more ...)

5. Perspectives: Which markers to be used for resolving the arvicoline radiation?

The resolved GHR tree is perhaps a first step in challenging the hypothesis of a hard-polytomy for Arvicolinae genera. As shown for other placental taxa – spiral-horn antelopes [32], hares [57], or bears [59] – our results suggest that phylogenetic signal has been progressively obliterated by higher mtDNA rates whereas it has persisted in at least one slower-evolving nuclear gene. However, the GHR gene tree is based on a single nuclear locus (GHR), and it does not necessarily reflect the organismal history of arvicoline rodents. Actually, the resolution of the polytomy is further complicated by the eventual incongruence of individual gene trees with species trees due to incomplete lineage sorting [60]. Consequently, the comparison of multiple gene trees is the only way to get a central tendency which could be interpreted as the species tree [61-63]. If multiple gene trees agree on a phylogenetic structure, we could then definitively dismiss the star-phylogeny hypothesis for the Arvicolinae. In addition to allowing us to make predictions on the amount of CYB/GHR data required to decipher the phylogeny of the Arvicolinae, variable length bootstrap analyses help us to identify which kind of markers should be considered for further sequencing.

Figure 4 Plots of the variable-length bootstrap percentages for CYB (filled circles), GHR (filled triangles), and combined data (open squares) for four nodes recovered in Bayesian and ML analyses (see Fig. 2). The X-axis is the number of sites resampled (in 250-bp (more ...)

The lack of phylogenetic signal as well as strong saturation and large among-sites heterogeneity questions the utility of CYB sequences for investigating phylogenetic relationships among taxa, like voles and lemmings, which have experienced relatively recent and rapid radiations. Actually, we concur with others [64,65] on the usefulness of CYB sequences – including third codon positions – when they are combined with less saturated sequences. Firstly, we have shown that CYB sometimes holds phylogenetic signal for specific nodes where GHR is non-informative (node O, see Results & Discussion, paragraph 1). Although our study mainly focuses on intergeneric phylogenetic relationships, we could expect that CYB contribution would increase with nodes younger than 2 Myr. Secondly, VLB curves clearly indicate that data combination (CYB + GHR) provides slightly higher bootstrap support than does GHR alone, especially when the number of resampled sites is higher than 2,000 (e.g., nodes L and P: Figure 4). Results of VLB analysis and comparisons of GHR and combined topologies show that the CYB signal does not contradict the GHR one. Moreover, the incorporation of CYB data to nuclear data does not lead to the degradation of phylogenetic signal but slightly improves it, revealing the existence of weak and initially hidden information within CYB. Because complete mitochondrial genomes include sequences undergoing various evolutionary pressures (e.g. CYB, but also COX1, or rRNAs), and are available for an increasing number of taxa, particularly among mammals, we suggest the evaluation of their ability to resolve rapid radiations at low taxonomic level [66], and a comparison of their efficiency to that of nuclear loci.

The nuclear exon 10 of the GHR performs better than the mitochondrial CYB for resolving Arvicolinae phylogenetic relationships. Support is found for a sister-group relationship between red-backed voles (Myodini) and a clade including water (Arvicola), snow (Chionomys), and meadow voles (Microtus and allies). Lemmings (Lemmus and Dicrostonyx) are found polyphyletic while the Caucasian long-clawed vole (Prometheomys) is among the basalmost arvicoline genera. Contrary to recent taxonomic suggestions, we do not obtain support for splitting Blanfordimys, Lasiopodomys, Neodon, and Phaiomys from Microtus. We concur with others [68-70] that the higher quality of nuclear genes resides in higher values of gamma parameter, uniform and stationary base compositions, and more uniform nucleotide substitution probabilities. The usefulness of nuclear exons for investigating 2–10 Myr old evolutionary histories is probably due to highly informative third codon positions, which keep a good compromise in their evolutionary rate: rapid enough to accumulate synapomorphies, yet slow enough to be not affected by substitutional saturation. The sequencing of nuclear DNA data in an order of magnitude three to five times greater than the size of GHR exon 10 might resolve most nodes of the Arvicolinae phylogeny provided that multiple genetically independent gene trees agree on the phylogenetic structure.

All ethanol-preserved arvicoline samples were stored in the mammalian tissue collection of the Institut des Sciences de l'Evolution de Montpellier [72]. Total DNA was extracted using the QIAamp DNA mini kit (Qiagen). Part of the GHR exon 10 was amplified and sequenced using the primers GHR5 forward (5' GGCRTTCATGAYAACTACAAACCTGACYTC 3') and GHR6 reverse (5' GAGGAGAGGAACCTTCTTTTTWTCAGGC 3'), and GHR3 forward (5' GACTTTATGCYCARGTRAG 3') and GHR4 reverse (5'-CTYACYTGRGCATAAAAGTC 3'). PCR conditions were 95°C 5 min, followed by 95°C 30 sec, 61°C 30 sec, 72°C 1 min (5 times), then 95°C 30 sec, 59°C 30 sec, 72°C 1 min (5 times), followed by 95°C 30 sec, 57°C 30 sec, 72°C 1 min (5 times), then 95°C 30 sec, 55°C 30 sec, 72°C 1 min (5 times), and then 95°C 30 sec, 53°C 30 sec, 72°C 1 min (20 times), with a final extension at 72°C 5 min.

The amplification and sequencing of the CYB were conducted using primers MVZ05 and MVZ14 [73] and additional internal ones MVZ16 [73] and H8 [74]. PCR products for GHR and CYB were purified from 1% agarose gels using Amicon Ultrafree-DNA columns (Millipore) and sequenced on both strands using automatic sequencing (Big Dye Terminator cycle kit) on an ABI 310 (PE Applied Biosystems).

All taxa included in our study were represented by both CYB and GHR sequences. A 921 bp segment of the exon 10 of GHR was sequenced for 22 specimens and the complete CYB (1140 bp) for 7 specimens. The arvicoline sequences new to this study have been deposited in the EMBL data bank, and we also used previously published sequences when available (see Table 1).

Phylogenetic analyses of CYB and GHR alignments were conducted under the maximum likelihood (ML) and Bayesian methods, using PAUP* [76] version 4b10, PHYML [77] version 2.4.4, and MrBayes version 3.04 [78]. Moreover, MrBayes provided the opportunity to run analyses assuming different models of sequence evolution for each predefined partition, thus permitting the parameters estimated for each model of sequence evolution to be directly compared among codon positions and between genes. When CYB and GHR sequences were combined, sequences from different specimens of the same species were sometimes used. We assumed that the phylogeographic structure detected in some arvicoline species [45,79] was negligible as compared to the level of genetic divergence between the distinct Microtus subgenera and even arvicoline tribes here compared.

The program Modeltest [80] version 3.06, was used to determine the sequence evolution model that best fits our data using the Akaike Information Criterion (AIC). This program examined the fit of 56 models, with either a proportion of invariable sites (I), a gamma distribution of among-sites variation of substitution rates (Γ), or both (I + Γ). The best-fitting substitution models were TrN93 + Γ + I [81] for the GHR data set, and GTR + Γ + I [82] for the CYB data set. However, to run the same model of sequence evolution under PHYML and MrBayes, GTR was chosen for all phylogenetic analyses (the optimal GHR topology was not model-dependent, as it appeared that both GTR and TrN93 topologies were identical). To allow a fair comparison of α estimates of Γ-shape among genes and partitions, we did not use a proportion of invariable sites, but rather assigned eight discrete Γ categories (Γ₈). The Γ distribution allows some sites to evolve at a very low rate, and the incorporation of a fraction of invariable sites does not necessarily lead to a significant increase in likelihood [83].

To avoid excessive calculation times, our PAUP* ML analyses were conducted in two steps. First, we estimated ML parameters on a neighbor-joining (NJ) starting tree. Second, a ML heuristic search was conducted by Tree Bisection Reconnection (TBR) branch swapping to identify the optimal tree under these constrained GTR + Γ₈ parameter estimates. This tree was re-used for a new round of parameter estimation/branch swapping, and the procedure was repeated until there was a stabilization of both topologies and parameters. The stability of nodes was estimated by ML bootstrap percentages (BP_ML) [84], computed by PHYML after 100 replicates of site resampling, with BioNJ starting trees. Because of its rapidity, PHYML was preferred over PAUP* for bootstrap analyses. To assess the amount of phylogenetic signal contained within an individual partition (each codon position of each gene), 500 replications of ML bootstrapping were also independently performed for each partition.

Bayesian analyses were performed with one distinct GTR + Γ₈ model per gene and codon position, with unlinking base frequencies, GTR, and Γ parameters. Metropolis-Coupled Markov chain Monte Carlo (MCMCMC) sampling was conducted during 10,000,000 generations with five incrementally heated chains. We used Dirichlet priors for base frequencies (1,1,1,1) and for GTR parameters (1,1,1,1,1) scaled to the G-T rate, a Uniform(0.05,50.00) prior for the Γ shape, and an Exponential(10.0) prior for branch lengths. Bayesian posterior probabilities (PP) were computed from trees sampled every 100 generations, after removing the 50,000 first trees as the "burn-in" stage. In order to discriminate between moderately and strongly supported nodes – for which initial PP were superior to 0.95 – we also calculated bootstrapped Bayesian posterior probabilities (BP_bay) as suggested by [85] and [86]. Due to computing time limitations, BP_bay were only computed for the combined data set (GHR + CYB). First, 100 bootstrap pseudo-replicates were independently generated from each of the six partitions (the three CYB and the three GHR codon positions) using the SEQBOOT program 3.6a2.1 [87] of the PHYLIP package. Second, for each of the 100 concatenated bootstrap data sets, MCMCMC sampling of trees was performed as previously described for the original data under the six GTR + Γ₈ partitioned model, except that trees were sampled every 100 generations for only 500,000 generations. To maximize the probability that the chains reached stationarity in each bootstrap replicate, one-half of the 5,000 trees sampled from the posterior probability distribution was systematically removed as the burnin [86]. BP_bay resulted from the overall 50% majority rule consensus of the 500,000 saved trees.

Following [34], sequences were also analyzed at the amino-acid level. We respectively used the JTT + Γ + I and mtREV + Γ + I ML models of protein evolution for GHR and CYB sequences as implemented in PHYML.

TG and EJPD initiated the study, assembled the data, designed and ran the calculations. TG contributed to collect specimens in the field, obtained DNA sequences, and wrote the manuscript under the supervision of EJPD. MKT intensively helped TG in DNA sequencing and sequence alignment. EP initiated the project on arvicoline phylogenetics and obtained financial support. He contributed with TG to collect specimens and assisted with drafting. PC and SS provided CYB sequences and helped to improve the manuscript. All authors read and approved the final manuscript.

Table 4 Genera and species of Arvicolinae, with cytochrome b (CYB) and Growth Hormon Receptor (GHR) accession numbers and references.

The authors wish to thank Francois Catzeflis who provided most of tissue samples for this study, Pierre-André Crochet, David Duneau and Johan R. Michaux who also contributed to increase our taxa set, Patrick Giraudoux, Jean-Pierre Quéré, and the European "Human cystic and alveolar echinococcosis in northwest China" STD TS3-CT94-0270 and National Institute of Health and National Science Fondation (USA) "Parasitic zoonosis transmission in China"IR01TW01565-01 programs for some precious samples of Chinese voles. TG acknowledges the financial support of a MENRT Doctoral grant. This work has also been supported by the "Institut Français de la Biodiversité", the "ACI Informatique-Mathématique-Physique en Biologie Moléculaire [ACI IMP-Bio]", the IFR119 "Biodiversité Continentale Méditerranéenne et Tropicale" (Montpellier), and benefited from discussions around the European Marie Curie project HOTSPOTS "Understanding and Conserving Earth Biodiversity" (contract MEST-CT-2005-020561). The authors wish to thank Andrew Rodrigues for his help, and two anonymous reviewers whose comments and suggestions greatly improved the paper. This publication is the contribution N° 2006-057 of the Institut des Sciences de l'Evolution de Montpellier (UMR 5554 – CNRS).