High-coverage sequencing and annotated assembly of the genome of the Australian dragon lizard Pogona vitticeps

Background The lizards of the family Agamidae are one of the most prominent elements of the Australian reptile fauna. Here, we present a genomic resource built on the basis of a wild-caught male ZZ central bearded dragon Pogona vitticeps. Findings The genomic sequence for P. vitticeps, generated on the Illumina HiSeq 2000 platform, comprised 317 Gbp (179X raw read depth) from 13 insert libraries ranging from 250 bp to 40 kbp. After filtering for low-quality and duplicated reads, 146 Gbp of data (83X) was available for assembly. Exceptionally high levels of heterozygosity (0.85 % of single nucleotide polymorphisms plus sequence insertions or deletions) complicated assembly; nevertheless, 96.4 % of reads mapped back to the assembled scaffolds, indicating that the assembly included most of the sequenced genome. Length of the assembly was 1.8 Gbp in 545,310 scaffolds (69,852 longer than 300 bp), the longest being 14.68 Mbp. N50 was 2.29 Mbp. Genes were annotated on the basis of de novo prediction, similarity to the green anole Anolis carolinensis, Gallus gallus and Homo sapiens proteins, and P. vitticeps transcriptome sequence assemblies, to yield 19,406 protein-coding genes in the assembly, 63 % of which had intact open reading frames. Our assembly captured 99 % (246 of 248) of core CEGMA genes, with 93 % (231) being complete. Conclusions The quality of the P. vitticeps assembly is comparable or superior to that of other published squamate genomes, and the annotated P. vitticeps genome can be accessed through a genome browser available at https://genomics.canberra.edu.au. Electronic supplementary material The online version of this article (doi:10.1186/s13742-015-0085-2) contains supplementary material, which is available to authorized users.

accordance with the Australian Capital Territory Animal Welfare Act 1992 and the approval of the University of Canberra Animal Ethics Committee. DNA was extracted and purified using standard protocols and transported to BGI-Shenzhen, China for sequencing. 13 insert libraries were constructed with insert sizes of 250 bp, 500 bp, 800 bp, 2 kbp (x2), 5 kbp, 6 kbp, 10 kbp (x2), 20 kbp (x2) and 40 kbp (x2) and subjected to paired-end sequencing on an Illumina HiSeq 2000 platform to generate 317 Gbp of raw sequence (Table 1). After filtering for low-quality reads and duplicated reads arising from Read depth was calculated on the basis of a genome size of 1.77 Gbp. Average read depth, number of times on average a particular base is included in a read. Physical coverage, the number of times on average a particular base is spanned by a paired read  35 Gbases of clean data in the form of 125 bp reads, which generated 76.89x109 17-mer sequences. The solid line shows the k-mer spectrum (percentage frequency against k-mer copy number). The second mode (copy number 48.5) represents homozygous single copy sequence, whereas the first mode (24.5), half the copy number of the first, represents heterozygous single copy sequence. Heterozygosity is high, which complicated assembly Gbp of data were retained for genome assembly. This amount of data represents an average read depth of 82.7 (Table 1), assuming a genome size of 1.81 pg, as estimated for a female P. vitticeps by flow cytometry [5]. This mass converts to a genome size of 1.77 Gbp [6]. Reads from the short-insert libraries (250, 500 and 800 bp) were decomposed into short sequences of length k (k-mers, with k = 17) using Jellyfish version 1.1 [7]. The histogram of k-mer copy number ( Fig. 1) was strongly bi-modal, the first mode with a copy number that was half that of the second, which reflects the high level of heterozygosity in this wild-caught lizard (0.85 % of single nucleotide polymorphisms [SNPs] plus sequence insertions or deletions [indels]). The second mode in the k-mer graph was used to obtain an estimate of the genome size using the formula: where L is the read length (125 bp), k is the k-mer length (17 bp), there are 98. 35 Gbp of sequence data, and the mode is taken from the k-mer graph (48.5, Fig. 1). Our sequence-based estimate of 1.768 Gbp agrees well with the estimate of 1.77 Gbp for the ZW genome that was previously made using flow cytometry data [4].

