The European sea bass Dicentrarchus labrax genome puzzle: comparative BAC-mapping and low coverage shotgun sequencing

BAC-end sequencing

For comparative mapping, a subset of 10,000 BAC-ES was chosen to perform BLASTN searches with an e-value cut-off of 1e^-5 against the genomes of D. rerio [10], T. nigroviridis [7]], O. latipes [9] and G. aculeatus http://www.ensembl.org. The genomic sequence of T. rubripes [8]] was not used as the genome assembly has not been assigned to chromosomes. The highest number of matches was obtained against G. aculeatus (4,359 ES matches), followed by O. latipes (2,702), T. nigroviridis (2,536) and D. rerio (1,128). The results reflect known phylogenies, with D. rerio (superorder Ostariophysii) distantly related to the other candidates (all from the superorder Acanthopterygii) [39].

Whole genome shotgun sequencing

Sequencing of whole genome shotgun libraries yielded >2 × 10⁶ reads with an average Q20 read length of 673 bp comprising ~1.4 Gbp and approximately twofold coverage of the D. labrax genome. Assembly of WGS-reads and BAC-ES yielded 273,453 contigs and 217,926 singlets covering ~580 Mbp. The N50 contig size was 2,891 bp and the largest contig was 15,629 bp. A part of this dataset, namely 36,166 contigs were useful to anchor additional BAC-ES to the stickleback genome (see below). These contigs have been submitted to EMBL nucleotide database [EMBL:CABK01000001 - EMBL:CABK01036166].

Comparative mapping

The whole BAC-ES dataset was aligned with the fully assembled stickleback genome. Further sorting and screening yielded 25,845 BAC-ES where only one end was sequenced or matched the stickleback genome and 13,996 BACs matched both ends to the same chromosome in stickleback. 18,013 BACs were matching weakly and were excluded (mostly due to repetitive motifs or possible chimeric BACs). BACs with both ES aligned were essential for comparative mapping and could be subdivided into 12,076 BACs with correct orientation and distance of aligned ES and 1,920 BACs not matching these consistency criteria due to possible rearrangements, miss-alignments or assembly failures in the stickleback genome. Plotting the frequency distribution of insert size of consistently mapped BACs resulted in a Gaussian-like distribution with a maximum at 115 kbp. This reflects a compression of the stickleback genome compared to the D. labrax genome, as the average insert size published for D. labrax is about 164 kbp [26] (see Fig. 1).

Comparative mapping was improved by aligning BAC-ES containing contigs from the WGS and BAC-ES assembly to stickleback chromosomes. This strategy yielded 20,635 BACs matching consistency criteria, an improvement of about 71% compared to comparative mapping using only BAC-ES data. The re-calculated minimal tiling path reduced total contig number to 588 and N50 contig number to 83 while increasing N50 BAC-contig size to 1.2 Mbp and coverage of stickleback chromosomes to 87% (see Table 2, Table 3, Table 4 and Fig. 2, a complete list of ordered paired-end aligned BAC clones on stickleback chromosomes may be downloaded [Additional file 1]).

Moreover the higher coverage with BACs enabled the identification of potential intra-chromosomal rearrangements between sea bass and stickleback (or failures in the stickleback assembly). A number of 214 potential chromosomal breakpoints spanned by BAC-clones were identified [Additional file 2]. To check if rearrangements were artefacts, the order of calculated BAC-contigs was cross-checked by alignment to the medaka genome (see Fig. 3). In total 139 cases (65%), had a neighbouring position at that site and thus confirmed the consistency of the identified BAC-clone on the second reference genome.

Fig. 4 shows PCR results that support the bioinformatic data on rearrangements between sea bass and stickleback. All of the seven rearrangements that were checked by PCR have been confirmed. For each of these rearrangements we found at least 2 BAC clones that gave positive results in the PCR, the average number of BAC clones spanning a rearrangement was 4.6 and the maximum number was 7 clones.

The consistently mapped BACs were also uploaded to the Ensembl genome browser and may be viewed in a user friendly format alongside the annotated stickleback chromosomes (see Fig. 5)

Comparative mapping compared to other mapping strategies

Since publication of the first BAC-vector [40] several strategies for the construction of physical genome maps from BAC-libraries have been published. Among these methods BAC-filter hybridization [41], BAC-fingerprinting [42] and PCR screening [43] have been applied most frequently. Comparative mapping approaches are likely to replace these methods because many genomes of higher eukaryotes have been published. Comparative maps are built by aligning paired end sequences of large insert clones (e.g. BACs) to a reference genome and thus detecting possible overlaps of clones that subsequently can be combined into contigs. This strategy has been successfully applied to closely related organisms such as chimpanzee and human [44] and also to more distantly related organisms like cattle and human [45]. Comparative mapping has some advantages for automated analysis over the methods mentioned above, as established pipelines for high-throughput sequencing and bioinformatics can be used.

