|Home | About | Journals | Submit | Contact Us | Français|
Multiple unlinked genetic loci often provide a more comprehensive picture of evolutionary history than any single gene can, but analyzing multigene data presents particular challenges. Differing rates and patterns of nucleotide substitution, combined with the limited information available in any data set, can make it difficult to specify a model of evolution. In addition, conflict among loci can be the result of real differences in evolutionary process or of stochastic variance and errors in reconstruction. We used 6 presumably unlinked nuclear loci to investigate relationships within the mammalian family Tupaiidae (Scandentia), containing all but one of the extant tupaiid genera. We used a phylogenetic mixture model to analyze the concatenated data and compared this with results using partitioned models. We found that more complex models were not necessarily preferred under tests using Bayes factors and that model complexity affected both tree length and parameter variance. We also compared the results of single-gene and multigene analyses and used splits networks to analyze the source and degree of conflict among genes. Networks can show specific relationships that are inconsistent with each other; these conflicting and minority relationships, which are implicitly ignored or collapsed by traditional consensus methods, can be useful in identifying the underlying causes of topological uncertainty. In our data, conflict is concentrated around particular relationships, not widespread throughout the tree. This pattern is further clarified by considering conflict surrounding the root separately from conflict within the ingroup. Uncertainty in rooting may be because of the apparent evolutionary distance separating these genera and our outgroup, the tupaiid genus Dendrogale. Unlike a previous mitochondrial study, these nuclear data strongly suggest that the genus Tupaia is not monophyletic with respect to the monotypic Urogale, even when uncertainty about rooting is taken into account. These data concur with mitochondrial DNA on other relationships, including the close affinity of Tupaia tana with the enigmatic Tupaia splendidula and of Tupaia belangeri with Tupaia glis. We also discuss the taxonomic and biogeographic implications of these results.
As multilocus genetic data proliferate in systematic studies, the question of how best to analyze them becomes more important than ever. Early analyses almost always involved simply concatenating data under a single model of evolution or optimality criterion. More recently, it has become easy to partition data—by locus, codon position, or any other criterion that appeals to the investigator's knowledge of the data—and assign a different model to each partition. However, it can be difficult to know what the best set of partitions is, and increasing the number of partitions means that each contains fewer data from which to estimate parameters. An alternative is mixture modeling, in which data are not partitioned but in which multiple models are averaged across all sites (Gelman et al. 2004; Pagel and Meade 2004). Other nonconcatenation methods treat multiple gene trees as having possibly separate evolutionary histories but use the information in each gene tree to infer a species tree through consensus methods (e.g., Holland et al. 2006; Huson and Bryant 2006) or using the properties of the coalescent (e.g., Degnan and Salter 2005; Edwards et al. 2007). Because these methods make different uses of the same data, require different assumptions, and have different strengths and weaknesses, using them together can allow complementary interpretations of a given data set. Within a general framework of Bayesian phylogenetics, we use 3 methods to examine phylogenetic hypotheses in a poorly studied mammalian group: mixture modeling, partitioned analyses, and splits-network consensus analysis. We use this combination of methods to help us determine where individual methods may be suggesting a solution with inflated confidence and to identify potential sources of conflict in our data.
The mammalian order Scandentia (treeshrews) is notable both for its close affinity to primates and for the lack of recent attention that has been paid to its evolutionary history. Treeshrews have long been considered to be among the closest living relatives of primates and have been well represented in recent large-scale studies of mammalian interordinal relationships (e.g., Murphy et al. 2001; Scally et al. 2001; Reyes et al. 2004; Janecka et al. 2007). As a result, some species have become important model organisms in biomedical research (e.g., Czeh et al. 2001; Bahr et al. 2003; von Weizsacker et al. 2004). However, knowledge of the evolutionary history of this group is woefully incomplete, and the last thorough revision of treeshrew taxonomy was published nearly a century ago (Lyon 1913). In many ways, treeshrews epitomize the legacy left by Victorian era typological taxonomy—of the more than 120 species and subspecies described between 1820 and 1920, only 20 species are currently recognized, although taxonomists readily acknowledge this as an underestimate (Corbet and Hill 1992; Helgen 2005). The order is generally Southeast Asian (Fig. 1), with a distribution almost perfectly delimited by Wallace's line (Wallace 1860, 1876) and extending west into India (Helgen 2005). This region encompasses several major conservation hot spots (Myers et al. 2000; Sodhi et al. 2004) with high biodiversity and endemism, and understanding evolutionary and biogeographic processes in the area is critical. The broad distribution of treeshrews makes them excellent candidates for regional-scale research in Southeast Asia, but questions about the biogeography, ecology, ecogeography, behavior, and morphological evolution of this group cannot be answered in a phylogenetic context without some knowledge of systematic relationships within the order.
One contentious question within treeshrews is the status of Urogale, a monotypic genus endemic to the Mindanao Faunal Region of the Philippines (Heaney et al. 1998). Urogale everetti was originally described as a species of the more widely distributed and speciose genus Tupaia (Thomas 1892) but was later elevated (Mearns 1905) on the basis of its unique morphology; modern evidence for its generic status is weak at best. DNA–DNA hybridization data (Han et al. 2000) suggest that Urogale is nested within Tupaia, as do some morphological and morphometric analyses (Steele 1973; Olson et al. 2004; Sargis 2004). Other analyses have supported the distinction between Tupaia and Urogale (Butler 1980; Luckett 1980), although results have been shown to vary depending on analytical methods and assumptions in character coding (Olson et al. 2004). Recently, Olson et al. (2005) used mitochondrial DNA to infer phylogenetic relationships within Scandentia but found little resolution for most intergeneric relationships in the family Tupaiidae, which includes Tupaia, Dendrogale, Anathana, and Urogale. Urogale has the easternmost distribution of any treeshrew, and Anathana has the westernmost distribution (Fig. 1); understanding the geographic basis of diversification in the family depends on the relationships among these genera. Dendrogale, from mainland Southeast Asia and Borneo, is the most basally divergent tupaiid based on 12S data (Olson et al. 2005). Ptilocercidae, the only other treeshrew family, contains the monotypic genus Ptilocercus (Helgen 2005) and has been recovered as a distant sister taxon to Tupaiidae (Olson et al. 2005; Janecka et al. 2007). Here, we use additional DNA sequences from 6 nuclear genes to investigate phylogenetic relationships in Tupaiidae and to assess the genetic evidence regarding the monophyly of Tupaia with respect to Urogale. Multiple genes provide a more complete view of evolutionary history than any single locus can, but multilocus data sets demand particular attention to phylogenetic models and conflicting signals. We use a combination of Bayesian systematics and splits-network methods to analyze our data. We also use a close examination of tree topology to investigate support for a variety of roots in these data and discuss the implications of conflict among genes.
We included all treeshrew species for which fresh DNA or tissue was available (see Appendix Table A1 for voucher information). This includes U. everetti, 9 species of Tupaia, and Dendrogale murina, which we used as the outgroup for all analyses (Olson et al. 2005); we lack Ptilocercus lowii, Dendrogale melanura, Anathana elliotti, and 6 species of Tupaia. We extracted DNA using the PureGene protocol for animal tissue (Gentra Systems, Valencia, CA). We sequenced DNA from 6 nuclear genes: brain-derived neurotrophic factor (BDNF, 566 bp; coding), the 3′ untranslated region (UTR) of cAMP responsive element modulator (CREM, 469 bp; noncoding), the 3′ UTR of phospholipase C beta 4 (PLCB4, 334 bp; noncoding), tyrosinase (Tyr, 440 bp; coding), recombination activating gene 2 (RAG2, 450 bp; coding), and Exon 28 of the von Willebrand factor (vWF, 574 bp; coding), for a total of 2833 bp. Although we have no way to be sure where in any given treeshrew genome these genes fall, they do not appear to be adjacent on the current shotgun assembly of the Tupaia belangeri genome (National Center for Biotechnology Information locus AAPY01000000), and we therefore have good reason to believe that none of them are closely linked.
Amplification followed standard polymerase chain reaction (PCR) methods, in general 35 cycles of denaturation, annealing, and extension, with a magnesium concentration of 1.5–4.0 mM, an annealing temperature of 50–65°C, and the primers shown in Table 1. PCR products were purified with either GeneClean (Bio101, Vista, CA), following the manufacturer's instructions, or Exo/SAP (USB Corp., Cleveland, OH) using one- quarter of the recommended volume of each reagent. We sequenced PCR products using 0.5–1.0 μL of ABI BigDyes v3.1 dye termination chemistry and an ABI 3130xl or 3100 sequencer in the Core Sequencing Facility of the University of Alaska Fairbanks Institute of Arctic Biology (Applied Biosystems, Foster City, CA). We used Sequencher 4.2–4.7 (GeneCodes Corp., Ann Arbor, MI) to inspect sequence chromatograms and initially align sequences; alignments were then adjusted by hand. For protein-coding genes, gaps were placed to eliminate frameshifts; for all genes, single large gaps were preferred over multiple small gaps. All sequences have been deposited in GenBank under accession numbers FJ554887–FJ555018, and alignments are available from TreeBASE (study accession number S2248; matrix accession number M4268) and from DRYAD.
We used BayesPhylogenies (Pagel and Meade 2004) to implement a mixture model. Particularly in data sets containing multiple genes or regions with substantial rate or pattern heterogeneity, a single likelihood model may be inappropriate. Indeed, recent studies have found that for large data sets, underpartitioning can result in misleading posterior probabilities (Brown and Lemmon 2007). Advances in computer power and Bayesian systematic methods have made it easy to partition multigene data sets, allowing genes with different underlying histories, or parts of genes (such as first, second, and third codon positions), different models of evolution, and different overall rates of evolution. It would be possible to model our data set with 6 partitions (1 per gene), 4 partitions (for noncoding regions and codon positions), or any other conceivable arrangement, and to assign a different model to each partition. However, the amount of data in any given partition may be insufficient to resolve a large number of parameters, making it easy to overparameterize systematic models (e.g., Rannala 2002; Sullivan and Joyce 2005). From a purely practical standpoint, using a model with many partitions and parameters can increase the number of Markov chain Monte Carlo (MCMC) generations required to converge on a stationary distribution and to obtain an adequate posterior sample for each parameter. Speed and performance can often be improved by eliminating some partitions; genes or parts of genes may have similar enough characteristics that a single data partition is sufficient to model more than 1, but there is usually no a priori way to know what the best arrangement of partitions is. An alternative to partitioning is mixture modeling, in which partitions are not specified but in which multiple sets of model parameters are simultaneously estimated; this can allow unsuspected shared characteristics of the data to be modeled together (Pagel and Meade 2004). Mixture models allow multiple data patterns to have different model parameters, without the need for defining data partitions; the likelihood for each site is the sum across all models, weighted by their fit. This approach often reduces the number of parameters over a highly partitioned model, and it can illuminate patterns of variation or similarity in the data that do not correspond to genes, codons, or other a priori partitions.
We therefore tested possible phylogenetic mixture models in BayesPhylogenies, using 1–6, 10, 15, and 20 independent GTR rate matrices, or patterns, which we notate “Q” following Pagel and Meade (2004), with and without gamma rate variation (notated “Γ”) and separate base frequencies (“π”). Rate matrices, gamma variation, and base-frequency variation interact and can sometimes compensate for each other—adding a rate matrix, for example, can allow “fast” sites an overall faster rate without explicit gamma variation—and it is possible for quite different parameterizations to fit the data similarly. We ran each of these 34 models for 2 runs of 5 million generations with a burn-in of 100 000 generations. We then used the harmonic mean of the model log likelihoods (post-burn-in) to approximate the marginal likelihood of the model (Raftery 1996) and compared model likelihoods using Bayes factors (Kass and Raftery 1995; Raftery 1996) to choose a model that balanced model fit with parameterization. We used the chosen model in 2 independent runs, 1 of 51 million generations and 1 of 26 million generations, each with 3 heated chains, and eliminated the first million generations as burn-in. We used MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003) to compute consensus trees and bipartition frequencies from post-burn-in tree sets.
To compare mixture modeling with partitioned analyses, we also analyzed the combined data in MrBayes. We used 5 partitioning schemes: unpartitioned (1 total “partition”); partitioned into coding and noncoding regions (2 partitions); partitioned by codon position, plus 1 partition for noncoding regions (4 partitions); partitioned by gene (6 partitions); and partitioned by gene and codon position (14 partitions). In all cases, for the best comparison with BayesPhylogenies, we assigned each partition a separate GTR rate matrix and separate base frequencies but included a single gamma parameter linked across all partitions. This makes, for example, the 6-partition model partitioned by gene as similar as possible to a 6Q + Γ + π model in BayesPhylogenies. In all partitioned analyses, we assigned the rate prior (“ratepr”) a Dirichlet(1,…,1) prior distribution, allowing relative rates to vary among partitions. We ran 2 runs of 20 million generations, sampling every 1000, using 4 chains with the default heating scheme.
For both MrBayes and BayesPhylogenies, we assessed convergence, mixing, and sampling using R 2.4.1 and 2.5.1 (R Development Core Team 2004). Our determination of convergence included comparing log-likelihood series for the 2 runs for each model, comparing means and distributions of other model parameters for the 2 runs for each model, and checking the variance of each parameter. Likelihoods and parameters were considered to have converged if the 2 runs appeared to be sampling the same stationary distribution. In 1 case (in BayesPhylogenies) in which 1 run had clearly gotten “stuck” in a different stationary distribution, with a poorer likelihood, we discarded and repeated a run. We also checked the effective sample size of parameters, which takes autocorrelation between successive samples into account and which can help indicate when 1 or more parameters are not mixing well (this can be particularly problematic with complex, highly parameterized models). When comparing the variance of models, we rescaled BayesPhylogenies rate matrix samples as proportions of the rate sum so that variances are roughly on the same scale as those from MrBayes. Calculations and manipulations in R used our own scripts or functions in the packages MCMCpack (Martin and Quinn 2007) and coda (Plummer et al. 2006).
Phylogenies based on different loci can vary both because of stochastic differences (the same evolutionary history can result in different trees) or because of evolutionary differences (different loci may have different histories). In either case, a multigene data set can contain evolutionary information that is not visible in a standard consensus tree, especially one limited to a strictly bifurcating branching pattern. Similarly, conflicting relationships can be present in trees from a single locus if more than 1 phylogenetic pattern is supported by the data. The strength of minority signals can be important in evaluating support for phylogenetic hypotheses. In order to investigate these minority signals and check the potential effect of individual-gene histories and unique trees, we constructed trees for each individual gene and computed bipartition frequencies in MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003; Altekar et al. 2004). We chose models for each gene using the Akaike Information Criterion by scoring each model on a neighbor-joining tree; the models chosen were HKY + I (BDNF, Tyr, and RAG2), HKY + G (CREM and vWF), and GTR (PLCB4). All runs used the default 4-chain heating scheme for 20 million generations with a 1 million generation burn-in. Convergence and sampling were assessed in R using the same methods as for the combined analyses.
Network methods can help identify conflict among multiple genes or among signals in multiple trees when support is ambiguous (e.g., Holland et al. 2004; Gadagkar et al. 2005; Huson and Bryant 2006; Gauthier and Lapointe 2007). Bayesian methods are well suited for the investigation of topological conflict among many trees, as they are fundamentally based on finding sets of trees rather than a single optimal tree. We used R (R Development Core Team 2004) to calculate the mean frequency of all bipartitions (splits) present in the 6 tree sets, based on the summary of bipartitions calculated by MrBayes. We then used SplitsTree (Huson and Bryant 2006) to create a network of all bipartitions with a mean frequency greater than 10%, some of which conflict with each other. This allowed us to see, simultaneously, bipartitions with low to medium support from several genes or high support from one, not just the partitions supported simultaneously by all the genes. This is a consensus method; we are making a network directly from taxon bipartitions—essentially from trees—and not from the sequence data. Our approach is slightly different than standard consensus network methods (e.g., Holland et al. 2004), as we select bipartitions by their mean frequency in a Markov chain–derived tree set rather than combining those that were present in a set of single trees found under some optimality criterion. It is equivalent to constructing a splits network directly from the Bayesian tree set, as suggested by Holland et al. (2006), but with a large number of trees, we found network construction to be substantially faster using our approach. In order to separate conflict involving only in-group taxa from that involving root placement, we also excluded the outgroup, D. murina, and created a second network. Both approaches enabled us to identify sources of conflict and places where individual genes “disagree” with each other and with the combined data topology.
In order to calculate the posterior probability of specific taxon bipartitions and tripartitions, including root placements, we used PAUP*4.0b10 (Swofford 2002) to define a constraint tree based on each partition and to filter posterior tree sets against these constraint trees, generating a list of partitions and their frequencies. To determine the frequency of bipartitions without the root, we pruned the outgroup from sets of trees before defining constraints or filtering. We then summarized these partition frequencies in R. The NEXUS network files are available from DRYAD.
We found that more complex, parameter-rich models did not always outperform less complex models according to Bayes factor model comparison. The single-pattern model with no rate variation had by far the worst marginal log likelihood, but most other models had similar scores (Fig. 2a). For each number of patterns (Q), the model with gamma rate variation (Γ) and independent base frequencies (π) had the best marginal likelihood. The maximum marginal log likelihood belonged to the 5Q + Γ + π model, but this model was only slightly better than several other models. The log-Bayes factor of this model compared with each of the 6Q + Γ + π, 4Q + Γ + π, 3Q + Γ + π, and 2Q + Γ + π models was between 0 and 1, implying that the support for the best model over these was “barely worth mentioning” (Raftery 1996, p. 165). Our results suggest that although 1 rate matrix, even with among-site rate variation, is insufficient for these data, adding parameters beyond an intermediate model is unnecessary and sometimes detrimental. With many more patterns, model likelihoods began to decline, especially for non–Γ+π models. These models also had the worst sampling and highest levels of autocorrelation among parameters; the effective sample size of likelihoods in the 20Q model was only 107 in the total sample of 9800 (compared with more than 1000 for most models). Overall, more highly parameterized models mixed more poorly, ran more slowly, and resulted in higher parameter variances (Fig. 2b). For our final tree search, we therefore chose the 2Q + Γ + π model, the simplest of the models with nearly identical support to the best model, giving us a compromise between model fit and number of parameters. The average variance in rate parameters of this model was also somewhat less than that for models with slightly higher marginal likelihoods (Fig. 2b). The final effective sample size was at least 100 per run and 450 total for every parameter. The 2 patterns had mean weights of 0.54 (95% highest posterior density interval [HPDI] 0.26–0.87) and 0.46 (95% HPDI 0.13–0.74), indicating that both are contributing substantially to the final weighted likelihood score.
In MrBayes, marginal likelihoods increased sharply with the number of parameters (Fig. 2a), as did parameter variance (Fig. 2b), and no leveling off was apparent in either. Model likelihoods were generally higher than those for mixture models in BayesPhylogenies.
The 50% consensus topology was identical in BayesPhylogenies and all MrBayes analyses. The combined data strongly support a single topology for the genera Tupaia and Urogale (Fig. 3). Of the taxa included here, sister relationships between Tupaia glis and T. belangeri and between Tupaia tana and Tupaia splendidula were strongly supported in all analyses. Likewise, the grouping of Tupaia longipes with T. glis/T. belangeri and Tupaia minor with T. splendidula/T. tana was supported. Relationships within the ingroup were generally well supported, with posterior probabilities greater than or equal to 95% for all but 2 branches shown in Figure 3. One exception is the support for the monophyly of T. belangeri, which has a posterior probability of 91%; the deepest split within T. belangeri is nearly as deep as the divergence between T. belangeri and T. glis. The other is a 66% posterior probability within T. tana; other intraspecific relationships in T. tana and T. splendidula had support less than 50% and are shown as multifurcations in Figure 3. The 90% posterior probability for the node subtending Urogale through T. tana reflects the location of the outgroup rather than uncertainty about relationships among ingroup taxa.
The combined nuclear genes also support a single root position with moderate posterior probability (89.9%), separating the T. belangeri/T. glis/T. longipes clade from the other included species. There were 13 other root positions in the posterior tree set (Fig. 3). Only 1 other root, on the internal branch separating U. everetti and Tupaia gracilis, was found in more than 5% of trees; 3 roots had probabilities between 1% and 5%, and the rest were <1%. The posterior probability of rooting on the branch to Urogale, making Tupaia monophyletic, was only 1.7%. Thus, the total posterior probability for Urogale nested within Tupaia, with the root at any position, is 98.3%.
With the most comparable model (using gamma variation and separate base frequency matrices), BayesPhylogenies inferred consistently longer total tree lengths than MrBayes (Fig. 4). Trees in MrBayes got slightly longer as the number of partitions increased, then much longer under our final partitioning scheme (with 14 separate partitions). Total tree lengths in the other 4 MrBayes analyses were very consistent but always slightly longer than comparable Γ + π models in BayesPhylogenies. There is no clear trend in total tree length in BayesPhylogenies as the number of patterns increases and no big increase with highly parameterized models.
Individual genes showed more topological variation and, overall, less support for most groupings than the combined data. In several cases, individual-gene consensus trees have strong support for 1 or a few relationships (Fig. 3 and Table 2) but little or no support for those elsewhere in the tree, which collapse into multifurcations. Some such well-supported relationships are present in the resolved combined data tree, but others are not. Conversely, some relationships in the combined data tree have high support from all 6 genes (e.g., the monophyly of our 2 Tupaia dorsalis samples), but some are supported by only 1 or 2 genes (e.g., the monophyly of our 5 individuals of T. belangeri). The single-gene trees did have substantial variability in root position—more than in the tree set from the combined data—and no single gene had strong support for any root (Table 3). The maximum root probabilities for individual genes ranged from 7.2% (BDNF) to 32.2% (RAG2), and 24 different root positions had posterior probabilities of at least 5% in at least 1 individual gene tree. The majority root in the combined analysis was not the most common root in any individual gene tree; its posterior probability ranged from 0.0% to 21.1% (Table 3). However, as with the combined analysis, support for a root on U. everetti was not found in any of the individual gene trees. Individual gene posterior probabilities for this root ranged from 0.6%(vWF) to 8.3% (Tyr) (mean: 5.1%; median: 5.7%).
In order to visualize the areas in which the 6 genes agreed and disagreed about topology, we created a splits network of all the bipartitions (splits) with a mean frequency of at least 10% (Fig. 5). In this type of network, any given split is represented by a set of 1 or more parallel lines; conflicting splits are perpendicular to each other. The length of a split is not a measure of evolutionary distance, as is usual for branch lengths in trees, but the amount of support for that split, in this case the mean posterior probability in the 6 tree sets. The number of relationships in conflict in a given area of the network is thus the number of perpendicular splits. A 2-dimensional rectangle suggests 2 conflicting relationships, whereas a 3-dimensional rectangular prism suggests 3; more complex structures are formed from more relationships. The shape of the rectangle or prism indicates relative support for conflicting relationships: the more even the edges, the more even the support.
The overall topology of our network (Fig. 5a) is similar to the consensus combined data tree. However, the network shows the conflicting signals present within and among data sets, which the consensus tree does not. There are several regions of conflict in these data. First, vWF and PLCB4 have moderate support for a clade containing T. tana, T. splendidula, and Tupaia palawanensis, conflicting with the clade uniting T. tana, T. splendidula, and T. minor that was supported by the other 4 genes and the combined analysis. In the network, the presence of conflict between these 2 topologies, and the relative support for each, is shown by the set of rectangular prisms connecting T. minor and T. palawanensis. This is even clearer when the outgroup is removed from the network (Fig. 5b); without the additional uncertainty about the root, these conflicting splits can be shown with a single rectangle.
Second, the complexity of the network around Urogale, T. gracilis, and T. dorsalis shows that many different topologies are in conflict with each other in this part of the tree. Again, this is clearest when Dendrogale is excluded. Without the outgroup, several conflicting relationships are present. There is some support for bipartitions of each of these 3 taxa with T. belangeri, T. glis, and T. longipes, and for either T. dorsalis or T. gracilis with T. tana, T. splendidula, T. minor, and T. palawanensis. There is also some support for a bipartition of T. gracilis and U. everetti. Because of support from 1 gene (CREM), there is also a visible partition of T. dorsalis with 2 individuals of T. belangeri. When the outgroup is included, the number of relationships present increases sharply, reflecting the lukewarm support these individual gene trees have for any root position.
The other clear areas of conflict involve relationships within species (T. tana, T. splendidula, T. belangeri) and within the T. glis/T. belangeri clade. Each of these species was recovered as monophyletic in the combined analysis, and this is not contradicted by individual genes for T. tana and T. splendidula. However, vWF groups 1 individual of T. belangeri with T. glis, and only 2 genes of these 6 strongly supported the monophyly of T. belangeri (Table 2).
Mixture and partitioned models offer fundamentally different ways of approaching multilocus data and of handling heterogeneous evolutionary processes. We found that model likelihoods were generally higher for partitioned models. This may reflect several of the ways in which these classes of model differ. First, in a mixture model, the likelihood at each site is determined by a weighted sum of likelihoods under each pattern. This will never be higher than the likelihood under the best-fitting pattern and may be substantially lower if some patterns fit that site poorly but have a relatively high weight. Second, the pattern weights are the same for all sites. This means that a pattern that fits a few sites well but the rest very poorly will have a relatively low weight and will not contribute much to the total likelihood, even at the sites that it fits, and those sites may have a poor overall likelihood.
It is misleading, although tempting, to think of mixture models as dynamically assigning sites to partitions. Sites under a mixture model do not belong to 1 partition; their likelihoods are averaged across all patterns. It is important to note that all sites sharing a single unique pattern in the data (e.g. “AAT” or “CCC” for 3 taxa) will have the same likelihood under a mixture model; in a partitioned model, the likelihood of a given site depends on which partition it is in. Thus, in a mixture model, all sites with a given site pattern will have their highest likelihood under the same pattern. This makes it very unlikely for sets of sites matching genes, codon positions, or other common a priori partitions to have their highest likelihood under 1 pattern, unless their evolution has been so heterogeneous that even site patterns are not shared across those partitions. In our data set, the dominant site patterns are the invariant ones—all A, all C, all G, and all T. These sites, unsurprisingly, are distributed throughout all the typical a priori partitions; it would be unusual for an entire gene or codon partition, for example, to have no invariant sites. As a result, sites that have their highest likelihoods under Patterns 1 and 2 in our analysis are not grouped in any clear way. Most of the 267 site patterns present here are found across the various likely a priori partitions, making it impossible for the identities of sites with their highest likelihood under a single pattern to conform to those partitions. It is also worth noting that in a mixture model, all the data are being used to estimate all the parameters. This is in stark contrast to partitioned models, in which different portions of the data, and different amounts of data, are used for each set of parameters. This approach may increase the variance of some parameter estimates but also keeps mixture modeling from suffering from one of the problems of partitioned models—partitions without sufficient variable data may be unable to estimate even a few parameters with precision.
We used Bayes factors to compare models. Bayes factors are not a hypothesis test and do not result in a P value for rejection of 1 model or another, but they can provide guidance in choosing among complex and nonnested models. Unlike many other popular model choice criteria, Bayes factors, which use an estimate of the total (marginal) likelihood of a model rather than the likelihood at a single combination of parameters, include information from across a wide range of parameter values. Because the marginal likelihood depends on both the data and the prior distribution for each parameter, Bayes factor tests include an inherent penalty for overparameterization. Each additional parameter adds a prior, and each additional prior (which is, at any given point, a number much smaller than 1) tends to reduce the total probability, possibly offsetting the gain in likelihood from better model fit. This is also evident in the behavior of MCMC samples from more complex models. As the number of parameters in a model goes up, the variance of samples along an MCMC run often increases, and the marginal likelihood can decrease—even for nested models, in which the maximum likelihood will always be better for more complex models. It is important to note that because this inherent penalty is related to variance, “uninformative” or diffuse priors provide more of a penalty than more focused ones. Choice of prior, not just the fit of the models to the data, can therefore affect the results of a Bayes factor comparison.
Phylogenetic use of Bayes factors has in general focused on the harmonic mean as an estimator of model likelihood, as we do here. The harmonic mean estimator is known to be unstable and to often have an infinite variance (e.g. Raftery 1996). This is manifested by occasional replicate runs having very different model likelihoods, even though by all typical convergence diagnostics they are sampling the same stationary distribution—the strong effect of rare low values on harmonic means makes the estimate sharply different. It is perfectly acceptable, and even desirable, for an MCMC chain to find an occasional low likelihood value, as chains that are mixing well should encounter such values frequently and accept them occasionally. Interestingly, we found that the difference in model likelihood between multiple runs was indeed related to how well the chain mixed. Poorly mixing chains, with high autocorrelation between successive samples, tend to result in nearly identical model likelihoods for multiple runs. Chains that mixed better had a greater chance of encountering an occasional low likelihood, which could have a big effect on the chain's marginal likelihood. However, this may cause the marginal likelihood of a poorly mixing model to be overestimated (if it is finding low values more rarely than it should be) and that of a well-mixing model to be underestimated (if it happens to find low values more frequently than their true probability). In a Bayes factor test, this could effectively penalize a model for mixing well, which would be undesirable.
We found that parameter variance increased with model complexity for both partitioned and mixture models. In the former, the mean variance of rate matrix parameters increased sharply as the number of partitions increased and did not level off. In the mixture models, the mean variance rapidly leveled off and appeared to be converging on the variance generated by randomly choosing values from the prior distribution. In both cases, it seems that adding parameters, while increasing model fit, does tend to increase variance as well. However, a Bayes factor test of mixture models suggests that these 2 attributes are best balanced by an intermediate model, whereas a similar test of the partitioned models we used would choose the most complex model. Interestingly, the most complex partitioned model also resulted in extremely long trees and a broad distribution of tree length (Fig. 4), despite its much higher harmonic mean likelihood. Highly partitioned models are known to cause problems estimating tree length under some conditions (e.g. Marshall et al. 2006), and in our case, the jump in tree length may reflect inadequate adjustment of branch length priors or (as with the other high parameter variances) the increasing effect of the priors on posterior branch length estimates.
For the taxon sample included here, the question of the monophyly of Tupaia with respect to Urogale is solely one of rooting: if the root falls on Urogale, Tupaia is reconstructed as monophyletic; if the root falls anywhere else, it is not. We found little evidence for the monophyly of Tupaia with respect to Urogale either in the combined analysis or in individual gene trees. The combined analysis strongly supports Urogale as nested within Tupaia, with posterior probability of 89.9% for the majority root position (Fig. 3) and 98.3% for the paraphyly of Tupaia; the posterior probability of a sister relationship between Tupaia and Urogale is only 1.7%. Support for this sister relationship varies among individual genes; this topology is occasionally recovered in some trees but was never the highest probability root and never had a posterior probability greater than 8.3%. Thus, these nuclear data reject the monophyly of Tupaia with respect to Urogale.
We are by no means the first to question the generic status of Urogale. Indeed, U. everetti was initially described as Tupaia everetti (Thomas 1892). Urogale was distinguished from Tupaia due to its unique morphology by E. A. Mearns in 1905, when he described Urogale cylindrura (which was later synonymized with U. everetti). Lyon (1913) upheld the status of Urogale in his seminal monograph on the order Scandentia, which remains the most thorough work on the taxonomy of the order. There is no question that Urogale is morphologically distinct from all treeshrews in the genus Tupaia, but autapomorphic morphological features cannot provide any evidence for its phylogenetic position. Han et al. (2000) concluded that Urogale was nested within Tupaia using DNA–DNA hybridization and some morphometric data, as did Sargis (2004) with postcranial evidence and Olson et al. (2004) in their reanalysis of Steele's (1973) morphological data. In a previous attempt with DNA sequence data, Olson et al. (2005) found that the mitochondrial gene encoding 12S ribosomal RNA (rRNA) did not resolve the position of Urogale. Their optimal maximum likelihood tree showed a polytomy of Urogale, T. gracilis, and the remainder of Tupaia, with poor support for most basal nodes in the Urogale/Tupaia radiation, leaving the monophyly of Tupaia an open question. However, we consider it significant that no phylogenetic study with good taxon sampling has recovered Tupaia as a monophyletic genus with strong support; neither genetic nor morphological data, analyzed with modern phylogenetic methods, supports the monophyly of the genus. We provisionally recommend synonymizing Urogale with Tupaia, although we acknowledge that broader taxon and character sampling could change our understanding of the history of this group.
Urogale is geographically as well as morphologically unique among treeshrews (Fig. 1). Its distribution is the farthest east in the group, and it is the only treeshrew species found on oceanic islands in the Philippines. It is thus the only treeshrew found east of Huxley's Line (Huxley 1868), which separates the oceanic region of the Philippines from the landbridge islands of the Sunda Shelf, including the Palawan island group. Urogale is endemic to the Mindanao Faunal Region of the southern Philippines, where it has been recorded only from the islands of Mindanao, Dinagat, and Siargao—that is, Mindanao and 2 small islands very close to Mindanao (Heaney et al. 1998; Helgen 2005). Several Philippine mammals have similar distributions, although many Mindanao Faunal Region endemics are somewhat more widespread (Heaney et al. 1998). Urogale‘s limited distribution suggests that it disperses poorly, that it has been extirpated in part of a previously more widespread range, or that it reached the Philippines only recently and has not yet had time to expand its distribution. On the other hand, its presence in the Philippines presumably reflects overwater dispersal, as these islands have never been connected to the Sunda Shelf (including Borneo) or to mainland Southeast Asia (Heaney 1985, 1986). The only other treeshrew whose current distribution can only be explained by overwater dispersal is Tupaia nicobarica, which is endemic to the oceanic, isolated Nicobar Islands (Fig. 1; Lyon 1913; Helgen 2005). These exceptions are notable given the apparent inability of treeshrews to cross seemingly insignificant channels (e.g., Lombok Strait, which separates Bali from Lombok Island and is only approximately 30 km at its narrowest point; Kitchener et al. 1990). The current position of Urogale in our phylogeny suggests that this treeshrew may have colonized the Philippines via overwater dispersal early in the Tupaia radiation. Our analysis of alternate rooting positions suggests that this is true not just for the majority root, but for the bulk of roots, as they are concentrated around T. gracilis and U. everetti on the backbone of the tree.
Our nuclear data do not support a close relationship between the 2 Philippine treeshrews in our sample, T. palawanensis and U. everetti. Although both of these are Philippine endemics, they are not sympatric and are endemic to regions that do not share close biogeographic affinities. Tupaia palawanensis occurs only in the Palawan Faunal Region, which is typically more similar to the Sunda Shelf than to other biogeographic areas in the Philippines (Dickerson 1928; Heaney 1986; Esselstyn et al. 2004). Our phylogenetic results suggest that each may be related to Bornean or Sundaland taxa and that treeshrews have colonized the Philippines at least twice, possibly through 2 different routes. However, the lack of some other Sundaic treeshrews in our data set (including the Bornean endemics Tupaia picta and Tupaia montana and the Indonesian Tupaia javanica) makes it impossible to draw strong conclusions about specific sister relationships. We also currently lack samples of Tupaia moellendorffi, the Calamian treeshrew, long considered a subspecies of T. palawanensis but recently restored to species status (Helgen 2005).
We continue to accumulate support for the close relationship of T. belangeri and T. glis, which together form a problematic species complex. Even the limited number of individuals of T. belangeri and T. glis included in this study demonstrate that current species taxonomy may not reflect monophyletic groups. Of the 5 individuals of T. belangeri we sequenced, 1 is deeply divergent from the rest and often groups with T. glis, as is clear in the network (Fig. 5). These data also confirm the association of T. longipes with the T. belangeri/T. glis species complex, an association further supported by a synapomorphic 6-bp deletion in vWF exon 28. Likewise, our nuclear data strongly support the close relationship between T. splendidula and T. tana. In addition to the phylogenetic evidence, this clade is supported by a 69-bp insertion/duplication event in the CREM 3′ UTR, found in all 4 T. tana and 3 T. splendidula in our data set but in no other taxon. Without complete taxon sampling for the family Tupaiidae (which contains several species for which, at present, no fresh tissue is available, and which are codistributed with these taxa), we cannot say for certain that these 2 are sister species, but there seems to be no question that they are each other's closest relatives in the taxon sample we currently have.
The key to understanding conflict in these data is 2-fold. First, it is useful to separate rooting questions from the topology of the ingroup. With Dendrogale, as in Figure 5a, the conflict appears to be spread throughout the tree and to be continuous along the backbone. The myriad root positions supported by 1 or more genes, in this case, overwhelm the conflict in the ingroup, which is actually much less. Separating the root from the in-group, as in Figure 5b, makes it clear that within the ingroup, conflict is concentrated in several separate regions. Indeed, when uncertainty about the root position is removed by pruning the outgroup, some posterior probabilities are considerably higher. The bipartition of T. belangeri/T. glis/T. longipes/U. everetti versus the rest of the ingroup, for example, has a posterior probability of 87.4% in the trees from the gene PLCB4 if the outgroup is removed, but is much lower when the outgroup is included because the position of the outgroup is so uncertain. The apparent low support for this relationship among ingroup taxa is really low support for the root.
When the root is removed, conflicting evidence for interspecific relationships (other than the T. belangeri/T. glis species complex) is limited to 2 areas, and the uncertainty of relationships among Urogale, T. gracilis, and T. dorsalis is independent of the lesser uncertainty involving T. palawanensis and T. minor. Indeed, without the confusion added by the outgroup, the latter is clearly a simple conflict between 2 opposing relationships, one with somewhat greater probability than the other. This is the benefit of using networks—in a standard consensus, this node would collapse to a polytomy and there would be no way of telling which possible relationships had some support. The more significant conflict in the tree involves U. everetti, T. gracilis, and T. dorsalis. In both networks, it is clear that several alternate relationships among these taxa have support from some genes, often at less than 50% posterior probability (and often with no relationships greater than 50%). The internodes separating these taxa appear to be very short, suggesting that successive divergences were rapid. In this situation, different genes may actually have different histories if polymorphisms present at the time of divergence sort differently. They may also appear to have different histories because of stochasticity in the mutation process; with only a few mutations defining a branch, it is harder to separate signal from noise than with many.
Second, it is helpful to pay attention to the difference between relationships that are supported by only a few genes and contradicted by others, and those that are supported by only a few genes but not contradicted. The monophyly of both T. splendidula and T. tana, for example, is strongly supported by 2 genes (RAG2 and vWF, and BDNF and vWF, respectively). In the network, both these are shown as relationships without conflict—no contradictory bipartitions exist at greater than 10% posterior probability. The clade combining T. splendidula, T. tana, and T. minor is also strongly supported by 2 genes (CREM and RAG2), but in this case, there is conflict: 2 other genes (Tyr and vWF) have strong support for T. splendidula+T. tana+T. palawanensis, and this is clear in the network. There is more support for the former relationship, but the network makes it possible to see and identify the latter one as well.
The 90% posterior probability we recovered for the majority root is hardly overwhelming evidence in Bayesian systematics, in which posterior probabilities for robust clades are often considerably higher. The evident difficulty in rooting the treeshrew phylogeny stems from an unfortunate combination of factors: the short branches separating nodes along the backbone of the tree, and the long branch separating the ingroup from our outgroup (D. murina). In the individual gene trees, some of the shorter internodes collapse, resulting in polytomies. This indicates that successive speciations may have been quite close together and is precisely the situation in which multiple genes are most likely to disagree with each other—the stochastic nature of both mutation and lineage sorting means that different genes may well support different evolutionary patterns. At the same time, no phylogenetic model is perfect, and this combination of a long outgroup branch and short internodes makes it particularly difficult for any phylogenetic method to reconstruct the root either accurately or consistently. We used a combination of phylogenetics and network methods to examine conflict surrounding both the root and ingroup relationships in more detail. Conflict in the ingroup was worst along the backbone of the tree, among T. gracilis, T. dorsalis, and U. everetti, and nearly every possible relationship among these taxa (including branching order and sister relationships among these taxa) was recovered with moderate posterior probability by at least 1 gene. These minority relationships are implicitly ignored by traditional consensus methods but can hold important information about sources of phylogenetic uncertainty. It is particularly interesting to note that among the ingroup taxa, a bipartition of T. gracilis+Urogale was found at a frequency greater than 15% in 4 of the 6 genes we used (Table 2); this is the ingroup relationship recovered by Olson et al. (2005) in their single maximum likelihood (ML) tree. This is, of course, not high posterior probability, but taken together, the evidence suggests that relationships among these taxa may be particularly sensitive to choice of gene, method, and parameters.
Examining conflict among individual genes allows us to see points that are not evident in the combined analysis. Likewise, examining minority signals allows us to get a better idea of all the relationships supported by the data. In this case, some high posterior probabilities in the combined analysis seem to be driven by single genes rather than by support across all the individual-gene data. This is worrisome because the sample of genes we have here is only a small proportion of the genome, and it is reasonable to assume that the distribution of support would change with a larger (or simply different) sample. We therefore have only limited confidence in relationships among T. gracilis, T. dorsalis, and U. everetti. On the other hand, relationships for which there is no substantial conflict in these genes, such as the relationship between T. tana and T. splendidula, are likely to be confirmed by further sampling.
This multigene nuclear data set confirms several of the results of Olson et al. (2005), including the close association between T. glis and T. belangeri, the association of these taxa with T. longipes, and the relationships among T. tana, T. splendidula, and T. minor. Indeed, with the exception of the root and the relationship between T. gracilis and U. everetti, our tree is identical to the single optimal ML tree reported by Olson et al. (2005) using the mitochondrial 12S rRNA gene. As we note above, the root and these relationships are precisely where there is the most uncertainty in the nuclear data. Comparing these trees suggests that 12S alone is in substantial agreement with the multigene nuclear data set and was probably doing an accurate job of reconstructing topology at this taxonomic level, although posterior probabilities in the mitochondrial tree are generally lower than those in the nuclear tree.
It is possible that additional taxa will, in the future, help solidify the position of the root in this group. In particular, data from Anathana, 1 of 2 genera missing from our analysis, and from D. melanura might help by breaking up the long branch between Tupaia/Urogale and Dendrogale. The species of Tupaia we were unable to include might also help, especially if their addition stabilized the unclear relationships among T. gracilis, T. dorsalis, and Urogale. In some cases, adding further outgroups, as well as taxon sampling in the ingroup, can help stabilize rooting. In treeshrews, however, the distance from the family Tupaiidae to the next available outgroups (Ptilocercus, Dermoptera, and Primates; Janecka et al., 2007) makes it nearly useless to add outgroups even farther from the ingroup; alignment uncertainty and the presence of multiple substitutions at segregating sites rapidly overwhelm the useful data that could be gained by this approach. Similarly, adding more data (nuclear or mitochondrial) is only marginally helpful in this situation. Faster-evolving genes result not only in longer internodes, but also in a longer branch to the outgroup; slower genes make the outgroup more tractable, but reduce short internodes even more, often to zero. Simultaneously resolving ingroup nodes and the root in this topology will probably require both types of character data, making the multigene approach even more necessary.
This research was supported by National Science Foundation Grant DEB-0542725/0542532 and an Alaska EPSCoR YIFA grant to L.E.O. and E.J.S., and startup funds from the University of Alaska Museum to L.E.O. DNA sequencing was conducted in the University of Alaska Fairbanks Institute of Arctic Biology Core Facility for Nucleic Acid Analysis with support from NSF EPSCoR Grant EPS-0346770. Analyses were performed on 2 parallel computing clusters administered by UAF Life Science Informatics, a core research resource supported by Grant RR016466 from the National Center for Research Resources, a component of the National Institutes of Health, with assistance from S. Houston. Fresh tissues of D. murina were the result of an expedition to Cambodia funded by a grant from the National Geographic Society's Committee on Research and Exploration to E.J.S. and L.E.O.
We are grateful to the following people and institutions for making tissues available for this study: W. T. Stanley and L. R. Heaney (Field Museum of Natural History); L. K. Gordon, J. F. Jacobs, and R. W. Thorington, Jr. (U.S. National Museum of Natural History, Smithsonian Institution); J. L. Patton, C. J. Conroy, and J. P. Boubli (Museum of Vertebrate Zoology, University of California, Berkeley); A. J. Gorog and P. Myers (University of Michigan Museum of Zoology); and J. South (University of Maryland). Fresh tissues provided by J. F. Jacobs resulted from fieldwork supported by the Smithsonian's National Museum of Natural History's Biodiversity Surveys and Inventory Program and facilitated by the Myanmar Nature and Wildlife Conservation Division. We thank J. Walston, E. Pollard, and M. Soryun of the Wildlife Conservation Society–Cambodia and Khiev Rithy Phion of the Cambodian Forestry Administration for their invaluable assistance with work in Cambodia. E. Humphries, H. Lanier, D. Sikes, N. Takebayashi, and 2 anonymous reviewers provided helpful comments, discussion, and technical assistance.
|Dendrogale murina||UAM 103000||Cambodia: Mondulkiri; Seima Biodiversity Conservation Area|
|Tupaia belangeri||FMNH 1654121||Zoo animal; no locality data available. ISIS M1194|
|USNM 5838573||Myanmar: Mon; 5 km northeast of Kinmun, at Yetagon Myaung in Kyaikhtiyo Wildlife Sanctuary|
|USNM 5837935||Myanmar: Bago; Dawe, in Bago Yoma mounts|
|MVZ 1864074, 1864082||Vietnam: Vinh Phu; Vinh Yen District, Tam Dao|
|Tupaia dorsalis||UMMZ 174651, 174427||Indonesia, Borneo: West Kalimantan; Ketapang Regency, Gunung Palung National Park, Cabang Panti Research Site|
|Tupaia glis||MVZ 192180||Indonesia, Sumatra: North Sumatra; 3 km northwest of Bukit Lawang|
|MVZ 192184||Indonesia, Sumatra: Aceh; Ketambe Research Station|
|Tupaia gracilis||USNZ 109023||Malaysia, Borneo: Sabah; Sepilok|
|Tupaia longipes||JS M02b||Malaysia, Borneo: Sabah; Danum Valley Field Center|
|Tupaia minor||USNZ 109751||Zoo animal; no locality data available|
|Tupaia palawanensis||FMNH 168969||Philippines, Palawan: Palawan; Puerto Princesa Municipality, Tarabanan River|
|Tupaia splendidula||UMMZ 1744282, 1744291, 1746523||Indonesia, Borneo: West Kalimantan; Ketapang Regency, Gunung Palung National Park, Cabang Panti Research Site|
|Tupaia tana||MVZ 1921931||Indonesia, Sumatra: Ketambe Research Station|
|JS M112||Malaysia, Borneo: Sabah; Danum Valley Field Center|
|USNZ 1090463, 1098504||Zoo animals; no locality data available|
|Urogale everetti||FMNH 147781||Philippines, Mindanao: Bukidnon; Mount Katanglad Range|
Notes: Voucher information for specimens sequenced in this study is given, arranged alphabetically by genus and species. Superscripts for species represented by 3 or more specimens correspond to Figures 3 and and4.4. Museum or zoo catalog numbers are given in the second column. Collecting localities are as provided by the loaning institution (more precise descriptions are available from the individual museums). Abbreviations: FMNH=Field Museum of Natural History; UAM =University of Alaska Museum; USNM=United States National Museum of Natural History [Smithsonian Institution]; USNZ =United States National Zoo; MVZ =Museum of Vertebrate Zoology, University of California, Berkeley; UMMZ =University of Michigan Museum of Zoology; JS =specimens collected by J. South and deposited at the Universiti Malaysia Sabah.