Assembly
Assembly was performed with SOAPdenovo (version 2.03) [8,9]. Briefly, the sequences derived from the shortinsert libraries were decomposed into k-mers to construct the de Bruijn graph, which was simplified to allow connection of the remaining k-mers into contiguous sequence (contigs). We then aligned the paired-end reads from small and large-insert library sets to the contigs, calculated the support for relationships between contigs, assessed the consistent and conflicting relationships, and constructed scaffolds. Finally, we retrieved paired reads that mapped to a unique contig but had the other member of the pair located in a gap region. Reads falling in the same gap region were then assembled locally.   genome sequence is contained in contigs, or scaffolds, equal to or greater than this length), with unclosed gap regions representing only 3.78 % of the assembly ( Table 2). Reads from small-insert libraries that satisfied our filtering criteria were aligned to the assembly using the Burrows-Wheeler Aligner (BWA, version 0.5.9-R16) [10], allowing for eight mismatches and one indel. Of the total number of reads (797.4 M), 96.4 % could be mapped back to the assembled genome and they covered 98.4 % of the assembly excluding gaps. Bases in the assembled scaffolds had, on average, reads mapped with 55X read depth. These data suggest that we have assembled most of the P. vitticeps genome. In addition, we used the CEGMA package (version 2.4) [11] to map 248 core eukaryotic genes to our P. vitticeps assembly. Our assembly captured 99 % (246 of 248) of the core CEGMA genes, with 93 % (231) being complete. This is a higher assembly rate than that estimated for the green anole Anolis carolinensis assembly (AnoCar2.0), which captured 93.6 % (232) of the core genes, with 85.9 % (213) being complete.

Transcriptomes
We generated transcriptome data from the brain, heart, lung, liver, kidney, skeletal muscle and gonads of male and female P. vitticeps ( Table 3). None of the seven animals from which we collected RNA was used in generating the genome sequence. Two sets of sequencing runs on two different male and female individuals were performed by BGI-Shenzen, producing 309,436,077 90 bp paired-end reads (ENA accession numbers ERR753524-ERR753530 and ERR413064-ERR413076). A third set of samples was sequenced by The Ramaciotti Centre, University of New South Wales, Australia, including a sex-reversed ZZ female, producing 89,687,526 101 bp paired-end reads (ENA accession numbers ERR413077-ERR413082). We assembled these datasets (from all seven individuals) into 595,564 contigs using Trinity (release r2013_08_14) [12] with default parameters (ENA accession number ERZ097159). Only the first set of RNA-seq reads was available for genome annotation (ENA accession numbers ERR753524-ERR753530) but we make the entire dataset, including our de novo assembly, available with this article (see ' Availability of supporting data' section).

Annotation
Transposable elements and other repetitive elements were identified using a combination of homology, at both the DNA and protein levels, and de novo prediction. In the homology-based approach, we searched Repbase [13] for known transposable elements, used RepeatMasker [14] for DNA homology search against the Repbase database, and used WuBlastX to search against the transposable element protein database provided within RepeatProteinMask (bundled in RepeatMasker). In the de novo approach, we used RepeatModeler [15] and LTR_FINDER [16] to predict repetitive elements. Tandem repeats were identified using Tandem Repeats Finder [17]. The relative success of the different approaches is shown in Table 4. Overall, we identified about 690 Mbp of repetitive sequences accounting for 39.47 % of the genome, of which the predominant elements were long interspersed nuclear elements (LINEs, which accounted for 33 % of repetitive sequences representing 12.2 % of the genome) ( Table 5).
We combined homology-based, de novo and transcriptome-based methods to predict gene content of the assembly. In the homology-based prediction, the assembly was annotated by generating reference sets of A. carolinensis, Gallus gallus and Homo sapiens proteins, and aligning the reference sets to the assembly using TBLASTN (version 2.2.23; E-value ≤ 1 × 10 −5 ). The resultant homologous genome sequences were then aligned against matching proteins using Genewise (version wise2-2-0) [18] to define gene models. In the de novo prediction, we randomly selected 1,000 genes with intact open reading frames (ORFs) as predicted by the homology-based approach to train the Augustus gene prediction tool (version 2.5.5) [19] with the parameters appropriate to P. vitticeps. The de novo gene prediction was then performed with Augustus applied to the genome after repeat sequences were masked as described above. In the transcriptomebased approach, we mapped transcriptome reads to the assembly using TopHat (version 1.3.1) [20], which can align reads across splice junctions. These mapped reads were assembled into transcripts using Cufflinks (version 1.3.0) [21] and then merged across samples (n = 7, Table 3) into a single transcriptome annotation using the Cuffmerge option.
The results of the three approaches were combined into a non-redundant gene set of 19,406 protein-encoding genes, 63 % of which included intact ORFs ( Table 6). Most of the predicted genes were supported by RNA-seq signals ( Table 3).
To assign gene names to each predicted protein-coding locus, we mapped the 19,406 genes to an Ensembl library collated from A. carolinensis, chicken G. gallus, human H. sapiens, western clawed frog Xenopus tropicalis and zebrafish Danio rerio. The name associated with the best hit for each P. vitticeps gene was assigned to each of 19,083 genes. Most of these genes (16,510) mapped to a homolog even at high stringency (>80 % of protein length aligned).

