A chromosome-level genome assembly of the Asian arowana, Scleropages formosus

Asian arowana (Scleropages formosus), an ancient teleost belonging to the Order Osteoglossomorpha, has been a valuable ornamental fish with some varieties. However, its biological studies and breeding germplasm have been remarkably limited by the lack of a reference genome. To solve these problems, here we report high-quality genome sequences of three common varieties of Asian arowana (the golden, red and green arowana). We firstly generated a chromosome-level genome assembly of the golden arowana, on basis of the genetic linkage map constructed with the restriction site-associated DNA sequencing (RAD-seq). In addition, we obtained draft genome assemblies of the red and green varieties. Finally, we annotated 22,016, 21,256 and 21,524 protein-coding genes in the genome assemblies of golden, red and green varieties respectively. Our data were deposited in publicly accessible repositories to promote biological research and molecular breeding of Asian arowana.

Asian arowana (Scleropages formosus), an ancient teleost belonging to the Order Osteoglossomorpha, has been a valuable ornamental fish with some varieties. However, its biological studies and breeding germplasm have been remarkably limited by the lack of a reference genome. To solve these problems, here we report high-quality genome sequences of three common varieties of Asian arowana (the golden, red and green arowana). We firstly generated a chromosome-level genome assembly of the golden arowana, on basis of the genetic linkage map constructed with the restriction site-associated DNA sequencing (RAD-seq). In addition, we obtained draft genome assemblies of the red and green varieties. Finally, we annotated 22,016, 21,256 and 21,524 protein-coding genes in the genome assemblies of golden, red and green varieties respectively. Our data were deposited in publicly accessible repositories to promote biological research and molecular breeding of Asian arowana.

Background & Summary
Asian arowana or Asian dragonfish (Scleropages formosus; Osteoglossidae, Osteoglossiformes, Osteoglossomorpha) represents an ancient lineage of teleost fish. The phylum Osteoglossomorpha is one of the basal branching lineages of teleosts with fossil records from the late Jurassic 1 . Because of its bright and shiny body colour, arowana has been considered as one of the most expensive ornamental fishes around the world 2 . At least five naturally occurring colour varieties of arowana have been identified. Among them, the red, golden and green varieties are the most popular. Within the three varieties, the red arowana has the highest commercial value, whilst the green one has the lowest price. In recent times the population size of arowana in its native habitat has declined due to overfishing. As a result, Asian arowana has been listed as an endangered species by the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) Appendix I since 1980 (ref. 3). It therefore raises people's awareness about this endangered ornamental fish to some extent. The Asian arowana, as a well-known representative of the Osteoglossiformes, possesses many primitive characters such as toothed tongue bone. On the other hand, arowana has evolved several derived traits such as mouth brooding and air-breathing function of the swim bladder. It is therefore considered as an important model for exploration of the teleost diversity. Furthermore, the genetic basis of body colour variations and the sex-determining mechanisms of arowana remain largely unknown.
Our previous study published in Scientific Reports 4 has reported a primary genome assembly of a female golden arowana and the draft genomes of the green and red varieties. The genomic and transcriptome comparisons among these three varieties have interpreted the possible molecular mechanism of colour variations. In addition, a potential sex chromosome was identified, revealing a solid clue for the sex-determination of arowana.
To construct these genomes, we extracted genomic DNAs from the three varieties of arowana and subsequently sequenced up to 100-fold coverage using the Illumina HiSeq2000 platform. After filtering low quality reads, we applied the SOAPdenovo2.2 (ref. 5) to assemble the clean reads of the three varieties separately. We obtained scaffold N50 values of 5.96 Mb (golden arowana), 1.63 Mb (red arowana) and 1.85 Mb (green arowana) respectively. The assembled genome sizes are 779, 753 and 759 Mb respectively, which are consistent with the estimates by k-mer analyses (Table 1). By using de novo-assembled transcipts and the CEGMA method 6 , we confirmed good completeness and accuracy of the three assembled genomes.
To construct chromosome-level assembly of the golden arowana, we performed restriction siteassociated DNA sequencing (RAD-seq), on basis of 94 F2 individuals from red grad 1 and Malaysian golden crossback, to generate a high-density genetic linkage map. Ultimately, we identified a total of 22,881 single-nucleotide polymorphisms (SNPs), in which 5,740 SNPs were clustered into 25 linkage groups (Fig. 1). The genetic linkage map spanned a genetic distance of 3,240 cM. According to the refined SNPs and their corresponding scaffolds, we calculated the size of the chromosomes-level assembly up to 683 Mb, which accounts for 87.7% of the achieved genome assembly. The high-quality chromosome-level assembly of golden arowana, along with the high-density genetic linkage map, will provide a valuable resource for further comparative genomic studies.
Transposable elements (TEs) in the three assembled genomes were examined. We observed that TEs account for 27, 27 and 28% of the genome assemblies in golden, red and green varieties, respectively (Table 1). Multiple methods for gene annotation, including ab initio, transciptome and homology based gene prediction, were employed to generate refined gene sets, which contain 22,016 (golden arowana), 21,256 (red arowana) and 21,524 (green arowana) protein-coding genes, respectively.
The availability of a high-quality reference genome of the golden arowana provides a good opportunity for biological characterization and molecular breeding of arowana.

