BACKGROUND
subterranean mammals have received extensive attention for their specialized adaptations to life underground from scientists in different fields  <cit> .
three orders of mammals include subterranean species: the rodents, the insectivores, and the marsupials, which independently display adaptive convergence and divergence, are good subjects for comparative studies in morphology, ecology, physiology, ethology, and most recently, genetics  <cit> .
all extant species of zokors , blind mole rats , and bamboo rats  are typical subterranean mammals  <cit> .
because of their cryptic lifestyle and the morphological convergence associated with the subterranean life in tubular burrows  <cit> , their systematic placement has been a long-standing puzzle  <cit> .
members of these subterranean rodents had been assigned to muridae , cricetidae , and rhizomyidae   <cit> .
recently, overwhelming molecular evidence has suggested the cladistic reunion of zokors, bamboo rats and blind mole rats, and the monophyletic origin of the three groups that form the family spalacidae has been generally accepted after a century of debate  <cit> .
within the family spalacidae, however, the phylogenetic affinities among the three subfamilies  remain unclear  <cit> .
mitochondrial genes  recovered a sister clade consisting of spalacinae and rhizomyinae  with moderate bootstrap support values of 63% or 71% depending on phylogenetic approaches, while myospalacinae is basal relative to the two additional subfamilies with a strong support values of 100%  <cit> .
by contrast, nuclear gene irbp suggested that myospalacinae is the sister to rhizomyinae , although the support  is fairly weak  <cit> .
by adding another nuclear gene, the concatenated genes of irbp and ghr united myospalacinae and spalacinae to be a monophyletic group  with a low bootstrap value, and rhizomyinae is the sister to the monophyletic group  <cit> .
consistent with the phylogeny inferred from irbp sequences , paleontological evidence has shown that zokors share more morphological features with bamboo rats than with blind mole rats  <cit> , supporting that myospalacinae and rhizomyinae are most closely related, while spalacinae is more basal than the others .
these controversies prompt us to revisit the phylogeny of the three groups within the family spalacidae using genomic data.
in this study, we employed the illumina hiseq  <dig> platform and sequenced the transcriptomes of brain and liver of a plateau zokor  and a hoary bamboo rat .
we compiled newly obtained transcriptome data from the zokor and the bamboo rat alongside the published transcriptome assembly from the middle east mole rat  <cit> , and provided a conclusive resolution of interrelationships of the three subfamilies within the family spalacidae using a phylogenomic approach.
RESULTS
illumina sequencing and assembly
total rnas were isolated from brain and liver of an adult male plateau zokor.
mrna was purified by passing total rna through a column of beads with oligo , and was subsequently fragmented into short sequences of 200 bp.
the fragmented mrna was ligated with the adaptor mix and was transcribed to cdna by reverse transcriptase before sequencing.
the sequenced raw data were filtered to obtain clean reads by removal of adaptors, low-quality reads, and ambiguous reads.
a total of  <dig> , <dig> and  <dig> , <dig> clean reads were obtained for the brain and liver of the zokor, respectively .
as for the bamboo rat,  <dig> , <dig> and  <dig> , <dig> clean reads were obtained for the brain and liver, respectively  .
the reads of brain and liver of each species were pooled and were assembled into unigenes.
eventually, a total of  <dig>  filtered transcripts  of the zokor were obtained, and from which,  <dig>  orfs  with a mean length of 975 bp were identified, while  <dig>  unigenes were generated and  <dig>  orfs with a mean length of  <dig>  bp were detected from the bamboo rat .
details of the number of reads, mean length and the n <dig> values are given in table  <dig>
functional annotation and classification
in total,  <dig>   and  <dig>   unigenes of the zokor and the bamboo rat, respectively, were matched at least once by either of the public databases: swiss-prot, trembl, ncbi non-redundant protein , ncbi nucleotide sequence , gene ontology , cluster of orthologous groups , and kyoto encyclopedia of genes and genomes  .
for example, we used blastx to search all possible translations of unigenes of the zokor against the swiss-prot protein database with an e-value of  <dig> e- <dig>  and obtained  <dig>  hits, covering  <dig> % of all unigenes .
specifically, both go  and cog clustering showed similar distributions of the unigenes between the two species, suggesting the two transcriptomes are comparable .
orthologous gene identification and multiple sequence alignment
the inparanoid program was employed to identify orthologs based on two-way best genome-wide pairwise matches for each species pair  <cit> .
the automatic clustering method implemented in the program can distinguish between orthologs and paralogs without using multiple alignments and phylogenetic trees  <cit> .
all possible pairwise orthologs identified from the inparanoid program were introduced into multiparanid program to generate multiple-species orthologs  <cit> .
all examined species include the plateau zokor, the hoary bamboo rat, the middle east blind mole rat, and mouse which was used as an outgroup.
multiple sequence alignments of each ortholog were performed with the macse  <cit>  and the prank  <cit>  programs.
each alignment was check manually in mega  <dig>  <cit> , alignments with poor quality or aligned region <100 bp were discarded.
the resulting multiple sequence alignments include:   <dig>  individual alignments for each nuclear gene ,   <dig> individual alignments for each protein-coding mitochondrial gene,  one alignment  for the concatenated  <dig> nuclear genes,  one alignment  for the concatenated fourfold degenerated sites  of the  <dig>  nuclear genes,  one alignment  for the concatenated  <dig> mitochondrial genes, and  one alignment  for the concatenated 4d-sites of the  <dig> mitochondrial genes.
note that 4d-sites refer to the four-fold degenerate third codon positions of protein-coding genes at which any of the four nucleotide substitutions cannot result in an amino acid replacement, thus 4d-sites are a subset of synonymous sites that are supposed to be under no or very weak selection  <cit> .
phylogenetic tree reconstruction
to avoid the confounding effects of sequence saturation in phylogenetic analysis, we assessed the substitution saturation with dambe <dig>  <cit> .
results showed that all the  <dig> randomly selected nuclear genes, the concatenated 4d-sites of the  <dig>  nuclear genes and the concatenated  <dig> protein-coding mitochondrial genes obviously have experienced little substitution saturation .
in contrast, the concatenated 4d-sites of  <dig> mitochondrial genes showed substantial substitution saturation  and were excluded for further analysis.
we used both maximum likelihood  and bayesian inference  approaches for phylogenetic tree reconstruction.
jmodeltest <dig> program  <cit>  was used to select the best-fitting substitution model for each alignment according to the akaike information criterion and bayesian information criterion  <cit> .
for the  <dig>  individual nuclear genes, totaling  <dig> models were selected as the best-fitting models for ml approach, the top five selected models were trn + g, trn + i, hky + i, tim3 + g, and tim3 + i, covering  <dig> % of the  <dig>  cases.
totaling  <dig> models were chosen as the best-fitting models for bayesian inference , the top five selected models were hky + i, gtr + g, gtr + i, hky + g, and hky, accounting for  <dig> % of the  <dig>  cases .
for the  <dig> individual mitochondrial genes, the tim <dig> models  were most common  for ml analysis, while the gtr models  were most popular  for bi approach .
in addition, the best-fitting models of the concatenated 4d-sites of the  <dig>  nuclear genes were gtr + i + g for both ml and bi approaches.
among the phylogenetic trees based on the  <dig>  individual nuclear genes,  <dig>  and  <dig>  were supported with bootstrap values >50% for ml and bi methods, respectively.
moreover, the majority of phylogenetic trees derived from  <dig> individual mitochondrial genes were shown with bootstrap values >50% for both ml  and bi  approaches .
after counting the genes that support a tree with bootstrap values >50%, we found that ~39% of nuclear genes recovered a sister clade consisting of the zokor and the bamboo rat, and ~34% and ~27% of nuclear genes supported two other topologies .
chi-square tests showed that the distributions of ml  and bayesian  tree topologies were significantly deviated from random distribution, meaning that the majority of nuclear genes support a sister group of myospalacinae and rhizomyinae.
by contrast, no significance was detected for either ml  or bayesian trees  inferred from  <dig> individual mitochondrial genes, possibly because of small number of genes used .
we also used a more stringent cutoff by counting the gene trees with boostrap values >70%.
we found that 40% of nuclear genes supported a sister relationship between the zokor and the bamboo rat , which is significantly deviated from random distribution .
in accordance with the individual gene analyses, the totality and the concatenated 4d-sites of the  <dig>  nuclear genes revealed a same topology as the majority of nuclear genes and received 100% bootstrap support in both analyses .
although the concatenated  <dig> protein-coding mitochondrial genes recovered the same phylogeny, the bootstrap values were relatively low  .
strikingly, we observed a very short branch leading to myospalacinae and rhizomyinae , suggesting that the time between the divergence of myospalacinae and rhizomyinae and the first speciation before the divergence of spalacidae is very short.
the short period of time for the divergence of spalacidae indicates a burst of rapid speciation occurred.
this observation was previously revealed by a phylogenetic analysis inferred from the nuclear gene irbp, although the monophyletic group comprising myospalacinae and rhizomyinae was weakly supported  <cit> .
discussion
when the genome sequence is unavailable, transcriptome sequencing is an effective way to obtain large numbers of transcripts.
in this study, we present the first transcriptome data in zokors and bamboo rats using massively parallel mrna sequencing.
we perform illumina sequencing of the brain which is the most complex organ in the body and the liver which plays a major role in metabolism.
although the raw data of the bamboo rat were much larger than that of the zokor, the results of sequence assembly, annotation, and classification of the two datasets were comparable.
furthermore, we provide a conclusive resolution of phylogenetic placement of the three subfamilies in the family spalacidae.
phylogenetic analyses using various data coding schemes and analytical approaches overwhelmingly support that the subfamily myospalacinae is more closely related to the subfamily rhizomyinae than to the subfamily spalacinae.
although we did not examine multiple individuals of each species, the genetic variations within species should not affect our phylogenetic resolution at subfamily levels.
traditional sanger sequencing technology produces longer est  sequences with low throughput, which becomes little advantages when new assembly programs can handle large throughput short sequences  <cit> .
further, traditional transcriptome studies generally use a cdna library of one tissue, while transcriptome sequencing based on ngs methods can work on a normalized cdna library comprising multiple tissues and individuals, which would considerably reduce the sequencing cost.
hence, the cost- and time-effective ngs technologies will provide more comprehensive transcriptome information and facilitate the transcriptome studies, particularly in less-studied species.
although the large-scale transcriptome sequencing was used to sequence cdna pools derived from spalax galili <cit> , a member of the subfamily spalacinae, the transcriptome surveys on the members of the other two subfamilies within the family spalacidae are virtually lacking.
our work provides a valuable genomic resource and will stimulate the further analysis of these taxa with exceptional scientific interest.
for example, subterranean rodents are believed to have greater hypoxia tolerance than their aboveground relatives  <cit> .
upon examining  <dig> candidate genes involving in adaptation to high-altitude hypoxia  <cit> , we identified  <dig>  and  <dig>  genes in the zokor and bamboo rat transcriptomes, respectively.
this finding will facilitate future comparative analysis of these functional genes.
before the emergence of molecular identification approaches, one determines taxonomic status and phylogenetic relationships of animals mainly rely on morphological characteristics.
however, the evolutionary change of morphological characters is extremely complicated , the phylogenetic trees derived from morphological data frequently remain controversial  <cit> .
in contrast, since the evolutionary change of dna follows a traceable pattern, it is possible to use a mathematical model to formulate the change and compare dna sequences among different organisms  <cit> .
as a result, molecular phylogenetics is expected to clarify the tree of life that has been difficult to be resolved by the classical approaches.
taking zokors as an example, these animals have been allied to several different muroid subfamilies including rhizomyinae, spalacinae, arvicolinae, and cricetinae based on morphological characteristics  <cit> , while overwhelming studies using molecular phylogenetic approaches concordantly placed zokors into spalacidae  <cit> .
it should be noticed that random noise will generally lead to poorly resolved phylogenetic trees, primarily because the number of nucleotide substitutions of a few genes is small.
this may be the major reason why the previous molecular phylogeny studies showed inconsistence of topology relationships among the three subfamilies in spalacidae  <cit> .
in this study, we used a large number of genes including  <dig>  nuclear genes and all  <dig> mitochondrial protein coding genes to examine the phylogenetic affinities of the three subfamilies within spalacidae.
our analyses based on various data sets overwhelmingly support that the spalacinae had a highest probability to be a basal clade relative to others within spalacidae, while rhizomyinae and myospalacinae form a sister group.
however, the most important challenges of phylogenomics involve different tree reconstruction methods and substitution saturation that can influence the validity of phylogenetic placements  <cit> .
to overcome the challenges, we used both ml and bi approaches based on both nuclear and mitochondrial genes, and found no incongruence between the two tree reconstruction approaches .
further, we performed saturation tests and excluded the saturated 4d-site concatenation of mitochondrial genes in our analyses.
interestingly, our proposed phylogeny of the three subfamilies within spalacidae is accordant with that inferred from massive morphological analyses on fossils of spalacidae  <cit> .
the known early miocene rhizomyines are closer to the stem zokor morphotype, suggesting that myospalacinae is more closely related to rhizomyinae than to spalacinae.
although we didn’t examine african mole rats  that were previously considered to be a separate subfamily in spalacidae, multiple lines of molecular evidence consistently supported a sister relationship between rhizomys and tachyoryctes, suggesting that tachyoryctes should be incorporated into the rhizomyinae  <cit> .
together, our phylogenomic analyses unambiguously resolve the phylogeny of the three subfamilies comprising the family spalacidae.
geographically, members of the subfamily spalacinae are found in east europe, west asia, near east, and north africa, members of myospalacinae inhabit siberia and northern china, while the rhizomyinae include two geographically distant relatives  that are from southeastern asia and east africa  <cit> .
our molecular phylogenetic evidence suggests that the common ancestor of spalacidae firstly split into spalacinae and rhizomyinae + myospalacinae clades in the northern asia adjacent to the middle east, and the latter clade subsequently split into rhizomyinae and myospalacinae in mongolia and northern china  <cit> .
interestingly, although rhizomyinae and myospalacinae are more phylogenetically related, myospalacinae appears to be more biologically similar with spalacinae than with rhizomyinae for certain traits.
for example, myospalacinae and spalacinae are totally blinded, while rhizomyinae retain their sight  <cit> .
the resolution of phylogenetic relationships among the three subfamilies will provide a framework for evolutionary studies on convergence and divergence of these animals.
CONCLUSIONS
to summarize, this work is the first report of transcriptome sequencing in zokors and bamboo rats, representing a valuable resource for future studies of comparative genomics in subterranean mammals.
phylogenomic analysis provides a conclusive resolution of interrelationships of the three subfamilies within the family spalacidae, and highlights the power of phylogenomic approach to dissect the evolutionary history of rapid radiations in the tree of life.

 subterranean mammals have received extensive attention for their specialized adaptations to life underground from scientists in different fields. three orders of mammals include subterranean species: the rodents, the insectivores, and the marsupials, which independently display adaptive convergence and divergence, are good subjects for comparative studies in morphology, ecology, physiology, ethology, and most recently, genetics. all extant species of zokors, blind mole rats, and bamboo rats are typical subterranean mammals. because of their cryptic lifestyle and the morphological convergence associated with the subterranean life in tubular burrows, their systematic placement has been a long-standing puzzle. members of these subterranean rodents had been assigned to muridae, cricetidae, and rhizomyidae. recently, overwhelming molecular evidence has suggested the cladistic reunion of zokors, bamboo rats and blind mole rats, and the monophyletic origin of the three groups that form the family spalacidae has been generally accepted after a century of debate. within the family spalacidae, however, the phylogenetic affinities among the three subfamilies remain unclear. mitochondrial genes recovered a sister clade consisting of spalacinae and rhizomyinae with moderate bootstrap support values of 63% or 71% depending on phylogenetic approaches, while myospalacinae is basal relative to the two additional subfamilies with a strong support values of 100%. by contrast, nuclear gene irbp suggested that myospalacinae is the sister to rhizomyinae, although the support is fairly weak. by adding another nuclear gene, the concatenated genes of irbp and ghr united myospalacinae and spalacinae to be a monophyletic group with a low bootstrap value, and rhizomyinae is the sister to the monophyletic group. consistent with the phylogeny inferred from irbp sequences, paleontological evidence has shown that zokors share more morphological features with bamboo rats than with blind mole rats, supporting that myospalacinae and rhizomyinae are most closely related, while spalacinae is more basal than the others. these controversies prompt us to revisit the phylogeny of the three groups within the family spalacidae using genomic data. in this study, we employed the illumina hiseq platform and sequenced the transcriptomes of brain and liver of a plateau zokor and a hoary bamboo rat. we compiled newly obtained transcriptome data from the zokor and the bamboo rat alongside the published transcriptome assembly from the middle east mole rat, and provided a conclusive resolution of interrelationships of the three subfamilies within the family spalacidae using a phylogenomic approach.  illumina sequencing and assembly total rnas were isolated from brain and liver of an adult male plateau zokor. mrna was purified by passing total rna through a column of beads with oligo, and was subsequently fragmented into short sequences of 200 bp. the fragmented mrna was ligated with the adaptor mix and was transcribed to cdna by reverse transcriptase before sequencing. the sequenced raw data were filtered to obtain clean reads by removal of adaptors, low-quality reads, and ambiguous reads. a total of, and, clean reads were obtained for the brain and liver of the zokor, respectively. as for the bamboo rat,, and, clean reads were obtained for the brain and liver, respectively. the reads of brain and liver of each species were pooled and were assembled into unigenes. eventually, a total of filtered transcripts of the zokor were obtained, and from which, orfs with a mean length of 975 bp were identified, while unigenes were generated and orfs with a mean length of  bp were detected from the bamboo rat. details of the number of reads, mean length and the n values are given in table  functional annotation and classification in total, and unigenes of the zokor and the bamboo rat, respectively, were matched at least once by either of the public databases: swiss-prot, trembl, ncbi non-redundant protein, ncbi nucleotide sequence, gene ontology, cluster of orthologous groups, and kyoto encyclopedia of genes and genomes. for example, we used blastx to search all possible translations of unigenes of the zokor against the swiss-prot protein database with an e-value of e- and obtained hits, covering % of all unigenes. specifically, both go and cog clustering showed similar distributions of the unigenes between the two species, suggesting the two transcriptomes are comparable. orthologous gene identification and multiple sequence alignment the inparanoid program was employed to identify orthologs based on two-way best genome-wide pairwise matches for each species pair. the automatic clustering method implemented in the program can distinguish between orthologs and paralogs without using multiple alignments and phylogenetic trees. all possible pairwise orthologs identified from the inparanoid program were introduced into multiparanid program to generate multiple-species orthologs. all examined species include the plateau zokor, the hoary bamboo rat, the middle east blind mole rat, and mouse which was used as an outgroup. multiple sequence alignments of each ortholog were performed with the macse and the prank programs. each alignment was check manually in mega, alignments with poor quality or aligned region <100 bp were discarded. the resulting multiple sequence alignments include: individual alignments for each nuclear gene, individual alignments for each protein-coding mitochondrial gene, one alignment for the concatenated nuclear genes, one alignment for the concatenated fourfold degenerated sites of the nuclear genes, one alignment for the concatenated mitochondrial genes, and one alignment for the concatenated 4d-sites of the mitochondrial genes. note that 4d-sites refer to the four-fold degenerate third codon positions of protein-coding genes at which any of the four nucleotide substitutions cannot result in an amino acid replacement, thus 4d-sites are a subset of synonymous sites that are supposed to be under no or very weak selection. phylogenetic tree reconstruction to avoid the confounding effects of sequence saturation in phylogenetic analysis, we assessed the substitution saturation with dambe. results showed that all the randomly selected nuclear genes, the concatenated 4d-sites of the nuclear genes and the concatenated protein-coding mitochondrial genes obviously have experienced little substitution saturation. in contrast, the concatenated 4d-sites of mitochondrial genes showed substantial substitution saturation and were excluded for further analysis. we used both maximum likelihood and bayesian inference approaches for phylogenetic tree reconstruction. jmodeltest program was used to select the best-fitting substitution model for each alignment according to the akaike information criterion and bayesian information criterion. for the individual nuclear genes, totaling models were selected as the best-fitting models for ml approach, the top five selected models were trn + g, trn + i, hky + i, tim3 + g, and tim3 + i, covering % of the cases. totaling models were chosen as the best-fitting models for bayesian inference, the top five selected models were hky + i, gtr + g, gtr + i, hky + g, and hky, accounting for % of the cases. for the individual mitochondrial genes, the tim models were most common for ml analysis, while the gtr models were most popular for bi approach. in addition, the best-fitting models of the concatenated 4d-sites of the nuclear genes were gtr + i + g for both ml and bi approaches. among the phylogenetic trees based on the individual nuclear genes, and were supported with bootstrap values >50% for ml and bi methods, respectively. moreover, the majority of phylogenetic trees derived from individual mitochondrial genes were shown with bootstrap values >50% for both ml and bi approaches. after counting the genes that support a tree with bootstrap values >50%, we found that ~39% of nuclear genes recovered a sister clade consisting of the zokor and the bamboo rat, and ~34% and ~27% of nuclear genes supported two other topologies. chi-square tests showed that the distributions of ml and bayesian tree topologies were significantly deviated from random distribution, meaning that the majority of nuclear genes support a sister group of myospalacinae and rhizomyinae. by contrast, no significance was detected for either ml or bayesian trees inferred from individual mitochondrial genes, possibly because of small number of genes used. we also used a more stringent cutoff by counting the gene trees with boostrap values >70%. we found that 40% of nuclear genes supported a sister relationship between the zokor and the bamboo rat, which is significantly deviated from random distribution. in accordance with the individual gene analyses, the totality and the concatenated 4d-sites of the nuclear genes revealed a same topology as the majority of nuclear genes and received 100% bootstrap support in both analyses. although the concatenated protein-coding mitochondrial genes recovered the same phylogeny, the bootstrap values were relatively low. strikingly, we observed a very short branch leading to myospalacinae and rhizomyinae, suggesting that the time between the divergence of myospalacinae and rhizomyinae and the first speciation before the divergence of spalacidae is very short. the short period of time for the divergence of spalacidae indicates a burst of rapid speciation occurred. this observation was previously revealed by a phylogenetic analysis inferred from the nuclear gene irbp, although the monophyletic group comprising myospalacinae and rhizomyinae was weakly supported. discussion when the genome sequence is unavailable, transcriptome sequencing is an effective way to obtain large numbers of transcripts. in this study, we present the first transcriptome data in zokors and bamboo rats using massively parallel mrna sequencing. we perform illumina sequencing of the brain which is the most complex organ in the body and the liver which plays a major role in metabolism. although the raw data of the bamboo rat were much larger than that of the zokor, the results of sequence assembly, annotation, and classification of the two datasets were comparable. furthermore, we provide a conclusive resolution of phylogenetic placement of the three subfamilies in the family spalacidae. phylogenetic analyses using various data coding schemes and analytical approaches overwhelmingly support that the subfamily myospalacinae is more closely related to the subfamily rhizomyinae than to the subfamily spalacinae. although we did not examine multiple individuals of each species, the genetic variations within species should not affect our phylogenetic resolution at subfamily levels. traditional sanger sequencing technology produces longer est sequences with low throughput, which becomes little advantages when new assembly programs can handle large throughput short sequences. further, traditional transcriptome studies generally use a cdna library of one tissue, while transcriptome sequencing based on ngs methods can work on a normalized cdna library comprising multiple tissues and individuals, which would considerably reduce the sequencing cost. hence, the cost- and time-effective ngs technologies will provide more comprehensive transcriptome information and facilitate the transcriptome studies, particularly in less-studied species. although the large-scale transcriptome sequencing was used to sequence cdna pools derived from spalax galili, a member of the subfamily spalacinae, the transcriptome surveys on the members of the other two subfamilies within the family spalacidae are virtually lacking. our work provides a valuable genomic resource and will stimulate the further analysis of these taxa with exceptional scientific interest. for example, subterranean rodents are believed to have greater hypoxia tolerance than their aboveground relatives. upon examining candidate genes involving in adaptation to high-altitude hypoxia, we identified and genes in the zokor and bamboo rat transcriptomes, respectively. this finding will facilitate future comparative analysis of these functional genes. before the emergence of molecular identification approaches, one determines taxonomic status and phylogenetic relationships of animals mainly rely on morphological characteristics. however, the evolutionary change of morphological characters is extremely complicated, the phylogenetic trees derived from morphological data frequently remain controversial. in contrast, since the evolutionary change of dna follows a traceable pattern, it is possible to use a mathematical model to formulate the change and compare dna sequences among different organisms. as a result, molecular phylogenetics is expected to clarify the tree of life that has been difficult to be resolved by the classical approaches. taking zokors as an example, these animals have been allied to several different muroid subfamilies including rhizomyinae, spalacinae, arvicolinae, and cricetinae based on morphological characteristics, while overwhelming studies using molecular phylogenetic approaches concordantly placed zokors into spalacidae. it should be noticed that random noise will generally lead to poorly resolved phylogenetic trees, primarily because the number of nucleotide substitutions of a few genes is small. this may be the major reason why the previous molecular phylogeny studies showed inconsistence of topology relationships among the three subfamilies in spalacidae. in this study, we used a large number of genes including nuclear genes and all mitochondrial protein coding genes to examine the phylogenetic affinities of the three subfamilies within spalacidae. our analyses based on various data sets overwhelmingly support that the spalacinae had a highest probability to be a basal clade relative to others within spalacidae, while rhizomyinae and myospalacinae form a sister group. however, the most important challenges of phylogenomics involve different tree reconstruction methods and substitution saturation that can influence the validity of phylogenetic placements. to overcome the challenges, we used both ml and bi approaches based on both nuclear and mitochondrial genes, and found no incongruence between the two tree reconstruction approaches. further, we performed saturation tests and excluded the saturated 4d-site concatenation of mitochondrial genes in our analyses. interestingly, our proposed phylogeny of the three subfamilies within spalacidae is accordant with that inferred from massive morphological analyses on fossils of spalacidae. the known early miocene rhizomyines are closer to the stem zokor morphotype, suggesting that myospalacinae is more closely related to rhizomyinae than to spalacinae. although we didn’t examine african mole rats that were previously considered to be a separate subfamily in spalacidae, multiple lines of molecular evidence consistently supported a sister relationship between rhizomys and tachyoryctes, suggesting that tachyoryctes should be incorporated into the rhizomyinae. together, our phylogenomic analyses unambiguously resolve the phylogeny of the three subfamilies comprising the family spalacidae. geographically, members of the subfamily spalacinae are found in east europe, west asia, near east, and north africa, members of myospalacinae inhabit siberia and northern china, while the rhizomyinae include two geographically distant relatives that are from southeastern asia and east africa. our molecular phylogenetic evidence suggests that the common ancestor of spalacidae firstly split into spalacinae and rhizomyinae + myospalacinae clades in the northern asia adjacent to the middle east, and the latter clade subsequently split into rhizomyinae and myospalacinae in mongolia and northern china. interestingly, although rhizomyinae and myospalacinae are more phylogenetically related, myospalacinae appears to be more biologically similar with spalacinae than with rhizomyinae for certain traits. for example, myospalacinae and spalacinae are totally blinded, while rhizomyinae retain their sight. the resolution of phylogenetic relationships among the three subfamilies will provide a framework for evolutionary studies on convergence and divergence of these animals.  to summarize, this work is the first report of transcriptome sequencing in zokors and bamboo rats, representing a valuable resource for future studies of comparative genomics in subterranean mammals. phylogenomic analysis provides a conclusive resolution of interrelationships of the three subfamilies within the family spalacidae, and highlights the power of phylogenomic approach to dissect the evolutionary history of rapid radiations in the tree of life. 
