Degenerate Tetraploidy Was Established Before Bdelloid Rotifer Families Diverged

Rotifers of Class Bdelloidea are abundant freshwater invertebrates known for their remarkable ability to survive desiccation and their lack of males and meiosis. Sequencing and annotation of approximately 50-kb regions containing the four hsp82 heat shock genes of the bdelloid Philodina roseola , each located on a separate chromosome, have suggested that its genome is that of a degenerate tetraploid. In order to determine whether a similar structure exists in a bdelloid distantly related to P. roseola and if degenerate tetraploidy was established before the two species separated, we sequenced regions containing the hsp82 genes of a bdelloid belonging to a different family, Adineta vaga , and the histone gene clusters of P. roseola and A. vaga . Our ﬁndings are entirely consistent with degenerate tetraploidy and show that it was established before the two bdelloid families diverged and therefore probably before the bdelloid radiation.


Introduction
Rotifers of Class Bdelloidea are common freshwater invertebrates thought to reproduce solely by ameiotic parthenogenesis (Hsu 1956a(Hsu , 1956b;;Normark et al. 2003) and are also unusual in being able to survive and reproduce after desiccation at any life stage (Ricci 1987) and after exposure to extremely high doses of ionizing radiation (Gladyshev and Meselson 2008).The class comprises some 460 described species (Segers 2007) classified among four families (Adinetidae, Habrotrochidae, Philodinavidae, and Philodinidae).The high synonymous nucleotide divergence (K s ) between orthologous genes in bdelloids of different families and the presence of bdelloid remains in 35-40-My-old amber indicate that bdelloids arose tens of millions of years ago (Waggoner and Poinar 1993;Mark Welch and Meselson 2000).
Sequencing of four 45-70-kb regions, each containing one of the four copies of the hsp82 heat shock gene in the bdelloid Philodina roseola (family Philodinidae), has indicated that P. roseola is a degenerate tetraploid (Mark Welch et al. 2008).The four regions constitute a quartet comprised of two colinear gene-rich pairs with only a minority of genes common to both pairs, in the same order and orientation and separated by gene-rich segments present in only one or the other pair.The average nucleotide difference between aligned members of a colinear pair is 4% and the nonweighted average K s between genes within a pair is 0.06 (range 0-0.15).There are tracts of sequence identity or near identity up to several kilobases in length between members of a pair, as could result from recent conversion.There are no conspicuous tracts of identity between members of different colinear pairs, and they can be aligned only within the genes they share, for which the nonweighted average K s is 0.93 (range 0.67-1.29).
Consistent with the above evidence that P. roseola is a degenerate tetraploid, its karyotype, ignoring dot chromosomes, consists of 11 equal-size chromosomes and one chromosome of about twice their length that appears to be an isochromosome (Mark Welch et al. 2004).Taking its two arms as a colinear pair, the P. roseola genome may therefore consist of three quartets, each made up of two colinear pairs.The karyotype of the distantly related bdelloid Adineta vaga (family Adinetidae) consists of 12 equal-size chromosomes, suggesting that it too may be a degenerate tetraploid (Mark Welch and Meselson 1998).
Here, we ask if degenerate tetraploidy is a general characteristic of bdelloid rotifers and, if so, whether it was established before the separation of bdelloid families or independently in different bdelloid families or taxa of lower rank.To address this question, we sequenced 38-67-kb regions containing the hsp82 genes of A. vaga and 25-39-kb regions containing the histone gene clusters of A. vaga and P. roseola.Our results are entirely consistent with degenerate tetraploidy.Moreover, the order and orientation of genes in sequenced regions are largely conserved between species, indicating that degenerate tetraploidy was established before the two families separated.Our findings are consistent with the occurrence of a whole genome duplication followed by numerous segmental deletions before the separation of the two bdelloid families but after the separation of Bdelloidea from rotifers of their sister Class, the Monogononta (Melone et al. 1998;Mark Welch 2000), which are facultatively sexual diploids (Robotti 1975;King and Snell 1977;Mark Welch and Meselson 2000).

Library Construction and Screening
Adineta vaga was cultured and DNA extracted from purified eggs as previously described except CsCl densitygradient purification was replaced by phenol:chloroform extraction (Mark Welch et al. 2004;Gladyshev and Meselson 2008).The fosmid library was constructed with an Epicentre CopyControl Fosmid Library Production Kit (Epicentre, Madison, WI), giving a primary fosmid library in Escherichia coli.The library was arranged onto GeneScreen Plus positively charged nylon membranes (PerkinElmer, Waltham, MA) using a Genetix Q-Bot (Genetix, New Milton, Hampshire, United Kingdom). 32P-labeled probes were made by nick translation of 732-bp polymerase chain reaction (PCR) amplicons (primers 5# GTGAATTGATTTCAAATGCATC 3# and 5# CATCTTTTTTCTTGTCATCATC 3#) of the previously cloned 927-bp region of A. vaga hsp82 (Mark Welch and Meselson 2000).Fosmids that hybridized with the probe were screened by PCR with the above primers, and the resulting amplicons were purified and sequenced.ClustalW (Higgins et al. 1994;Larkin et al. 2007) alignments and Blast searches (Altschul et al. 1990) were used to determine which copy of hsp82 was present in each fosmid.As a further screen, hsp82 regions were sequenced directly from fosmids.

Fosmid Sequencing and Contig Construction
Fosmid DNA was prepared with an Eppendorf Perfectprep Plasmid Midi kit (Eppendorf, Westbury, NY), and approximately 10 lg of each fosmid was sheared to 3-4-kb fragments (Genemachines Hydroshear, Genomic Solutions, Ann Arbor, MI), subcloned into pCR 4-TOPO vectors (TOPO TA Cloning Kit for Sequencing, Invitrogen, Carlsbad, CA), and sequenced as described (Mark Welch et al. 2008).Approximately 800 bp at both ends of approximately 300 plasmid subclones from each fosmid were sequenced (ca., 9-fold coverage of a 50-kb-fosmid), and remaining gaps were closed by direct sequencing from fosmid templates.Plasmid sequences and any additional gapspanning sequences from each fosmid were assembled into contigs using STRAW (Josephine Bay Paul Center Computational Facility, Marine Biological Laboratories, Woods Hole, MA), a script that combines PHRED, PHRAP, and CROSS_MATCH (Ewing and Green 1998;Ewing et al. 1998) to call bases from sequence chromatograms, align overlapping sequences, and organize them into contigs.There was complete identity within the approximately 8-kb overlap of the two sequenced fosmids containing hsp82-1, and these were combined into one 67-kb-contig designated Avhsp-1.All other contigs presented in this study were single fosmids.

Southern Blot
The 732-bp regions of hsp82-1 and hsp82-3 used to screen the library were also used to probe a Southern blot of A. vaga genomic DNA (gDNA) (Sambrook et al. 1989).Probes were separately labeled and hybridization to a lane containing equal amounts of separate Nsil digests of three fosmids containing hsp82-1, hsp82-3, and hsp82-4 showed them to have equal specific activity (fig.1A, lane 2).The intensities of the bands in each lane were read with a Storm 860 Imager (GE Healthcare, Piscataway, NJ), measured with IMAGEJ (Abramoff et al. 2004), and corrected for background.

Annotation
Putative genes were identified by BlastX searches of the GenBank nonredundant database (Altschul et al. 1990;Benson et al. 2005); by de novo prediction (GeneMark.hmmES-3.0, self-trained on Caenorhabditis elegans) (Lomsadze et al. 2005); and by TblastX searches for conserved regions between contigs within and between species.Candidate genes were then screened against the Pfam (Bateman et al. 2004) and NCBI Conserved Domain (Marchler-Bauer et al. 2005) databases for homology to known protein domains.Intron-exon boundaries were identified by comparison with homologous amino acid sequences using GENEWISE2 (mode-local, organism-worm, no intronic bias, splice site-GT/AG only, model-flat, and algorithm-GeneWise 623) (Birney et al. 2004) and, in a few cases, were identified manually as regions with GT/AG boundaries containing stop codons and joining open reading frames.The results for P. roseola hsp82 contigs differed slightly from the previous annotation in the number of putative genes and the assignment of gene boundaries (Mark Welch et al. 2008).

Divergence Calculations
Contigs were aligned with ClustalW, and total divergence was measured with DNASP 4.0 (Rozas et al. 2003).Genes were aligned based on translation with REVTRANS 1.4 (Wernersson and Pedersen 2003), and synonymous (K s ) and nonsynonymous (K a ) divergence was measured with DIVERGE (GCG Wisconsin Package, Accelrys, San Diego, CA), which corrects for multiple hits and transversiontransition bias (Li et al. 1985;Li 1993;Pamilo and Bianchi 1993).Where the ratio of transversions to transitions was 1, adjustments were made only for multiple hits (Jukes and Cantor 1969) using DNASP 4.0.No corrections were made in the two cases in which raw divergence was 0.75.

Results hsp82 Regions in A. vaga
We constructed a primary genomic fosmid library of A. vaga and probed it with a mixture of approximately 700-bp segments from the same region of each of the three different previously described copies of hsp82 that had been found in this species (Mark Welch and Meselson 2000).Of the approximately 110,000 fosmids screened, 102 exhibited hybridization signals.These were screened by PCR and direct sequencing with hsp82 primers, yielding 67 that contained an hsp82 probe sequence.The hsp82 genes represented in these fosmids were of three kinds, designated hsp82-1 (35 fosmids), hsp82-3 (17 fosmids), and hsp82-4 (15 fosmids), with sequences identical to the three copies found previously and designated copies 1, 2, and 3, respectively (Mark Welch and Meselson 2000).No other hsp82 sequences were found.Two fosmids containing hsp82-1 and one each containing hsp82-3 and hsp82-4 were fully sequenced, giving contigs Avhsp-1 (EU637017), Avhsp-3 (EU637018), and Avhsp-4 (EU831279), respectively.
The finding of approximately twice as many fosmids containing hsp82-1 as those containing hsp82-3 or hsp82-4 suggested that there are two copies of the hsp82-1 gene in the A. vaga genome that are identical within the sequenced region.This possibility was tested by Southern blotting and by FISH.
When probed with a mixture of the 32 P-labeled hsp82-1 and hsp82-3 PCR products used to screen the fosmid library, an NsiI digest of A. vaga gDNA revealed two bands of equal intensity, corresponding to the expected 3.5-kb fragment from hsp82-1 and the expected 2.6-kb fragment from both hsp82-3 and hsp82-4 (fig.1A, lane 1).That the specific activities and hybridization efficiencies of the two probes were the same was verified by the observation that the restriction fragment from hsp82-3 and hsp82-4 was labeled with twice the intensity as that from hsp82-1 when an equal mixture of NsiI digests of the three fosmids, each containing a different copy of hsp82, was simultaneously hybridized on the same membrane with both probes (fig.1A, lane 2).The observation that the 2.6-and 3.5-kb bands are of equal intensity therefore shows that the A. vaga genome contains twice as many copies of hsp82-1 as of hsp82-3 or hsp82-4.
The chromosomal distribution of the various hsp82containing regions was visualized by hybridizing squash preparations of A. vaga embryos with one of the hsp82-1 containing fosmids labeled with a green fluorophore together with a fosmid containing hsp82-3 and a fosmid FIG.1.-Adineta vaga contains four hsp82-containing regions, two of which are identical, on four separate chromosomes.(A) Southern blots of NsiI-digested A. vaga gDNA hybridized to a mixture of hsp82-1 and hsp82-3 probes of equal specific activity.The hsp82-3 probe diverges by 1.6% from the corresponding region of hsp82-4 and hybridizes equally with hsp82-3 and hsp82-4 (data not shown).The 3.5-kb fragment derives from hsp82-1, and the 2.6-kb fragment derives from hsp82-3 and hsp82-4.Lane M, size marker.Lane 1, 7-lg NsiI-digested gDNA; the ratio of corrected intensity of the 2.5-kb band to that of the 3.5 kb-band is 0.92.Lane 2, 5 ng each of digested fosmids containing hsp82-1, hsp82-3, and hsp82-4 mixed with 7 lg of NsiI-digested gDNA; the ratio of background-corrected intensity of the 2.5-kb band to that of the 3.5-kb band is 0.52.(B) Left, A. vaga chromosomes hybridized to a mixture of an hsp82-1 containing fosmid labeled with green fluorophore and hsp82-3 and hsp82-4 containing fosmids labeled with a red fluorophore.Arrows point to hybridization signals.One of the two red signals is photographed in a different focal plane (inset).Right, the same squash stained with 4#,6-diamidino-2-phenylindole.Scale bar, 5 lm.(C) Alignment of the termini of 32 fosmids containing hsp82-1 on Avhsp-1.Each line below Avhsp-1 represents a fosmid that contains hsp82-1.Thickened lines represent the sequenced regions at the termini and in the hsp82 gene of each fosmid insert and more extensive regions sequenced in one of the fosmids, all of which were found to be identical to the corresponding sequence in Avhsp-1.Scale bar, 5 kb.
Degenerate Tetraploidy of Bdelloid Rotifers 377 containing hsp82-4, both labeled with a red fluorophore.As shown in figure 1B, there is a green signal on each of two chromosomes and a red signal on each of two others, consistent with the presence of two copies of Avhsp-1 and one each of Avhsp-3 and Avhsp-4, all on separate chromosomes.
Seeking to determine how far the region of identity extends beyond the hsp82 genes in the two regions containing hsp82-1, we sequenced the terminal approximately 500 bp of the inserts in each of 32 fosmids containing the hsp82-1 sequence from the library screen.These 64 end sequences fall randomly along Avhsp-1 and should be equally likely to derive from one or the other of the two regions (fig.1C).In no case was any difference found between a fosmid insert end sequence and the corresponding sequence in Avhsp-1.It may be concluded that the two regions containing hsp82-1 have sequences that are identical or nearly identical over a region of at least 67 kb containing the gene.
Contigs Avhsp-3 and Avhsp-4 differ by 3% in the approximately 32-kb region of their overlap (gene 23 through gene 43) and contain the same genes in the same order and orientation, with nonweighted average K s 5 0.04 (range 0-0.15) (fig.2, table 1).The only exception involves a region of Avhsp-4 containing four tandem copies of a putative 7 transmembrane receptor (genes 20-23).An approximately 2.5-kb segment of Avhsp-4 containing one of these copies is not present in Avhsp-3.In contrast to the colinearity of Avhsp-3 and Avhsp-4 and the apparent identity of the two regions containing hsp82-1, only 5 of the 19 different putative genes present in Avhsp-3 or Avhsp-4 are also present in Avhsp-1.These five genes are in the same order and orientation, separated by regions up to 12 kb in length that are present in only one or the other of the two colinear pairs.The nonweighted average K s between genes in different colinear pairs is 0.63 (range 0.32-1.16).Thus, the organization of the hsp82 regions of A. vaga, like the hsp82 regions of P. roseola, is that of a degenerate tetraploid.
Histone Gene Cluster Regions in A. vaga and P. roseola Further evidence for bdelloid degenerate tetraploidy is seen in the regions surrounding the clusters of canonical H4, H3, and H2B histone genes in A. vaga and P. roseola.The genomes of both bdelloid species contain four such clusters comprising two colinear pairs with nonweighted average K s 5 0.004 (range 0-0.01) between corresponding histone genes within a colinear pair in A. vaga and K s 5 0.04 (range 0.02-0.06) in P. roseola.The nonweighted average K s between corresponding histone genes in different colinear pairs is 0.38 (range 0.24-0.62) in A. vaga and 0.51 (range 0.43-0.59) in P. roseola (table 2).In order to investigate larger regions containing these histone gene clusters, fosmids covering each cluster of each species were fully sequenced and annotated (Avhis-1, EU652315; Avhis-2, EU850438; Avhis-3, EU652316; Avhis-4, EU850439; Prhis-1, EU850440; Prhis-2, EU652317; Prhis-3, EU652318; and Prhis-4, EU850441), revealing a pattern of degenerate tetraploidy like that of the hsp82 regions of these species (fig.3, table 2).For example, between the most distantly separated genes common to both colinear Degenerate Tetraploidy of Bdelloid Rotifers 379 pairs in P. roseola (genes 11 and 25), there are eight genes present in both pairs, in the same order and orientation, and four genes in only one or the other pair.As in the hsp82 regions of A. vaga and P. roseola, colinear regions containing the histone gene clusters are closely similar, differing by 2% in A. vaga and 6% in P. roseola and contigs belonging to different colinear pairs can be aligned only within the genes they share.The nonweighted average K s for genes in members of a colinear pair is 0.02 (range 0-0.03) in A. vaga and 0.07 (range 0.02-0.18) in P. roseola, whereas the nonweighted average K s for genes in different colinear pairs is 0.68 (range 0.24-1.46) in A. vaga and 1.02 (range 0.43-1.87) in P. roseola (tables 1 and 2).

Introns and Transposable Elements
The hsp82 and histone gene cluster contigs in both species are gene rich with exons constituting approximately 50% of the sequenced 800 kb.The lack of transposable genetic elements and the abundance of small introns in these regions are consistent with findings reported for P. roseola hsp82 regions and other gene-rich regions of bdelloid genomes (Arkhipova andMeselson 2000, 2005;Mark Welch and Meselson 2001;Mark Welch et al. 2008;Gladyshev, Meselson, and Arkhipova 2008).There are introns even in the canonical histone genes contrary to the general lack of introns in histone genes of other species.

Discussion
As has been found in the hsp82 region of P. roseola, we find that this region in the distantly related bdelloid species A. vaga and the histone gene cluster regions of A. vaga and P. roseola each comprise two colinear pairs with only a minority of genes common to both pairs, in the same order and orientation.Thus, our findings are entirely consistent with degenerate tetraploidy in both species.
Within regions in which both colinear pairs can be compared in both A. vaga and P. roseola (hsp82 region genes 19-44 and histone gene cluster region genes 11-23), we find that 30 of 39 genes are present in both species (figs. 2 and 3, tables 1 and 2).Further, each colinear pair in one species resembles a colinear pair in the other species, as seen in the sequence of genes and the pattern of presence and absence of gene-containing segments.Moreover, there is no case in which a gene present in only one colinear pair in a given species is found only in the dissimilar pair in the other species.This interspecies conservation almost certainly means that degenerate tetraploidy was established before the separation of A. vaga and P. roseola.Although it remains to be seen if the same pattern is characteristic of the two other bdelloid families, it appears likely that degenerate tetraploidy was established before the bdelloid radiation and that, consistent with other evidence, Class Bdelloidea is monophyletic (Melone et al. 1998;Mark Welch 2000).
In sexually reproducing species, the colinearity of homologous chromosomes is maintained by segregation and haplotype drift.The presence of colinear chromosome pairs in bdelloids, however, may have an altogether different explanation, being required in G1 nuclei for faithful templatedependent repair of DNA double strand breaks (DSBs) as may be particularly frequent as a consequence of the anhydrobiotic bdelloid life style (Gladyshev and Meselson 2008).As in the hsp82 regions of P. roseola (Mark Welch et al. 2008), there are numerous tracts of sequence identity between members of a colinear pair, possibly resulting from gene conversion associated with DSB repair.The longest such tract entirely within sequenced regions is a 1.3 kb tract between Avhsp-3 and Avhsp-4.Further, the entire 67 kb covered by Avhsp-1, which is present as two identical or nearly identical regions, may be an unusually long conversion tract, as could result from break-induced replication, or, if the identity extends to one or both ends of the colinear chromosomes involved, may be the result of mitotic crossing-over in G2 or chromosome loss and reduplication, respectively.The possibility is not excluded, however, that some of the regions of identity we see are the result of some form of homologous DNA exchange between closely related individuals.
In contrast to the tracts of identity or near identity seen within colinear pairs, the longest stretch of identity between contigs belonging to different colinear pairs is only 99 bp, found within the hsp82-coding region of Avhsp-1 and Avhsp-4.Considering the possibility that conversion tracts are associated with DSB repair, it is interesting that even ignoring the two apparently identical hsp82 regions in A. vaga, the average K s as well as the total divergence within colinear pairs in A. vaga is consistently lower than the corresponding values for P. roseola (tables 1 and 2), possibly reflecting a higher frequency of DSB repair in the more ephemeral aquatic habitats characteristic of A. vaga (Ricci 1998).
The apparent lack of conversion between different colinear pairs also suggests that, all else being equal and assuming no homologous exchange between individuals, K s between corresponding genes in different colinear pairs in one species would be the same as that in the other species.Instead, the nonweighted average K s in the 12 genes present in both colinear pairs in both species (tables 1 and 2) is 0.96 (range 0.48-1.76) in P. roseola and 0.66 (range 0.24-1.37) in A. vaga and is greater in the former species in 11 of these 12 genes (Spearman's rank order correlation; P 5 0.001).The observed difference, if it is a general characteristic of the two species, could reflect conversion or mitotic crossing-over between colinear pairs early in the ancestry of A. vaga or a lower average mutation rate in A. vaga, as might result if it has spent more time than P. roseola in a desiccated state (Ricci 1998).
There are two genes present in both colinear pairs of one species that, although surrounded by genes that are present in both colinear pairs, are present in only one pair of the other.Gene 34 in the hsp82 region is present in all A. vaga hsp82 contigs but absent in Prhsp-1 and Prhsp-2 (fig.2).Gene 19 in the histone gene cluster regions is present in all P. roseola histone gene cluster contigs but missing in Avhis-1 and Avhis-2 (fig.3).The presence of these genes in only one colinear pair in only one species suggests that although most deletions occurred before the two species separated, a few have occurred subsequently.
FIG. 2.-Distribution of genes in hsp82 regions of Adineta vaga and Philodina roseola.Contigs are shown as lines with putative genes shown as boxes numbered as in table 1. Boxes on opposite sides of a line represent genes transcribed in opposite directions.Genes common to all colinear pairs are connected by purple lines.Colinear pairs of one species that are similar to those in the other are designated by letters A and B. Philodina roseola, red; A. vaga, blue.Scale bar, 5 kb.
FIG. 3.-Distribution of genes in histone gene cluster regions of Adineta vaga and Philodina roseola.Contigs are shown as lines with putative genes shown as boxes numbered as in table 2. Other designations as in figure 1. Scale bar, 5 kb.

Table 1
Putative Genes Identified in hsp82 Regions NOTE.-AverageK s values corrected for multiple hits and transition/transversion bias.O-sequenced in only one contig; -present in; *-present in one colinear pair of one species; [ ]-only corrected for multiple hits using Jukes-Cantor corrections.

Table 2
Putative Genes Identified in Histone Gene Cluster Regions NOTE.-Corrections for multiple hits and symbols as described in table 1.