BACKGROUND
bdelloid rotifers are microscopic aquatic invertebrates of particular interest for two features of their lifestyles.
first, they have persisted for tens of millions of years  <cit>  and diversified into nearly  <dig> recognized species  <cit> , despite the apparent absence of meiosis and sex  <cit> .
second, adult bdelloids can survive essentially complete desiccation in a state known as anhydrobiosis  <cit> .
moreover, bdelloid genomes show a number of unusual features that might relate to their long-term asexuality and ability to survive desiccation  <cit> .
understanding these features can shed light not only on bdelloids’ unusual lifestyle, but also on the range of mechanisms by which animal genomes can evolve.
in efforts to understand bdelloid genome organization, much interest has focused on gene copy numbers.
the first work on bdelloid genomes found that bdelloids possessed divergent copies of genes present as a single diploid locus in their closest sexual relatives, the monogonont rotifers  <cit> .
it was hypothesized that these were former alleles that had diverged following the loss of meiosis and recombination  <cit> , a process called the meselson effect  <cit> .
later work suggested that putative former alleles of lea  protein had diverged so much as to adopt different biochemical roles protecting bdelloid cells from desiccation  <cit> .
although functional divergence of gene copies is not restricted to asexual lineages, it was argued that functional divergence of former alleles might have contributed towards adaptation to desiccating environments.
subsequent work has shown that bdelloid genome complexity cannot be explained simply through divergence of former diploid alleles: bdelloids are degenerate tetraploids  <cit> .
sequencing the genomic regions surrounding four copies of hsp <dig> revealed that the chromosomes comprise two collinear pairs.
between pairs, divergence is so high that they only share a handful of genes, whereas within collinear pairs divergence is much lower and there are regions of complete identity.
the same pattern was found in other genomic regions and in different bdelloid families  <cit> .
it was proposed that bdelloids use pairs of homologous chromosomes as templates to repair dna breaks caused during desiccation, and in the process sequences are homogenised at broken regions by gene conversion  <cit> .
this provides an explanation for the patchy nature of divergence within collinear pairs and for bdelloids’ uncommon ability to repair dna breaks caused by desiccation and ionizing radiation  <cit> .
the original cause of tetraploidy is unknown but could have involved either a genome-wide duplication event or hybridization between sexual ancestors – the latter being commonly associated with asexuality in many organisms  <cit> .
together, these findings make bdelloids an interesting case study for understanding opposing forces on gene copy evolution.
although divergent copies of bdelloid genes can no longer be interpreted simply as former alleles, it remains of interest to know how often divergent copies of genes have been maintained and diverged in function.
in principle, the maintenance of functionally divergent copies of genes could provide bdelloids with an enhanced genetic repertoire and, perhaps, increased phenotypic plasticity  <cit> .
here, we explore copy evolution in alpha tubulin and ask whether divergent copies have been maintained in bdelloid genomes.
alpha tubulin forms microtubules with beta tubulin and provides a range of functions in cellular movement and cell division.
previous studies of copy variation of alpha tubulin in other organisms revealed contrasting patterns making it an interesting gene to study in bdelloids  <cit> .
in many animals, it is present in low copy number and highly conserved, which has led to its use as a marker for deep phylogenetic relationships  <cit> .
however, the different roles of microtubules have led to functional divergence in many cases.
four distinct copies are present in drosophila melanogaster and expressed in different tissues: including ovary and testis-specific forms  <cit> .
if there were male-function or meiosis specific copies in sexual rotifers  <cit> , bdelloids might be expected to have lost those copies or for them to have evolved new functions.
duplication and diversification of tubulin genes has also been implicated in adaptation to extreme environments in the antarctic fish notothenia coriiceps <cit> , which expresses at least nine different alpha tubulin isotypes that vary in stability and polymerization rates at extremes of temperature  <cit> .
by expressing copies stable at different temperatures, function is maintained despite variation in environmental conditions.
in plant seedlings alpha tubulin is up-regulated during both dehydration and rehydration  <cit> .
therefore, bdelloids might make use of their ancestral tetraploid state and maintain alpha tubulin copies to perform specialist functions during cycles of desiccation  <cit> .
we used cloning and preliminary whole transcriptome data to characterize alpha tubulin copies in several monogonont and bdelloid rotifer species.
monogononts contained either one or two copies with very similar amino acid sequence .
contrary to the four or eight copies expected if bdelloids retained copies from their tetraploid ancestor, we found between  <dig> and  <dig> copies in each of four species of bdelloids.
phylogenetic relationships indicate that multiple copies were inherited from the common ancestor of bdelloids and some of them were lost or homogenized in some species.
approximately half of the copies retained amino acid sequences similar to the monogonont copies, whereas the rest have diverged into four distinct classes that differ in intron structure and biochemical properties predicted by computational analysis.
distinct classes have diverged significantly in amino acid sequence relative to the purifying selection observed within classes.
the findings do not contradict current understanding of bdelloid genome organization, but show how further duplication events in the ancestor of bdelloids led to proliferation in both conserved and functionally divergent copies of this gene.
the possible functional significance of alpha tubulin diversity is discussed in relation to sequence characteristics.
methods
study species, dna extraction, cloning and sequencing
four species of bdelloid rotifer belonging to two bdelloid families were compared to four species of monogonont rotifers also from two families .
sequences for one of the monogononts, brachionus plicatilis, were taken from the cdna library of suga et al.
<cit> .
for the other monogonont species dna was extracted from wild-caught individuals using chelex instagene matrix  following the manufacturers instructions.
genomic dna from bdelloid species was extracted using the dneasy blood and tissue kit .
for each bdelloid species dna was extracted from ~ <dig> individuals from populations grown from a single female.
extractions were performed according to manufacturers instructions, and samples were heated at 56°c for 1 h for protein digestion.
dna quality was checked on a 1% agarose gel.
pcr amplification was performed using 2 μl template dna and illustra pure-taq ready-to-go pcr beads.
alpha tubulin was amplified using the primers toa1f 5′-rgtnggnaaygcntgytggga-3′ and toa1r 5′-ccatnccytcnccnacrtacca-3′ <cit> , which amplify >90% of the coding region of each gene  <cit> .
the amplification profile consisted of an initial denaturation of 5 min at 94°c, followed by 35 cycles of 45 s at 94°c, 60s at 58°c and 60s at 72°c and a final extension of 5 min at 72°c.
products of expected size  were cloned using the topo ta cloning kit  and alpha tubulin inserts were re-amplified from individual clones.
amplicons were sequenced in both directions using the bigdye terminator v <dig>  cycle sequencing kit  and an abi  <dig> genetic analyzer .
sequence data are available in genbank .
sequence alignment
sequences were confirmed to be alpha tubulin using blastn, and aligned in geneious pro  <dig> .
<dig>  <cit>  using muscle  <cit>  followed by manual adjustments.
wise <dig>  <cit>  was used to define introns, which were manually checked by comparison to alpha tubulin cdna from the monogonont rotifer brachionus plicatilis .
sequences used in our analyses contain no premature stop codons.
taq polymerase has a baseline error rate of  <dig> in  <dig> base pairs  <cit> .
normally, when sequencing from uncloned pcr products this error rate is negligible, as errors occur at random sites across the amplified products and are averaged out upon sequencing.
however, cloning to separate different copies of alpha tubulin within each sample introduces greater errors into the sequencing process  <cit> .
based on the estimated error rate of taq of  <dig> in  <dig> and the 35 cycles in our pcr reaction, we estimated an overall error rate of  <dig> % per sequence , which we rounded to  <dig> % to allow for any additional error in sequencing these long products.
consequently, sequences obtained from cloning were considered to be the same if they were <1% different, and collapsed to a single majority rule consensus sequence, thus providing us with a conservative estimate of copy number.
plots of the number of distinct copies against the number of clones sequenced were used to assess whether most copies present in the genome had been discovered.
the chao estimator was calculated to estimate the potential number of missing copies  <cit> , based on the assumption that if a few copies have been sampled only once or twice, then it is likely that some copies remain unsampled.
phylogenetic analyses
relationships among copies were inferred from exon sequences using bayesian markov chain monte carlo  analyses implemented in p <dig> .
we used p <dig> because preliminary investigations revealed compositional heterogeneity in our sequences [see additional file  <dig>  introns were excluded from phylogenetic analysis because their extreme variability in sequence and insertion site confounded the evolutionary signal.
jmodeltest  <cit>  selected the gtr  model  <cit>  with gamma-distributed rate variation and a proportion of invariable sites  as the most likely substitution model using the akaike information criterion.
a chi-squared test revealed significant compositional homogeneity over the tree , and comparison of the observed compositional heterogeneity among sequences to a null distribution simulated using the fitted gtr + i + g model, confirmed that composition of the observed sequences was significantly heterogeneous .
the gtr + i + g model was therefore compared to node-discrete composition heterogeneity models ndch and ndch  <cit> , which incorporate two and five tree-heterogeneous base-composition vectors respectively [see additional file  <dig>  the preferred model was chosen using bayes factors: a loge of  <dig> or more was considered as very strong support for improved fit of a more complex model  <cit> .
bayesian mcmc analysis was run with the chosen model and two independent runs of four chains each for  <dig> generations.
posterior probabilities were calculated for support.
a second tree was reconstructed, using maximum likelihood implemented in phyml  <dig> .
<dig>  <cit> , from the translated amino acid alignment.
prottest  <dig>   <cit>  selected the jtt model of amino acid substitution  <cit>  with gamma-distributed rate variation.
statistical support was assessed using alrt non-parametric branch support  <cit> .
trees were rooted on the branch connecting bdelloid and monogonont sequences.
tests for selection
we estimated the ratio of the rates of non-synonymous and synonymous substitutions  using likelihood methods implemented in paml version  <dig>  <cit> .
first, we fitted a model with a single dn/ds ratio across all bdelloid copies of alpha tubulin.
second, to test for divergent selection acting on different classes of alpha tubulin, we fitted a model with a different dn/ds ratio on branches connecting different classes  versus on branches connecting sequences of the same class  <cit> .
third, we fitted a branch-site model to test for the presence of codons showing purifying selection on within-class branches  but evolving either neutrally or by positive selection on between-class branches   <cit> .
bayes empirical bayes  inference was used to identify codons belonging to this set.
we predict that the branch-site model should be preferred if amino acid divergence between the different classes has evolved through selection for functional specialization of some parts but not all of the protein.
the best model was selected using the akaike information criterion.
models were run using the average nucleotide frequencies at each codon position to specify equilibrium codon frequencies .
to check for sensitivity to the observed compositional heterogeneity, we repeated the analyses using free parameters for each codon frequency .
all inferences from model results were qualitatively unchanged over the simpler versions.
adineta ricciae transcriptome searches
to further assess the completeness of sampling of copies by cloning and sequencing, and to establish if multiple copies were expressed under experimental conditions, we searched an expressed sequence tag  library  <cit>  and draft transcriptome  <cit>  sequences of a. ricciae for alpha tubulin copies.
the transcriptome was generated by paired-end sequencing of ~200bp cdna fragments using illumina parallel sequencing.
the resulting  <dig>  million 76 bp reads were assembled with the clc-bio  assembler to give a library of  <dig>  transcripts of size range 118-3674 bp.
more details are provided in  <cit> .
the databases were searched by blastn  <cit>  for matches to the sequences obtained by cloning .
protein chemical analyses
three parameters indicative of predicted protein chemistry, and hence potentially of changes in protein function, were calculated for each alpha tubulin amino acid sequence using protparam  <cit> .
i) the grand average of hydropathy  index indicates the hydropathy of proteins: a positive gravy score means the protein is hydrophobic and a negative score hydrophilic  <cit> .
detrich  <cit>  hypothesized that increasing contact surface hydrophobicity might enhance polymer stabilizing interactions in microtubules, so a more hydrophobic  gravy score could indicate a more stable microtubule.
ii) isoelectric point  is the ph at which a protein carries no net electrical charge, which can affect its solubility at a given ph.
zambito et al.
<cit>  found a wide variation in the pi of rat brain tubulin isotypes, and showed that some variants are associated with particular cellular compartments.
variation in pi of our alpha tubulin copies may therefore imply variation in cellular function.
iii) the instability index provides an estimate of a protein’s folding stability in vitro, based on estimated instabilities of the protein’s dipeptides.
a protein with an instability index of less than  <dig> is predicted as stable, a value above  <dig> predicts that the protein may be unstable, and therefore less viable as a functional protein.
predictive discriminant analysis was performed using the predict.lda function of the mass library version  <dig> - <dig>  <cit>  in the r statistical programming language r  <dig> .
<dig>  <cit>  to see if classes of copies based on intron structure are predicted by these three protein chemical properties, i.e.
whether copies with different intron structure have divergent chemical properties.
multivariate analysis of protein polymorphism  scores of physiochemical impact were estimated for all codons  <cit> , providing a measure of constraint at each site .
to determine if sites that diverged significantly between classes had significantly different mapp scores to the rest of the protein, the mean mapp score at codons identified by the beb analysis  to be significantly divergent between classes  was compared to  <dig> replicates of the mean mapp score at the same number of random sites drawn from a set of all sites in the protein.
to determine if sites identified by the beb analysis had significantly different mapp scores to polymorphic sites not evolving faster between classes than within them, the mean mapp score at the codons identified by the beb analysis  was compared to  <dig> replicates of the mean mapp score at the same number of random sites drawn from a set of all polymorphic sites.
protein structural analyses
translated protein sequences were mapped to the structural and functional regions of alpha tubulin identified by  <cit> .
regions of alpha helices and beta sheets were classified as structural.
regions falling into the following four types were classified as functional: the interface with beta tubulin; the regions interacting with gtp ; the tubulin nonexchangeable site ; and the c-terminal helices  where maps  and motor proteins usually bind.
to determine whether amino acid divergence was concentrated in particular structural regions, amino acid positions identified by the beb analysis to be significantly divergent between classes were mapped onto the identified functional and structural regions.
monte carlo tests were used to test for significant concentration of divergent amino acids either in structural or in functional regions of the genome.
RESULTS
number of copies of alpha tubulin genes from cloning and sequencing
the b. plicatilis est database of suga et al.
<cit>  contained two alpha tubulin copies  that were 13% divergent at the nucleotide level but only differed in three amino acids.
amplifying and sequencing alpha tubulin in the three selected monogonont species  revealed just one copy of alpha tubulin in each species.
pairwise amino acid divergence between all monogonont alpha tubulin sequences was only  <dig> %.
both b. plicatilis copies contained our primer binding sites and should therefore have been amplified using our primers, if they had been present in the other monogonont species.
cloning and sequencing of bdelloid alpha tubulin yielded between  <dig> and  <dig> different copies in each species .
this is a conservative estimate because sequences were only considered unique when they varied by >1% to account for  <dig> % pcr fidelity of the cloned copies.
including all sequences obtained by cloning would increase the number of copies per species to almost as many as the number sequenced, but would not alter our broad findings beyond increasing the number of very similar copies within species [see additional file  <dig>  saturation curves showed that the number of copies reached a plateau against the number of clones sequenced in each species .
the chao <dig> estimator indicated that all copies were likely to have been sampled in all species.
intron and exon structure
no introns were present in the monogonont sequences.
in contrast, each bdelloid alpha tubulin copy contained between  <dig> and  <dig> introns located across  <dig> different sites .
the copies fall into five major classes based on differences in intron content and location, with some additional losses or gains of single introns in particular species leading to a total of eight different intron patterns .
most introns were between  <dig> and  <dig> nucleotides in length , but two introns were consistently longer across copies and species .
all introns present in two or more classes were always found in class  <dig> and there were no unique introns in class  <dig>  class  <dig> only differed from class  <dig> by the loss of one intron  whereas classes  <dig>   <dig> and  <dig> each contained one or more unique introns.
all seven introns present in the class  <dig> copies shared a ggt motif at the 5’ splice site across all the classes.
one intron unique to class  <dig> had the same 5’ splice site as the introns in class  <dig> copies and represents a putative translocation of an adjacent intron that is present in all the other classes but lacking in class  <dig>  all other introns had different 5′ splice site motifs.
only minor changes in length of coding regions were observed among classes.
phylogenetic relationships among classes based on exon nucleotide sequences
bdelloid alpha tubulin copies are monophyletic with respect to monogonont copies .
the five classes of bdelloid copy types based on intron structure also form groups in the phylogenetic analysis of exon nucleotide sequences.
four classes  are monophyletic and supported by posterior probabilities of  <dig> .
classes  <dig> and  <dig> are sister clades , and together form a monophyletic clade sister to all other bdelloid alpha tubulin sequences .
classes  <dig> and  <dig> are also sister clades .
their position relative to other sequences is unstable among different analyses: in the p <dig> analysis they apparently derive from within the paraphyletic grade of class  <dig> copies , whereas in the amino acid tree they are sister clade to the class  <dig> and  <dig> clade .
class  <dig> copies are monophyletic  in the amino acid tree but paraphyletic in the nucleotide tree.
all four species have class  <dig> and class  <dig> copies, whereas a. ricciae lacks class  <dig>  m. quadricornifera lacks class  <dig>  p. flaviceps lacks class  <dig> copies .
sequence divergence among copies
across all classes, there were fifteen cases of pairs of very similar copies  belonging to the same species .
this is consistent with there being pairs of alpha tubulin copies homogenized by gene conversion as found for other bdelloid genes.
one case of three similar sequences in a. vaga class  <dig> was also found and in this cluster, one copy has gained an intron lacking in the other two.
class  <dig> lacked any clusters of very similar sequences.
the groups of similar copies within species were always separated from their next closest sequence by >10% pairwise sequence divergence.
we refer to copies differing by more than 10% sequence divergence as divergent copy types .
classes  <dig>   <dig> and  <dig> only contained one copy type of alpha tubulin per species .
class  <dig> copies included three  or four  divergent copy types in each species and class  <dig> contained two copy types per species.
the nearest divergent copy type to each sequence  was always from a different species except for two pairs of divergent p. flaviceps class  <dig> copies.
in all cases where both a. vaga and a. ricciae are present in a mini-species tree of copy types, contrary to classical taxonomy, a. vaga appears more closely related to m. quadricornifera than to a. ricciae.
amino acid variation and tests for selection
comparison of the nucleotide and amino acid trees revealed high amino acid divergence between classes compared to nucleotide divergence between them : classes are separated by long amino acid branches relative to branch lengths within classes.
this pattern is particularly striking in class  <dig>  in which amino acid sequences differ by less than  <dig> % whereas the nucleotide sequences have diverged to yield highly supported clades.
class  <dig> copies are similar to monogonont alpha tubulin copies .
the branch-class model in paml optimized a higher average dn/ds ratio between the five classes of alpha tubulin copies than within them, and a significantly higher dn/ds ratio within classes  <dig>   <dig>   <dig> and  <dig> than within class  <dig> .
this rejects the null model that all copies are under the same functional constraints, i.e.
targeted by purifying selection towards the same amino acid sequence.
instead, it is consistent with each class being subjected to purifying selection, but that different classes have diverged in amino acid sequence more than expected based on the strength of purifying selection within each class.
within classes  <dig>   <dig> and  <dig>  amino acid divergence occurs primarily between species, whereas in class  <dig> there are multiple copies within species that differ by  <dig>  to 7% in amino acid sequence.
the branch-site model confirmed these findings.
an estimated 26% of codons belonged to a set having average dn/ds ratios of  <dig> ± <dig>  on within class branches but of  <dig> ± <dig>  on between class branches.
the improved aic of  <dig> compared to the best branch-class model  indicates significant heterogeneity among codons in how they differ between branch classes.
there is no evidence for positive selection at these sites according to the standard interpretation  <cit> , because their average dn/ds ratio is optimized as  <dig> , which equates to neutral evolution.
there is evidence, however, that a quarter of codons that are under strong purifying selection within classes have diverged in amino acid sequence between classes at an average rate equivalent to synonymous sites.
twenty-two codons were shown by bayes empirical bayes analysis to have > <dig>  posterior probability of belonging to this class .
adineta ricciae transcriptome search
searching the a. ricciae transcriptome using blast yielded  <dig> hits for alpha tubulin that varied in length from  <dig> to  <dig> bases .
transcriptome sequences were all less than 3% divergent in nucleotide sequence from their nearest cloned and sequenced copy .
a distinct clade of four copies of class  <dig> alpha tubulin retrieved by cloning and sequencing was not represented by similar transcriptome sequence, possibly because those copies were not expressed in the experimental conditions used for isolating rna for the transcriptome sequencing.
one other copy, from class  <dig>  also lacked a transcriptome match, but this copy was similar in sequence to other copies for which transcriptome matches were obtained.
together, we conclude that  <dig> out of  <dig> divergent copy types of a. ricciae alpha tubulin were expressed in laboratory conditions and  <dig> out of all  <dig> copies  were found.
predicted protein chemistry of alpha tubulin classes
classes based on intron structure differed significantly in predicted chemistry.
the predictive discriminant analysis correctly identified copies of classes  <dig>   <dig> and  <dig> of alpha tubulin based on the three indices of protein chemistry .
class  <dig> and  <dig> could be distinguished from the other classes but not from each other using the three indices of protein chemistry we considered .
class  <dig> and  <dig> copies are most stable and have the lowest isoelectric point; class  <dig> copies are significantly more hydrophobic  than copies of any other class ; class  <dig> copies have the highest instability index ; and class  <dig> copies are most hydrophilic , have the highest isoelectric point  and high instability .
the same qualitative differences are found for the region of protein at the beta tubulin interface.
mean mapp scores of the  <dig> amino acid positions identified as significantly more divergent between classes than within them  were significantly higher than mean mapp scores of a random selection of  <dig> codons across the alignment , and also significantly higher than mean mapp scores of other polymorphic sites .
amino acid divergence in relation to protein structure
of the  <dig> amino acid positions evolving faster between classes than within them identified by the beb analysis, none were in regions identified as functional ; there were significantly fewer of these codons in functional regions than would be expected by chance .
the gtp binding nonexchangeable site , notably, was completely conserved across all copies in all species, and the conservation of the amino acids that interact with the nucleotides in gtp was 96%, 100% and 100% respectively.
by contrast, the number of amino acid positions identified by the beb analysis that fell within structural regions of the protein was not significantly different from random expectation .
discussion
assuming that bdelloid rotifers are tetraploid and inherited alpha tubulin copies from a diploid common ancestor with monogononts, we expected bdelloids to have twice as many copies of alpha tubulin as diploid monogononts, or fewer than this if monogononts have copies with sex-specific functions that were lost in bdelloids.
instead, we found that alpha tubulin has undergone proliferation in bdelloid rotifers relative to the one or two copies found in monogonont rotifers, and beyond numbers expected solely based on degenerate tetraploidy.
each species contains between  <dig> and  <dig> copies that can be divided into five main classes supported both by exon phylogeny and intron structure.
roughly half the copies belong to a single class  and have amino acid sequences conserved both within the class and relative to copies present in monogonont rotifers.
within class  <dig>  copies fall into four distinct groups  in all four species.
the remaining alpha tubulin copies belong to four classes that have diverged significantly in intron structure, amino acid sequence and, in the case of classes  <dig> to  <dig>  in predicted biochemical properties.
the four copy types within class  <dig> and the four divergent classes  <dig> to  <dig> together make eight divergent copy types per species.
because all five classes  are present in species from both sampled families of bdelloids, we propose that bdelloids inherited these eight divergent copy types from their common ancestor followed by the loss of some copy types in some species.
the repeated observation of pairs of similar sequences  of a given copy type in the same species is consistent with there being pairs of copies maintained as templates for dna repair through gene conversion, as proposed to occur between collinear pairs in other bdelloid genome regions  <cit> .
an alternative explanation would be recent duplication of copies within species, but this would be unlikely to generate a widespread pattern of recent divergence between similar copies across classes, unless a similar event had happened independently in all four species.
cases in which only a single copy of a divergent copy type is present could be due to very recent conversion generating identical copies, or to the loss of some duplicate copies.
only in one case were more than two similar copies present within the same species , and two of these sequences were just outside our error threshold, so perhaps reflect sequencing error.
given the presence of paired copies within species’ copy types, we therefore propose that bdelloids inherited sixteen copies of alpha tubulin from their common ancestor: eight divergent copy types, each present as a pair of similar sequences.
some copies were later lost or homogenized in particular lineages leading to the observed numbers of copies in the sampled species.
alternative scenarios invoking duplication events since the species diverged cannot readily explain why all four species share these patterns, unless copies originating independently in separate species were either transferred horizontally to other species - which seem extremely unlikely, even with bdelloids’ known ability to take up foreign dna - or had converged in both amino acid sequence and intron structure, again unlikely.
our results provide robust evidence for numbers of copies and likely inheritance from their common ancestor, yet reconstructing the ancestry of alpha tubulin copies using phylogenetics proved challenging for two reasons.
first, we detected major differences in base composition among sequences, which we attempted to deal with using models of composition heterogeneity in p <dig>  second, functional divergence between different classes has led to extreme rate heterogeneity and departures from the assumptions of neutral models of sequence evolution.
for example, the relationships of class  <dig> copies relative to classes  <dig> to  <dig> were weakly supported and unstable among alternative reconstructions.
the preferred tree using five base composition vectors had classes  <dig> and  <dig> as sister to one a. vaga class  <dig> copy type, which cannot reflect the true ancestry of the genes for reasons argued above .
the amino acid tree does have class  <dig> copies as monophyletic and classes  <dig> to  <dig> copies as sister clade to them, which is consistent with inheritance of multiple copies from a common ancestor.
however, the amino acid tree lacks resolution to infer ancestry among class  <dig> copies.
together, these problems prevented the use of formal genealogical methods for reconstructing patterns of gene duplication and loss .
also, the homogenizing effect of gene conversion is not taken into account in current approaches, and therefore the existence of recent pairs of copies within species would have been inferred as multiple duplications occurring in each species, rather than ancestrally inherited copies maintained by gene conversion .
our dataset presents a useful test case for developers of more complex models of sequence evolution.
most monogonont species had one copy of alpha tubulin, which likely reflects a pair of identical alleles in these sexually reproducing diploids.
this low copy number eliminates the possibility of sex-specific roles for alpha tubulin in most monogononts or that multiple copies in bdelloids arose from formerly sex-specific copies in their sexual ancestor.
assuming that the background level of ploidy in bdelloids conforms to the tetraploidy as demonstrated for two separate genome regions and in both bdelloid families we studied here  <cit> , then if the ancestor contained four copies of alpha tubulin , a further two rounds of duplication of all genes are required to explain our findings.
one possible mechanism would be tandem duplication, possibly facilitated by dna breakage and repair during desiccation.
tandem duplication has been observed for alpha tubulin in species ranging from trypanosoma brucei <cit>  to zea mays <cit> .
if the ancestral sexual progenitor of bdelloids had two distinct copies of alpha tubulin, as found in the b. plicatilis transcriptome, then only one round of additional duplication events would be needed to explain our findings.
why might bdelloids have so many copies of alpha tubulin?
first, multiple identical copies might provide redundancy and increase the rate of production of proteins.
gene duplicates have been shown to provide robustness against null mutations in saccharomyces cerevisiae <cit> .
second, different copies might be expressed in different tissues or at different developmental stages.
for example, bombyx mori and drosophila melanogaster alpha tubulin isotypes fall into classes that are ubiquitously expressed throughout the body at all life stages, classes that are expressed only at certain stages of the life-cycle, and classes expressed only in certain tissues, such as the testes or ovaries  <cit> .
third, different copies might have specialist functions in different cellular processes or locations.
modifications in amino acid sequence can result in changes in the stability and kinetics of microtubule assemblies and in the binding of different microtubule associated proteins  <cit> .
fourth, different copies might provide robustness against environmental variation, by maintaining microtubule function under extreme environmental conditions.
for example nine copies of alpha tubulin expressed in brain tissue of the antarctic fish notothenia coriiceps vary in their stability and polymerization rates at extreme temperatures  <cit> .
which of these mechanisms might explain alpha tubulin proliferation in bdelloids?
class  <dig> copies have conserved protein structure and predicted biochemical properties shared with alpha tubulin copies in monogononts.
yet, each bdelloid species has four divergent types of this class of copies, plus additional pairs of similar copies.
interestingly, only two of the six class  <dig> copies were recovered in the transcriptome of a. ricciae pooled from several laboratory cultures both in hydrated and desiccated conditions.
the missing copies might be expressed in small amounts or only in life-stages  or particular environments  that were not sampled in the experiments.
because no biochemical differences were detected among class  <dig> copies, these copies might be truly functionally redundant.
functional redundancy of housekeeping genes might provide extra protection in the face of dna breakage by increasing the probability that at least one copy still functions after desiccation.
alternatively, multiple copies might provide the machinery to enhance rapid production of protein during rehydration.
in arabidopsis thaliana alpha tubulin is up-regulated upon rehydration after a period of dehydration  <cit> .
tubulin has also been shown to form aggregates and dense networks in brassica napus seedlings during dehydration  <cit> .
it is therefore possible that having multiple functionally redundant tubulin copies might enable rapid up-regulation of tubulin expression at important times.
in contrast, the remaining classes are divergent in amino acid sequence, and those present in a. ricciae were all recovered in the transcriptome and hence had been expressed under laboratory conditions.
there was significant evidence for excess amino acid divergence between classes relative to the strong purifying selection within them.
this is consistent with a standard model of functional specialization after gene duplication in which neo- or sub- functionalized copies evolve, perhaps following a period of relaxed selection, and are then subject to purifying selection.
none of the amino acid positions identified as significantly divergent between classes are within regions identified a priori as functionally important for alpha tubulin, and very few are in structural regions, indicating that selection has not acted to alter fundamentally the functioning of the copies as microtubule proteins.
these positions did, however, have significantly higher mapp scores than random selections of other polymorphic codons, which indicates that they are predicted to have larger consequences for protein function than randomly selected amino acid positions.
this finding adds further evidence that divergence at these sites is associated with functional divergence between classes.
class  <dig> alpha tubulin copies differed in amino acid sequence but have the same predicted chemistry as class  <dig> for the parameters examined here.
however, this does not preclude a different function, as small amino acid changes can alter the function of tubulin  <cit> .
classes  <dig> to  <dig> have diverged significantly both in amino acid sequence and in predicted protein chemistry, and are therefore likely to play different roles.
increased hydrophobicity of the domain that interacts with beta tubulin is predicted to increase microtubule stability  <cit> .
although there was no evidence of excess amino acid divergence in this region, contact surface hydrophobicity did differ significantly between classes and class  <dig> copies are predicted to form the most stable microtubule arrays and class  <dig> the least.
amino acid variation can also be involved in sub-cellular localization and cellular function.
variation in isoelectric point has been noted in tubulin copies with different sub-cellular distributions  <cit> , so the significantly higher isoelectric point in class  <dig> relative to the others might point to a variation in sub-cellular localization or function between classes.
future experiments are needed to test these alternative hypotheses for the function of divergent copies.
another striking difference between classes was in intron structure.
in contrast to the lack of introns in monogonont alpha tubulin, many introns are present and their number and location varies markedly between classes.
this finding matches broader evidence that intron loss and gain occurs at a significantly faster rate between paralogs than between orthologs  <cit> .
genetic manipulation experiments in rice and yeast have shown that introns influence the timing and cellular location of expression of alpha tubulin  <cit> .
such experiments are not feasible at present in bdelloid rotifers, for which no method of genetic modification is currently available.
however, the significant pattern of divergence of intron structure between classes, and conservation within classes, adds further evidence for functional divergence among classes.
CONCLUSIONS
alpha tubulin has expanded into multiple copies in bdelloid rotifers compared to bdelloids’ nearest relatives, the monogonont rotifers.
the copies differ in predicted functions: some provide redundancy for core functions similar to those in monogononts, whereas others have diverged in amino acid sequence and in predicted chemistry.
functional specialization of multiple copies might provide phenotypic flexibility to allow bdelloids to thrive in a range of different environments, which has been argued as a potential mechanism to compensate for their presumed slow rates of genetic adaptation  <cit> .
however, expansion of gene families potentially reduces the potential efficiency of selection in bdelloid rotifers even further, by increasing the number of linked genes and the potential for interference among clones harboring alternative beneficial mutations  <cit> .
future work needs to link characterization of gene families like our study with understanding of the functional consequences of genetic variants for bdelloids in natural environments.
abbreviations
aic, akaike information criterion; alrt, approximate likelihood ratio test; beb, bayes empirical bayes; blast, basic local alignment search tool; dn/ds, ratio ; est, expressed sequence tag; g, gamma-distributed rate variation; gravy index, the grand average of hydropathy index; gtr, general time-reversible; i, proportion of invariable sites; lea, late embryogenesis abundant; mapp, multivariate analysis of protein polymorphism; maps, microtubule associated proteins; mcmc, markov chain monte carlo; ndch, node-discrete composition heterogeneity; pi, isoelectric point; pp, posterior probability.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
ie, tb, df, eh, cb and at wrote manuscript.
ie, ef, df and eh performed the labwork.
ie and tb analysed the data.
cb, ac, gm and at provided the draft transcriptome.
all authors read and approved the final manuscript.
supplementary material
additional file 1
average gc content at the third codon position for sequences belonging to each species.
bp=brachionus plicatilis ; pf =p.
flaviceps; av=a.
vaga; ar=a.
ricciae; mq=m.
quadricornifera.
kruskal-wallis test for variation among species: chi-squared =  <dig> , df =  <dig>  p< <dig> .
click here for file
additional file 2
model fit comparisons.
model fit comparisons.
log marginal likelihood estimated from posterior samples using equation  in newton and raftery  implemented using func.newtonraftery94_eqn <dig>  <cit>  in p <dig> .
click here for file
additional file 3
alpha tubulin gene tree based on exon nucleotides including all cloned copies from all bdelloid species sequenced.
alrt support values are shown and the scale bar is in units of substitutions per site.
click here for file
additional file 4
alpha tubulin gene tree foradineta ricciaeincluding cloned and sequenced copies and copies found in the transcriptome.
alpha tubulin gene tree for adineta ricciae including both cloned and sequenced copies  and the copies found by searching the transcriptome using blastn .
click here for file
additional file 5
average dn/ds  obtained by sliding window analysis of the alpha tubulin alignment across classes.
window size =  <dig>  step size =  <dig>  colours along the x-axis show functional and structural regions: red = alpha helix; blue = beta sheet; green = functional regions .
arrows = codons identified using bayes empirical bayes as experiencing divergent selection .
click here for file
acknowledgements
this work was supported by bbsrc grant bb/f020562/ <dig> to at and tgb, a european research council advanced investigator grant  <dig> to at, a nerc studentship to ie, a marie curie fellowship to df, and a royal society dorothy hodgkin fellowship and european research council grant  <dig> ‘genovir’ to eah.

 bdelloid rotifers are microscopic aquatic invertebrates of particular interest for two features of their lifestyles. first, they have persisted for tens of millions of years and diversified into nearly recognized species, despite the apparent absence of meiosis and sex. second, adult bdelloids can survive essentially complete desiccation in a state known as anhydrobiosis. moreover, bdelloid genomes show a number of unusual features that might relate to their long-term asexuality and ability to survive desiccation. understanding these features can shed light not only on bdelloids’ unusual lifestyle, but also on the range of mechanisms by which animal genomes can evolve. in efforts to understand bdelloid genome organization, much interest has focused on gene copy numbers. the first work on bdelloid genomes found that bdelloids possessed divergent copies of genes present as a single diploid locus in their closest sexual relatives, the monogonont rotifers. it was hypothesized that these were former alleles that had diverged following the loss of meiosis and recombination, a process called the meselson effect. later work suggested that putative former alleles of lea protein had diverged so much as to adopt different biochemical roles protecting bdelloid cells from desiccation. although functional divergence of gene copies is not restricted to asexual lineages, it was argued that functional divergence of former alleles might have contributed towards adaptation to desiccating environments. subsequent work has shown that bdelloid genome complexity cannot be explained simply through divergence of former diploid alleles: bdelloids are degenerate tetraploids. sequencing the genomic regions surrounding four copies of hsp revealed that the chromosomes comprise two collinear pairs. between pairs, divergence is so high that they only share a handful of genes, whereas within collinear pairs divergence is much lower and there are regions of complete identity. the same pattern was found in other genomic regions and in different bdelloid families. it was proposed that bdelloids use pairs of homologous chromosomes as templates to repair dna breaks caused during desiccation, and in the process sequences are homogenised at broken regions by gene conversion. this provides an explanation for the patchy nature of divergence within collinear pairs and for bdelloids’ uncommon ability to repair dna breaks caused by desiccation and ionizing radiation. the original cause of tetraploidy is unknown but could have involved either a genome-wide duplication event or hybridization between sexual ancestors – the latter being commonly associated with asexuality in many organisms. together, these findings make bdelloids an interesting case study for understanding opposing forces on gene copy evolution. although divergent copies of bdelloid genes can no longer be interpreted simply as former alleles, it remains of interest to know how often divergent copies of genes have been maintained and diverged in function. in principle, the maintenance of functionally divergent copies of genes could provide bdelloids with an enhanced genetic repertoire and, perhaps, increased phenotypic plasticity. here, we explore copy evolution in alpha tubulin and ask whether divergent copies have been maintained in bdelloid genomes. alpha tubulin forms microtubules with beta tubulin and provides a range of functions in cellular movement and cell division. previous studies of copy variation of alpha tubulin in other organisms revealed contrasting patterns making it an interesting gene to study in bdelloids. in many animals, it is present in low copy number and highly conserved, which has led to its use as a marker for deep phylogenetic relationships. however, the different roles of microtubules have led to functional divergence in many cases. four distinct copies are present in drosophila melanogaster and expressed in different tissues: including ovary and testis-specific forms. if there were male-function or meiosis specific copies in sexual rotifers, bdelloids might be expected to have lost those copies or for them to have evolved new functions. duplication and diversification of tubulin genes has also been implicated in adaptation to extreme environments in the antarctic fish notothenia coriiceps, which expresses at least nine different alpha tubulin isotypes that vary in stability and polymerization rates at extremes of temperature. by expressing copies stable at different temperatures, function is maintained despite variation in environmental conditions. in plant seedlings alpha tubulin is up-regulated during both dehydration and rehydration. therefore, bdelloids might make use of their ancestral tetraploid state and maintain alpha tubulin copies to perform specialist functions during cycles of desiccation. we used cloning and preliminary whole transcriptome data to characterize alpha tubulin copies in several monogonont and bdelloid rotifer species. monogononts contained either one or two copies with very similar amino acid sequence. contrary to the four or eight copies expected if bdelloids retained copies from their tetraploid ancestor, we found between and copies in each of four species of bdelloids. phylogenetic relationships indicate that multiple copies were inherited from the common ancestor of bdelloids and some of them were lost or homogenized in some species. approximately half of the copies retained amino acid sequences similar to the monogonont copies, whereas the rest have diverged into four distinct classes that differ in intron structure and biochemical properties predicted by computational analysis. distinct classes have diverged significantly in amino acid sequence relative to the purifying selection observed within classes. the findings do not contradict current understanding of bdelloid genome organization, but show how further duplication events in the ancestor of bdelloids led to proliferation in both conserved and functionally divergent copies of this gene. the possible functional significance of alpha tubulin diversity is discussed in relation to sequence characteristics. methods study species, dna extraction, cloning and sequencing four species of bdelloid rotifer belonging to two bdelloid families were compared to four species of monogonont rotifers also from two families. sequences for one of the monogononts, brachionus plicatilis, were taken from the cdna library of suga et al.. for the other monogonont species dna was extracted from wild-caught individuals using chelex instagene matrix following the manufacturers instructions. genomic dna from bdelloid species was extracted using the dneasy blood and tissue kit. for each bdelloid species dna was extracted from ~ individuals from populations grown from a single female. extractions were performed according to manufacturers instructions, and samples were heated at 56°c for 1 h for protein digestion. dna quality was checked on a 1% agarose gel. pcr amplification was performed using 2 μl template dna and illustra pure-taq ready-to-go pcr beads. alpha tubulin was amplified using the primers toa1f 5′-rgtnggnaaygcntgytggga-3′ and toa1r 5′-ccatnccytcnccnacrtacca-3′, which amplify >90% of the coding region of each gene. the amplification profile consisted of an initial denaturation of 5 min at 94°c, followed by 35 cycles of 45 s at 94°c, 60s at 58°c and 60s at 72°c and a final extension of 5 min at 72°c. products of expected size were cloned using the topo ta cloning kit and alpha tubulin inserts were re-amplified from individual clones. amplicons were sequenced in both directions using the bigdye terminator v cycle sequencing kit and an abi genetic analyzer. sequence data are available in genbank. sequence alignment sequences were confirmed to be alpha tubulin using blastn, and aligned in geneious pro. using muscle followed by manual adjustments. wise was used to define introns, which were manually checked by comparison to alpha tubulin cdna from the monogonont rotifer brachionus plicatilis. sequences used in our analyses contain no premature stop codons. taq polymerase has a baseline error rate of in base pairs. normally, when sequencing from uncloned pcr products this error rate is negligible, as errors occur at random sites across the amplified products and are averaged out upon sequencing. however, cloning to separate different copies of alpha tubulin within each sample introduces greater errors into the sequencing process. based on the estimated error rate of taq of in and the 35 cycles in our pcr reaction, we estimated an overall error rate of % per sequence, which we rounded to % to allow for any additional error in sequencing these long products. consequently, sequences obtained from cloning were considered to be the same if they were <1% different, and collapsed to a single majority rule consensus sequence, thus providing us with a conservative estimate of copy number. plots of the number of distinct copies against the number of clones sequenced were used to assess whether most copies present in the genome had been discovered. the chao estimator was calculated to estimate the potential number of missing copies, based on the assumption that if a few copies have been sampled only once or twice, then it is likely that some copies remain unsampled. phylogenetic analyses relationships among copies were inferred from exon sequences using bayesian markov chain monte carlo analyses implemented in p. we used p because preliminary investigations revealed compositional heterogeneity in our sequences [see additional file introns were excluded from phylogenetic analysis because their extreme variability in sequence and insertion site confounded the evolutionary signal. jmodeltest selected the gtr model with gamma-distributed rate variation and a proportion of invariable sites as the most likely substitution model using the akaike information criterion. a chi-squared test revealed significant compositional homogeneity over the tree, and comparison of the observed compositional heterogeneity among sequences to a null distribution simulated using the fitted gtr + i + g model, confirmed that composition of the observed sequences was significantly heterogeneous. the gtr + i + g model was therefore compared to node-discrete composition heterogeneity models ndch and ndch, which incorporate two and five tree-heterogeneous base-composition vectors respectively [see additional file the preferred model was chosen using bayes factors: a loge of or more was considered as very strong support for improved fit of a more complex model. bayesian mcmc analysis was run with the chosen model and two independent runs of four chains each for generations. posterior probabilities were calculated for support. a second tree was reconstructed, using maximum likelihood implemented in phyml., from the translated amino acid alignment. prottest selected the jtt model of amino acid substitution with gamma-distributed rate variation. statistical support was assessed using alrt non-parametric branch support. trees were rooted on the branch connecting bdelloid and monogonont sequences. tests for selection we estimated the ratio of the rates of non-synonymous and synonymous substitutions using likelihood methods implemented in paml version. first, we fitted a model with a single dn/ds ratio across all bdelloid copies of alpha tubulin. second, to test for divergent selection acting on different classes of alpha tubulin, we fitted a model with a different dn/ds ratio on branches connecting different classes versus on branches connecting sequences of the same class. third, we fitted a branch-site model to test for the presence of codons showing purifying selection on within-class branches but evolving either neutrally or by positive selection on between-class branches. bayes empirical bayes inference was used to identify codons belonging to this set. we predict that the branch-site model should be preferred if amino acid divergence between the different classes has evolved through selection for functional specialization of some parts but not all of the protein. the best model was selected using the akaike information criterion. models were run using the average nucleotide frequencies at each codon position to specify equilibrium codon frequencies. to check for sensitivity to the observed compositional heterogeneity, we repeated the analyses using free parameters for each codon frequency. all inferences from model results were qualitatively unchanged over the simpler versions. adineta ricciae transcriptome searches to further assess the completeness of sampling of copies by cloning and sequencing, and to establish if multiple copies were expressed under experimental conditions, we searched an expressed sequence tag library and draft transcriptome sequences of a. ricciae for alpha tubulin copies. the transcriptome was generated by paired-end sequencing of ~200bp cdna fragments using illumina parallel sequencing. the resulting million 76 bp reads were assembled with the clc-bio assembler to give a library of transcripts of size range 118-3674 bp. more details are provided in. the databases were searched by blastn for matches to the sequences obtained by cloning. protein chemical analyses three parameters indicative of predicted protein chemistry, and hence potentially of changes in protein function, were calculated for each alpha tubulin amino acid sequence using protparam. i) the grand average of hydropathy index indicates the hydropathy of proteins: a positive gravy score means the protein is hydrophobic and a negative score hydrophilic. detrich hypothesized that increasing contact surface hydrophobicity might enhance polymer stabilizing interactions in microtubules, so a more hydrophobic gravy score could indicate a more stable microtubule. ii) isoelectric point is the ph at which a protein carries no net electrical charge, which can affect its solubility at a given ph. zambito et al. found a wide variation in the pi of rat brain tubulin isotypes, and showed that some variants are associated with particular cellular compartments. variation in pi of our alpha tubulin copies may therefore imply variation in cellular function. iii) the instability index provides an estimate of a protein’s folding stability in vitro, based on estimated instabilities of the protein’s dipeptides. a protein with an instability index of less than is predicted as stable, a value above predicts that the protein may be unstable, and therefore less viable as a functional protein. predictive discriminant analysis was performed using the predict.lda function of the mass library version - in the r statistical programming language r. to see if classes of copies based on intron structure are predicted by these three protein chemical properties, i.e. whether copies with different intron structure have divergent chemical properties. multivariate analysis of protein polymorphism scores of physiochemical impact were estimated for all codons, providing a measure of constraint at each site. to determine if sites that diverged significantly between classes had significantly different mapp scores to the rest of the protein, the mean mapp score at codons identified by the beb analysis to be significantly divergent between classes was compared to replicates of the mean mapp score at the same number of random sites drawn from a set of all sites in the protein. to determine if sites identified by the beb analysis had significantly different mapp scores to polymorphic sites not evolving faster between classes than within them, the mean mapp score at the codons identified by the beb analysis was compared to replicates of the mean mapp score at the same number of random sites drawn from a set of all polymorphic sites. protein structural analyses translated protein sequences were mapped to the structural and functional regions of alpha tubulin identified by. regions of alpha helices and beta sheets were classified as structural. regions falling into the following four types were classified as functional: the interface with beta tubulin; the regions interacting with gtp; the tubulin nonexchangeable site; and the c-terminal helices where maps and motor proteins usually bind. to determine whether amino acid divergence was concentrated in particular structural regions, amino acid positions identified by the beb analysis to be significantly divergent between classes were mapped onto the identified functional and structural regions. monte carlo tests were used to test for significant concentration of divergent amino acids either in structural or in functional regions of the genome.  number of copies of alpha tubulin genes from cloning and sequencing the b. plicatilis est database of suga et al. contained two alpha tubulin copies that were 13% divergent at the nucleotide level but only differed in three amino acids. amplifying and sequencing alpha tubulin in the three selected monogonont species revealed just one copy of alpha tubulin in each species. pairwise amino acid divergence between all monogonont alpha tubulin sequences was only %. both b. plicatilis copies contained our primer binding sites and should therefore have been amplified using our primers, if they had been present in the other monogonont species. cloning and sequencing of bdelloid alpha tubulin yielded between and different copies in each species. this is a conservative estimate because sequences were only considered unique when they varied by >1% to account for % pcr fidelity of the cloned copies. including all sequences obtained by cloning would increase the number of copies per species to almost as many as the number sequenced, but would not alter our broad findings beyond increasing the number of very similar copies within species [see additional file saturation curves showed that the number of copies reached a plateau against the number of clones sequenced in each species. the chao estimator indicated that all copies were likely to have been sampled in all species. intron and exon structure no introns were present in the monogonont sequences. in contrast, each bdelloid alpha tubulin copy contained between and introns located across different sites. the copies fall into five major classes based on differences in intron content and location, with some additional losses or gains of single introns in particular species leading to a total of eight different intron patterns. most introns were between and nucleotides in length, but two introns were consistently longer across copies and species. all introns present in two or more classes were always found in class and there were no unique introns in class class only differed from class by the loss of one intron whereas classes and each contained one or more unique introns. all seven introns present in the class copies shared a ggt motif at the 5’ splice site across all the classes. one intron unique to class had the same 5’ splice site as the introns in class copies and represents a putative translocation of an adjacent intron that is present in all the other classes but lacking in class all other introns had different 5′ splice site motifs. only minor changes in length of coding regions were observed among classes. phylogenetic relationships among classes based on exon nucleotide sequences bdelloid alpha tubulin copies are monophyletic with respect to monogonont copies. the five classes of bdelloid copy types based on intron structure also form groups in the phylogenetic analysis of exon nucleotide sequences. four classes are monophyletic and supported by posterior probabilities of. classes and are sister clades, and together form a monophyletic clade sister to all other bdelloid alpha tubulin sequences. classes and are also sister clades. their position relative to other sequences is unstable among different analyses: in the p analysis they apparently derive from within the paraphyletic grade of class copies, whereas in the amino acid tree they are sister clade to the class and clade. class copies are monophyletic in the amino acid tree but paraphyletic in the nucleotide tree. all four species have class and class copies, whereas a. ricciae lacks class m. quadricornifera lacks class p. flaviceps lacks class copies. sequence divergence among copies across all classes, there were fifteen cases of pairs of very similar copies belonging to the same species. this is consistent with there being pairs of alpha tubulin copies homogenized by gene conversion as found for other bdelloid genes. one case of three similar sequences in a. vaga class was also found and in this cluster, one copy has gained an intron lacking in the other two. class lacked any clusters of very similar sequences. the groups of similar copies within species were always separated from their next closest sequence by >10% pairwise sequence divergence. we refer to copies differing by more than 10% sequence divergence as divergent copy types. classes and only contained one copy type of alpha tubulin per species. class copies included three or four divergent copy types in each species and class contained two copy types per species. the nearest divergent copy type to each sequence was always from a different species except for two pairs of divergent p. flaviceps class copies. in all cases where both a. vaga and a. ricciae are present in a mini-species tree of copy types, contrary to classical taxonomy, a. vaga appears more closely related to m. quadricornifera than to a. ricciae. amino acid variation and tests for selection comparison of the nucleotide and amino acid trees revealed high amino acid divergence between classes compared to nucleotide divergence between them: classes are separated by long amino acid branches relative to branch lengths within classes. this pattern is particularly striking in class in which amino acid sequences differ by less than % whereas the nucleotide sequences have diverged to yield highly supported clades. class copies are similar to monogonont alpha tubulin copies. the branch-class model in paml optimized a higher average dn/ds ratio between the five classes of alpha tubulin copies than within them, and a significantly higher dn/ds ratio within classes and than within class. this rejects the null model that all copies are under the same functional constraints, i.e. targeted by purifying selection towards the same amino acid sequence. instead, it is consistent with each class being subjected to purifying selection, but that different classes have diverged in amino acid sequence more than expected based on the strength of purifying selection within each class. within classes and amino acid divergence occurs primarily between species, whereas in class there are multiple copies within species that differ by to 7% in amino acid sequence. the branch-site model confirmed these findings. an estimated 26% of codons belonged to a set having average dn/ds ratios of ± on within class branches but of ± on between class branches. the improved aic of compared to the best branch-class model indicates significant heterogeneity among codons in how they differ between branch classes. there is no evidence for positive selection at these sites according to the standard interpretation, because their average dn/ds ratio is optimized as, which equates to neutral evolution. there is evidence, however, that a quarter of codons that are under strong purifying selection within classes have diverged in amino acid sequence between classes at an average rate equivalent to synonymous sites. twenty-two codons were shown by bayes empirical bayes analysis to have > posterior probability of belonging to this class. adineta ricciae transcriptome search searching the a. ricciae transcriptome using blast yielded hits for alpha tubulin that varied in length from to bases. transcriptome sequences were all less than 3% divergent in nucleotide sequence from their nearest cloned and sequenced copy. a distinct clade of four copies of class alpha tubulin retrieved by cloning and sequencing was not represented by similar transcriptome sequence, possibly because those copies were not expressed in the experimental conditions used for isolating rna for the transcriptome sequencing. one other copy, from class also lacked a transcriptome match, but this copy was similar in sequence to other copies for which transcriptome matches were obtained. together, we conclude that out of divergent copy types of a. ricciae alpha tubulin were expressed in laboratory conditions and out of all copies were found. predicted protein chemistry of alpha tubulin classes classes based on intron structure differed significantly in predicted chemistry. the predictive discriminant analysis correctly identified copies of classes and of alpha tubulin based on the three indices of protein chemistry. class and could be distinguished from the other classes but not from each other using the three indices of protein chemistry we considered. class and copies are most stable and have the lowest isoelectric point; class copies are significantly more hydrophobic than copies of any other class; class copies have the highest instability index; and class copies are most hydrophilic, have the highest isoelectric point and high instability. the same qualitative differences are found for the region of protein at the beta tubulin interface. mean mapp scores of the amino acid positions identified as significantly more divergent between classes than within them were significantly higher than mean mapp scores of a random selection of codons across the alignment, and also significantly higher than mean mapp scores of other polymorphic sites. amino acid divergence in relation to protein structure of the amino acid positions evolving faster between classes than within them identified by the beb analysis, none were in regions identified as functional; there were significantly fewer of these codons in functional regions than would be expected by chance. the gtp binding nonexchangeable site, notably, was completely conserved across all copies in all species, and the conservation of the amino acids that interact with the nucleotides in gtp was 96%, 100% and 100% respectively. by contrast, the number of amino acid positions identified by the beb analysis that fell within structural regions of the protein was not significantly different from random expectation. discussion assuming that bdelloid rotifers are tetraploid and inherited alpha tubulin copies from a diploid common ancestor with monogononts, we expected bdelloids to have twice as many copies of alpha tubulin as diploid monogononts, or fewer than this if monogononts have copies with sex-specific functions that were lost in bdelloids. instead, we found that alpha tubulin has undergone proliferation in bdelloid rotifers relative to the one or two copies found in monogonont rotifers, and beyond numbers expected solely based on degenerate tetraploidy. each species contains between and copies that can be divided into five main classes supported both by exon phylogeny and intron structure. roughly half the copies belong to a single class and have amino acid sequences conserved both within the class and relative to copies present in monogonont rotifers. within class copies fall into four distinct groups in all four species. the remaining alpha tubulin copies belong to four classes that have diverged significantly in intron structure, amino acid sequence and, in the case of classes to in predicted biochemical properties. the four copy types within class and the four divergent classes to together make eight divergent copy types per species. because all five classes are present in species from both sampled families of bdelloids, we propose that bdelloids inherited these eight divergent copy types from their common ancestor followed by the loss of some copy types in some species. the repeated observation of pairs of similar sequences of a given copy type in the same species is consistent with there being pairs of copies maintained as templates for dna repair through gene conversion, as proposed to occur between collinear pairs in other bdelloid genome regions. an alternative explanation would be recent duplication of copies within species, but this would be unlikely to generate a widespread pattern of recent divergence between similar copies across classes, unless a similar event had happened independently in all four species. cases in which only a single copy of a divergent copy type is present could be due to very recent conversion generating identical copies, or to the loss of some duplicate copies. only in one case were more than two similar copies present within the same species, and two of these sequences were just outside our error threshold, so perhaps reflect sequencing error. given the presence of paired copies within species’ copy types, we therefore propose that bdelloids inherited sixteen copies of alpha tubulin from their common ancestor: eight divergent copy types, each present as a pair of similar sequences. some copies were later lost or homogenized in particular lineages leading to the observed numbers of copies in the sampled species. alternative scenarios invoking duplication events since the species diverged cannot readily explain why all four species share these patterns, unless copies originating independently in separate species were either transferred horizontally to other species - which seem extremely unlikely, even with bdelloids’ known ability to take up foreign dna - or had converged in both amino acid sequence and intron structure, again unlikely. our results provide robust evidence for numbers of copies and likely inheritance from their common ancestor, yet reconstructing the ancestry of alpha tubulin copies using phylogenetics proved challenging for two reasons. first, we detected major differences in base composition among sequences, which we attempted to deal with using models of composition heterogeneity in p second, functional divergence between different classes has led to extreme rate heterogeneity and departures from the assumptions of neutral models of sequence evolution. for example, the relationships of class copies relative to classes to were weakly supported and unstable among alternative reconstructions. the preferred tree using five base composition vectors had classes and as sister to one a. vaga class copy type, which cannot reflect the true ancestry of the genes for reasons argued above. the amino acid tree does have class copies as monophyletic and classes to copies as sister clade to them, which is consistent with inheritance of multiple copies from a common ancestor. however, the amino acid tree lacks resolution to infer ancestry among class copies. together, these problems prevented the use of formal genealogical methods for reconstructing patterns of gene duplication and loss. also, the homogenizing effect of gene conversion is not taken into account in current approaches, and therefore the existence of recent pairs of copies within species would have been inferred as multiple duplications occurring in each species, rather than ancestrally inherited copies maintained by gene conversion. our dataset presents a useful test case for developers of more complex models of sequence evolution. most monogonont species had one copy of alpha tubulin, which likely reflects a pair of identical alleles in these sexually reproducing diploids. this low copy number eliminates the possibility of sex-specific roles for alpha tubulin in most monogononts or that multiple copies in bdelloids arose from formerly sex-specific copies in their sexual ancestor. assuming that the background level of ploidy in bdelloids conforms to the tetraploidy as demonstrated for two separate genome regions and in both bdelloid families we studied here, then if the ancestor contained four copies of alpha tubulin, a further two rounds of duplication of all genes are required to explain our findings. one possible mechanism would be tandem duplication, possibly facilitated by dna breakage and repair during desiccation. tandem duplication has been observed for alpha tubulin in species ranging from trypanosoma brucei to zea mays. if the ancestral sexual progenitor of bdelloids had two distinct copies of alpha tubulin, as found in the b. plicatilis transcriptome, then only one round of additional duplication events would be needed to explain our findings. why might bdelloids have so many copies of alpha tubulin? first, multiple identical copies might provide redundancy and increase the rate of production of proteins. gene duplicates have been shown to provide robustness against null mutations in saccharomyces cerevisiae. second, different copies might be expressed in different tissues or at different developmental stages. for example, bombyx mori and drosophila melanogaster alpha tubulin isotypes fall into classes that are ubiquitously expressed throughout the body at all life stages, classes that are expressed only at certain stages of the life-cycle, and classes expressed only in certain tissues, such as the testes or ovaries. third, different copies might have specialist functions in different cellular processes or locations. modifications in amino acid sequence can result in changes in the stability and kinetics of microtubule assemblies and in the binding of different microtubule associated proteins. fourth, different copies might provide robustness against environmental variation, by maintaining microtubule function under extreme environmental conditions. for example nine copies of alpha tubulin expressed in brain tissue of the antarctic fish notothenia coriiceps vary in their stability and polymerization rates at extreme temperatures. which of these mechanisms might explain alpha tubulin proliferation in bdelloids? class copies have conserved protein structure and predicted biochemical properties shared with alpha tubulin copies in monogononts. yet, each bdelloid species has four divergent types of this class of copies, plus additional pairs of similar copies. interestingly, only two of the six class copies were recovered in the transcriptome of a. ricciae pooled from several laboratory cultures both in hydrated and desiccated conditions. the missing copies might be expressed in small amounts or only in life-stages or particular environments that were not sampled in the experiments. because no biochemical differences were detected among class copies, these copies might be truly functionally redundant. functional redundancy of housekeeping genes might provide extra protection in the face of dna breakage by increasing the probability that at least one copy still functions after desiccation. alternatively, multiple copies might provide the machinery to enhance rapid production of protein during rehydration. in arabidopsis thaliana alpha tubulin is up-regulated upon rehydration after a period of dehydration. tubulin has also been shown to form aggregates and dense networks in brassica napus seedlings during dehydration. it is therefore possible that having multiple functionally redundant tubulin copies might enable rapid up-regulation of tubulin expression at important times. in contrast, the remaining classes are divergent in amino acid sequence, and those present in a. ricciae were all recovered in the transcriptome and hence had been expressed under laboratory conditions. there was significant evidence for excess amino acid divergence between classes relative to the strong purifying selection within them. this is consistent with a standard model of functional specialization after gene duplication in which neo- or sub- functionalized copies evolve, perhaps following a period of relaxed selection, and are then subject to purifying selection. none of the amino acid positions identified as significantly divergent between classes are within regions identified a priori as functionally important for alpha tubulin, and very few are in structural regions, indicating that selection has not acted to alter fundamentally the functioning of the copies as microtubule proteins. these positions did, however, have significantly higher mapp scores than random selections of other polymorphic codons, which indicates that they are predicted to have larger consequences for protein function than randomly selected amino acid positions. this finding adds further evidence that divergence at these sites is associated with functional divergence between classes. class alpha tubulin copies differed in amino acid sequence but have the same predicted chemistry as class for the parameters examined here. however, this does not preclude a different function, as small amino acid changes can alter the function of tubulin. classes to have diverged significantly both in amino acid sequence and in predicted protein chemistry, and are therefore likely to play different roles. increased hydrophobicity of the domain that interacts with beta tubulin is predicted to increase microtubule stability. although there was no evidence of excess amino acid divergence in this region, contact surface hydrophobicity did differ significantly between classes and class copies are predicted to form the most stable microtubule arrays and class the least. amino acid variation can also be involved in sub-cellular localization and cellular function. variation in isoelectric point has been noted in tubulin copies with different sub-cellular distributions, so the significantly higher isoelectric point in class relative to the others might point to a variation in sub-cellular localization or function between classes. future experiments are needed to test these alternative hypotheses for the function of divergent copies. another striking difference between classes was in intron structure. in contrast to the lack of introns in monogonont alpha tubulin, many introns are present and their number and location varies markedly between classes. this finding matches broader evidence that intron loss and gain occurs at a significantly faster rate between paralogs than between orthologs. genetic manipulation experiments in rice and yeast have shown that introns influence the timing and cellular location of expression of alpha tubulin. such experiments are not feasible at present in bdelloid rotifers, for which no method of genetic modification is currently available. however, the significant pattern of divergence of intron structure between classes, and conservation within classes, adds further evidence for functional divergence among classes.  alpha tubulin has expanded into multiple copies in bdelloid rotifers compared to bdelloids’ nearest relatives, the monogonont rotifers. the copies differ in predicted functions: some provide redundancy for core functions similar to those in monogononts, whereas others have diverged in amino acid sequence and in predicted chemistry. functional specialization of multiple copies might provide phenotypic flexibility to allow bdelloids to thrive in a range of different environments, which has been argued as a potential mechanism to compensate for their presumed slow rates of genetic adaptation. however, expansion of gene families potentially reduces the potential efficiency of selection in bdelloid rotifers even further, by increasing the number of linked genes and the potential for interference among clones harboring alternative beneficial mutations. future work needs to link characterization of gene families like our study with understanding of the functional consequences of genetic variants for bdelloids in natural environments. abbreviations aic, akaike information criterion; alrt, approximate likelihood ratio test; beb, bayes empirical bayes; blast, basic local alignment search tool; dn/ds, ratio; est, expressed sequence tag; g, gamma-distributed rate variation; gravy index, the grand average of hydropathy index; gtr, general time-reversible; i, proportion of invariable sites; lea, late embryogenesis abundant; mapp, multivariate analysis of protein polymorphism; maps, microtubule associated proteins; mcmc, markov chain monte carlo; ndch, node-discrete composition heterogeneity; pi, isoelectric point; pp, posterior probability. competing interests the authors declare that they have no competing interests. authors’ contributions ie, tb, df, eh, cb and at wrote manuscript. ie, ef, df and eh performed the labwork. ie and tb analysed the data. cb, ac, gm and at provided the draft transcriptome. all authors read and approved the final manuscript. supplementary material additional file 1 average gc content at the third codon position for sequences belonging to each species. bp=brachionus plicatilis; pf =p. flaviceps; av=a. vaga; ar=a. ricciae; mq=m. quadricornifera. kruskal-wallis test for variation among species: chi-squared =, df = p<. click here for file additional file 2 model fit comparisons. model fit comparisons. log marginal likelihood estimated from posterior samples using equation in newton and raftery implemented using func.newtonraftery94_eqn in p. click here for file additional file 3 alpha tubulin gene tree based on exon nucleotides including all cloned copies from all bdelloid species sequenced. alrt support values are shown and the scale bar is in units of substitutions per site. click here for file additional file 4 alpha tubulin gene tree foradineta ricciaeincluding cloned and sequenced copies and copies found in the transcriptome. alpha tubulin gene tree for adineta ricciae including both cloned and sequenced copies and the copies found by searching the transcriptome using blastn. click here for file additional file 5 average dn/ds obtained by sliding window analysis of the alpha tubulin alignment across classes. window size = step size = colours along the x-axis show functional and structural regions: red = alpha helix; blue = beta sheet; green = functional regions. arrows = codons identified using bayes empirical bayes as experiencing divergent selection. click here for file acknowledgements this work was supported by bbsrc grant bb/f020562/ to at and tgb, a european research council advanced investigator grant to at, a nerc studentship to ie, a marie curie fellowship to df, and a royal society dorothy hodgkin fellowship and european research council grant ‘genovir’ to eah.
alpha tubulin is conserved and normally present in low copy numbers in animals, but multiplication of alpha tubulin copies has occurred in animals adapted to extreme environments, such as cold-adapted antarctic fish. studies into the impact of asexuality and anhydrobiosis on bdelloid genomes have focused on understanding gene copy number. here we investigate copy number and sequence divergence in alpha tubulin. results were verified using transcriptome data from one bdelloid species, adineta ricciae. our findings are consistent with recent evidence that bdelloids are degenerate tetraploids and that functional divergence of ancestral copies of genes has occurred, but show how further duplication events in the ancestor of bdelloids led to proliferation in both conserved and functionally divergent copies of this gene. copies of each class were expressed under the laboratory conditions used to construct the transcriptome. in common with the typical pattern for animals, monogonont rotifers contain either one or two copies of alpha tubulin, but bdelloid species contain between and different copies, distributed across five classes. bdelloid rotifers are microscopic animals that have apparently survived without sex for millions of years and are able to survive desiccation at all life stages through a process called anhydrobiosis. approximately half of the copies form a highly conserved group that vary by only % amino acid pairwise divergence with each other and with the monogonont copies. using cloning and sequencing we compared alpha tubulin copy variation in four species of bdelloid rotifers and four species of monogonont rotifers, which are facultatively sexual and cannot survive desiccation as adults. both of these characteristics are believed to have played a role in shaping several unusual features of bdelloid genomes discovered in recent years. the other copies have divergent amino acid sequences that evolved significantly faster between classes than within them, relative to synonymous changes, and vary in predicted biochemical properties.
 bdelloid rotifers are microscopic animals that have apparently survived without sex for millions of years and are able to survive desiccation at all life stages through a process called anhydrobiosis. both of these characteristics are believed to have played a role in shaping several unusual features of bdelloid genomes discovered in recent years. studies into the impact of asexuality and anhydrobiosis on bdelloid genomes have focused on understanding gene copy number. here we investigate copy number and sequence divergence in alpha tubulin. alpha tubulin is conserved and normally present in low copy numbers in animals, but multiplication of alpha tubulin copies has occurred in animals adapted to extreme environments, such as cold-adapted antarctic fish. using cloning and sequencing we compared alpha tubulin copy variation in four species of bdelloid rotifers and four species of monogonont rotifers, which are facultatively sexual and cannot survive desiccation as adults. results were verified using transcriptome data from one bdelloid species, adineta ricciae.  in common with the typical pattern for animals, monogonont rotifers contain either one or two copies of alpha tubulin, but bdelloid species contain between and different copies, distributed across five classes. approximately half of the copies form a highly conserved group that vary by only % amino acid pairwise divergence with each other and with the monogonont copies. the other copies have divergent amino acid sequences that evolved significantly faster between classes than within them, relative to synonymous changes, and vary in predicted biochemical properties. copies of each class were expressed under the laboratory conditions used to construct the transcriptome.  our findings are consistent with recent evidence that bdelloids are degenerate tetraploids and that functional divergence of ancestral copies of genes has occurred, but show how further duplication events in the ancestor of bdelloids led to proliferation in both conserved and functionally divergent copies of this gene. bdelloid rotifersgene copiestubulinevolution 