A duck RH panel and its potential for assisting NGS genome assembly

The development and commercialization of next-generation massively parallel DNA sequencing approaches, by dramatically decreasing the cost of sequencing, have revolutionized genomic research. The main innovation of NGS, as compared to Sanger sequencing, is the parallelisation of the process, allowing between a few thousands and up to millions of sequencing reactions to be processed simultaneously. The three main NGS technologies available on the market use different approaches for library construction, template immobilisation and sequencing reaction, but the basic principles remain the same. NGS approaches also have some drawbacks compared with Sanger sequencing: (1) sequence reads produced currently by NGS (100 bp for Illumina, 500 bp for 454) are shorter than Sanger sequencing reads (1000 bp) and have a higher error rate, making the sequence assembly more problematic; (2) the pairing of reads in Sanger sequencing is limited by the size of the DNA fragments that can be inserted in cloning vectors, ranging from 1–2 kb or less up to 100–200 kb (plasmids, fosmids, BAC), whereas pairing of reads is limited to 40 kb with NGS, limiting the average assembled scaffold length and leading to more difficulty in segmental duplication and copy number variation detection; (3) as a consequence, a higher sequencing depth is required for assembly and a very high number of small scaffolds are produced. Nevertheless, the sequencing and de novo assembly of a Chinese individual [1] proved the feasibility of sequencing and assembling whole genomes by NGS. Many species have been sequenced and/or resequenced by NGS, such as the giant panda Ailuropoda melanoleura[2], the silk worm Bombyx mori[3], the cucumber Cucumis sativus[4], the chicken Gallus gallus domesticus[5, 6], the turkey Meleagris gallopavo[7] and the Mallard duck Anas platyrhybchos domesticus (Huang et al., in prep).

The Pekin duck (Anas platyrhynchos, APL) is an obvious target for detailed genomic studies due to its agricultural importance [8–10] as well as for its role as a natural reservoir of all influenza A viruses. It can usually carry the infection with no sign of disease and thus propagates the virus to other bird species and potentially to mammals such as pigs or humans [11–15]. The duck genome presents most of the characteristics encountered in birds, which are: (i) a more compact genome, one third the size of a typical mammalian one, (ii) a large number of chromosomes (2n = 80), (iii) the presence of macrochromosomes and microchromosomes, the latter being as small as a few Mb [16] and (iv) the females are the heterogametic sex (ZW) and males the homogametic one (ZZ). Due to its importance both in the economic and scientific fields, the sequencing of the Pekin duck genome was initiated in 2008 using the same strategy recently published for the giant panda [2] at the Beijing Genomics Institute (BGI). The sequence reads provided a depth of 65X and a total of 78,487 scaffolds were assembled in which the N50 scaffold was 1.2 Mb and the largest was 5.9 Mb in length (Huang et al., in prep). However, owing to the lack of a clone-based physical map and other supplementary mapping data, apart from a first generation genetic map composed of 155 microsatellite markers, 115 of which located in only 19 linkage groups spanning 1353.3 cM [17], it is possible to assign only very few assembled scaffolds to chromosomes.

Several studies have shown that birds seem to have a slower rate of chromosome rearrangements than mammals, with only very little inter-chromosomal rearrangements [18–23]. Between chicken and zebra finch, whole genome comparisons revealed 114 tentative intrachromosomal rearrangements (56 inversions and 58 translocations) in which some were confirmed by FISH (Fluorescent In Situ Hybridization) experiments [24, 25]. Recently, Zhang et al. (2011) provided confirmed evidence for 20–27 major rearrangements between turkey and chicken, almost all of which are inversions [26]. The mean reported phylogenetic distance between chicken and turkey is 47 million years, whereas it is 81 between chicken and duck [27], so the number or rearrangements reported between chicken and turkey provide the minimum level of difference expected between chicken and duck. To date, only one interchromosomal difference has been reported between the chicken and the duck karyotypes, with APL4 and APL10 corresponding to GGA4q and GGA4p respectively [21]. This interchromosomal rearrangement explains the difference in diploid chromosomal number between the two species, which is 2n = 78 in chicken and 2n = 80 in duck and therefore the nomenclature for numbering the duck chromosomes follows mainly that of chicken. Macrochromosomes APL1 to APL9 correspond to GGA1 to GGA9, then APL10 corresponds to GGA4p and finally, the rest of the karyotype is offset by one, with GGA10 corresponding to APL11 and so on. Cross-species fluorescent in-situ hybridization (FISH) studies using chicken BAC clones on duck metaphase spreads showed only a few large scale intrachromosomal rearrangements concerning the largest chromosomes [23, 28]. All this demonstrates a high karyotype stability despite 80 million years of divergence between the two species [29, 30].

