BACKGROUND
caecilians  together with frogs  and salamanders  constitute the three living orders of the class amphibia.
caecilians live in the tropical regions of southeast asia, africa, the seychelles islands and central and south america.
they primarily have a terrestrial, surface-cryptic or burrowing lifestyle as adults.
because of their secretive habits, caecilians are not frequently observed in the wild and are relatively unknown.
there are currently  <dig> recognized caecilian species, which is far fewer than the numbers of salamanders  and frogs   <cit> .
although the group is small, caecilians are intriguing to scientists because they have many traits that readily distinguish them from frogs and salamanders.
for example, consistent with their fossorial lifestyle, caecilians’ body is greatly elongated and segmented by annular grooves  <cit> .
they also lack sternal elements, girdles and limbs, and their tail is either very short or absent.
these characteristics make them look like snakes.
however, as a non-model cryptic organism, many aspects of caecilian’s biology are poorly investigated.
phylogenetically, caecilians are placed in a key position on the vertebrate tree.
they are the basal lineage of modern amphibians , and are therefore the most primitive group of land vertebrates.
the genomic information of caecilian will not only aid the biological studies on these animals but also will be important in addressing issues concerning the evolution of early land vertebrates, such as the water-to-land transition.
currently, however, only a sparse number of caecilian nucleotide sequences are available in public databases, and most of them are mitochondrial genes from phylogenetic studies.
obtaining more genomic information on caecilians is urgently important.
hox genes are usually among the first genes investigated in a non-model organism.
they are an old family of transcription factors characterized by a highly conserved 180 bp motif, the homeobox  <cit> .
they govern the timing and route of segmental development along the bilaterian animal body axis and are also involved in the patterning of the limbs and organogenesis .
in vertebrates, hox genes are organized into several clusters that exhibit colinearity between the gene order in the cluster and the temporal and spatial expression order of the hox genes during embryonic development: genes at the 3′ ends of the clusters are expressed earlier in development and more anterior than genes at the 5′ ends of the clusters  <cit> .
this colinearity is thought to be regulated by the many cis-regulator elements within and outside the hox clusters.
thus, both the hox gene coding sequences and the noncoding regions in the hox clusters are important in maintaining proper hox functions.
as a result of their significant involvement in early development, changes in expression of hox genes may cause severe morphological and/or physiological alterations, which are expected to play important roles in driving evolutionary changes.
hox genes have been extensively studied in a variety of animals, not only for their importance during development, but also for their genomic organization, such as number of genes, number of clusters and elements in the non-coding regions.
often, the genomic organization of hox clusters provides insights about the animal’s genome.
for example, a pattern from one single hox cluster in the cephalochordate amphioxus  <cit>  to four paralogous hox clusters  in tetrapods was one of the empirical findings that support the two round genome duplication hypothesis  <cit> .
the subsequent discovery of extra hox clusters  in teleost fishes was the first indication of an additional teleost-specific genome duplication  <cit> .
the comparatively slow evolution of the hox clusters in coelacanth was among some of the earliest molecular evidence supporting the slow evolution of its genome  <cit> .
in anole lizard, a massive enrichment of transposable elements at the hox loci was revealed, reflecting a genomic character that a large number of transposable elements accumulated in the development-related-gene-containing regions  <cit> .
thus, information on hox clusters is valuable, especially for animals of which the genome has not been sequenced.
for the cryptic caecilians, we and other labs have used comprehensive pcr surveys of the conserved homeodomains to make inferences on the number of hox genes  <cit> .
sequencing of the posterior fragment of the hoxd cluster  in the caecilian typhlonectes natans was recently achieved based on a bac library and revealed an expansion and multiple repeats in the intergenic region between hoxd <dig> and hoxd <dig>  <cit> .
to better characterize the caecilian hox clusters and better understand the caecilian genome, we aimed to sequence the four hox clusters of the banna caecilian.
with no access to a bac library, we used the results from our previous pcr survey and an alternative strategy that combined long-range pcr and genome walking.
we successfully obtained the majority of the hox cluster sequences, which greatly increased our knowledge of the caecilian hox clusters.
based on these findings, we determined most of the hox genes’ genomic linkages; screened for repetitive sequences, identified conserved non-coding elements using a comparative genomics approach, and measured the relative rates of evolution of the caecilian sequences in comparison with those of other sarcopterygians.
in our obtained caecilian hox cluster sequences, there were fewer repetitive sequences and more conserved noncoding elements than found in the frog homolog sequences, and the caecilian hox genes appear to be evolving relatively slowly.
this information from the caecilian hox clusters hinted at a slowly-evolved caecilian genome and we further tested it using a large dataset of  <dig> orthologous protein genes.
methods
sequencing and annotation of the banna caecilian hox clusters
this study was performed in strict accordance with the guidelines developed by the china council on animal care and use.
all animal processing procedures were approved by the institutional animal care and use committee of sun yat-sen university .
genomic dna was extracted from ethanol-preserved tissues  of banna caecilian  using the standard salt extraction protocol.
we amplified fragments of the hox clusters with a combination of la pcr and gw.
primers with high annealing temperature  that were suitable for la pcr and gw were designed in various anchors, including the banna caecilian hox coding regions and the highly conserved noncoding regions within the hox clusters.
the highly conserved noncoding regions were identified based on the multiz alignments of human, mouse, chicken, anole lizard, frog and zebrafish hox clusters constructed at the ucsc genome bioinformatics  <cit> .
la pcr was conducted to amplify the regions between every two adjacent anchors using takara la taq dna polymerase  according to standard protocol.
if there were weak or no bands on the agarose gel electrophoresis after the first-round of la pcr, nested pcr amplification was performed.
for regions that still failed to amplify, gw was conducted using a genome walking kit based on the tail-pcr technique  to shorten the distance between the two adjacent anchors.
the physical linkages between the gw acquired fragments and the anchors were further verified by pcr amplification with primers designed in the anchors and the distal ends of the gw acquired fragments.
after verification, the gw acquired fragments served as new anchors in which new primers were designed, and la-pcr was performed again with the new primers to amplify the regions.
finally, we obtained  <dig> la pcr fragments and a number of gw fragments.
fragments less than 5 kb long were sequenced on the abi3730xl sequencing platform.
longer fragments  were sequenced using ion torrent at life technologies corporation.
to prepare the sample, long pcr products  were pooled in proportion to their length and digested with neb next ds dna fragmentase according to the manufacturer’s instructions to generate a sequencing library with 300–600 bp dna fragments.
sequence reads obtained from ion pgm were de novo assembled using the mira program   <cit>  in its ‘genome’ assembly type with quality grades set as ‘accurate’.
together, the assembled contigs from mira and sequence fragments of  <dig> banna caecilian hox genes were reassembled using seqman in dnastar lasergene software  and manually edited.
the ambiguous contigs were verified by both blastn search against indonesian coelacanth or human hox clusters and further pcr surveys.
the remaining gaps were filled through pcr amplification and primer walking.
to annotate the newly obtained caecilian hox cluster sequences, we used genomescan  <cit>  and blastx search against indonesian coelacanth or human hox clusters to identify the hox coding sequences.
exon-intron boundaries were determined manually based on the consensus splice motifs.
screening of repetitive elements
identification and classification of repetitive elements was conducted using censor  <cit>  with the default parameters against the repbase library of vertebrate repeat sequences as well as using repeatmasker  <dig> .
<dig>  <cit> .
simple repeats and low-complexity sequences were excluded.
in addition to our newly obtained caecilian hox cluster sequences, genomic sequences of the four hox clusters in human , tammar wallaby , chicken , anole lizard , western clawed frog , indonesian coelacanth , spotted gar , and elephant shark  were extracted from the ensembl  <cit>  and ncbi nucleotide databases  <cit>  for comparison.
when comparing the repeat content among the human, wallaby, chicken, lizard, frog, caecilian and coelacanth hox clusters, only regions homologous to the obtained hox sequences in caecilian were analyzed, and regions homologous to the gaps in caecilian were excluded.
and “density of repetitive sequences”, which was calculated as the ratio of the total length of repetitive sequences to the length of the clusters in each species, was used.
when comparing the distribution of the repetitive sequences in caecilian and frog, all the available hox cluster sequences of the two amphibians were analyzed.
identification and analyses of cnes
hox clusters of human, tammar wallaby, chicken, anole lizard, western clawed frog, banna caecilian, indonesian coelacanth, spotted gar and elephant shark were aligned using the global alignment program lagan  <cit>  available on the vista website  <cit>  and screened for the presence of cnes.
because of the poor sequence coverage , the chicken hoxc cluster was not included in the analysis.
banna caecilian and the western clawed frog were used as reference sequence, respectively.
the cnes were restricted to the sequences with a cut-off of ≥65 % identity across windows ≥50 bp.
to determine the evolutionary origin of the cnes, all the cnes were classified into  <dig> phylogenetic groups, “gnathostome”, “osteichthyan”, “sarcopterygian” and “tetrapod”.
the “gnathostome” group was defined as cnes conserved in the elephant shark and caecilian  or frog .
the “osteichthyan” group referred to cnes shared by the spotted gar and caecilian or frog.
the “sarcopterygian” group referred to cnes found in coelacanth and caecilian or frog.
the remaining cnes were all classified in the “tetrapod” group.
and cnes already included in one phylogenetic group were not counted in the next group.
relative rate test of gene evolution
to determine the rate of evolution of caecilian relative to other species, we performed the tajima relative rate test  <cit>  on the hox family.
first, each gene-set was separately aligned using muscle  <dig>   <cit>  and ambiguous alignment regions were removed using gblocks   <cit>  with no gaps allowed.
sites with gaps or unknown data were excluded from the rrt.
each comparison included two ingroups and one outgroup.
the rrt on each gene-set was performed using in-house python scripts.
for most of the hox genes, the rrts were performed using banna caecilian, western clawed frog, puerto rican worm lizard  /chinese softshell turtle and human as the ingroups and indonesian coelacanth as the outgroup.
the rrts on hoxc <dig> were conducted with elephant shark as the outgroup and caecilian, african lungfish, eel and zebrafish as the ingroups because the western clawed frog, puerto rican worm lizard and human do not have hoxc <dig> gene.
for hoxc <dig>  the ingroups were changed to caecilian, frog, alpine stream salamander and african lungfish and the outgroup was indonesian coelacanth.
multiple rrts can be combined to determine the relative evolutionary speeds of several species.
a hasse diagram, in which the slower-evolving genes were placed below the faster-evolving ones, was used to show the combined rrt results.
the significance test corresponding to each rrt was drawn as a solid line  or a dotted line .
when comparing the rates between caecilian and coelacanth, elephant shark was used as the outgroup.
we also performed an rrt on a carefully curated dataset consisting of  <dig> orthologous genes, in which caecilian, frog and human were used as the ingroups and coelacanth was used as the outgroup.
for both hox family and the  <dig> orthologous genes, both nucleotide sequences and protein sequences were analyzed.
phylogenetic analysis
phylogenetic trees were constructed using hox cluster sequences and the  <dig> orthologous genes.
the four hox clusters from various gnathostomes were aligned with lagan.
and the  <dig> orthologous genes were aligned with muscle  <dig> .
the ambiguous positions in each alignment were eliminated by gblocks  with half gaps allowed.
maximum likelihood trees were constructed with raxml   <cit>  under the gtr + gamma + i model for nucleotide sequences and the lg + gamma model for protein sequences.
node support was estimated from  <dig> rapid bootstrap replicates.
RESULTS
hox clusters in banna caecilian
based on the hox gene fragments from our previous pcr survey, we used a strategy that combined long-range pcr  and genome walking  to determine the sequences of the four hox clusters in banna caecilian.
primers were designed in the coding regions of banna caecilian hox genes and the highly conserved noncoding regions in the sarcopterygian hox cluster alignments.
after assembly, all the physical linkages between the newly acquired fragments were determined, which were then verified by further pcr amplification.
in total, we retrieved  <dig>  kb sequences of banna caecilian hox clusters, which included the nearly complete hoxa cluster , most of the hoxb cluster , part of the hoxc cluster  and most of the hoxd cluster  .
the ncbi accession numbers of the four hox clusters are from kf <dig> to kf <dig>  we were unable to sequence seven intergenic regions most likely because of the limitations of our strategy.
they were mainly located in the hoxc cluster and the extremities of the hoxa, hoxb and hoxd clusters.
four of these regions  were putatively very large as estimated by previous studies in other animals  <cit> , and therefore, beyond the scope of long range pcr.
the remaining three gaps  may have been due to either enlarged gene distances or complex structures that hindered the progression of the taq dna polymerase.fig.
<dig> genomic organization of the hox clusters in banna caecilian and western clawed frog.
schematic representations of the hoxa, hoxb, hoxc and hoxd clusters in caecilian are shown, with the frog hox clusters drawn for comparison.
for each hox cluster, location of the genes, repetitive sequences and conserved noncoding elements are separately illustrated from top to bottom.
total numbers of the repetitive sequences and conserved noncoding elements are indicated at the end of each cluster.
exons of the genes are characterized by black bars and pseudogenes are denoted by the symbol Ψ. gaps in the intergenic regions are indicated by dotted lines.
repetitive sequences identified by censor are represented by grey bars.
conserved noncoding sequences predicted using mvista are classified into  <dig> phylogenetic groups : “gnathostome ”, “osteichthyan ”, “sarcopterygian ” and “tetrapod ”.
for clarity, cnes conserved in caecilian and frog only are not shown in the figure
estimated from the sequences we obtained, the sizes of the caecilian hox a, hoxb and hoxd clusters are comparable to their counterparts in another amphibian, the western clawed frog.
in both species, the hoxa , hoxb  and hoxd clusters  are ~ <dig>  ~ <dig> and ~95 kb in size, respectively.
for the less complete hoxc cluster, three fragments were compared.
the hoxc11-hoxc <dig> and hoxc6-hoxc <dig> fragments in the two species were ~40 kb in size; while the hoxc13-hoxc <dig> fragment in caecilian  was approximately only one-third of the size of its counterpart in frog .
for the hoxa, hoxb and hoxd clusters, gene content and gene organization were largely determined.
for hoxc cluster, gene linkages were determined except for hoxc12-hoxc <dig>  hoxc8-hoxc <dig>  hoxc4-hoxc <dig> and hoxc3-hoxc <dig>  the gene content and organization of the caecilian hox clusters were similar to those in other amphibians, with slight differences.
caecilian retained hoxc <dig>  which was lost in many other vertebrates, including frog and salamander, but detected in coelacanth and lungfish.
the complete coding sequence of hoxc <dig> that we obtained demonstrates that the caecilian hoxc <dig> is fully functional.
hoxd <dig> was completely absent in frog and salamander, but remnants of the caecilian hoxd <dig> was observed in region 5′ to hoxd <dig> after performing a global alignment with other gnathostomes .
in our caecilian hox cluster sequences, we found five microrna genes, two belong to family mir- <dig> and three belong to family mir- <dig>  at the same genomic location observed in other vertebrates.
mirna 196a- <dig> which is usually located approximately 6 kb 5′ to hoxb <dig> in other vertebrates was not identified even though  <dig>  kb of sequence 5′ to the caecilian hoxb <dig> was screened .
we could not detect mirna 196a- <dig> in the corresponding genomic region in western claw frog, either.
thus, the mirna 196a- <dig> microrna in amphibians may be either lost or located further away from hoxb <dig>
in summary, we found in banna caecilian  <dig> hox genes  and  <dig> pseudogene  organized into four clusters,  <dig> evx paralogs associated with the hoxa  and hoxd clusters, respectively, and  <dig> micrornas .
distribution of repetitive sequences
obtaining the majority of the caecilian hox clusters allowed us to investigate the details of the hox cluster structure and elements in the intergenic regions.
first, we evaluated the repetitive sequences which are usually strongly excluded from vertebrate hox clusters.
we used censor and the repbase database to detect repetitive sequences in the obtained banna caecilian hox clusters and the corresponding regions in the human, wallaby, chicken, lizard, frog and coelacanth hox clusters.
the density of repetitive sequences  in the caecilian hox clusters was  <dig>  %, among the lowest densities of all the species studied .
it is possible that the density of repetitive sequences in banna caecilian may have been underestimated without a specific repeat library available.
wallaby, lizard and coelacanth, however, do not have specific repeat libraries either, and they did not show a general bias towards low densities of repetitive sequences; and all of these species had a higher density of repetitive sequences than banna caecilian.
thus, the effect of databases in identifying repeats does not appear to be predominant and the low density of repetitive sequences was considered primarily as a characteristic of the caecilian hox cluster sequences that we obtained.
when the two amphibians were compared, the repeat content of the hox clusters in caecilian was approximately only a quarter of that in frog .
to avoid any algorithm bias, repeatmasker was also used to detect repetitive sequences.
the density of repetitive sequences in the caecilian hox clusters  was again among the lowest densities and much lower than that of frog .
because the general trend was consistent in the two analyses, only the results of censor are shown for the following comparisons.fig.
<dig> densities of repetitive sequences of the gnathostome hox clusters.
density of repetitive sequences for a specific hox cluster was calculated by dividing the total length of the repetitive sequences by the total length of the cluster.
only repetitive sequences identified by censor with default settings were calculated.
the chicken hoxc cluster was not included in the analysis because it is currently missing many sequences
the attributes and distribution of the repeats detected in the caecilian and frog hox cluster sequences  were further investigated .
note that because all the unsequenced gaps might contain repeats, the number of repeats in each caecilian hox cluster was likely underestimated.
the most abundant type of repeats in caecilian was non-ltr retrotransposons , while in frog approximately  <dig>  % of the repeats were dna transposons, such as hats, harbingers and koloboks.
among the four hox clusters, the caecilian hoxa cluster contained the least number of repetitive sequences  interspersed in the cluster; the frog hoxa cluster contained  <dig> repeats, 59 % of which were concentrated in the evx1-hoxa <dig> region.
in the caecilian hoxb cluster,  <dig> repeats were found, and over half of them were located between hoxb <dig> and hoxb1; while in the frog hoxb cluster,  <dig> repeats were identified and a large portion of them were located in regions 5′ to hoxb <dig>  hoxb8-hoxb <dig> and hoxb2-hoxb <dig>  for both species, the hoxc cluster contained the most abundant amount of repetitive sequences, with  <dig> and  <dig> repeats in caecilian and frog, respectively.
while the caecilian repeats were largely interspersed in the cluster,  <dig>  % of the frog repeats were concentrated in the two terminal regions, i.e., upstream of hoxc <dig> and downstream of hoxc <dig>  for the hoxd cluster,  <dig>  % of the  <dig> repeats in caecilian were located in the hoxd13-hoxd <dig> and hoxd8-hoxd <dig> regions, while approximately 50 % of the  <dig> repeats in frog were in the hoxd3-hoxd <dig> region.
thus, the distribution patterns of the repeats had “hotspots” that were in different locations in the two amphibians.
and most of these “hotspots” were enlarged compared with the orthologous regions in other vertebrates.
two of the gapped intergenic regions in caecilian, hoxb2-hoxb <dig> and hoxd13-hoxd <dig>  were among the repeat “hotspots”.
they might be also enlarged, thus were hard to amplify.
to identify direct and inverted repeats that were not in the known repeat libraries, self-self blastn analyses  were performed.
a direct/inverted repeat consists of two repeat copies  that are approximately identical/complementary to each other.
this analysis revealed the presence of  <dig> inverted repeats and  <dig> direct repeats in the obtained caecilian hox clusters, over half of which were identified in the hoxd cluster.
in particular, we found two long inverted repeats in the hoxd8-hoxd <dig> intergenic region.
the arms of these two repeats were both longer than 100 bp , and the arm identity was almost 100 %.
the spacers between the two arms of both repeats were small .
in frog,  <dig> inverted repeats and  <dig> direct repeats were found in the self-self blastn analysis in corresponding regions of the hox clusters, but no long inverted repeat was identified.
indeed, inverted repeats with an arm size longer than 100 bp are rare in genomes.
there are only  <dig> such repeats found in the human genome  <cit> , and currently, no long inverted repeat has ever been reported in the hox clusters of other gnathostomes species.
determining whether there is any functional significance of these two long inverted repeats in the caecilian hox cluster is worthy of further investigation.
conserved noncoding elements in caecilian hox clusters
next, we screened the caecilian hox clusters for conserved noncoding sequences, which may be indicative of potential regulatory elements.
multiple global alignments of the four caecilian hox clusters and the homologous regions of other gnathostomes  were conducted using lagan.
from these alignments, conserved noncoding elements were predicted with vista.
this method has been shown to be effective at identifying and visualizing overtly conserved non-coding elements  <cit> .
because of the incompleteness of the caecilian hoxc cluster, the number of cnes detected in the caecilian hoxc clusters may be an underestimation.
furthermore, the chicken hoxc cluster was not included in this analysis because it is currently missing many sequences .
the vista plots of the four hox clusters are shown in additional file  <dig>
the comparisons between caecilian and other gnathostomes identified  <dig> cnes in intergenic and some intronic regions of the caecilian hox clusters.
the distribution of these cnes is shown in fig.
<dig> and table  <dig>  of the four clusters, the hoxa cluster contained the highest number of cnes, and the hoxc cluster had the highest number of cnes in introns.
the cne densities in the 3′ part of the caecilian hoxa, hoxb and hoxd clusters  tended to be higher than those in the 5′ part, consistent with the previous observations in many other gnathostomes .
this density difference was not readily observed in the caecilian hoxc cluster because of the incompleteness especially in its 3′ part.
overall, cnes tended to be located around the coding regions of the hox genes, and the regions containing more repeats  tended to have lower cne densities.table  <dig> cnes in the four hox clusters of banna caecilian
for comparison, cnes in the homologous regions of the frog hox clusters were also identified using a similar method, and  <dig> cnes were found, which was fewer than those identified in caecilian.
comparing the cnes in the two amphibians,  <dig> of the caecilian cnes  corresponded to  <dig> of the frog cnes  as several long cnes in caecilian were broken into shorter cnes in frog.
in addition to the cnes common to caecilian and frog,  <dig> cnes were identified by comparing caecilian with other non-frog gnathostomes , but only  <dig> cnes were identified by comparing frog with other non-caecilian gnathostomes .
we next classified the cnes in caecilian and frog into  <dig> phylogenetic groups, “gnathostome”, “osteichthyan”, “sarcopterygian” and “tetrapod”.
the “gnathostome” group was defined as cnes at the same positions in the hox clusters of elephant shark and caecilian or frog.
the “osteichthyan” group referred to cnes shared by spotted gar and caecilian or frog.
the “sarcopterygian” group referred to cnes found in both indonesian coelacanth and caecilian or frog.
the remaining cnes were all placed in the “tetrapod” group.
note that the cnes already included in a previous phylogenetic group were not counted in the next groups.
as shown in table  <dig>  a large portion of the cnes found in both amphibians was categorized into the “tetrapod” and “gnathostome” groups.
most of the cnes in the “gnathostome” group were conserved in almost all of the species studied and most of the cnes in the “tetrapod” group were also found in one or two other species .
because similarities in sequences from highly divergent organisms often imply functional constraints, the cnes we identified, especially those in the “gnathostome” group, were most likely indicative of potential important functional elements.
we analyzed the caecilian “gnathostome” cnes by blast searching the corresponding human sequences against the ncbi database of human ests,  <dig> of which  had est matches  and were therefore putative non-coding rna genes.
the remaining “gnathostome” cnes likely contained many potential cis-regulatory elements;  <dig> of them were on the list of functionally-verified hox enhancers summarized in a study by ravi et al.
<cit> .
furthermore, cnes in the other three groups, especially those present in multiple organisms, may also be putative non-coding rnas or cis-regulatory elements and deserve further functional analysis.table  <dig> comparing the caecilian or frog cnes in the four phylogenetic groups
differences were observed between caecilian and frog in the retention of cnes at different evolutionary stages.
among them, the most prominent were in the “sarcopterygian” group in which there was  <dig>  fold greater difference in the cnes in caecilian compared with frog  and this difference increased to more than  <dig> -fold if the cnes common to caecilian and frog were excluded .
this difference was mostly observed in the hoxa, hoxb and hoxc clusters, especially in the hoxb cluster where caecilian had  <dig> “sarcopterygian” cnes while frog had none.
this variation in the “sarcopterygian” group was the primary reason for the different number of cnes found in caecilian and frog.
slow evolutionary rates of caecilian genes
the evolutionary rates of the banna caecilian hox genes were estimated using the tajima’s relative rate test  in which banna caecilian, western clawed frog, puerto rican worm lizard/chinese softshell turtle and human were first used as the ingroups and indonesian coelacanth as the outgroup.
because of the missing hoxc <dig> and hoxc <dig> in certain tetrapods,  <dig> of the  <dig> hox genes were analyzed.
as shown in fig.
<dig>  over 62 % of the analyzed hox coding sequences showed significantly  slower evolutionary rates in caecilian than in frog, whereas the reverse was true for no hox genes.
nearly 84 % of the analyzed hox coding sequences were identified evolving significantly  more slowly in caecilian than in the mammal representative, human, whereas the reverse was true for no hox genes.
because anole lizard’s hox genes generally evolved fast, reptiles were represented by slow worm  and chinese softshell turtle .
over 45 % of the analyzed hox genes revealed significantly slower evolutionary rates in caecilian than in the reptile representatives.
thus, most of the hox genes analyzed evolved significantly more slowly in banna caecilian than in other species,  <dig> of which evolved the most slowly in banna caecilian compared with all other species studied.
because indonesian coelacanth hox clusters have been reported to evolve comparatively slowly, we next compared the relative evolutionary rates of caecilian and coelacanth using the elephant shark as the outgroup.
all of the hox genes showed no significant difference in evolutionary rate, but  <dig> appeared to evolve faster in caecilian than in coelacanth.
in summary, the hox gene sequences evolved much more slowly in caecilian than in other tetrapod species, except for a few genes located in the hoxc cluster, and approximately 81 % of the hox genes demonstrated similar evolutionary rates to those in the slowly evolving indonesian coelacanth.
rrts were also performed using the hox protein sequences  to avoid inaccurate estimation due to the issue of saturation when using nucleotide sequences.
twenty one of the caecilian hox proteins showed slower evolutionary rates than at least one other species, while only  <dig> frog hox proteins were slower than at least one other species.
thirteen caecilian hox proteins evolved more slowly than the frog ones, but only one frog hox protein evolved more slowly than its caecilian counterpart.
thus, though the evolutionary rate differences of the hox protein sequences were not as significant as the nucleotide sequences, caecilian still appears a slow evolving species.fig.
<dig> rrts of the gnathostome hox gene sequences.
a evolutionary relationship of the gnathostome species used in the rrts of the hox genes.
b summary of the rrts conducted on the hox gene nucleotide sequences.
for most of the hox genes, rrts were performed using banna caecilian, western clawed frog, puerto rican worm lizard/chinese softshell turtle and human as the ingroups and indonesian coelacanth as the outgroup.
rrts of hoxc <dig> were carried out using elephant shark as the outgroup and caecilian, african lungfish, eel and zebrafish as the ingroups to replace the species that do not have hoxc <dig> gene.
for hoxc <dig>  the ingroups included caecilian, frog, alpine stream salamander and african lungfish and the outgroup was indonesian coelacanth.
results of rrts for each gene are shown in a hasse diagram, in which the slower-evolving genes are placed below the faster-evolving ones, with statistical significance denoted as a solid line  or a dotted line .
c summary of the significant rrts among banna caecilian, western clawed frog, puerto rican worm lizard and human.
for each pair of species, the significant rrts are denoted by arrows pointing to the slower-evolving one.
the number of highly significant  tests for each hox cluster is indicated on the side of the arrow
additionally, phylogenetic trees of the four hox clusters were constructed using all of the alignable sequences .
branch lengths reflect the differences in evolutionary rates among and within lineages with longer branch lengths representing faster evolutionary rates.
as observed in fig.
<dig>  in the hoxa, hoxb and hoxd cluster trees branch lengths of frog were  <dig> -fold ,  <dig> -fold  and 2-fold  longer than those of caecilian, respectively.
in the hoxc cluster tree, the branch length of frog was also longer than that of caecilian, though this difference appeared to be smaller than the differences observed in the other three clusters.
all of the phylogenetic trees further demonstrate that the caecilian hox clusters evolved more slowly than the frog hox clusters.fig.
<dig> phylogenetic trees of the gnathostome hox gene clusters.
all the alignable sequences  for each hox cluster were used.
the trees were constructed using raxml under the gtr + gamma + i model .
elephant shark was used as the outgroup.
values above branches denote maximum likelihood bootstrap support
next, we wondered whether the slow evolution of the caecilian hox clusters reflected a genome-wide phenomenon, as in the case of coelacanth  <cit> .
to test this idea, we performed an rrt analysis with both nucleotide and protein sequences using  <dig> orthologous genes from human, frog and caecilian  .
for the nucleotide dataset, as many as  <dig>  % genes in caecilian were significantly slower than the orthologs from either frog or human.
approximately  <dig>  % caecilian genes showed the slowest evolutionary rate, i.e., slower than both frog and human.
in frog, the percentage of genes slower than their caecilian and/or human orthologs was only  <dig>  %.
the rrt result of the protein sequences showed a similar trend .
up to  <dig>  % genes in caecilian were significantly slower than the orthologs from either frog or human.
in frog, this percentage was  <dig>  %.
the percentages of genes showed the slowest evolutionary rate in caecilian and frog were  <dig>  and  <dig>  %, respectively.
additionally, we constructed a phylogenetic tree with the  <dig> concatenated orthologous genes from  <dig> vertebrate species .
both nucleotide sequences and protein sequences were analyzed.
in both trees, the branch lengths of caecilian were short which suggested that it evolved slowly.
among the three distinct living amphibian lineages, the evolutionary rates of frog and axolotl were both greater than that of caecilian.
thus, not only the hox genes, but also a large portion of the protein-coding genes are evolving significantly more slowly in caecilian.fig.
<dig> genome-wide rrts of banna caecilian, western clawed frog and human using nucleotide sequences.
rrts were performed on the nucleotide sequences of  <dig> orthologous genes, with banna caecilian, western clawed frog and human as the ingroups and indonesian coelacanth as the outgroup.
the number in a cell represents the number of statistical significance  for rrts of a species compared with the others; hence, there are six possibilities for each ingroup species.
increasing warm color intensity indicates faster relative evolutionary rate, whereas increasing cold color intensity indicates slower relative evolutionary rate
discussion
in this study, we successfully obtained most of the banna caecilian hox cluster sequences using la pcr and gw.
we identified  <dig> hox genes,  <dig> evx paralogs,  <dig> microrna genes and  <dig> pseudogene  in four clusters.
most of their physical linkages were determined which was not accomplished in previous pcr surveys.
the overall organization of the caecilian hox clusters was conserved, similar to those from other species with four clusters.
the presence of hoxc <dig> in caecilian is unique compared with all the tetrapods studied so far, which suggests that the tetrapod ancestor had hoxc <dig> but lost it in different lineages.
loss or pseudogenization of hoxd <dig> appears to be a distinct feature of the amphibians’ hox clusters and has been proposed to be related to digit reduction as frogs and salamanders normally only have four fingers and caecilians are limbless  <cit> .
other surveys have shown that hoxd <dig> is also lost in snakes, which are limbless, and the african lungfish, which has only thread-like fins, but it is present in many limbless lizards  <cit> .
therefore, there may be different mechanisms responsible for digit reduction and hoxd <dig> may not be directly involved in all of them.
alternatively, it could be the insertion of repetitive sequences that first disrupted the regulatory network of the posterior hoxd genes and altered the developmental program for fingers and even limbs.
these effects of repeat insertion may then account for the pseudogenization of hoxd <dig>  it would be interesting to characterize the posterior hoxd cluster in limbless lizards and see whether there are overtly multiple repeats in the region though the hox genes remain intact.
from the obtained sequences, the caecilian hox clusters appear compact and have fewer repetitive sequences than the hox clusters in the other organisms used in this study.
though both are amphibians, caecilian and frog have much different repeat contents, which may reflect the different evolutionary constraints on the hox clusters of the two amphibians.
there were, however, gaps in several intergenic regions that we could not fill in caecilian.
based on the hox cluster information from other amphibians, one possible reason is that these regions have been greatly enlarged, such as the expanded hoxd11-hoxd <dig> in the caecilian typhlonectes natans  <cit>  and the red spotted newt  <cit>  and the expanded hoxc13-hoxc <dig>  hoxc5-hoxc <dig> and hoxd3-hoxd <dig> in the western claw frog  <cit> .
repetitive sequences are likely abundant in these possibly expanded regions in banna caecilian as they have accumulated in the known large intergenic regions in the hox clusters of many organisms  <cit> .
two of them  have actually shown signs of repeat accumulation .
thus, these gapped intergenic regions were possibly repeat “hotspots” and the amount of repetitive sequences in the banna caecilian hox clusters may currently be underestimated.
in the future, obtaining the full sequences of these gaps and determining the structural features will elucidate whether they are indeed expanded and have multiple repeat insertions.
if this were the case, it would also help to understand why repeat hotspots exist.
expansion of intergenic regions may affect gene regulation by long-range regulatory elements.
insertion of repetitive dna has the potential to alter local regulatory structures by undermining the conserved motifs and/or creating new ones, further influencing gene expression  <cit> .
alternatively, the presence of repetitive dna may just indicate weakened coherence or relaxed constraint within the clusters.
why repeats are enriched in certain intergenic regions and whether they are related to developmental variations of caecilians or just a reflection of the relaxation of evolutionary constraints remains unknown.
vista plots demonstrated a high level of conservation in the non-protein-coding regions of the caecilian hox clusters.
in contrast to repetitive sequences, cnes were abundant in the obtained caecilian hox cluster sequences.
and cnes tended to be located in the regions containing few repeats.
this finding may reflect the preservation of regulatory modules that formed the complex regulatory network of the hox genes.
the phylogenetic group classifications suggested that the cnes in the hox clusters had different conservation depths.
the ancient  and well conserved  cnes may be fundamental for the regulatory network and proper hox functions.
other cnes may then be gradually added to or removed from the regulatory network at different evolutionary stages and in different lineages, thus complicating and rewiring the regulatory network in different species.
comparative analyses of the hox cluster sequences from more vertebrate species will improve the accuracy of this classification and elucidate the functions of the cne.
that caecilian hox clusters contained more cnes than frog, with the major difference in the “sarcopterygian” group, implies that the frog hox cluster sequences evolved faster than the caecilian sequences and that the caecilian hox cluster sequences retained more characteristics of the early tetrapods than the frog sequences, such as hoxc <dig> and the hoxd <dig> remnant.
the slow evolution of the caecilian hox clusters was supported by the rrts performed for the hox protein sequences and the phylogenetic tree analyses of all the alignable hox cluster sequences.
the majority of the caecilian hox genes had evolutionary rates similar to the indonesian coelacanth hox genes which were known to have evolved slowly.
as in the case of coelacanth, the slow evolution of the caecilian hox clusters may be indicative of its slowly evolving genome.
it was supported by the rrts and the phylogenetic tree analysis of  <dig> protein homologs which showed a general pattern of slow evolutionary rate of the caecilian sequences.
thus from the evolutionary rate of the hox cluster and the  <dig> protein dataset, it was quite likely that the caecilian genome was evolving slowly.
however, this still require support from more genomic information.
for this rarely found organism, deciphering more of its genomic information and studying it from the inside out would greatly complement the studies on its morphology, physiology and function.
furthermore, to better understand the genome biology of living amphibians, a caecilian genome project is especially of practical importance, because salamanders possess genome sizes that are intractably large for routine genomic analyses  <cit> .
and as the most basal tetrapod lineage, the caecilian genome is also highly favorable for comparative genomic studies.
for example, comparative analyses between the caecilian and coelacanth genomes would provide important clues for reconstructing genomic perspectives of ancestral tetrapods and understanding the water-to-land transition.
currently, there are fewer than  <dig> nucleotide sequences for caecilian in genbank, and most of these sequences are mitochondrial gene markers from phylogenetic studies.
our study sequencing the majority of the four caecilian hox clusters, is an important data resource and provides an amount of genomic information for the cryptic caecilian.
in the future, caecilian shall be added to the list of organisms in demand of genome projects, and characterizing the caecilian genome will not only help to understand its biology, but also will provide insights into the genome biology and evolution of early tetrapods.
CONCLUSIONS
in this study, we have cloned and sequenced most of the banna caecilian hox cluster.
thirty nine hox genes,  <dig> microrna genes and  <dig> pseudogene  were identified and most of their physical linkages in the four clusters were determined.
the presence of hoxc <dig> in caecilian suggests that the tetrapod ancestor had hoxc <dig> but then lost it in different lineages.
loss of hoxd <dig> function appears to be a feature of the amphibians hox clusters.
whether and how it was related to digit reduction requires further investigation.
from the obtained sequences, the caecilian hox clusters appear compact and have fewer repetitive sequences than the hox clusters in the other organisms studied in this work.
however, to fully understand the repetitive sequences in the caecilian hox clusters, the intergenic gaps are to be sequenced and their structural features be determined in future studies.
the analysis on cnes, the rrts on the hox genes and the phylogenetic trees constructed from the hox clusters all revealed stronger constraints on the caecilian hox clusters than their frog counterparts.
these information may be indicative of a slowly evolving genome, which was supported by rrt analyses and the phylogenetic tree construction using a large orthologous protein dataset.
therefore, characterizing more of the caecilian genomic information is in demand, which will not only help to understand the basic biology of this clade, but also be highly valuable for the comparative genomics of amphibians and the evolutionary studies of early tetrapods.
additional files
additional file 1:
global genomic alignments of the posterior hoxd cluster .
exons  and transcription direction of the genes are indicated.
sequence comparisons were conducted using the human sequence as the reference.
nucleotide identities relative to the human sequence are given by histogram peaks.
cnes are depicted by red peaks.
additional file 2:
global genomic alignments of part of the hoxb cluster .
exons  and transcription direction of the genes are indicated.
sequence comparisons were conducted using the human sequence as the reference.
nucleotide identities relative to the human sequence are given by histogram peaks.
cnes are depicted by red peaks.
additional file 3:
global genomic alignments of the four hox clusters.
exons  and transcription direction of the genes are indicated.
sequence comparisons were conducted using the caecilian hox clusters as the reference sequence.
nucleotide identities relative to the caecilian sequences are given by histogram peaks.
cnes are depicted by red peaks.
additional file 4:
conservation of the cnes in banna caecilian hox clusters belonging to the  <dig> phylogenetic groups.
nine species were used for the global genomic alignments of the hoxa, hoxb and hoxd clusters; for the hoxc cluster, only  <dig> species were used due to the low sequence coverage of the chicken hoxc cluster.
numbers indicate the count of species that retained an identical cne.
from the “tetrapod” group to the “gnathostome” group, the color changes from light red to dark red gradually.
additional file 5:
rrts of the gnathostome hox gene using protein sequences.
evolutionary relationship of the gnathostome species used in the rrts of the hox protein sequences.
summary of the rrts conducted on the hox gene protein sequences.
for most of the hox genes, rrts were performed using banna caecilian, western clawed frog, puerto rican worm lizard/chinese softshell turtle and human as the ingroups and indonesian coelacanth as the outgroup.
rrts of hoxc <dig> were carried out using elephant shark as the outgroup and caecilian, african lungfish, eel and zebrafish as the ingroups to replace the species that do not have hoxc <dig> gene.
for hoxc <dig>  the ingroups included caecilian, frog, alpine stream salamander and african lungfish and the outgroup was indonesian coelacanth.
results of rrts for each gene are shown in a hasse diagram, in which the slower-evolving genes are placed below the faster-evolving ones, with statistical significance denoted as a solid line  or a dotted line .
summary of the significant rrts among banna caecilian, western clawed frog, puerto rican worm lizard and human.
for each pair of species, the significant rrts are denoted by arrows pointing to the slower-evolving one.
the number of highly significant  tests for each hox cluster is indicated on the side of the arrow.
additional file 6:
genome-wide rrts of banna caecilian, western clawed frog and human using protein sequences.
rrts were performed on the protein sequences of  <dig> orthologous genes, with banna caecilian, western clawed frog and human as the ingroups and indonesian coelacanth as the outgroup.
the number in a cell represents the number of statistical significance  for rrts of a species compared with the others; hence, there are six possibilities for each ingroup species.
increasing warm color intensity indicates faster relative evolutionary rate, whereas increasing cold color intensity indicates slower relative evolutionary rate.
additional file 7:
phylogenetic tree inferred from a dataset of  <dig> nuclear protein-coding genes.
nucleotide tree .
protein tree .
putative orthologous genes in two species were identified using the mutual best hit  in basic local alignment search tool .
the tree was constructed using the maximum likelihood method with raxml under the gtr + gamma + i model for nucleotide sequences and the lg + gamma model for protein sequence .
elephant shark was used as the outgroup.
values above branches denote maximum likelihood bootstrap support.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
dl, pz and rw conceived the study.
sm provided the species sample.
rw and ql designed and carried out the experiments.
rw and ql analyzed the data.
dl, rw and ql wrote the manuscript.
all authors read and approved the final manuscript.
we thank rongli mao for helping in experiments and all the lab members for insightful discussions on the manuscript.
this work was supported by the national natural science foundation of china grants  to pz and sm, respectively, the china national fund for excellent young scientists  and the cultivation project for young teachers of sysu  to pz and the cultivation project of the department of education of guangdong province  to dl.

 caecilians together with frogs and salamanders constitute the three living orders of the class amphibia. caecilians live in the tropical regions of southeast asia, africa, the seychelles islands and central and south america. they primarily have a terrestrial, surface-cryptic or burrowing lifestyle as adults. because of their secretive habits, caecilians are not frequently observed in the wild and are relatively unknown. there are currently recognized caecilian species, which is far fewer than the numbers of salamanders and frogs. although the group is small, caecilians are intriguing to scientists because they have many traits that readily distinguish them from frogs and salamanders. for example, consistent with their fossorial lifestyle, caecilians’ body is greatly elongated and segmented by annular grooves. they also lack sternal elements, girdles and limbs, and their tail is either very short or absent. these characteristics make them look like snakes. however, as a non-model cryptic organism, many aspects of caecilian’s biology are poorly investigated. phylogenetically, caecilians are placed in a key position on the vertebrate tree. they are the basal lineage of modern amphibians, and are therefore the most primitive group of land vertebrates. the genomic information of caecilian will not only aid the biological studies on these animals but also will be important in addressing issues concerning the evolution of early land vertebrates, such as the water-to-land transition. currently, however, only a sparse number of caecilian nucleotide sequences are available in public databases, and most of them are mitochondrial genes from phylogenetic studies. obtaining more genomic information on caecilians is urgently important. hox genes are usually among the first genes investigated in a non-model organism. they are an old family of transcription factors characterized by a highly conserved 180 bp motif, the homeobox. they govern the timing and route of segmental development along the bilaterian animal body axis and are also involved in the patterning of the limbs and organogenesis. in vertebrates, hox genes are organized into several clusters that exhibit colinearity between the gene order in the cluster and the temporal and spatial expression order of the hox genes during embryonic development: genes at the 3′ ends of the clusters are expressed earlier in development and more anterior than genes at the 5′ ends of the clusters. this colinearity is thought to be regulated by the many cis-regulator elements within and outside the hox clusters. thus, both the hox gene coding sequences and the noncoding regions in the hox clusters are important in maintaining proper hox functions. as a result of their significant involvement in early development, changes in expression of hox genes may cause severe morphological and/or physiological alterations, which are expected to play important roles in driving evolutionary changes. hox genes have been extensively studied in a variety of animals, not only for their importance during development, but also for their genomic organization, such as number of genes, number of clusters and elements in the non-coding regions. often, the genomic organization of hox clusters provides insights about the animal’s genome. for example, a pattern from one single hox cluster in the cephalochordate amphioxus to four paralogous hox clusters in tetrapods was one of the empirical findings that support the two round genome duplication hypothesis. the subsequent discovery of extra hox clusters in teleost fishes was the first indication of an additional teleost-specific genome duplication. the comparatively slow evolution of the hox clusters in coelacanth was among some of the earliest molecular evidence supporting the slow evolution of its genome. in anole lizard, a massive enrichment of transposable elements at the hox loci was revealed, reflecting a genomic character that a large number of transposable elements accumulated in the development-related-gene-containing regions. thus, information on hox clusters is valuable, especially for animals of which the genome has not been sequenced. for the cryptic caecilians, we and other labs have used comprehensive pcr surveys of the conserved homeodomains to make inferences on the number of hox genes. sequencing of the posterior fragment of the hoxd cluster in the caecilian typhlonectes natans was recently achieved based on a bac library and revealed an expansion and multiple repeats in the intergenic region between hoxd and hoxd. to better characterize the caecilian hox clusters and better understand the caecilian genome, we aimed to sequence the four hox clusters of the banna caecilian. with no access to a bac library, we used the results from our previous pcr survey and an alternative strategy that combined long-range pcr and genome walking. we successfully obtained the majority of the hox cluster sequences, which greatly increased our knowledge of the caecilian hox clusters. based on these findings, we determined most of the hox genes’ genomic linkages; screened for repetitive sequences, identified conserved non-coding elements using a comparative genomics approach, and measured the relative rates of evolution of the caecilian sequences in comparison with those of other sarcopterygians. in our obtained caecilian hox cluster sequences, there were fewer repetitive sequences and more conserved noncoding elements than found in the frog homolog sequences, and the caecilian hox genes appear to be evolving relatively slowly. this information from the caecilian hox clusters hinted at a slowly-evolved caecilian genome and we further tested it using a large dataset of orthologous protein genes. methods sequencing and annotation of the banna caecilian hox clusters this study was performed in strict accordance with the guidelines developed by the china council on animal care and use. all animal processing procedures were approved by the institutional animal care and use committee of sun yat-sen university. genomic dna was extracted from ethanol-preserved tissues of banna caecilian using the standard salt extraction protocol. we amplified fragments of the hox clusters with a combination of la pcr and gw. primers with high annealing temperature that were suitable for la pcr and gw were designed in various anchors, including the banna caecilian hox coding regions and the highly conserved noncoding regions within the hox clusters. the highly conserved noncoding regions were identified based on the multiz alignments of human, mouse, chicken, anole lizard, frog and zebrafish hox clusters constructed at the ucsc genome bioinformatics. la pcr was conducted to amplify the regions between every two adjacent anchors using takara la taq dna polymerase according to standard protocol. if there were weak or no bands on the agarose gel electrophoresis after the first-round of la pcr, nested pcr amplification was performed. for regions that still failed to amplify, gw was conducted using a genome walking kit based on the tail-pcr technique to shorten the distance between the two adjacent anchors. the physical linkages between the gw acquired fragments and the anchors were further verified by pcr amplification with primers designed in the anchors and the distal ends of the gw acquired fragments. after verification, the gw acquired fragments served as new anchors in which new primers were designed, and la-pcr was performed again with the new primers to amplify the regions. finally, we obtained la pcr fragments and a number of gw fragments. fragments less than 5 kb long were sequenced on the abi3730xl sequencing platform. longer fragments were sequenced using ion torrent at life technologies corporation. to prepare the sample, long pcr products were pooled in proportion to their length and digested with neb next ds dna fragmentase according to the manufacturer’s instructions to generate a sequencing library with 300–600 bp dna fragments. sequence reads obtained from ion pgm were de novo assembled using the mira program in its ‘genome’ assembly type with quality grades set as ‘accurate’. together, the assembled contigs from mira and sequence fragments of banna caecilian hox genes were reassembled using seqman in dnastar lasergene software and manually edited. the ambiguous contigs were verified by both blastn search against indonesian coelacanth or human hox clusters and further pcr surveys. the remaining gaps were filled through pcr amplification and primer walking. to annotate the newly obtained caecilian hox cluster sequences, we used genomescan and blastx search against indonesian coelacanth or human hox clusters to identify the hox coding sequences. exon-intron boundaries were determined manually based on the consensus splice motifs. screening of repetitive elements identification and classification of repetitive elements was conducted using censor with the default parameters against the repbase library of vertebrate repeat sequences as well as using repeatmasker.. simple repeats and low-complexity sequences were excluded. in addition to our newly obtained caecilian hox cluster sequences, genomic sequences of the four hox clusters in human, tammar wallaby, chicken, anole lizard, western clawed frog, indonesian coelacanth, spotted gar, and elephant shark were extracted from the ensembl and ncbi nucleotide databases for comparison. when comparing the repeat content among the human, wallaby, chicken, lizard, frog, caecilian and coelacanth hox clusters, only regions homologous to the obtained hox sequences in caecilian were analyzed, and regions homologous to the gaps in caecilian were excluded. and “density of repetitive sequences”, which was calculated as the ratio of the total length of repetitive sequences to the length of the clusters in each species, was used. when comparing the distribution of the repetitive sequences in caecilian and frog, all the available hox cluster sequences of the two amphibians were analyzed. identification and analyses of cnes hox clusters of human, tammar wallaby, chicken, anole lizard, western clawed frog, banna caecilian, indonesian coelacanth, spotted gar and elephant shark were aligned using the global alignment program lagan available on the vista website and screened for the presence of cnes. because of the poor sequence coverage, the chicken hoxc cluster was not included in the analysis. banna caecilian and the western clawed frog were used as reference sequence, respectively. the cnes were restricted to the sequences with a cut-off of ≥65 % identity across windows ≥50 bp. to determine the evolutionary origin of the cnes, all the cnes were classified into phylogenetic groups, “gnathostome”, “osteichthyan”, “sarcopterygian” and “tetrapod”. the “gnathostome” group was defined as cnes conserved in the elephant shark and caecilian or frog. the “osteichthyan” group referred to cnes shared by the spotted gar and caecilian or frog. the “sarcopterygian” group referred to cnes found in coelacanth and caecilian or frog. the remaining cnes were all classified in the “tetrapod” group. and cnes already included in one phylogenetic group were not counted in the next group. relative rate test of gene evolution to determine the rate of evolution of caecilian relative to other species, we performed the tajima relative rate test on the hox family. first, each gene-set was separately aligned using muscle and ambiguous alignment regions were removed using gblocks with no gaps allowed. sites with gaps or unknown data were excluded from the rrt. each comparison included two ingroups and one outgroup. the rrt on each gene-set was performed using in-house python scripts. for most of the hox genes, the rrts were performed using banna caecilian, western clawed frog, puerto rican worm lizard /chinese softshell turtle and human as the ingroups and indonesian coelacanth as the outgroup. the rrts on hoxc were conducted with elephant shark as the outgroup and caecilian, african lungfish, eel and zebrafish as the ingroups because the western clawed frog, puerto rican worm lizard and human do not have hoxc gene. for hoxc the ingroups were changed to caecilian, frog, alpine stream salamander and african lungfish and the outgroup was indonesian coelacanth. multiple rrts can be combined to determine the relative evolutionary speeds of several species. a hasse diagram, in which the slower-evolving genes were placed below the faster-evolving ones, was used to show the combined rrt results. the significance test corresponding to each rrt was drawn as a solid line or a dotted line. when comparing the rates between caecilian and coelacanth, elephant shark was used as the outgroup. we also performed an rrt on a carefully curated dataset consisting of orthologous genes, in which caecilian, frog and human were used as the ingroups and coelacanth was used as the outgroup. for both hox family and the orthologous genes, both nucleotide sequences and protein sequences were analyzed. phylogenetic analysis phylogenetic trees were constructed using hox cluster sequences and the orthologous genes. the four hox clusters from various gnathostomes were aligned with lagan. and the orthologous genes were aligned with muscle. the ambiguous positions in each alignment were eliminated by gblocks with half gaps allowed. maximum likelihood trees were constructed with raxml under the gtr + gamma + i model for nucleotide sequences and the lg + gamma model for protein sequences. node support was estimated from rapid bootstrap replicates.  hox clusters in banna caecilian based on the hox gene fragments from our previous pcr survey, we used a strategy that combined long-range pcr and genome walking to determine the sequences of the four hox clusters in banna caecilian. primers were designed in the coding regions of banna caecilian hox genes and the highly conserved noncoding regions in the sarcopterygian hox cluster alignments. after assembly, all the physical linkages between the newly acquired fragments were determined, which were then verified by further pcr amplification. in total, we retrieved  kb sequences of banna caecilian hox clusters, which included the nearly complete hoxa cluster, most of the hoxb cluster, part of the hoxc cluster and most of the hoxd cluster. the ncbi accession numbers of the four hox clusters are from kf to kf we were unable to sequence seven intergenic regions most likely because of the limitations of our strategy. they were mainly located in the hoxc cluster and the extremities of the hoxa, hoxb and hoxd clusters. four of these regions were putatively very large as estimated by previous studies in other animals, and therefore, beyond the scope of long range pcr. the remaining three gaps may have been due to either enlarged gene distances or complex structures that hindered the progression of the taq dna polymerase.fig. genomic organization of the hox clusters in banna caecilian and western clawed frog. schematic representations of the hoxa, hoxb, hoxc and hoxd clusters in caecilian are shown, with the frog hox clusters drawn for comparison. for each hox cluster, location of the genes, repetitive sequences and conserved noncoding elements are separately illustrated from top to bottom. total numbers of the repetitive sequences and conserved noncoding elements are indicated at the end of each cluster. exons of the genes are characterized by black bars and pseudogenes are denoted by the symbol Ψ. gaps in the intergenic regions are indicated by dotted lines. repetitive sequences identified by censor are represented by grey bars. conserved noncoding sequences predicted using mvista are classified into phylogenetic groups: “gnathostome ”, “osteichthyan ”, “sarcopterygian ” and “tetrapod ”. for clarity, cnes conserved in caecilian and frog only are not shown in the figure estimated from the sequences we obtained, the sizes of the caecilian hox a, hoxb and hoxd clusters are comparable to their counterparts in another amphibian, the western clawed frog. in both species, the hoxa, hoxb and hoxd clusters are ~ ~ and ~95 kb in size, respectively. for the less complete hoxc cluster, three fragments were compared. the hoxc11-hoxc and hoxc6-hoxc fragments in the two species were ~40 kb in size; while the hoxc13-hoxc fragment in caecilian was approximately only one-third of the size of its counterpart in frog. for the hoxa, hoxb and hoxd clusters, gene content and gene organization were largely determined. for hoxc cluster, gene linkages were determined except for hoxc12-hoxc hoxc8-hoxc hoxc4-hoxc and hoxc3-hoxc the gene content and organization of the caecilian hox clusters were similar to those in other amphibians, with slight differences. caecilian retained hoxc which was lost in many other vertebrates, including frog and salamander, but detected in coelacanth and lungfish. the complete coding sequence of hoxc that we obtained demonstrates that the caecilian hoxc is fully functional. hoxd was completely absent in frog and salamander, but remnants of the caecilian hoxd was observed in region 5′ to hoxd after performing a global alignment with other gnathostomes. in our caecilian hox cluster sequences, we found five microrna genes, two belong to family mir- and three belong to family mir- at the same genomic location observed in other vertebrates. mirna 196a- which is usually located approximately 6 kb 5′ to hoxb in other vertebrates was not identified even though  kb of sequence 5′ to the caecilian hoxb was screened. we could not detect mirna 196a- in the corresponding genomic region in western claw frog, either. thus, the mirna 196a- microrna in amphibians may be either lost or located further away from hoxb in summary, we found in banna caecilian hox genes and pseudogene organized into four clusters, evx paralogs associated with the hoxa and hoxd clusters, respectively, and micrornas. distribution of repetitive sequences obtaining the majority of the caecilian hox clusters allowed us to investigate the details of the hox cluster structure and elements in the intergenic regions. first, we evaluated the repetitive sequences which are usually strongly excluded from vertebrate hox clusters. we used censor and the repbase database to detect repetitive sequences in the obtained banna caecilian hox clusters and the corresponding regions in the human, wallaby, chicken, lizard, frog and coelacanth hox clusters. the density of repetitive sequences in the caecilian hox clusters was  %, among the lowest densities of all the species studied. it is possible that the density of repetitive sequences in banna caecilian may have been underestimated without a specific repeat library available. wallaby, lizard and coelacanth, however, do not have specific repeat libraries either, and they did not show a general bias towards low densities of repetitive sequences; and all of these species had a higher density of repetitive sequences than banna caecilian. thus, the effect of databases in identifying repeats does not appear to be predominant and the low density of repetitive sequences was considered primarily as a characteristic of the caecilian hox cluster sequences that we obtained. when the two amphibians were compared, the repeat content of the hox clusters in caecilian was approximately only a quarter of that in frog. to avoid any algorithm bias, repeatmasker was also used to detect repetitive sequences. the density of repetitive sequences in the caecilian hox clusters was again among the lowest densities and much lower than that of frog. because the general trend was consistent in the two analyses, only the results of censor are shown for the following comparisons.fig. densities of repetitive sequences of the gnathostome hox clusters. density of repetitive sequences for a specific hox cluster was calculated by dividing the total length of the repetitive sequences by the total length of the cluster. only repetitive sequences identified by censor with default settings were calculated. the chicken hoxc cluster was not included in the analysis because it is currently missing many sequences the attributes and distribution of the repeats detected in the caecilian and frog hox cluster sequences were further investigated. note that because all the unsequenced gaps might contain repeats, the number of repeats in each caecilian hox cluster was likely underestimated. the most abundant type of repeats in caecilian was non-ltr retrotransposons, while in frog approximately  % of the repeats were dna transposons, such as hats, harbingers and koloboks. among the four hox clusters, the caecilian hoxa cluster contained the least number of repetitive sequences interspersed in the cluster; the frog hoxa cluster contained repeats, 59 % of which were concentrated in the evx1-hoxa region. in the caecilian hoxb cluster, repeats were found, and over half of them were located between hoxb and hoxb1; while in the frog hoxb cluster, repeats were identified and a large portion of them were located in regions 5′ to hoxb hoxb8-hoxb and hoxb2-hoxb for both species, the hoxc cluster contained the most abundant amount of repetitive sequences, with and repeats in caecilian and frog, respectively. while the caecilian repeats were largely interspersed in the cluster,  % of the frog repeats were concentrated in the two terminal regions, i.e., upstream of hoxc and downstream of hoxc for the hoxd cluster,  % of the repeats in caecilian were located in the hoxd13-hoxd and hoxd8-hoxd regions, while approximately 50 % of the repeats in frog were in the hoxd3-hoxd region. thus, the distribution patterns of the repeats had “hotspots” that were in different locations in the two amphibians. and most of these “hotspots” were enlarged compared with the orthologous regions in other vertebrates. two of the gapped intergenic regions in caecilian, hoxb2-hoxb and hoxd13-hoxd were among the repeat “hotspots”. they might be also enlarged, thus were hard to amplify. to identify direct and inverted repeats that were not in the known repeat libraries, self-self blastn analyses were performed. a direct/inverted repeat consists of two repeat copies that are approximately identical/complementary to each other. this analysis revealed the presence of inverted repeats and direct repeats in the obtained caecilian hox clusters, over half of which were identified in the hoxd cluster. in particular, we found two long inverted repeats in the hoxd8-hoxd intergenic region. the arms of these two repeats were both longer than 100 bp, and the arm identity was almost 100 %. the spacers between the two arms of both repeats were small. in frog, inverted repeats and direct repeats were found in the self-self blastn analysis in corresponding regions of the hox clusters, but no long inverted repeat was identified. indeed, inverted repeats with an arm size longer than 100 bp are rare in genomes. there are only such repeats found in the human genome, and currently, no long inverted repeat has ever been reported in the hox clusters of other gnathostomes species. determining whether there is any functional significance of these two long inverted repeats in the caecilian hox cluster is worthy of further investigation. conserved noncoding elements in caecilian hox clusters next, we screened the caecilian hox clusters for conserved noncoding sequences, which may be indicative of potential regulatory elements. multiple global alignments of the four caecilian hox clusters and the homologous regions of other gnathostomes were conducted using lagan. from these alignments, conserved noncoding elements were predicted with vista. this method has been shown to be effective at identifying and visualizing overtly conserved non-coding elements. because of the incompleteness of the caecilian hoxc cluster, the number of cnes detected in the caecilian hoxc clusters may be an underestimation. furthermore, the chicken hoxc cluster was not included in this analysis because it is currently missing many sequences. the vista plots of the four hox clusters are shown in additional file the comparisons between caecilian and other gnathostomes identified cnes in intergenic and some intronic regions of the caecilian hox clusters. the distribution of these cnes is shown in fig.  and table  of the four clusters, the hoxa cluster contained the highest number of cnes, and the hoxc cluster had the highest number of cnes in introns. the cne densities in the 3′ part of the caecilian hoxa, hoxb and hoxd clusters tended to be higher than those in the 5′ part, consistent with the previous observations in many other gnathostomes. this density difference was not readily observed in the caecilian hoxc cluster because of the incompleteness especially in its 3′ part. overall, cnes tended to be located around the coding regions of the hox genes, and the regions containing more repeats tended to have lower cne densities.table cnes in the four hox clusters of banna caecilian for comparison, cnes in the homologous regions of the frog hox clusters were also identified using a similar method, and cnes were found, which was fewer than those identified in caecilian. comparing the cnes in the two amphibians, of the caecilian cnes corresponded to of the frog cnes as several long cnes in caecilian were broken into shorter cnes in frog. in addition to the cnes common to caecilian and frog, cnes were identified by comparing caecilian with other non-frog gnathostomes, but only cnes were identified by comparing frog with other non-caecilian gnathostomes. we next classified the cnes in caecilian and frog into phylogenetic groups, “gnathostome”, “osteichthyan”, “sarcopterygian” and “tetrapod”. the “gnathostome” group was defined as cnes at the same positions in the hox clusters of elephant shark and caecilian or frog. the “osteichthyan” group referred to cnes shared by spotted gar and caecilian or frog. the “sarcopterygian” group referred to cnes found in both indonesian coelacanth and caecilian or frog. the remaining cnes were all placed in the “tetrapod” group. note that the cnes already included in a previous phylogenetic group were not counted in the next groups. as shown in table  a large portion of the cnes found in both amphibians was categorized into the “tetrapod” and “gnathostome” groups. most of the cnes in the “gnathostome” group were conserved in almost all of the species studied and most of the cnes in the “tetrapod” group were also found in one or two other species. because similarities in sequences from highly divergent organisms often imply functional constraints, the cnes we identified, especially those in the “gnathostome” group, were most likely indicative of potential important functional elements. we analyzed the caecilian “gnathostome” cnes by blast searching the corresponding human sequences against the ncbi database of human ests, of which had est matches and were therefore putative non-coding rna genes. the remaining “gnathostome” cnes likely contained many potential cis-regulatory elements; of them were on the list of functionally-verified hox enhancers summarized in a study by ravi et al.. furthermore, cnes in the other three groups, especially those present in multiple organisms, may also be putative non-coding rnas or cis-regulatory elements and deserve further functional analysis.table comparing the caecilian or frog cnes in the four phylogenetic groups differences were observed between caecilian and frog in the retention of cnes at different evolutionary stages. among them, the most prominent were in the “sarcopterygian” group in which there was fold greater difference in the cnes in caecilian compared with frog and this difference increased to more than -fold if the cnes common to caecilian and frog were excluded. this difference was mostly observed in the hoxa, hoxb and hoxc clusters, especially in the hoxb cluster where caecilian had “sarcopterygian” cnes while frog had none. this variation in the “sarcopterygian” group was the primary reason for the different number of cnes found in caecilian and frog. slow evolutionary rates of caecilian genes the evolutionary rates of the banna caecilian hox genes were estimated using the tajima’s relative rate test in which banna caecilian, western clawed frog, puerto rican worm lizard/chinese softshell turtle and human were first used as the ingroups and indonesian coelacanth as the outgroup. because of the missing hoxc and hoxc in certain tetrapods, of the hox genes were analyzed. as shown in fig.  over 62 % of the analyzed hox coding sequences showed significantly slower evolutionary rates in caecilian than in frog, whereas the reverse was true for no hox genes. nearly 84 % of the analyzed hox coding sequences were identified evolving significantly more slowly in caecilian than in the mammal representative, human, whereas the reverse was true for no hox genes. because anole lizard’s hox genes generally evolved fast, reptiles were represented by slow worm and chinese softshell turtle. over 45 % of the analyzed hox genes revealed significantly slower evolutionary rates in caecilian than in the reptile representatives. thus, most of the hox genes analyzed evolved significantly more slowly in banna caecilian than in other species, of which evolved the most slowly in banna caecilian compared with all other species studied. because indonesian coelacanth hox clusters have been reported to evolve comparatively slowly, we next compared the relative evolutionary rates of caecilian and coelacanth using the elephant shark as the outgroup. all of the hox genes showed no significant difference in evolutionary rate, but appeared to evolve faster in caecilian than in coelacanth. in summary, the hox gene sequences evolved much more slowly in caecilian than in other tetrapod species, except for a few genes located in the hoxc cluster, and approximately 81 % of the hox genes demonstrated similar evolutionary rates to those in the slowly evolving indonesian coelacanth. rrts were also performed using the hox protein sequences to avoid inaccurate estimation due to the issue of saturation when using nucleotide sequences. twenty one of the caecilian hox proteins showed slower evolutionary rates than at least one other species, while only frog hox proteins were slower than at least one other species. thirteen caecilian hox proteins evolved more slowly than the frog ones, but only one frog hox protein evolved more slowly than its caecilian counterpart. thus, though the evolutionary rate differences of the hox protein sequences were not as significant as the nucleotide sequences, caecilian still appears a slow evolving species.fig. rrts of the gnathostome hox gene sequences. a evolutionary relationship of the gnathostome species used in the rrts of the hox genes. b summary of the rrts conducted on the hox gene nucleotide sequences. for most of the hox genes, rrts were performed using banna caecilian, western clawed frog, puerto rican worm lizard/chinese softshell turtle and human as the ingroups and indonesian coelacanth as the outgroup. rrts of hoxc were carried out using elephant shark as the outgroup and caecilian, african lungfish, eel and zebrafish as the ingroups to replace the species that do not have hoxc gene. for hoxc the ingroups included caecilian, frog, alpine stream salamander and african lungfish and the outgroup was indonesian coelacanth. results of rrts for each gene are shown in a hasse diagram, in which the slower-evolving genes are placed below the faster-evolving ones, with statistical significance denoted as a solid line or a dotted line. c summary of the significant rrts among banna caecilian, western clawed frog, puerto rican worm lizard and human. for each pair of species, the significant rrts are denoted by arrows pointing to the slower-evolving one. the number of highly significant tests for each hox cluster is indicated on the side of the arrow additionally, phylogenetic trees of the four hox clusters were constructed using all of the alignable sequences. branch lengths reflect the differences in evolutionary rates among and within lineages with longer branch lengths representing faster evolutionary rates. as observed in fig.  in the hoxa, hoxb and hoxd cluster trees branch lengths of frog were -fold, -fold and 2-fold longer than those of caecilian, respectively. in the hoxc cluster tree, the branch length of frog was also longer than that of caecilian, though this difference appeared to be smaller than the differences observed in the other three clusters. all of the phylogenetic trees further demonstrate that the caecilian hox clusters evolved more slowly than the frog hox clusters.fig. phylogenetic trees of the gnathostome hox gene clusters. all the alignable sequences for each hox cluster were used. the trees were constructed using raxml under the gtr + gamma + i model. elephant shark was used as the outgroup. values above branches denote maximum likelihood bootstrap support next, we wondered whether the slow evolution of the caecilian hox clusters reflected a genome-wide phenomenon, as in the case of coelacanth. to test this idea, we performed an rrt analysis with both nucleotide and protein sequences using orthologous genes from human, frog and caecilian. for the nucleotide dataset, as many as  % genes in caecilian were significantly slower than the orthologs from either frog or human. approximately  % caecilian genes showed the slowest evolutionary rate, i.e., slower than both frog and human. in frog, the percentage of genes slower than their caecilian and/or human orthologs was only  %. the rrt result of the protein sequences showed a similar trend. up to  % genes in caecilian were significantly slower than the orthologs from either frog or human. in frog, this percentage was  %. the percentages of genes showed the slowest evolutionary rate in caecilian and frog were and  %, respectively. additionally, we constructed a phylogenetic tree with the concatenated orthologous genes from vertebrate species. both nucleotide sequences and protein sequences were analyzed. in both trees, the branch lengths of caecilian were short which suggested that it evolved slowly. among the three distinct living amphibian lineages, the evolutionary rates of frog and axolotl were both greater than that of caecilian. thus, not only the hox genes, but also a large portion of the protein-coding genes are evolving significantly more slowly in caecilian.fig. genome-wide rrts of banna caecilian, western clawed frog and human using nucleotide sequences. rrts were performed on the nucleotide sequences of orthologous genes, with banna caecilian, western clawed frog and human as the ingroups and indonesian coelacanth as the outgroup. the number in a cell represents the number of statistical significance for rrts of a species compared with the others; hence, there are six possibilities for each ingroup species. increasing warm color intensity indicates faster relative evolutionary rate, whereas increasing cold color intensity indicates slower relative evolutionary rate discussion in this study, we successfully obtained most of the banna caecilian hox cluster sequences using la pcr and gw. we identified hox genes, evx paralogs, microrna genes and pseudogene in four clusters. most of their physical linkages were determined which was not accomplished in previous pcr surveys. the overall organization of the caecilian hox clusters was conserved, similar to those from other species with four clusters. the presence of hoxc in caecilian is unique compared with all the tetrapods studied so far, which suggests that the tetrapod ancestor had hoxc but lost it in different lineages. loss or pseudogenization of hoxd appears to be a distinct feature of the amphibians’ hox clusters and has been proposed to be related to digit reduction as frogs and salamanders normally only have four fingers and caecilians are limbless. other surveys have shown that hoxd is also lost in snakes, which are limbless, and the african lungfish, which has only thread-like fins, but it is present in many limbless lizards. therefore, there may be different mechanisms responsible for digit reduction and hoxd may not be directly involved in all of them. alternatively, it could be the insertion of repetitive sequences that first disrupted the regulatory network of the posterior hoxd genes and altered the developmental program for fingers and even limbs. these effects of repeat insertion may then account for the pseudogenization of hoxd it would be interesting to characterize the posterior hoxd cluster in limbless lizards and see whether there are overtly multiple repeats in the region though the hox genes remain intact. from the obtained sequences, the caecilian hox clusters appear compact and have fewer repetitive sequences than the hox clusters in the other organisms used in this study. though both are amphibians, caecilian and frog have much different repeat contents, which may reflect the different evolutionary constraints on the hox clusters of the two amphibians. there were, however, gaps in several intergenic regions that we could not fill in caecilian. based on the hox cluster information from other amphibians, one possible reason is that these regions have been greatly enlarged, such as the expanded hoxd11-hoxd in the caecilian typhlonectes natans and the red spotted newt and the expanded hoxc13-hoxc hoxc5-hoxc and hoxd3-hoxd in the western claw frog. repetitive sequences are likely abundant in these possibly expanded regions in banna caecilian as they have accumulated in the known large intergenic regions in the hox clusters of many organisms. two of them have actually shown signs of repeat accumulation. thus, these gapped intergenic regions were possibly repeat “hotspots” and the amount of repetitive sequences in the banna caecilian hox clusters may currently be underestimated. in the future, obtaining the full sequences of these gaps and determining the structural features will elucidate whether they are indeed expanded and have multiple repeat insertions. if this were the case, it would also help to understand why repeat hotspots exist. expansion of intergenic regions may affect gene regulation by long-range regulatory elements. insertion of repetitive dna has the potential to alter local regulatory structures by undermining the conserved motifs and/or creating new ones, further influencing gene expression. alternatively, the presence of repetitive dna may just indicate weakened coherence or relaxed constraint within the clusters. why repeats are enriched in certain intergenic regions and whether they are related to developmental variations of caecilians or just a reflection of the relaxation of evolutionary constraints remains unknown. vista plots demonstrated a high level of conservation in the non-protein-coding regions of the caecilian hox clusters. in contrast to repetitive sequences, cnes were abundant in the obtained caecilian hox cluster sequences. and cnes tended to be located in the regions containing few repeats. this finding may reflect the preservation of regulatory modules that formed the complex regulatory network of the hox genes. the phylogenetic group classifications suggested that the cnes in the hox clusters had different conservation depths. the ancient and well conserved cnes may be fundamental for the regulatory network and proper hox functions. other cnes may then be gradually added to or removed from the regulatory network at different evolutionary stages and in different lineages, thus complicating and rewiring the regulatory network in different species. comparative analyses of the hox cluster sequences from more vertebrate species will improve the accuracy of this classification and elucidate the functions of the cne. that caecilian hox clusters contained more cnes than frog, with the major difference in the “sarcopterygian” group, implies that the frog hox cluster sequences evolved faster than the caecilian sequences and that the caecilian hox cluster sequences retained more characteristics of the early tetrapods than the frog sequences, such as hoxc and the hoxd remnant. the slow evolution of the caecilian hox clusters was supported by the rrts performed for the hox protein sequences and the phylogenetic tree analyses of all the alignable hox cluster sequences. the majority of the caecilian hox genes had evolutionary rates similar to the indonesian coelacanth hox genes which were known to have evolved slowly. as in the case of coelacanth, the slow evolution of the caecilian hox clusters may be indicative of its slowly evolving genome. it was supported by the rrts and the phylogenetic tree analysis of protein homologs which showed a general pattern of slow evolutionary rate of the caecilian sequences. thus from the evolutionary rate of the hox cluster and the protein dataset, it was quite likely that the caecilian genome was evolving slowly. however, this still require support from more genomic information. for this rarely found organism, deciphering more of its genomic information and studying it from the inside out would greatly complement the studies on its morphology, physiology and function. furthermore, to better understand the genome biology of living amphibians, a caecilian genome project is especially of practical importance, because salamanders possess genome sizes that are intractably large for routine genomic analyses. and as the most basal tetrapod lineage, the caecilian genome is also highly favorable for comparative genomic studies. for example, comparative analyses between the caecilian and coelacanth genomes would provide important clues for reconstructing genomic perspectives of ancestral tetrapods and understanding the water-to-land transition. currently, there are fewer than nucleotide sequences for caecilian in genbank, and most of these sequences are mitochondrial gene markers from phylogenetic studies. our study sequencing the majority of the four caecilian hox clusters, is an important data resource and provides an amount of genomic information for the cryptic caecilian. in the future, caecilian shall be added to the list of organisms in demand of genome projects, and characterizing the caecilian genome will not only help to understand its biology, but also will provide insights into the genome biology and evolution of early tetrapods.  in this study, we have cloned and sequenced most of the banna caecilian hox cluster. thirty nine hox genes, microrna genes and pseudogene were identified and most of their physical linkages in the four clusters were determined. the presence of hoxc in caecilian suggests that the tetrapod ancestor had hoxc but then lost it in different lineages. loss of hoxd function appears to be a feature of the amphibians hox clusters. whether and how it was related to digit reduction requires further investigation. from the obtained sequences, the caecilian hox clusters appear compact and have fewer repetitive sequences than the hox clusters in the other organisms studied in this work. however, to fully understand the repetitive sequences in the caecilian hox clusters, the intergenic gaps are to be sequenced and their structural features be determined in future studies. the analysis on cnes, the rrts on the hox genes and the phylogenetic trees constructed from the hox clusters all revealed stronger constraints on the caecilian hox clusters than their frog counterparts. these information may be indicative of a slowly evolving genome, which was supported by rrt analyses and the phylogenetic tree construction using a large orthologous protein dataset. therefore, characterizing more of the caecilian genomic information is in demand, which will not only help to understand the basic biology of this clade, but also be highly valuable for the comparative genomics of amphibians and the evolutionary studies of early tetrapods. additional files additional file 1: global genomic alignments of the posterior hoxd cluster. exons and transcription direction of the genes are indicated. sequence comparisons were conducted using the human sequence as the reference. nucleotide identities relative to the human sequence are given by histogram peaks. cnes are depicted by red peaks. additional file 2: global genomic alignments of part of the hoxb cluster. exons and transcription direction of the genes are indicated. sequence comparisons were conducted using the human sequence as the reference. nucleotide identities relative to the human sequence are given by histogram peaks. cnes are depicted by red peaks. additional file 3: global genomic alignments of the four hox clusters. exons and transcription direction of the genes are indicated. sequence comparisons were conducted using the caecilian hox clusters as the reference sequence. nucleotide identities relative to the caecilian sequences are given by histogram peaks. cnes are depicted by red peaks. additional file 4: conservation of the cnes in banna caecilian hox clusters belonging to the phylogenetic groups. nine species were used for the global genomic alignments of the hoxa, hoxb and hoxd clusters; for the hoxc cluster, only species were used due to the low sequence coverage of the chicken hoxc cluster. numbers indicate the count of species that retained an identical cne. from the “tetrapod” group to the “gnathostome” group, the color changes from light red to dark red gradually. additional file 5: rrts of the gnathostome hox gene using protein sequences. evolutionary relationship of the gnathostome species used in the rrts of the hox protein sequences. summary of the rrts conducted on the hox gene protein sequences. for most of the hox genes, rrts were performed using banna caecilian, western clawed frog, puerto rican worm lizard/chinese softshell turtle and human as the ingroups and indonesian coelacanth as the outgroup. rrts of hoxc were carried out using elephant shark as the outgroup and caecilian, african lungfish, eel and zebrafish as the ingroups to replace the species that do not have hoxc gene. for hoxc the ingroups included caecilian, frog, alpine stream salamander and african lungfish and the outgroup was indonesian coelacanth. results of rrts for each gene are shown in a hasse diagram, in which the slower-evolving genes are placed below the faster-evolving ones, with statistical significance denoted as a solid line or a dotted line. summary of the significant rrts among banna caecilian, western clawed frog, puerto rican worm lizard and human. for each pair of species, the significant rrts are denoted by arrows pointing to the slower-evolving one. the number of highly significant tests for each hox cluster is indicated on the side of the arrow. additional file 6: genome-wide rrts of banna caecilian, western clawed frog and human using protein sequences. rrts were performed on the protein sequences of orthologous genes, with banna caecilian, western clawed frog and human as the ingroups and indonesian coelacanth as the outgroup. the number in a cell represents the number of statistical significance for rrts of a species compared with the others; hence, there are six possibilities for each ingroup species. increasing warm color intensity indicates faster relative evolutionary rate, whereas increasing cold color intensity indicates slower relative evolutionary rate. additional file 7: phylogenetic tree inferred from a dataset of nuclear protein-coding genes. nucleotide tree. protein tree. putative orthologous genes in two species were identified using the mutual best hit in basic local alignment search tool. the tree was constructed using the maximum likelihood method with raxml under the gtr + gamma + i model for nucleotide sequences and the lg + gamma model for protein sequence. elephant shark was used as the outgroup. values above branches denote maximum likelihood bootstrap support. competing interests the authors declare that they have no competing interests. authors’ contributions dl, pz and rw conceived the study. sm provided the species sample. rw and ql designed and carried out the experiments. rw and ql analyzed the data. dl, rw and ql wrote the manuscript. all authors read and approved the final manuscript. we thank rongli mao for helping in experiments and all the lab members for insightful discussions on the manuscript. this work was supported by the national natural science foundation of china grants to pz and sm, respectively, the china national fund for excellent young scientists and the cultivation project for young teachers of sysu to pz and the cultivation project of the department of education of guangdong province to dl.
