High-coverage sequencing and annotated assemblies of the budgerigar genome

Raw genome DNA sequence reads

DNA samples were obtained from a blood sample taken from a single male Melopsittacus undulatus, who we aptly named Mr. B. For Illumina sequencing, reads were generated at Duke University (16×), Illumina UK (54×), and BGI (219×) using Illumina’s TruSeq [1] version2 or version3 chemistries (Table 1 and GigaDB [2]). The version3 chemistry reads through GC-rich regions, which are often found in promoters, more evenly than does version2 [3]. The insert sizes for the BGI libraries ranged from 220 bp to 40 Kbp, and the insert sizes for the Duke libraries ranged from 400–600 bp, in order to assist assemblies. Fragment sizes for the mate pair libraries, based on genome mapping, and the per base sequence quality distribution for the libraries are shown in GigaDB [2]. The Duke University Illumina libraries were sequenced at two different cluster densities: 8× coverage reads at the normal 420 k clusters/mm density and 8× coverage at a lower 350 k clusters/mm. The lower cluster density was used to increase the number of GC-rich regions sequenced. For PacBio sequencing, 6.76 Gbp (~5.5× coverage) of PacBio RS reads [4] were generated at Pacific Biosciences from two insert size libraries (7.5 K bp at 1.93× and 13 Kbp at 3.56×; PacBio reads error-corrected with Illumina can be downloaded from the supplementary webpage associated with [5]). With all reads combined, the total coverage exceeds 300× (assuming a haploid genome size of 1.23 Gbp) (Table 1), perhaps making Mr. B one of the most sequenced individual vertebrate animals as of to date. The read length distributions of these different types of reads are shown in Figure 1.

Fosmid Library

To validate the assemblies in the Assemblathon 2 competition, a fosmid library was created from sheared genomic DNA (35–40 Kbp) of Mr. B [6]. Ten pools of clones were generated and sequenced using Illumina as described in [7]. Each pool of reads was individually assembled using Velvet [8]. The fosmid assemblies have been deposited at GigaDB [2]).

Transcriptome Reads

454 FLX transcriptome reads were generated from brain RNA isolated from two males, neither of whom was Mr. B. An initial set of sequencing runs of both males was conducted at Washington University at St. Louis, producing 89.2 Mb of transcriptome sequence as reported in [9] (NCBI accession numbers SRR029329–30) and were assembled using Newbler [10] into 19,198 contigs. An additional 21× coverage (run label GK0K2XF01) was generated at Duke University from one of the males.

Assemblies

We present three hybrid assemblies: 1) Budgerigar 454-illumina hybrid v6.3 using the CABOG assembler; 2) Budgerigar PBcR hybrid using the CABOG assembler; and 3) Budgerigar illumina-454 hybrid using the SOAPdenovo2 assembler. The first two assemblies were annotated, after which, optical-map assisted megascaffolds were constructed based on them. As of yet, the SOAPdenovo2 assemblies have not been annotated or aligned to optical maps. The quality statistics of these assemblies are in listed in Table 2, and brief descriptions of their construction and relative quality are provided in Additional file 1.

Validating sequence assemblies with optical maps

Optical Mapping is a single molecule system for the construction of ordered restriction maps of whole genomes [11], and it has been used to guide and validate sequence assemblies [12]. An optical map for the budgerirgar genome was created, using a method described in Additional file 1. The optical map contigs ranged in size from 2 Mbp to 74 Mbp and spanned over 900 Mbp with a resolution of 13.94 Kbp (i.e., one non-redundant SwaI every 13.94 Kbp). The contigs were then aligned to in silico restriction maps generated from Budgerigar_v6.3 and PBcR assembly scaffolds in order to validate the scaffolds. An approximate 859.21 Mb of the optical maps aligned to the Budgerigar_v6.3 assembly, in 146 scaffolds with 3 or more SwaI restriction fragments (excluding ends and fragments less than 0.4 Kbp). Of these 146 scaffolds, 43 appeared chimeric (i.e., aligned to two or more optical map contigs). For the PBcR assembly, 796.63 Mbp optical map contigs aligned, in 673 scaffolds. Of the 673 scaffolds, only 51 were chimeric. Thus, while the Budgerigar_v6.3 assembly has a higher N50 scaffold metric and hence longer scaffolds compared to the PBcR assembly, 30% the v6.3 scaffolds are chimeric, whereas only 7.6% of the PBcR assembly are chimeric.

Optical map assisted assemblies

We took both Budgerigar_v6.3 and PBcR assemblies and filtered out alignments that did not extend to the end of either the genomic sequence scaffold or the optical map. The remaining high-quality alignments were then used to identify optical map alignments that bridged scaffolds, such that a single optical map aligned to the ends of at least two sequence scaffolds. We then iteratively extended the megascaffolds beyond pairs of sequence scaffolds, using three heuristics: (1) we limited the overhangs (i.e., the portion of the scaffold sequence that does not align to the optical map) to 2 Mbp total; (2) we bridged two scaffolds together only if the size of the gap separating them is less than 2 Mbp of Ns; and (3) if a sequence scaffold aligned to more than one optical map, we placed it into the largest one it aligns with. The above procedure slightly reduced the number of scaffolds from 25,212 to 25,163 in the Budgerigar_v6.3 assembly, and from 54,668 to 54,138 in the PBcR assembly. This relatively small change in number is expected as our procedure tended to join only sequence scaffolds that were already fairly large into even larger megascaffolds, since it is only possible to confidently align an optical map to a fairly large sequence scaffold bearing numerous SwaI restriction sites. However, this analysis substantially improved the scaffold N50 sizes from 10.6 Mbp to 13.8 Mbp in the Budgerigar_v6.3, and 1.7 Mbp to 7.3 Mbp in the PBcR assemblies, respectively (Table 2). Without limiting the length of the overhangs and gap sizes to 2 Mbp, the increase in N50 scaffold sizes in the Budgerigar_v6.3 is 17.1 Mbp (which we think could be an artifact). We speculate that some of the large gaps in the optical map correspond to centromeres or highly repetitive DNA that are difficult to assemble.