As a first attempt to order the duck sequence scaffolds, we aligned the 7,205 ones larger than 1 kb to the current chicken assembly using the Narcisse database [31, 32] and successfully positioned 1,787 of them. This still leaves a large number of scaffolds to assign and also means that the ordering of the duck scaffolds and genes we obtained will follow the chicken genome and will be wrong whenever large- or small-scale rearrangements will have happened between the two species.

To assemble the scaffolds in an order corresponding to the real duck chromosomes, several approaches can be used. High density SNP genetic maps allow high precision in mapping. However, the SNP markers need first to be discovered by a sequencing approach, such as published by Kraus et al. (2011) [33] and must be informative in a reference population to be used for mapping. However, despite several thousand SNP discovered to date [33], only a small subset of 384 were genotyped [34], mainly due to the high cost that would have been required for additional markers. Finally, out of these, only a small subset was informative in mapping populations (R. Kraus, R. Crooijmans, personal communication), which will allow only for low resolution maps and poor marker ordering. Furthermore, for high precision mapping, very large populations counting several hundred individuals are required, yet again increasing cost and labour. Further sequencing for SNP detection and a consortium for generating a SNP chip would help improve genetic maps in duck and may happen in the future. Physical maps can be based on the mapping of BAC clones by FISH for chromosome assignment and large-scale ordering. The BAC clones from large libraries can be used for contig construction by fingerprinting or high throughput hybridization. End-sequencing of the clones allow linking sequence scaffolds together. BAC contig maps are thus usually a backbone to the sequence assembly. A BAC library has been made for Duck [35], but to the best of our knowledge, there are no plans yet to build physical maps. In this context, radiation hybrid mapping can be an excellent complementary mapping approach, as it does not require complex marker development and large-scale genotyping. Any STS can be placed on the map by simple PCR on as little as 90 hybrids. Thus with a minimal effort, maps with a resolution intermediate between the genetic and the BAC contig maps can be constructed to propose a correct chromosomal assignment and ordering of scaffolds. We report here the production of a duck whole genome radiation hybrid panel and demonstrate its utility to verify the quality of sequence scaffolds and for assigning and positioning scaffolds onto chromosomes. Large-scale culture of radiation hybrid clones is a time-consuming process and moreover causes the loss of donor fragments in the hybrids. To avoid the necessity of cultivating the radiation hybrid clones at a large scale, we tested three approaches. One involves whole genome amplification (WGA) and conventional genotyping by PCR and gel electrophoresis and the other two use minute amounts of DNA and Fluidigm BioMark^TM IFC Dynamic Array^TM genotyping by quantitative PCR. The advantage of the Fluidigm approaches is low cost combined with simple and rapid high-scale genotyping.

Overall, the pattern of retention for the broken duck chromosome fragments in the hamster cells obtained here is very similar to that observed for the chicken radiation hybrid panel, with higher retentions for microchromosomes than for macrochromosomes. However, whereas only 23% of the chicken-hamster hybrids produced had sufficient retention frequency values to be retained in the final panel, 50% of the duck-hamster clones did. Indeed, although the fusion efficiency for chicken-hamster hybrids was reported to be as high as 2–9 x 10^-6 by Kwok et al.[44], it was only approximately 1.4 x 10^-6 in our hands when we produced the 452 clones for the chicken whole genome RH panel. Here, the fusion efficiency is close to 3.5 x 10^-6 which is three times higher. Such differences could be due to variations in chromosome structure and/or content between the two bird species or to differences in culture conditions. For instance, the HPRT gene used as a selection marker for the clones is on the short arm of macrochromosome GGA4 in chicken [45] and thus very likely to be on microchromosome APL10 in duck. Microchromosomes being better retained than the macrochromosomes, having the selection gene on one of them could help recovering a higher number of clones in each fusion experiment. It is also very likely that these results are due to our change in culture conditions after the cell fusions: the chicken-hamster hybrids were cultivated at 40°C, the usual temperature for avian cells, whereas the duck-hamster ones at 37°C. Similarly, Ekker et al.[46] succeeded in producing zebrafish somatic hybrids at 37°C but not at 28°C, which is the normal temperature for the culture of zebrafish cells. More generally, the difference in optimal temperatures for the growth of donor and recipient cells may be one of the possible causes for the lower retention frequencies usually observed for somatic and radiation hybrids in non mammalian species.