all clean reads were assembled into and unigenes, which were annotated by the public databases: the swiss-prot, trembl, ncbi non-redundant protein, ncbi nucleotide sequence, gene ontology, cluster of orthologous groups, and kyoto encyclopedia of genes and genomes. phylogenomic analysis provides a conclusive resolution of interrelationships of the three subfamilies within the family spalacidae, and highlights the power of phylogenomic approach to dissect the evolutionary history of rapid radiations in the tree of life. phylogenetic analysis revealed a sister group relationship between the zokor and the bamboo rat, which is supported by the majority of gene trees inferred from individual orthologous genes, suggesting subfamily myospalacinae is more closely related to subfamily rhizomyinae. owing to the convergence of morphological traits and the incongruence of molecular evidence, the phylogenetic relationships among three subfamilies myospalacinae, spalacinae and rhizomyinae within the family spalacidae remain unresolved. here, we performed de novo transcriptome sequencing of four rna-seq libraries prepared from brain and liver tissues of a plateau zokor and a hoary bamboo rat, and analyzed the transcriptome sequences alongside a published transcriptome of the middle east blind mole rat. this is the first report of transcriptome sequencing in zokors and bamboo rats, representing a valuable resource for future studies of comparative genomics in subterranean mammals. the same topology was recovered from concatenated sequences of nuclear genes, fourfold degenerate sites of the nuclear genes and concatenated sequences of protein coding mitochondrial genes. subterranean mammals have been of great interest for evolutionary biologists because of their highly specialized traits for the life underground. approximately million clean reads from the zokor and million clean reads from the bamboo ratwere generated by illumina paired-end rna-seq technology. we characterize the transcriptome assemblies of the two spalacids, and recover the phylogeny of the three subfamilies using a phylogenomic approach.
 subterranean mammals have been of great interest for evolutionary biologists because of their highly specialized traits for the life underground. owing to the convergence of morphological traits and the incongruence of molecular evidence, the phylogenetic relationships among three subfamilies myospalacinae, spalacinae and rhizomyinae within the family spalacidae remain unresolved. here, we performed de novo transcriptome sequencing of four rna-seq libraries prepared from brain and liver tissues of a plateau zokor and a hoary bamboo rat, and analyzed the transcriptome sequences alongside a published transcriptome of the middle east blind mole rat. we characterize the transcriptome assemblies of the two spalacids, and recover the phylogeny of the three subfamilies using a phylogenomic approach.  approximately million clean reads from the zokor and million clean reads from the bamboo ratwere generated by illumina paired-end rna-seq technology. all clean reads were assembled into and unigenes, which were annotated by the public databases: the swiss-prot, trembl, ncbi non-redundant protein, ncbi nucleotide sequence, gene ontology, cluster of orthologous groups, and kyoto encyclopedia of genes and genomes. a total of nuclear orthologous genes were identified in the three spalacids and mouse, which was used as an outgroup. phylogenetic analysis revealed a sister group relationship between the zokor and the bamboo rat, which is supported by the majority of gene trees inferred from individual orthologous genes, suggesting subfamily myospalacinae is more closely related to subfamily rhizomyinae. the same topology was recovered from concatenated sequences of nuclear genes, fourfold degenerate sites of the nuclear genes and concatenated sequences of protein coding mitochondrial genes.  this is the first report of transcriptome sequencing in zokors and bamboo rats, representing a valuable resource for future studies of comparative genomics in subterranean mammals. phylogenomic analysis provides a conclusive resolution of interrelationships of the three subfamilies within the family spalacidae, and highlights the power of phylogenomic approach to dissect the evolutionary history of rapid radiations in the tree of life. spalacidaephylogenomicstranscriptomemitochondrial genomesubterranean rodents 