Annotations

The Budgerigar_v6.3 and PBcR assemblies were annotated at BGI for protein coding genes by first generating a reference set of human, chicken and zebra finch proteins, and then aligning the reference set to the assemblies, and propagating annotations to 30% coverage of the reference at TBlastN, E = 1e⁻⁵. For the Budgerigar_v6.3 assembly, the reference set comprised of human proteins from Ensembl 60 and a set of zebra finch and chicken proteins re-annotated based on these human proteins, using a custom BGI pipeline reported on separately (Jarvis et al. in preparation; Zhang et al., in preparation). For the PBcR assembly, the reference set comprised of the Ensembl 60 human, chicken and zebra finch proteins. The propagation of these reference sets to the budgerigar assemblies is described in more detail in Additional file 1. Further, in the PBcR assembly, UTRs were annotated for 6,203 genes using the GK0K2XF01 transcriptome runs with a pipeline similar to the one described in [13]. The assembly annotations were then propagated to the corresponding sets of megascaffolds. No de novo gene annotations were performed.

The annotated Budgerigar assemblies had fewer genes (15,470 and 16,204 genes in the Budgerigar_v6.3 and PBcR assemblies respectively) than the published Zebra Finch (18,618 genes) and Chicken genome assemblies (17,108 genes in the 2011 Galgal4 assembly [14]). We believe the lower number of annotated genes in budgerigar assemblies is due to the differences in annotation methods rather than assembly completeness, for two reasons: (1) These annotations were produced based on similarities to zebra finch, chicken and human proteins, and hence they cannot contain more genes than the source genome annotations; and (2) The independent GenScan annotation of the Budgerigar_v6.3 assembly at the UCSC Genome Browser contains more genes than in zebra finch and chicken, 24,095 in total.

Comparisons to other avian assemblies

Our budgerigar genome assemblies were compared with the zebra finch, chicken, and falcon genomes [15-17]. The other assemblies from the Assemblathon 2 competition are available from GigaDB [18]. The zebra finch and chicken had similar contig and scaffold N50 values (38.5 kb and 10.4 Mb for zebra finch, and 279.8 kb and 90.2 Mb for chicken, respectively). In addition, since the Peregrine Falcon is the closest relative to parrots [19], we also compared the budgerigar genome assemblies to this bird. However, it was not possible to do an in depth comparison of these genomes to the recently sequenced Scarlet Macaw and Puerto Rican Parrot genomes [20,21], because both bird genomes had N50 scaffold sizes under 20,000 and N50 contig sizes under 7,000. A summary of assemblies, including the Scarlet Macaw and Puerto Rican Parrot, are shown in Table 2. Apart from the standard genome assembly quality statistics, we assessed the quality of the budgerigar assemblies along two other dimensions: (1) the coverage of highly conserved avian exons, and (2) the number of gaps 10 Kbp upstream and downstream of each gene (gene territories), and conversely, the number gene territories assembled without gaps. Of 3,288 highly conserved exons (>86% coverage across >87% of their length) we identified between chicken and zebra finch, 3,165 (96.25%) and 3,134 (95.31%) were covered with >86% identity across >87% of their length in the Budgerigar_v6.3 and PBcR assemblies respectively, pointing to good coverage of coding regions in these assemblies. The budgerigar assemblies had fewer gaps within the coding sequences and gene territories than all other avian genomes examined, except the newer unpublished Galgal4 chicken assembly that is similar to the budgerigar in that it is a hybrid that includes both short and long sequences (Sanger and 454 FLX+) (Figure 2). This suggests that our budgerigar assemblies have very well assembled genes and promoter regions.

Using the online CoGe tool [22-24], we assessed the structural similarities between the various budgerigar assemblies and other avian assemblies [25-30], by computing the level of coding sequence synteny among assemblies. The highest numbers of genes in synteny were observed, as expected, between a budgerigar assembly and the optical map assisted version of the same assembly (Figure 3A). However, the number of genes in synteny between the Budgerigar_v6.3 and the PBcR assemblies was similar to the number of genes in synteny between budgerigar and falcon (Figure 3A, B). Further, the number of genes in synteny did not strictly reflect phylogenetic relationships, as the zebra finch and budgerigar, close relatives [19], had a lower level of synteny than budgerigar and chicken. In addition, a number of inversions were observed even in the syntenic dotplots between the original budgerigar assemblies and their optical map-assisted assemblies (88 inversions between Budgerigar_v6.3 and Budgerigar_v6.3_mega; 209 inversions between PBcR and PBcR_mega, plots shown in GigaDB [2]). This suggests that synteny based on CoGE syntenic maps is affected by the quality of the assemblies and the characteristics of the synteny algorithm. Thus, the number of genes in synteny computed using the available methods is only a rough measure of the actual structural similarity between the assemblies compared.

In summary, this study shows that the budgerigar genomic resource we have generated has provided [5,6] (and is still expected to provide more) valuable data and material for genome technology development and for further investigating complex behavioral traits at the genomics level.

All procedures on live animals were approved by the Institutional Animal Care and Use Committee of Duke University.