To obtain the DNA quantity required for building genome-wide maps, large-scale culture of the hybrid clones is necessary. However, in this process, donor DNA is lost. For instance, Karere et al.[47] reported a genome wide half-life of the donor DNA of 8.7 passages and when preparing the whole genome RH (WGRH) panel in chicken, we observed the loss of 10% of the chicken genome after large cell culture of the hybrids [38]. This problem, in addition to the fact that large-scale culture of a RH panel requires lots of labor, encouraged us to find an alternative, such as WGA or scaling down the reaction volumes. Since the 1990s three major whole genome amplification techniques including primer extension pre-amplification (PEP) [48], degenerate oligonucleotide primed (DOP) PCR [49] and multiple displacement amplification (MDA) have been developed to address the problem of limiting amounts of DNA samples. PEP and DOP are both PCR-based methods and are limited by features of the Taq polymerase: typical amplification fragment length of < 3 kb and an error rate of 3 × 10^-5. These methods also suffer from incomplete coverage and uneven amplification of the genomic loci of several orders of magnitude, with 10^-2 ~ 10^-4 and 10^-3 ~ 10^-6 fold amplification biases for PEP and DOP-PCR methods, respectively [50]. MDA is an isothermal amplification employing the high fidelity Phi29 phage DNA polymerase for DNA synthesis and strand displacement [51]. The genome coverage is much improved, with an estimation of only 2.2% missing after WGA by the MDA method in mammals [52]. Karere et al.[47] confirmed that MDA was suitable for RH mapping and reported a high concordance rate of 97.6% with data from genomic DNA. However, even if this is true for mammals, it might not be the case of microchromosomes in an avian genome.

When comparing retention frequencies before and after WGA in the 90 hybrids, with the 35 markers used for clone selection, only slight variations of retention, either gains or losses, were usually observed. However three markers, CAUD064, S618 and CAUD022, show an important loss of retention frequency after WGA while two others, CAUD013 and S2870, show a high increase, suggesting potential coverage problems by the WGA, either by lack of coverage (losses) or by the over-representation of a region (gains). Moreover, genotyping of eight no hit EST markers on WGA DNA, either using conventional PCR and Agarose or FLDMqPCR, demonstrated a very low retention which is not in accordance with the retention levels usually observed for microchromosomes. Therefore, we suggest that the genomic features in the smallest microchromosomes causing coverage problems in whole genome sequencing projects may also interfere with the efficiency of WGA. As we have already shown, RH mapping can allow building maps for non-sequenced chromosomes [42, 43], it is important that we produce genotyping results for them.

In this context, the Fluidigm BioMark^TM IFC Dynamic Array^TM genotyping method can be an alternative to WGA, as only minute amounts of DNA (as little as 70 pg) are required. High throughput gene expression analysis by real time PCR in a microfluidic dynamic array was first introduced by Spurgeon et al.[40], and has since been successfully applied to copy number variation studies [53] and quantitative miRNA expression analysis [54]. In our case, by performing qPCR with the Fluidigm BioMark^TM IFC Dynamic Array^TM genotyping, the additional benefit is high throughput, as the identification of bands on gel electrophoresis is replaced by monitoring the PCR with Ct (Cycle threshold) and end point Tm (melting temperature) values, allowing the distinction between specific and non-specific amplification profiles. The Tm value is mainly influenced by base composition of amplicons, making it a specifically interesting parameter to follow when using markers defined from coding regions, which are more prone to cross-amplifying the hamster DNA.

