On the phylogeny of Mustelidae subfamilies: analysis of seventeen nuclear non-coding loci and mitochondrial complete genomes

Background Mustelidae, as the largest and most-diverse family of order Carnivora, comprises eight subfamilies. Phylogenetic relationships among these Mustelidae subfamilies remain argumentative subjects in recent years. One of the main reasons is that the mustelids represent a typical example of rapid evolutionary radiation and recent speciation event. Prior investigation has been concentrated on the application of different mitochondrial (mt) sequence and nuclear protein-coding data, herein we employ 17 nuclear non-coding loci (>15 kb), in conjunction with mt complete genome data (>16 kb), to clarify these enigmatic problems. Results The combined nuclear intron and mt genome analyses both robustly support that Taxidiinae diverged first, followed by Melinae. Lutrinae and Mustelinae are grouped together in all analyses with strong supports. The position of Helictidinae, however, is enigmatic because the mt genome analysis places it to the clade uniting Lutrinae and Mustelinae, whereas the nuclear intron analysis favores a novel view supporting a closer relationship of Helictidinae to Martinae. This finding emphasizes a need to add more data and include more taxa to resolve this problem. In addition, the molecular dating provides insights into the time scale of the origin and diversification of the Mustelidae subfamilies. Finally, the phylogenetic performances and limits of nuclear introns and mt genes are discussed in the context of Mustelidae phylogeny. Conclusion Our study not only brings new perspectives on the previously obscured phylogenetic relationships among Mustelidae subfamilies, but also provides another example demonstrating the effectiveness of nuclear non-coding loci for reconstructing evolutionary histories in a group that has undergone rapid bursts of speciation.

Previously, molecular studies of the phylogenetic reconstruction of subfamilies within the Mustelidae were based on short fragments of nuclear and mt DNA. Recently, some efforts have been made to obtain the phylogenetic tree based on large datasets, including those using 12 mt protein-coding genes [15], 5 nuclear genes [4,12], 4 nuclear genes and 1 mt gene [10,14], 21 nuclear genes and 1 mt gene [5] and 25 nuclear genes and mt genes [13]. Despite numerous efforts, however, evolutionary relationships among the Mustelidae subfamilies remain controversial (see Figure 1).
Hence, it is necessary to exploit larger independent sources of phylogenetic characters to clarify these enigmatic problems. Several studies have shown that relative to the commonly used nuclear protein-coding and mt genes, the noncoding intron sequences can be an equally fruitful source of phylogenetic characters as they possess a number of traits that are desirable for molecular phylogenetics [19][20][21][22][23][24], for example, lack of functional constraints, a high substitution rate and less homoplasy [19,25,26]. In these studies, the nuclear introns have been shown to provide powerful complementary data to address the ambiguous relationships of different taxonomic levels, including the beaked whale species [23], the Asian pitvipers genus [19], the carnivoran families [24,27], and the eutherian orders [21].
In the present study, we aim to sequence 17 nuclear intron loci comprising a total of >15 kb from 17 mustelids. The mustelids examined here represent all subfamilies of Mustelidae, except for Mellivorinae and Galictinae. Of the 17 nuclear loci, 14 were first applied in the studies of Mustelidae phylogeny. In addition, we undertook the sequencing of the mt genome from these species and presented a phylogeny based on the mt genome data currently available for mustelids. Our objectives were to: (1) provide new insights into the relationships among the Mustelidae subfamilies, and (2) examine the utilities and evolutionary dynamics of the nuclear and mt genes in the context of Mustelidae phylogeny, with special attention to the previously unexplored nuclear intron genes.