Methods
These methods are expanded from detailed descriptions in our previous work 4 .  Fish samples for RAD sequencing. The Qian Hu fish farm (Singapore) obtained F1 hybrid individuals that originated from crossing two unrelated and genetically divergent founder (F0) Asian arowana grandparents (Red grade 1 × Malaysian golden varieties). Previously, we have generated two mapping families by crossing two pairs of these F1 hybrid brooders. The fostering methods of arowana are described in detail in our related work 4 . We used 94 F2 progenies and their parents for RAD sequencing to construct the high-density genetic map.
Library construction, sequencing and genome assembly Library construction and sequencing. High-quality genomic DNAs of three arowana varieties were extracted from the mixture of three tissues (muscle, skin and liver) independently using Puregene Tissue Core Kit A (Qiagen, MD, USA) for constructing libraries with different insert-sizes (170 bp to 40 kb; see more details in Table 2). Overall, we generated 16 pair-end libraries (8 for the golden arowana, 4 for the red arowana and 4 for the green variety, respectively) with the standard operating procedure provided by Illumina (San Diego, USA). Pair-end sequencing with the whole genome sequencing (WGS) strategy was performed on the Illumina HiSeq2000 platform using the standard operating procedure.
Raw data filtering. Raw reads were subjected to quality filtering (Raw Data Clean Procedure). Details for the disciplines of filtering were described in our related work 4 (Table 1).
RAD sequencing, genetic map construction, and chromosome-level assembly RAD sequencing. High-quality genomic DNAs were extracted from the scales and/or fin clips of the 94 progeny individuals and their parents for RAD sequencing by using Mag Attract HMW DNA Kit (Qiagen, MD, USA). The DNAs were subsequently digested with the restriction endonuclease EcoRI and divided into 3 RAD libraries 8 . Sequencing was performed by an Illumina HisSeq2000 platform. Finally, 72.8 Gb of reads with 101-bp length (the average depth is 1 × ) were obtained.
RAD SNP calling. After performing the above-mentioned Raw Data Clean Procedure to filter the adaptor sequences and remove low-quality reads, we re-aligned the clean reads onto the golden assembly (reference) using SOAP2 v2.21 (ref. 9) software with optimized parameters (-m 100 -× 888 -s 35 -l 32 -v 3 -p4). The SNPs were called using the Samtools 10 in each individual. Those low-quality SNPs with the missing rates higher than 30% were discarded.
Genetic map construction. SNPs showing significant Mendelian segregation distortion (χ2 test, P o0.01, d.f. = 1) were discarded. Then, the filtered SNPs were uploaded into JoinMap v4.1 (ref. 11) and analyzed with default parameters. The pairwise recombination estimations and logarithm of odds (LOD) scores of all SNPs were calculated, and SNPs pairs were then clustered into linkage groups at a LOD threshold of 10.0. Map distances (cM) were counted using the Regression mapping algorithm with the Kosambi function. All the SNP markers were clustered into 25 linkage groups (Fig. 1), which is consistent with our previously reported chromosome karyotype (2n = 48 and one additional W chromosome) 12 .
Anchoring the genome assembly to the linkage groups. A total of 5,740 SNP markers which located on 194 scaffolds were used for chromosomal-level assembly. The relative order between anchored scaffolds on each chromosome was determined by the following disciplines. For those scaffolds with the number of SNPs more than two, the two SNPs with the best quality on each scaffold were chosen to determine the order and the orientation. However, the orientation of those scaffolds with only one SNP was not ordered due to the lack of enough markers. Subsequently these ordered scaffolds were directly anchored to the chromosomes. Finally, 87.65% of the assembled genome sequences were able to be anchored onto the 25 pairs of chromosomes (Fig. 1, Table 3).