phylogenetic tree of the four hox clusters also revealed shorter branch length especially for the caecilian hoxa, hoxb and hoxd clusters. hox genes, an old family of transcription factors playing central roles in the establishment of metazoan body plan. we obtained the majority of the four caecilian hox clusters and identified hox genes, microrna genes and pseudogene. relative rate of sequence evolution showed that caecilian hox genes evolved significantly more slowly than the other tetrapod species used in this study and were comparable to the slowly evolving coelacanth hox genes. there remained seven intergenic gaps we were unable to fill. the caecilian genomic sequences are also important and practical in the comparative genomics of amphibians. from the obtained hox sequences and the orthologous protein dataset, the caecilian hox loci and its genome appear evolving comparatively slowly. from the obtained sequences, the caecilian hox clusters contained less repetitive sequences and more conserved noncoding elements than the frog counterparts. as the basal lineage of amphibians and land vertebrate, this characteristic of the caecilian genome is valuable in the study concerning the genome biology and evolution of amphibians and early tetrapods. we sequenced and characterized the hox clusters of banna caecilian with a strategy combining long range pcr and genome walking. these features of the caecilian hox clusters suggested a slowly evolving genome, which was supported by further analysis of a large orthologous protein dataset. our analyses greatly extended the knowledge about the caecilian hox clusters from previous pcr surveys. understanding their structure and genomic organization may provide insights into the animal’s genome, which is valuable for animals without a sequenced genome.
 caecilians, with a discrete lifestyle, are the least explored group of amphibians. though with distinct traits, many aspects of their biology are poorly investigated. obtaining the caecilian genomic sequences will offer new perspectives and aid the fundamental studies in caecilian biology. the caecilian genomic sequences are also important and practical in the comparative genomics of amphibians. currently, however, only sparse genomic sequences of caecilians are available. hox genes, an old family of transcription factors playing central roles in the establishment of metazoan body plan. understanding their structure and genomic organization may provide insights into the animal’s genome, which is valuable for animals without a sequenced genome.  we sequenced and characterized the hox clusters of banna caecilian with a strategy combining long range pcr and genome walking. we obtained the majority of the four caecilian hox clusters and identified hox genes, microrna genes and pseudogene. there remained seven intergenic gaps we were unable to fill. from the obtained sequences, the caecilian hox clusters contained less repetitive sequences and more conserved noncoding elements than the frog counterparts. we found that caecilian and coelacanth shared many more cnes than frog and coelacanth did. relative rate of sequence evolution showed that caecilian hox genes evolved significantly more slowly than the other tetrapod species used in this study and were comparable to the slowly evolving coelacanth hox genes. phylogenetic tree of the four hox clusters also revealed shorter branch length especially for the caecilian hoxa, hoxb and hoxd clusters. these features of the caecilian hox clusters suggested a slowly evolving genome, which was supported by further analysis of a large orthologous protein dataset.  our analyses greatly extended the knowledge about the caecilian hox clusters from previous pcr surveys. from the obtained hox sequences and the orthologous protein dataset, the caecilian hox loci and its genome appear evolving comparatively slowly. as the basal lineage of amphibians and land vertebrate, this characteristic of the caecilian genome is valuable in the study concerning the genome biology and evolution of amphibians and early tetrapods. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. issue-copyright-statement© the author 2015 