BAC end sequencing results

Sanger sequencing of BAC-ends remains restricted to a 96 well format in many sequencing centers, because of low template yields and large amounts of template used for the sequencing reactions. Thus the successful development of an automated DNA purification process to purify BAC-DNA from 384 well plates was a crucial step to enable the comparative genome mapping of D. labrax. With an average read length of 650 bp on 36 cm and 750 bp on 50 cm capillaries the read length of BAC-end sequences was substantially higher than reported in comparable projects [45]. Failed reactions were less than 11%.

Reference genomes

Besides read quality, the choice of a suitable reference genome is influencing mapping success and quality. Several sequenced genomes of model teleosts are available (e.g. D. rerio, T. nigroviridis, T. rubripes, O. latipes and G. aculeatus). With the highest number of mappable reads, the stickleback genome sequences shared the highest homology to D. labrax, making it the genome of choice for a comparative approach. The stickleback and the European sea bass belong to the superorder Percomorpha, and the evolutionary related orders of Gasterosteiformes and Perciformes, respectively [46]. Additional beneficial features of the stickleback genome sequence is the high sequencing coverage (~12 fold) and the mapping of most scaffolds to chromosomes.

The comparative map

After BAC-ES data for sea bass became available, a first comparative BAC-map was built. Results were already usable to render megabase sized contigs and to screen for BACs covering genes of interest. Subsequently with a WGS dataset of the sea bass genome available, BAC-ES sequences and WGS data were combined by assembly. If a BAC-ES alone could not be matched with the reference earlier on, the length extension of the BAC-ES by aligned WGS reads now increased the probability to find matches with good alignments to the reference genome. Final mapping (Fig. 2) shows in green that most of stickleback chromosomes are covered by D. labrax BACs with consistent orientation and distance (87% of reference genome), while red regions have a weak mapping, where no or only one BAC-end sequence could be matched. These regions may be either due to highly repetitive fragments, gaps and/or failures in the assembly of the sticklebackgenomic sequence or regions that are underrepresented in the BAC-library. It is obvious that especially centromeric and telomeric regions, known for the problems mentioned above, account for weakly mapped regions.

Calculating the genome size of D. labrax

When comparing the insert size distribution of consistently mapped BACs on the reference genome with the published insert size distribution of the D. labrax BAC-library (Fig. 1), a shift to lower insert sizes is observed. An explanation for this may be found in the evolution of genome size. It has been shown that teleost genomes tend to accumulate most indels in intergenic or intronic regions leading towards large differences in genome size, while synteny of genes is conserved [47]. Thus one may conclude that the ratio of the maxima in the insert size distributions of BAC-clones equals the ratio of genome sizes. From this calculation one may conclude that the D. labrax genome is about 1.3 fold larger than the 460 Mbp of the G. aculeatus genome. The calculated haploid genome size of 600 Mbp is smaller than the estimated haploid genome size of 763 Mbp derived from flow cytometric measurements of diploid nuclear DNA content [31]. A smaller genome size is also suggested by the first assembly of our twofold coverage WGS dataset (see WGS sequencing results). Nevertheless, genome size estimates from sequencing may be biased towards the euchromatic portion of the genomes and different results of the methods may be explained by underrepresentation or different size evolution of heterochromatic regions.

Comparing the BAC and the linkage map of D. labrax

BAC contigs represented by green regions in Fig. 2 are considered blocks with a high level of synteny between D. labrax and G. aculeatus. Nevertheless it is questionable whether neighbouring BAC-contigs on the reference genome are really neighbours in the D. labrax genome or whether chromosomes have undergone extensive inter-chromosomal rearrangements during evolution. To decide either whether a rearrangement has taken place or the order of BAC-contigs is consistent in both genomes, it is helpful to compare results from the D. labrax genetic linkage map and the radiation hybrid map of the closely related sparid Sparus aurata (gilthead sea bream). Such comparisons with stickleback have been done by Chistiakov et al. [30] and Sarropoulou et al. [16] and showed synteny of complete chromosomes between these species. Chromosome identity and re-shuffling are common features among closely related organisms. The different chromosome number of G. aculeatus (n = 21) and D. labrax (n = 24) is a common feature between related taxa and can be explained by fusions/fissions of complete orthologous groups. Thus it is unlikely that BAC-contigs mapped to one G. aculeatus chromosome are not located on a single D. labrax chromosome. These results also allow assigning the comparatively mapped BAC-contigs to D. labrax linkage groups (Table 2A, Table 3A and Table 4A).