Characteristics of the Nuclear Intron Data and Mt Genomes
The general characteristics of the nuclear intron data and mt genomes are summarized in Table 1. The 17 nuclear introns of 21 species varied in length from 642 (Fgb-7) to 1685 (Plod2-14) aligned positions. The removal of ambiguous areas resulted in length variation of the aligned sequences from 523 (Fgb-4) to 1332 (Coro1c-4) positions. The numbers of parsimony-informative sites range from 71 (13.58%) (Fgb-4) to 343 (30.25%) (Cidea-1). According to different gap selection criteria in Gblocks, the alignment of the combined dataset comprised 15688 (allowed gap positions = all), 15038 (allowed gap positions = with half) and 12570 (allowed Figure 1 Hypotheses of phylogenetic relationships among Mustelidae subfamilies. Trees were reconstructed based on (a) 46 morphological characters [8], (b) analyses of 21 nuclear genes and 1 mt gene [5], (c) analyses of 5 nuclear genes and 1 mt gene [10] (support values are indicated above the line) and 25 nuclear genes and mt genome [13] (support values are indicated below the line), (d) analyses of 5 nuclear genes [4] (support values are indicated above the line) [12] support values are indicated below the line), (e) supertree analyses of 5 nuclear genes [4]. gap positions = none) positions. The parsimony-informative sites in these three datasets are 2170 (13.83%), 2137(14.21%) and 1735 (13.80%), respectively. An A-T bias (average = 58.64%) and low transition (Ti)/transversion (Tv) rate ratio (average = 1.93) were observed in most introns. In addition, most introns showed gamma shape parameters (α) close to or larger than 1.0. The nuclear sequence divergence among ingroup taxa ranged from 3% (Wasf1-7) to 8.2% (Tbc1d7-6), and averaged 4.7%.
The complete mt genomes of 25 species ranged from 16388 to 16623 bp in size. Length differences are largely due to the variation in tandem repeats within the control region. All genomes shared the same 13 proteincoding genes, 22 tRNAs genes, 2 rRNAs, and a control region, and also the same gene order. These mt genomes are apparently AT-biased (average = 60.1%). The sequence divergence among ingroup taxa ranged from 16.3 (ND6) to 23% (ND2) for the protein-coding dataset (average 20.4%), from 9.2 (12S rRNA) to 10.5% (16S rRNA) for the rRNA dataset (average 9.85%), 7.3% for the tRNA dataset, 11.4% for the control region, and 16.6% for the complete dataset.

Occurrence of Transposable Elements (TEs)
In our intron datasets, pervasive transposable element (TE) insertions were discovered (Table 2), which are mainly non-long-terminal repeat retrotransposons (non-LTR), e.g., long interspersed elements (LINEs), short interspersed elements (SINEs), mammalian-wide interspersed repeats (MIRs), and DNA transposons. MIRs and DNA transposons integrated into the orthologous loci of all examined species, which suggested an ancient origin. The result was consistent with the earlier finding that these two classes of TEs represented remnants or "fossils" of TEs, predating the radiation of mammalian orders, and had long ago become inactive in mammalian lineages [28][29][30].
The great majority of LINEs and SINEs identified here were members of the L1_Canid (Fc) and CAN SINE groups that have been exclusively found in Carnivora [31][32][33][34][35][36][37][38]. Most of SINEs showed restricted taxonomic distributions and were characterized by sporadic locations in the intronic regions. This suggests that those SINEs emerged after species diversification, likely retaining their ability to retrotranspose.
The insertions of TEs at genomic sites are often considered irreversible and random [39], which suggests that they may be excellent homoplasy-free markers in phylogenetic analyses [40][41][42][43]. Here, we identified one SINE insertion shared by Martes flavigula, Martes zibellina, Martes foina, Martes pennanti, Martes amaricana, and Gulo gulo, supporting their close relationship and the monophyly of Martinae subfamily.

Occurrence of Intra-individual Allele Heterozygotes (IIAHs)
The overall incidence of intra-individual allele heterozygotes (IIAHs) in our 14 new introns appears universal (Table 3). There were 106 cases of IIAHs in total. Of the 21 species examined, 3 to 11 cases of IIAHs were detected from each intron. IIAHs were observed to be either of equal or variable length. 11 of 14 introns had allele length variant heterozygotes due to a 1-bp indel, with the other nucleotide sites either the same or distinct at 1-11 bp. IIAHs of identical length were discovered in all introns with 1-10 substitutional differences.
Generally, IIAHs formed monophyletic pairs on the phylogenetic trees as expected (see Additional file 1). Two cases of strongly-supported nonmonophyletic IIAHs were illustrated by the close relatedness of one allele of the least weasel Mustela nivalis to one allele of the European polecat Mustela putorius (Plod2-13 and Guca1b-3 genes; Figure 2), which are most likely to indicate the cases of incomplete lineage sorting.