Genome annotation
In brief, we applied three independent methods to predict gene sets.
Homology annotation. We used the protein sequences from eight reported genomes (Table 4) to detect homology-based genes. All the protein sequences were downloaded from Ensembl (release 75). An all-against-all TblastN search was performed with an e-value of 10 − 5 . Best hits in each of the analyzed genome were searched, and the potential gene structures were then predicted by using Genewise2.2.0 (ref. 17). Those genes with length less than 150 bp, or prematurely terminated or frame-shifted, were discarded.
De novo annotation. We implemented de novo similarity-based gene prediction using AUGUSTUS v.2.5 (ref. 18) with optimized parameters (pre-trained with 1,500 randomly selected genes from the homology annotation set). To avoid pseudogene annotation, the repetitive regions of three arowana varieties were masked as 'N' seqeunces. Then AUGUSTUS v2.5 and GENSCAN v1.0 (ref. 19) were performed on the three-draft repeat-mask genome sequences. Subsequently, three gene sets were filtered using the same method applied for the homology annotation. Integration of annotation results. After merging all results generated from the above three methods, we applied GLEAN 22 to obtain a comprehensive and non-redundant gene set. Expression levels of the GLEAN genes and the alignments of Tophat were calculated using the Cuffdiff package of Cufflink 21 software with optimized parameters (-FDR 0.05 --geometric-norm TRUE -compatible-hits-norm TURE).
To find the best hit of each deduced protein, we employed BlastP with an e-value of 10 − 5 to map the protein sequences of GLEAN results against the SwissPort and TrEMBL databases 23 (Uniprot release 2011.06). The predicted protein sequences of the three arowana varieties were then aligned against the public databases (Pfam, PRINTS, ProDom and SMART) for detection of function motifs and domains. Ultimately, the GLEAN gene sets were filtered by the following three steps: 1) discarded genes with the length less than 150 bp, 2) discarded genes recognized as TEs, and 3) removed genes that were only obtained from the de novo method but without functional assignments and with an expression value lower than 1.

Data Records
In our previous work 4 , annotated genome sequences of the golden, red and green arowana varieties were uploaded to GenBank under the GenBank assembly accession numbers GCA_001624265.1 (Data Citation 1), GCA_001624255.1 (Data Citation 2) and GCA_001624245.1 (Data Citation 3), respectively. The RAD-seq raw read files were stored at NCBI SRA under experiment accession numbers SRX1728941 to SRX1728946 (Data Citation 4). The RNA-seq raw read files can be found at NCBI SRA SRX1668426 to SRX1668432 (Data Citation 5). In this paper, the gene annotation information of three varieties of arowana are available from Dryad (Data Citation 6). The chromosome annotation of golden arowana are available from Dryad (Data Citation 6). Data are given in tabular, tab-delimited value format.

Technical Validation
We took the following two steps to assess the generated genome assemblies. The first approach is Transcriptome evaluation. We used the Trinity 24 software to de novo assemble the RNA sequences of skin tissues from three varieties, and then utilized the BLAT 25 (E-value = 10e − 6 , identity = 90% and coverage>90%) to align the genome assemblies ( Table 5). The second approach is CEGMA 6 (Core Eukaryotic Genes Mapping Approach; http://korflab.ucdavis.edu/Datasets/genome_completeness, version 2.3) assessment with 248 conserved Core Eukaryotic Genes (CEGs) subjected to evaluation of the gene space completeness within the three assemblies. The evaluation results revealed a high coverage of more than 90% of gene coding-regions and over 95% of core eukaryotic genes, confirming the high-level completeness of the three genome assemblies (Table 6).

Usage Notes
Except for the Joinmap v4.1 that was ran on a Windows system, the other softwares were run on a linux system. The optimized parameters are provided in the main text.