Comparison of the D. labrax linkage map with the G. aculeatus genome has suggested some intra-chromosomal rearrangements [30]. Due to the higher resolution of the comparative BAC-map, it is possible to pinpoint potential rearrangements by focussing on inconsistently mapped BACs that connect two BAC-contigs at their boundary regions. Since BAC-libraries are known to harbour some chimeric clones, the location of potentially neighbouring BAC-contigs was confirmed by cross-checking their position in the medaka genome. If BAC-contigs connected by a rearrangement spanning BAC were located next to each other in the medaka genome, a true rearrangement was considered (Fig. 3). In this way 139 BACs spanning rearrangements between the reference and D. labrax genome were identified. Seven rearrangements between chr III of stickleback and the corresponding sea bass linkage group 10 were also tested by means of PCR. All of them could be confirmed (Fig. 4).

Applications

The main advantage of BAC-maps over other mapping methods, like genetic linkage maps or radiation hybrid maps, is the possibility to access defined portions of a genome for subsequent analysis by common methods of molecular genetics. As the comparative BAC-map covers about 85.4% of predicted G. aculeatus genes, it is now possible to easily access orthologous D. labrax genes by selecting a BAC-clone that covers the genomic region of interest. As proof of principle, we have successfully identified and shotgun sequenced 10 overlapping BAC-clones that cover a 1.3 Mbp genomic region on sea bass linkage group 5 (Negrisolo et al. in preparation). The BAC map was also used to analyze two clones that contain a novel immune-type receptor (NITR) gene cluster [48] and to sequence the fatty acid delta-6 desaturase gene in European sea bass (Santigosa et al. in preparation).

BAC end-sequencing

The Dicentrarchus labrax BAC-library constructed by Whitaker et al. [26] was obtained from the German resources center for genome research (RZPD, Berlin, Germany). The library comprises pCC1BAC-clones arrayed in 180 × 384 well microtiter plates. The total genome coverage of the library is >12 fold with an average insert size of 164 kbp per BAC-clone. For end sequencing, BAC-clones were inoculated in 2 × 384 deep well plates containing 190 μl of 2YT media and 12.5 mg/l chloramphenicol and cultivated for 18 h at 37°C with rigorous shaking at 1100 rpm in Titramax 1000 incubators (Heidolph Instruments). BAC-DNA was purified by an automated process that was developed at the MPI for molecular genetics. The process applies size selective precipitation in polyethylene-glycol 6000/2-propanol mixtures and a final washing step with ethanol 70% (v/v).

BAC-templates were end sequenced using ABI BigDyeV3.1 Terminator chemistry and T7 or SP6 primers. After post-sequencing cleanup by ethanol/NaAcetate precipitation, sequence analysis was performed on ABI3730 × l capillary sequencers with either 36 cm or 50 cm capillary arrays. Processing of raw sequencing data was done by the PHRED basecaller [51], quality clipping and vector-clipping by LUCY [52].

Whole genome shotgun sequencing

For the construction of WGS plasmid libraries of Dicentrarchus labrax, we obtained genomic DNA isolated from the same specimen (male 57 originating from the Adriatic clade) that was used for BAC-library construction (kindly provided by A. Libertini, CNR, Venice, Italy through J. B. Taggart, University of Stirling, UK).

Genomic DNA was sheared by ultrasonic sound and size selected for fragment sizes of 0.9 - 1.5 kbp and 1.5 - 4 kbp. Fragments were polished by T4-DNA-polymerase/DNA-polymerase I (Klenow) and ligated with T4-DNA-Ligase into SmaI digested pUC19 sequencing vector. Competent E. coli DH10B cells (Invitrogen) were transformed by electroporation, plated on 22 × 22 cm agarplates (Nunc) containing LB media with 110 mg/l Ampicillin, X-GAL and IPTG. After 16 h of incubation at 37°C white colonies were arrayed into 384 well microtiter library plates by a picking robot (Q Bot, Genetix). These plates (media: LB+HMFM+Ampicillin) were again incubated for 16 h at 37°C and stored at -80°C. Plasmid DNA preparation was done as described for BAC-DNA with the difference that the final washing step with 70% (v/v) ethanol was not necessary and a single 384 deepwell microtiter plate filled with 190 μl of 2YT + 110 mg/l Ampicillin yielded enough template amounts for several sequencing reactions.