Phylogenetic Inference
Although individual nuclear gene analyses produced inconsistent topologies with low levels of support (Additional file 1), possibly due to limited phylogenetic information harbored in a single gene, the analyses of the combined nuclear data set using three gap selection criteria in Gblocks (allowed gap positions = none, with half, and all) and different tree-building methods (MP, ML and Bayesian methods) yielded nearly identical, well-resolved trees with strong support for all nodes, except for the relationships among Gulo gulo, Martes americana, and Martes pennanti ( Figure 2). In the tree, Taxidiinae diverged first (MP BS = 100%, ML BS = 100%, PP = 1.00), followed by Melinae (MP BS = 85%, ML BS = 98%, PP = 1.00). The remaining mustelids were divided into two clades, one consisting of Martinae and Helictidinae (MP BS = 99%, ML BS = 99%, PP = 1.00), and the other one consisting of Lutrinae and Mustelinae (MP BS = 82%, ML BS = 99%, PP = 1.00). In addition, both the inclusion of IIAHs in the combined nuclear analysis using software POFAD [44] and the Bayesian concordance analysis (BCA) analysis (Additional file 2) using software BUCKy [45] produced the same tree topologies as that in Figure 2. At the subfamily level, all nodes in the BCA analysis received high concordance factors (CF) value (1.000).
For individual mt gene analyses, the rRNA and tRNA data sets demonstrated reduced resolving power for phylogenetic inference compared to protein-coding gene analysis (Additional file 3). The complete mtDNA genome-based analyses, irrespective of the used tree-building methods and parameter sets, produced a well-resolved and well-supported tree ( Figure 3), with the tree topology and branch supports identical to that of the combined protein-coding gene analysis. At the subfamilial level, the single difference between mt genome tree ( Figure 3) and combined nuclear gene tree ( Figure 2) is the phylogenetic position of Helictidinae. In mt genome tree, Helictidinae is closer to the clade uniting Lutrinae and Mustelinae than to Martinae (MP = 60%; ML BS = 95%; PP = 1.00). The phylogenetic tree reconstructed from the combined nuclear and mt genome data set by using the BCA analysis [45] produced a tree topology (Additional file 2) identical to that from the combined nuclear gene analysis (Figure 2). At the subfamily level, all nodes received high CFs (1.000), except for that of Helictidinae (CF = 0.502), whose position is the single discrepancy between the combined nuclear and mt gene trees.
Shimodaira-Hasegawa (SH) test and the approximately unbiased (AU) test were carried out to examine the degree of significant difference between the nuclear and mt trees produced in the present study. Both tests indicated significant topological incongruence concerning the phylogenetic position of Helictidinae between the nuclear and mt trees. When using nuclear data, the mt genome tree topology, in which Helictidinae is closer to the clade uniting Lutrinae and Mustelinae, was rejected by the AU and SH tests (P < 0.05). When using mt genome data, the nuclear tree topology, in which Helictidinae and Martinae are grouped together, was rejected by the AU and SH tests (P < 0.05).

Divergence Time Estimation
Divergence time estimates for the origin and diversification of Mustelidae subfamilies yielded broadly consistent results between combined nuclear and mt genome data set (Table 4). Combined nuclear gene analysis placed the earliest branching Taxidiinae around 23.73 Mya (95% confidence intervals = 22.80-24.70 Mya). After

