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.  from 17 mustelids and 4 non-mustelid 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 . 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 protein-coding genes. Due to the presence of ambiguous areas in the nuclear alignments, we excluded the alignment ambiguities using Gblocks 0.91b  with default parameters (allowed gap positions = all).
Pairwise comparisons and sequence characterizations were estimated using MEGA 4.0 . Given that nuclear introns tend to be favorable chromosomal regions for integration of transposable elements (TEs) , they were screened for interspersed repeats by using the program RepeatMasker (Smit, Hubley and Green, RepeatMasker Open-3.0. 1996-2004, http://www.repeatmasker.org).
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  for maximum parsimony (MP) and maximum likelihood (ML) analyses, and using MrBayes 3.1.2  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 . 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) . 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 × 106 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  for MP analysis, RAxML online web server  for partitioned ML analysis and using MrBayes  for partitioned Bayesian analysis (pBI) . 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  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 × 106 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)  implemented in the program BUCKy , 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 × 106 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)  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* . These distances then served as the input for the calculation of standardized pairwise distances between species in POFAD . The standardized distances were then used as input for the neighbor-joining analysis conducted using PAUP*  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 . 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  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  and the approximately unbiased (AU) test , as implemented in the CONSELV0.1i program  with default scaling and replicate values. The site-wise log-likelihood values were estimated by PAUP* .