Sequencing, sequence analysis and sequence processing of plasmids was done as described above using ABI BigDyeV3.1 Terminator chemistry and M13(-40) or M13(-28) primers.

Alignment of BAC-ES to reference genome

BAC-ES were aligned by BLAST [53, 54] algorithms to genomic sequence of G. aculeatus (Assembly: BROAD S1, Feb 2006, http://www.ensembl.org). To minimize computational time BLAST searches were done incrementally beginning with stringent parameters (Megablast, word size 20, and nucleotide mismatch penalty -1). Results were filtered for alignments that matched with an e-value equal or lower than 10^-5. Additionally, alignments were only submitted to further analysis, if the second best alignment resulted in an e-value that was at least 10⁵-fold larger than the e-value of the best alignment. Sequences with alignments not matching these criteria were extracted by notseq [55] and subsequently aligned by BLAST searches with lower stringency. Stringency in the following rounds was adjusted by choosing word sizes of 15, 11 and 7.

The number of BAC-ES with alignments meeting our criteria was further improved by adding sequences from whole genome shotgun sequencing of D. labrax. All sequences available were assembled by the Celera Assembler [56]. Contigs that contained BAC-ES were filtered and again aligned to the G. aculeatus genome as described above. Match coordinates of contigs on G. aculeatus chromosomes were corrected by position of the BAC-ES in the contigs and assigned to the corresponding BAC-ES.

Calculating and visualizing the comparative map

Resulting BLAST-tables of both approaches were screened for BACs that were aligned with both ends. These BACs were further screened for matches to the same chromosome in G. aculeatus and then checked for consistent orientation and distance.

BACs that matched all consistency criteria were chosen for the calculation of a minimal tiling path. Starting with a first BAC-clone, BAC-contigs were constructed by choosing BACs that were overlapping and maximizing the contig in length. These analyses were done by common spreadsheet software and scripting language. BAC-contigs arranged on the 21 G. aculeatus chromosomes were visualized by passing coordinates to a vector graphics application (CorelDRAW version 11, Corel corp., Ottawa, Canada). To view the BAC map alongside the annotated stickleback genome the mapping coordinates were uploaded to the ensembl genome browser as a GFF formatted textfile.

Dealing with possible rearrangements and checking them on a second reference genome

The subset of BAC-ES that aligned to the same chromosome but did not match consistency criteria could be due to intra-chromosomal rearrangements between stickleback and sea bass. These clones were visualized as black arcs on the stickleback chromosomes. If these arcs were starting not at the edge of contigs they were manually removed. To check if the rearrangement spanning BAC-ES have a consistent order in the medaka genome, we exploited stickleback and medaka synteny of orthologous genes. Using the biomart tool a table was prepared that showed the genes annotated to stickleback with their orthologous position in medaka. Furthermore the coordinates of contig starts and ends from the sea bass BAC-map were implemented in the table. In this way the position of sea bass contig starts and ends on medaka could be mapped to stickleback chromosome coordinates. BAC-contig edges that are located next to each other in medaka were subsequently visualized by arcs, in many cases confirming a connection between contigs that was also found before by non-consistent matching sea bass BAC-ES.

Evaluation of several rearrangements by means of PCR

Seven potential rearrangements between stickleback chr III and the corresponding sea bass linkage group were checked by means of PCR. Primers for PCR were designed on BAC-ES representing the end of BAC-contigs that seem to be connected by a rearrangement spanning clone. Subsequently amplification of the chosen markers was carried out using the rearrangements spanning clones as templates. If both BAC-contig end markers can be amplified on a rearrangement spanning BAC, the overlap and therefore connection of the two BAC-contigs in sea bass is confirmed. Additionally markers were amplified on genomic DNA of sea bass to check that they are unique markers in the genome. The PCR was set up as 50 μl reactions. For amplification of BAC-templates we added 2 μl of overnight culture to the PCR, while amplification of genomic DNA was carried out by adding 2 μl DNA with a concentration of 45 ng/μl to the PCR. Composition of PCR was as follows: 0.3 μM for each primer, 300 μM dNTPs, 75 mM TRIS-HCl, pH 9, 20 mM (NH₄)₂SO₄, 0.01% Tween 20, 2.5 mM MgCl₂, 0.1 U/μl Taq-DNA-polymerase and 0.5 M Betaine. Thermocycler profile was: Step I: 5 min at 94°C. Step II: 30 s at 94°C. Step III: 30 s at 55°C. Step IV: 1 min at 72°C. Step V: 7 min at 72°C. Step VI: hold at 4°C. Steps II-IV were repeated 25 times. PCR products were analyzed on 1.5% agarose gels and stained with ethidium bromide.