We tested the Fluidigm genotyping method on WGA DNA and on standard DNA, with a pre-amplification step using a mix of all primers of the 96 markers analyzed together in a run [55]. In the WGA-FLDMqPCR runs, Ct values for the duck positive control was high with an average of 22 cycles (data not shown), as opposed to an average of 12 cycles, which is in the recommended scale, for the Pre-ampFLDMqPCR runs (Figure 3). These high Ct values suggested the quantity of DNA template was too low [55]. For a variety of reasons, WGA coupled with either FLDMqPCR or conventional PCR and agarose electrophoresis was unsuitable for genotyping on the smallest microchromosomes. Therefore, although the combination of WGA and FLDMqPCR would have allowed us to use less RH DNA, we decided the best genotyping method was to use standard DNA by FLDMqPCR genotyping, with a pre-amplification step performed using a mix of all primers for a set of 96 markers.

The drawbacks of genotyping by Fluidigm BioMark^TM IFC Dynamic Array^TM come from the fact that all 96 markers are genotyped with the same condition and therefore special care must be taken in the marker design. As a consequence, approximately 10% of the markers were discarded during the final analysis due to poor quality data.

Apart from improving the genome assembly by assigning and ordering scaffolds to chromosomes, the duck RH panel can be used to test the scaffold assemblies. We tested this by genotyping markers at Mb density on the 13 scaffolds larger than 4 Mb, spanning altogether 60 Mb and thus accounting for 5.5% of the current duck genome assembly. A total of 70 markers were genotyped and only one marker (sca504F) on the end of sca504 was not linked with other markers derived from the same scaffold (Additional File 2: Figure S1), suggesting an overall good quality of the final genome. To test further our capacity for detecting potentially misassembled scaffolds, we took advantage of previously published data indicating that on the whole, avian chromosomes are known to be well conserved throughout evolution and more specifically, that no inter-chromosomal rearrangements, apart from the well documented case of GGA4 = APL4 + APL10, have been discovered between chicken and duck by current comparative cytogenetic approaches [17, 23, 28]. The 41 scaffolds (including sca504) we detected as potential chimeras by comparative mapping had poor pair-end sequence support (BGI, personal communication), suggesting most of them could indeed be misassembled (Additional file 4: Table S2). We tested nineteen of them by genotyping markers flanking the potential breakpoints (Additional File 3: Figure S2). As a result, all but one scaffold (sca649) could be misassembled, and sca649 possibly suggesting the first detection of a small inter-chromosomal rearrangement between the duck and chicken genomes, or perhaps more likely a segmental duplication in duck or in the last common ancestor of the two species. This would need further confirmation by FISH mapping with chicken BAC clones. It can be noted that the pair-end sequence support for this scaffold was high, showing an agreement between sequencing and RH mapping data. When disagreements between assembly and our RH data are detected in large scaffolds, they tend to happen towards the end. To achieve better assembly accuracy, higher sequencing depth or more efforts on developing sequencing libraries with longer inserts are needed.