Bacterial artificial chromosome library
A large-insert genomic DNA bacterial artificial chromosome (BAC) library was constructed from DNA from a wild-caught female dragon lizard (TC1542) confirmed to have the ZW genotype using sex-linked PCR markers [3,4] and cytologically [3]. The library is estimated to represent 6.3× of genome coverage, and is comprised of 92,160 clones with an average insert size of 120 kbp. This resource is commercially available through Amplicon Express (Pullman, WA, USA; http://ampliconexpress.com).

Sex chromosome sequences
The sex of P. vitticeps is determined by a combination of chromosomal constitution and influence of environmental   temperature on the developing embryo. P. vitticeps has female heterogamety (with ZZ male and ZW female individuals), and the Z and W chromosomes are among the ten pairs of microchromosomes [2]. Sex chromosome heteromorphy is evident by C-banding, but the degree of differentiation of the Z and W chromosomes is slight [2]. The sex chromosomes of P. vitticeps are not homologous to the sex chromosomes of chicken (G. gallus) or other reptiles so far examined [22]. The ZZ genotype is reversed to a female phenotype at high incubation temperatures [3,4]. Our laboratory has previously identified a sex-linked sequence using amplified fragment length polymorphism screening and genome walking [4,23]. Five contiguous BAC clones containing sex-linked markers that map to the sex chromosome pair were sequenced to reveal 352 kbp of P. vitticeps sex chromosome sequence [24]. This region contained five protein-coding genes (oprd1, rcc1, znf91, znf131 and znf180) and several major families of repetitive sequences (long terminal repeat [LTR] and non-LTR retrotransposons, including chicken repeat 1 [CR1] and bovine B LINEs [Bov-B LINEs]) [1,24].
More recently, we amplified micro-dissected Wchromosome fragments to yield many sex chromosome sequence tags that were reciprocally mapped to their Z homologs (Matsubara et al., unpublished data). All putative sex chromosome scaffolds were confirmed to co-localize with the known ZW-BAC Pv3-L07 when physically mapped (Deakin et al., unpublished data). In this way we identified 12.8 Mbp of the Z chromosome (on three scaffolds) and increased the number of confirmed sex chromosome genes to 240 (Deakin et al., unpublished data).
The macrochromosomes of P. vitticeps are largely devoid of variation in GC content at small (5 kbp) spatial scales. In fact, P. vitticeps macrochromosomes are more uniform in terms of GC distribution than is the A. carolinensis genome (standard deviation 0.029 versus 0.032 respectively; Table 7). With the exception of the Z microchromosome, P. vitticeps microchromosomes possess a heterogeneous distribution of GC-rich sequences over 5 kbp windows (Fig. 4a). In this regard, P. vitticeps microchromosomes resemble those of birds but differ markedly from those of A. carolinensis, whose GC content more closely resembles that of the macrochromosomes [33]. Intriguingly, the Z microchromosome of P. vitticeps has an average GC content comparable to that of coding sequences and short interspersed nuclear elements (SINEs) (Fig. 4b), which suggests that this chromosome may be enriched in these GC-rich components of the genome. However, with the exception of LTR transposable elements, all components (CDS, introns, tandem repeats and transposable elements) showed greater GC content if they resided on the Z chromosome than elsewhere (Z chromosome > microchromosomes > macrochromosomes; Fig. 4c), suggesting that there are other, as yet unidentified, reasons for the observed variation in GC content across different chromosome classes.
When variation in GC distribution is considered over larger spatial scales (tens to hundreds of kbp, Fig. 3), the P. vitticeps macrochromosomes are similar to the A. carolinensis genome, which lacks substantial variation in GC composition, a striking departure from isochore patterns seen in mammals and birds [33]. The Z chromosome, too, lacks substantial heterogeneity over larger spatial scales, which perhaps reflects the uniform distribution of repetitive elements along its length [24]. Only the autosomal microchromosomes of P. vitticeps bear any similarity in GC distribution to the other sauropsid genomes examined. The P. vitticeps genome, therefore, has compositional patterns distinct from that of A. carolinensis, which indicates that different processes have shaped the genomes of the two lizards since they shared a common ancestor 144 million years ago.
Comparison with other assemblies P. vitticeps and A. carolinensis had similar scaffold N50 values (2.29 Mbp and 4.03 Mbp, respectively). These values for P. vitticeps are surprisingly good, as its genome was assembled from short read sequences, whereas that of A. carolinensis was generated using Sanger sequencing. Our assembly compares well to nine other sauropsid genomes, including those of two squamates, two turtles and three crocodilians ( Table 8).
The gene parameters listed in Table 6 compare well to those of other vertebrates (see also Fig. 5).

Concluding remarks
The quality of the P. vitticeps assembly is comparable to that of other published squamate genomes. This genome assembly, coupled with the availability of a BAC library and the development of a high-density physical map for each chromosome, provides an unparalleled resource for accelerating research on sex determination, major histocompatibility complex evolution, and the evolution of adaptive traits in squamates to complement the advances brought about by the sequencing of the A. carolinensis genome [25].

Availability of supporting data
The genomic and transcriptomic sequence reads and assemblies have been deposited in the ENA under the project accession number PRJEB5206 (see Additional file 1 for a complete list of accession numbers). The genome sequence has been submitted to GigaDB [34] along with other supporting resources, including: SoapDeNovo2 pvi1.1.Jan2013 genome assembly (ENA accession number ERZ094017) Trinity de novo transcriptome assembly (ENA accession number ERZ097159) Peptide and coding sequences for the pvi1.1.Jan2013 assembly Gene annotations and repeat annotations for the scaffolds Sequence Read Archive accession numbers for all sequencing runs.
The annotated P. vitticeps genome sequence can be accessed through a publicly available genome browser [35].

Additional file
Additional file 1: ENA accession numbers. (XLSX 18 kb) Abbreviations BAC: Bacterial artificial chromosome; bp: Base pair; ENA: European Nucleotide Archive; indel: Sequence insertion or deletion; k-mer: Short sequence of length k; LINE: Long interspersed nuclear element; LTR: Long terminal repeat; N50: 50 % of the genome sequence is contained in contigs (or scaffolds) equal to or greater than this length; ORF: Open reading frame; SINE: Short interspersed nuclear element; SNP: Single nucleotide polymorphism.
undertook the gene and repeat annotation. SS, JAMG and AG provided access to the P. vitticeps BAC library and associated resources. JD oversaw the cytological work of KM, who undertook the work necessary to identify the chromosomal sex of Fabian, complemented by work of CH and XZ using sex-linked markers. DO'M and HP handled the bioinformatics work required to deliver the genome sequence through the browser. PZ analysed the transcriptomes and construct transcript models that were used for gene annotation. MF and YZ undertook the CG and isochore analyses. All authors contributed to the preparation of the manuscript. All authors read and approved the final manuscript.