Phylogeny of Mustelidae Subfamilies
Among mammalian phylogenies, those characterized by rapid species radiations have long been one of the plaguing and challenging problems in species tree reconstruction [46]. This is the first study utilizing data from such large-scale nuclear non-coding loci from Mustelidae.
Both our combined nuclear intron and mt genome phylogenies not only strongly favor the prevailing view that Taxidiinae was the most basal member within family Mustelidae [4,5,[10][11][12][13][14], but also provide strong evidence that Melinae diverged between Taxidiinae and all the other mustelids examined as well. The latter is in contradiction to morphological investigations [8], but supports the nuclear gene results from Sato et al. [11,12] and Wolsan and Sato [13], and disagree those from Koepfli et al. [5] and Yu et al. [15].
Notably, the sister relationship between Lutrinae and Mustelinae was reinforced by consistent recovery from both our mt genome and combined nuclear intron analyses with high confidence, upholding and strengthening the hypothesis drawn by almost all sequence-based analyses in previous studies [4,5,7,[10][11][12][13][14][15]. In contrast, the position of Helictidinae varied between our nuclear and mt genome analyses. Nuclear data analysis placed it as sister to Martinae, whereas mt genome data indicated a sister-taxa association of it to the clade uniting Lutrinae and Mustelinae. Our mt genome result (Figure 3) is consistent with that inferred from most previous nuclear studies [3][4][5][10][11][12][13], but dissented from morphological view and karyological analyses [8,[47][48][49]. Interestingly, our combined nuclear analysis ( Figure 2) yields a result that is different from all previous hypotheses, suggesting for the first time Helictidinae and Martinae are more closely related to each other than any other taxa in Mustelidae.
Corresponding tests (AU and KH tests) have indicated significant topological incongruence between nuclear and mt trees. When using nuclear data, the mt genome tree topology was rejected by the AU and SH tests (P < 0.05), and vice versa. Phylogenetic incongruence between nuclear and mitochondrial genes has also been reported in Drosophila, Aves and bears [50][51][52][53]. Various elements may bear the responsibility for the presence of conflicting signal regarding the placement of Helictidinae, including different evolutionary histories and gene properties in gene regions from different genomes, sampling error and lineage sorting. The probability of their occurrence increased especially when separation time between different species is short [54][55][56][57], as in the present study. Although the BCA analysis of the combined nuclear intron and mt genome sequences, which is an approach that allows for gene tree discordance, retrieved  Node numbers correspond to those indicated in Figure 2 and 3. a the divergence times were estimated based on combined nuclear genes and tree topology in Figure 2. b the divergence times were estimated based on mt genomes and tree topology in Figure 3. c 95% confidence intervals. Mt gene nd5 8 1 9 * * * # * * * * * an identical tree topology to that of the nuclear intron gene, the position of Helictidinae received week support. Therefore, the analyses involving more characters in the future may help to confirm the precise position of Helictidinae.

Implications for Mustelidae Radiation
Mustelidae has one of the most extensive fossil records of extant Carnivora families [8], and molecular dating of the Mustelidae radiation has also been attempted in several studies previously. Our results, from an independent character source, provide important insights into the time scale of the origin and diversification of extant subfamilies of Mustelidae. It is interesting to draw a comparison of the dating results estimated from prior and present studies (Table 4). Among the subfamilies, Taxidiinae is the basal-most branching lineage in mustelid diversification dating to the late Oligocene, which is more concordant with the time estimate of Koepfli et al. [5], but earlier than those of Bininda-Emonds et al. [58], Yonezawa et al. [15] and Eizirik et al. [59]. The divergence time of the next branching lineage, i.e., Melinae clade, is dated to the early Miocene in our analysis, which is more broadly consistent with the paleontological data [60] and the sequence-based data from Yonezawa et al. [15] than those from Sato et al. [2], Bininda-Emonds et al. [54], Koepfli et al. [5], and Eizirik et al [59]. The estimates of the origin of the remaining mustelids are more recent than those fossil-based ages [61], and older than the other sequence-based dates [2,5,15,58,59].
Among the remaining mustelids, one important event of the mustelid diversification is the divergence between Lutrinae and Mustelinae. Our estimated dates are consistent with the corresponding fossil records [60,62] and that from Yonezawa et al. [15], but more recent than those from Hosoda et al. [63], Wayne [64], and Sato et al. [2], and older than those from Koepfli et al. [5] and Eizirik et al. [59]. Other important events are the origins of Mustelinae and Martinae. Our estimates for them are both much older than the existing fossil record [65,66]. The dates of Mustelinae origin are older than most of the other molecular estimates, but in good agreement with that of Koepfli et al. [5]. As regards the origin of Martinae, there is large difference between our nuclear and mt genome estimates. The nuclear estimate is more in agreement with those from Koepfli et al. [5] and Eizirik et al. [59], while the mt estimate is more concordant with those from Hosoda et al. [63] and Yonezawa et al. [15].
In addition, our analyses resulted in time estimates of divergence of Helictidinae that more agree with that from Koepfli et al. [5] than that from Bininda-Emonds et al. [58], which is younger than the present results.
More intensive taxonomic sampling will improve accuracy on the time estimation of mustelid diversification.

Utilities of the nuclear introns in phylogenetic study of Mustelidae subfamilies
Several recent studies have indicated that nuclear introns hold considerable signals for resolution of difficult phylogenies at both shallow and deeper species level hierarchies [21][22][23][24]27,67,68]. We are among the first to use large-scale nuclear intron genes in inferring phylogenies of Mustelidae. Our analysis not only brings new perspectives on the phylogenetic relationship of Mustelidae subfamilies, but provides another example demonstrating that the nuclear non-coding genes can be an effective data source for reconstructing evolutionary histories in a group that has undergone rapid bursts of speciation as well.
We assessed the phylogenetic utilities of individual introns and mt genes in resolution of the inter-subfamilial relationships of Mustelidae by counting the number of congruent nodes between the individual phylogenies and the combined gene trees (Table 5). In the individual nuclear gene analyses, the Plod2-13 gene recovered all 8 nodes of the combined nuclear gene tree. Anyway, the Plod2-13 and Plod2-14 genes recovered the highest number of congruent nodes of the combined nuclear gene tree, whereas the Wasf1-7, Guca1b-3, Ssr1-5, and Tbc1d7-6 genes showed the lowest phylogenetic performance. As regards the mt gene analyses, we observed that the ND5 and CYTB genes recovered all 9 nodes of the mt genome tree. Ranking the single mt gene shows that the ND5, CYTB, 16SrRNA and ND2 genes are better indicators of Mustelidae phylogeny at subfamilial level than are other genes, such as Atp8, Atp6 and ND4L genes. This result agrees broadly with previous conclusions about the rough classification of mt genes into good, medium, and poor performance categories [69][70][71][72][73][74] (Additional file 4). In summary, the assessment of phylogenetic utility and limits of these individual nuclear and mt genes makes it possible to preselect subsets of genes for future molecular studies of vertebrate phylogeny.
Although the use of nuclear introns as genetic markers has now been implemented as a powerful approach in recent phylogenetic studies, there are, however, numerous potential problems associated with this approach. Chief among these is the difficulties in sequence data acquisition, alignment, and analysis of the nuclear introns compared to the traditional mt and nuclear protein-coding genes, as a result of the higher rates of variation and frequencies of indels. The common presence of indels (including TEs, small gaps and tandem repeats) and IIAHs, as reported in this study, makes experimental work labor-intensive by virtue of the additional time and money required to isolate alleles and optimize PCR amplification and sequencing. In addition, they can create positional homology problems associated with areas of ambiguous alignment [75]. Several studies have also shown that the phylogenetic inference is sensitive to the various treatments of indels. These issues are central to the appropriate application of intron data in phylogenetic reconstruction and they should be comprehensively and explicitly addressed in the future studies [19].

Conclusions
The phylogenetic relationships among Mustelidae subfamilies have posed one of the major problems concerning Carnivora systematics. In this study, phylogenetic relationships among Mustelidae subfamilies are presented based on 17 nuclear intron loci and mt whole genomes. Our results resolve some of the ambiguous issues in Mustelidae phylogeny, whereas some phylogenetic relationships require confirmation by analyzing additional samples and character information, such as the precise position of Helictidinae. Our study not only brings new perspectives on the previously obscured phylogenetic relationships among Mustelidae subfamilies, but also provides another example demonstrating the effectiveness of nuclear non-coding loci for reconstructing evolutionary histories in a group that has undergone rapid bursts of speciation.

Sequence Data
Detailed information of the 14 nuclear intron loci first used in the Mustelidae phylogeny is shown in Additional file 5. These loci were amplified with primers as described in Yu et al. [27] from 17 mustelids and 4 nonmustelid carnivoran species. Three other nuclear introns (Ttr-1, Fgb-4 and Fgb-7), which were available from our previous published studies [14,76,77], were also included in the present nuclear data set.
The species examined in this study and their Genbank accession numbers are listed in Additional file 6. A "touch-down" PCR amplification was carried out using the following parameters: 95°C hot start (5 min), 10 cycles of 94°C denaturation (1 min), 60-50°C annealing (1 min), 72°C extension (1 min), and finally 25 cycles of 94°C denaturation (1 min), 50°C annealing (1 min), 72°C extension (1 min). The amplified DNA fragments were purified and sequenced in both directions with an ABI PRISM(tm) 3730 DNA sequencer following the manufacturer's protocol. In the case of poor performance of direct sequencing resulting from complex DNA structures, tandem repeats or intron heterozygotes, the amplified PCR products were gel-purified and cloned into the pMD18-T vector (TaKaRa Biotechnology Co., Ltd. Dalian, China) and transformed into ultracompetent E. coli cells (TaKaRa Biotechnology Co., Ltd. Dalian, China). Thirty positive clones per ligation reaction were sequenced. All sequences obtained were checked carefully and queried in BLAST searches of GenBank to assess homology. In a few cases, PCR attempts using different primer pairs and cloning methods failed to produce sequence data (see Additional file 6). These sequences were excluded from the independent gene analyses and treated as missing data in the combined analyses. The newly determined nuclear sequences have been deposited in GenBank with accession numbers HM063147-HM063412.
The mt complete genome sequences were amplified using LA PCR™ Kit (Takara Biotechnology Co., Ltd) and 9 universal long PCR primers from 14 mustelids and 2 non-mustelid carnivoran species (Additional file 6). In addition, 23 species-specific primers were designed when the universal PCR primers failed to produce successful PCR amplification. A "touch-down" long PCR amplification was carried out using the following parameters: 95°C hot start (2 min), 10 cycles of 98°C denaturation (10 sec), 67-58°C annealing (1 min 30 sec;°C /cycle), 72°C extension (5 min), and finally 25 cycles of 98°C denaturation (10 sec), 58°C annealing (1 min 30 sec), 72°C extension (5 min). At the end, a final 10-min extension at 72°C was performed. Long PCR products were sequenced in both directions using a primer walking strategy with a total of 306 primers. Sequencing was performed in an ABI PRISM(tm) 3700 DNA sequencer following the manufacturer's protocol. Primer sequence information is available upon request. Where necessary, PCR products were cloned into the pMD18-T vector and transformed into ultracompetent E. coli cells (TaKaRa Biotechnology Co., Ltd. Dalian, China) in order to resolve the difficulty of direct sequencing of control regions arising from long tandem repeats. Five positive clones per ligation reaction were sequenced. Mt sequences obtained were checked to ensure that they did not include nuclear copies of mtDNA-like pseudogenes. In addition to these new 16 mt genomes, 6 other mustelid and 3 other non-mustelid mt genomes available in public database were included in the mt analyses (Additional file 6).

Alignments and sequence Characterizations
Sequences were aligned using CLUSTAL X under the default settings [78]. The nuclear alignment was divided into two data sets: (1) each of the 17 intron loci and (2) combined sequences of all introns. The mt alignment was divided into four data sets: (1) each of the 13 protein-coding genes, (2) 22 tRNAs, (3) two rRNAs, and (4) combined sequences of the tRNAs, rRNAs and proteincoding genes. Due to the presence of ambiguous areas in the nuclear alignments, we excluded the alignment ambiguities using Gblocks 0.91b [79] with default parameters (allowed gap positions = all).

Phylogenetic Analyses
Phylogenetic analyses of the individual introns and mt genes, i.e., nuclear data set (1) and mt data set (1), (2), and (3), were performed using PAUP* 4.0b10 [82] for maximum parsimony (MP) and maximum likelihood (ML) analyses, and using MrBayes 3.1.2 [83] for the Bayesian inference. In MP analyses, a heuristic search was performed with tree-bisection-reconnection (TBR) branch swapping, random addition of taxa, and 1000 replicates per search. Only one of the best trees found during branch swapping was saved. In ML analysis, the best-fit models of sequence evolution were selected using the Akaike Information Criterion (AIC) [84,85] with Modeltest version 3.7 [86]. The chosen models and their parameters were used to infer ML trees with the heuristic algorithm, 10 random-addition sequence replicates, and TBR branch swapping. The tree reliability under ML analysis was assessed using a bootstrap resampling of 100 replicates (BP) [87]. In Bayesian inference, each Metropolis-coupled Markov chain Monte Carlo (MCMC) run for all individual genes employed the model selected by ModelTest for that gene, or the nearest model to that model that could be implemented in MrBayes. Three heated chains and a single cold chain were used in all MCMC analyses and run for 2 × 10 6 generations, Trees were sampled every 100 generation. The average standard deviation of split frequencies was close to 0.001 when the run was end. The first 25% were discarded as the burn-in. A 50% majority-rule consensus of post burn-in trees was constructed to summarize posterior probabilities (PP) for each branch.
In addition to individual analyses, phylogenetic reconstruction were performed based on the combined data sets, i.e., nuclear data set (2) and mt data set (4), using PAUP* 4.0b10 [82] for MP analysis, RAxML online web server [88] for partitioned ML analysis and using MrBayes [83] for partitioned Bayesian analysis (pBI) [89]. For the combined data sets, we identified model partitions based on partitioning matrices by locus. That is, in the analysis of combined nuclear data set, each nuclear intron gene was considered as a different partition, whereas in that of combined mt data set, each of the 13 individual protein-coding genes, all tRNAs, and each of the two rRNA genes were considered as different partitions. Based on the selected models using the AIC [84] as mentioned above for individual analyses, we assigned a separate substitution model for each of the data partitions. Three heated chains and a single cold chain were used in all MCMC analyses and run for 5 × 10 6 generations, sampling trees every 100 generations. The average standard deviation of split frequencies was close to 0.001 when the run was end. The first 25% were discarded as the burn-in. A 50% majority-rule consensus of post burn-in trees was constructed to summarize posterior probabilities (PP) for each branch.
In addition, given the heterogeneous gene trees observed among 17 nuclear intron gene analyses and between the combined nuclear intron and mt genome analyses, Bayesian concordance analysis (BCA) [45] implemented in the program BUCKy [45], which uses individually calculated gene trees to infer the species tree that maximizes the bipartition concordance among each gene tree, was also performed for both the nuclear intron gene datasets and the nuclear plus the mt genome datasets. Tree reliability was evaluated by sample concordance factors (CFs). The MCMCMC sampled with 2 × 10 6 generations was employed (4 runs and 4 chains) and a priori level of discordance α = 2.5 was used in BCA.
In all analyses, trees were rooted with the red panda Ailurus fulgens and the skunk Mephitis mephitis, based on the general consensus that they branched off earlier than the mustelids [2,4,11,12,76,77].

Intra-individual allele heterozygotes
For individual intron analyses, both copies of alleles from a species were included. For combined nuclear data sets, we performed phylogenetic analyses by (1) choosing randomly an allele per species for portioned ML analysis and Bayesian analysis (without the inclusion of IIAHs), and (2) using POFAD v1.03 algorithm (Phylogeny Analysis From Allelic Data) [90] to incorporate IIAHs. POFAD is a recently developed method of constructing phylogeny from multiple datasets that contain allelic information. It converts a distance matrix of alleles into a distance matrix of organisms so that individuals become the terminals of the analyses. First, we calculated the average uncorrected pairwise distances in PAUP* [82]. These distances then served as the input for the calculation of standardized pairwise distances between species in POFAD [90]. The standardized distances were then used as input for the neighbor-joining analysis conducted using PAUP* [82] to produce a phylogenetic tree.

Divergence time estimation
Divergence times based on combined nuclear intron and combined mt data sets were estimated using the Bayesian relaxed phylogenetic approach implemented in BEAST v1.5.4 [91]. We assumed a GTR+I+G model of DNA substitution with four rate categories. Uniform priors were employed for GTR substitution parameters (0, 100), gamma shape parameter (0, 100) and proportion of invariant sites parameter (0, 1). The uncorrelated lognormal relaxed molecular clock model was used to estimate substitution rates for all nodes in the tree, with uniform priors on the mean (0, 100) and standard deviation (0, 10) of this clock model. We employed the Yule process of speciation as the tree prior and a UPGMA tree to construct a starting tree, with the ingroup assumed to be monophyletic with respect to the outgroup.
Three calibration points from the fossil records were applied in the dating analyses. These calibration points are all implemented as minimum age constraints, including 27.6 Mya for the split between Procyonidae and Mustelidae [12,15,92], 24 Mya for the crown Mustelidae [2,5,17], and 5.3 Mya for the origin of the genus Mustela in Mustelinae [5,93]. All fossil constraint priors were set as means of a normal distribution, with a standard deviation of 1.0 MYA. Two independent MCMC runs of 30,000,000 generations were performed for each data set with parameters logged every 1,000 generations. The Auto Optimize Operators function was enabled to maximize efficiency of MCMC runs. Two independent MCMC runs for each analysis were combined to estimate the posterior distribution of the substitution model and tree model parameters, as well as node ages. Analyses of these parameters in Tracer 1.5 [94] suggested that the number of MCMC steps was more than adequate, with effective sample sizes of all parameters often exceeding 1,000 and Tracer plots showing strong equilibrium after discarding burn-in.

Testing Tree Incongruence
The incongruence among different tree topologies was evaluated using the Shimodaira-Hasegawa (SH) test [95] and the approximately unbiased (AU) test [96], as implemented in the CONSELV0.1i program [97] with default scaling and replicate values. The site-wise loglikelihood values were estimated by PAUP* [82].