Concerning the smallest duck microchromosomes, paralogous to those absent from the chicken assembly, we suspect similar problems will arise: lack of sequence information, difficulties in cloning, in genetic mapping, etc. RH mapping has proved useful for getting a grip on these regions and one striking example is the case of some regions of HSA19, to which no corresponding chicken genome data could be assigned by sequence similarity and to which many chicken no hit EST showed significant sequence similarity. RH mapping with these markers allowed building maps for GGA30 and GGA32 [43]. By developing markers targeted to this region, a small linkage group composed of 4 no hit markers (absent in the chicken genome assembly) orthologous to HSA19 was obtained. When aligned to HSA19, we found they spanned a 5 Mb region on HSA19p. Due to the lack of BAC clones for FISH or other supplementary information, we cannot identify the duck chromosome, but according to known data on synteny conservation between chicken and duck, we suggest that this small linkage group should be assigned to APL31. Of the 8 no hit to chicken markers we studied three have hits with small to medium-size scaffolds (between 23 and 96 kb) of the duck assembly, suggesting that more sequence from the smallest microchromosomes could be obtained in NGS (Table 2). Chicken/duck comparative mapping of GGA21/APL22 and GGA11/APL12 microchromosomes demonstrate several intrachromosomal rearrangements, the first described for microchromosomes in this pair of species. The maps obtained using the usual Carthagene mapping approach or the comparative approach are very similar, apart for a few markers, especially non-framework / non-robust ones, for which lower reliability in map position can be due to the limits of the possible resolution of the mapping or to genotyping errors. As the comparative approach starts with an ordering of markers corresponding to chicken, it is interesting to note that the major duck-chicken rearrangements found with the Carthagene approach are confirmed. A second advantage of the comparative mapping approach and the associated construction of robust maps is that the number of robust markers obtained is usually higher than the number of framework markers in the classical approach. The major inversion found between GGA11 and APL12 is confirmed by FISH mapping, but also by sequence alignment of duck scaffolds on the chicken assembly. Indeed, scaffold736, whose integrity is demonstrated by RH mapping, with markers sca736A and sca736B positioned close to one another at 153 and 154 cR on the CarthageneRH map and 402 and 441 cR on the ComparativeRH-Robust map, is separated in two locations when aligned to the chicken sequence (Figure 5). Likewise, although a more complex pattern of events accounts for the differences between GGA21 and APL22, one of them is supported by scaffold246, whose integrity is demonstrated by RH mapping with three markers on the robust map, each of which are positioned in different regions when aligned to the chicken sequence. Another is supported by markers sca871_1 and sca871_2, which are co-localized on the RH map and are 1.4 Mb apart in chicken (Figures 3 and 4). When comparing the turkey and chicken genomes, Zhang et al. also confirmed evidence for 20–27 major rearrangements between the two species and found one inversion between GGA11 and MGA13. However, they did not observe any rearrangement between GGA21 and MGA23. The mean estimated divergence time between chicken and turkey is 47 million years and 81 between chicken and duck [56]. A higher number of rearrangements are thus expected between the two latter pair of species. Only one major interchromosomal difference -with APL4 and APL10 corresponding to GGA4q and GGA4p respectively [21]- and very few intrachromosomal rearrangements have been reported between the chicken and the duck karyotypes [18–23]. The rearrangements observed with our data between GGA21 and APL22 seem more complex for example, Sca246B, Sca246C and Sca246D are split by Sca1885 in both RH maps. Likewise, Sca367B and Sca367C are split by Sca3327 in both RH maps, whereas they are adjacent in the chicken sequence, and Sca148 markers are widely split in the ComparativeRH map, while adjacent in the chicken. Further investigations and more precise maps using different techniques such as FISH or BAC contig maps will be needed to confirm these rearrangements. The increased resolution obtained by RH mapping as compared to the FISH mapping performed to date show that intrachromosomal rearrangements might happen on a finer scale than shown until now. This means that although the simple ordering of the duck scaffolds along the chicken genome by sequence similarity helps for chromosome assignment, the duck sequence thus obtained will be wrong whenever large or small-scale rearrangements will have happened between the two species. The whole duck assembled sequence will have to be ordered using the whole genome RH map which will be constructed in our laboratory, in conjunction with other mapping methods, such as genetic and/or BAC contig physical maps, the latter allowing finer mapping and orientation of small scaffolds.

The chicken WGRH panel has been used to construct chromosome RH maps and helped in the genome assembly or the mapping of some of the smallest micro chromosomes [42, 43, 57]. Similarly, the duck WGRH panel presented here will also be a major contribution to duck genomics. RH mapping can be a complementary approach to NGS by allowing the assignment of scaffolds to duck chromosomes and furthermore, detailed RH maps will allow a precise estimation of the intrachromosomal rearrangements that have occurred between chicken and duck.

Using the chicken genome as model and in combination with survey sequencing, the construction of dense RH maps of a less studied bird such as duck can be made. By taking advantage of the conservation of syntenies, optimal orders can be proposed [58, 59], thus maximizing the information obtained as first proposed by Hitte et al.[60, 61]. Indeed in duck, a dense RH map combined with scaffold sequencing and comparison to the chicken sequence, should lead to an improved genome assembly.