Chromosome genome assembly and annotation of the yellowbelly pufferfish with PacBio and Hi-C sequencing data

Pufferfish are ideal models for vertebrate chromosome evolution studies. The yellowbelly pufferfish, Takifugu flavidus, is an important marine fish species in the aquaculture industry and ecology of East Asia. The chromosome assembly of the species could facilitate the study of chromosome evolution and functional gene mapping. To this end, 44, 27 and 50 Gb reads were generated for genome assembly using Illumina, PacBio and Hi-C sequencing technologies, respectively. More than 13 Gb full-length transcripts were sequenced on the PacBio platform. A 366 Mb genome was obtained with the contig of 4.4 Mb and scaffold N50 length of 15.7 Mb. 266 contigs were reliably assembled into 22 chromosomes, representing 95.9% of the total genome. A total of 29,416 protein-coding genes were predicted and 28,071 genes were functionally annotated. More than 97.7% of the BUSCO genes were successfully detected in the genome. The genome resource in this work will be used for the conservation and population genetics of the yellowbelly pufferfish, as well as in vertebrate chromosome evolution studies.


Background & Summary
The yellowbelly pufferfish (FishBase ID: 24266), Takifugu flavidus, is an economically and ecologically important fish species in coastal regions of East Asia, including the East China Sea, Yellow Sea and Bohai Bay 1,2 . The yellowbelly pufferfish is also a temperate bottom fish that exhibits only short-distance seasonal migration 3 . The yellowbelly pufferfish is caught and cultivated as a delicious fish species with high market value 2,4,5 . However, due to environmental deterioration and overfishing, the wild populations of the species have declined in the last decade 6,7 . Additionally, a low survival rate in artificial breeding has greatly limited the development of the marine aquaculture of the yellowbelly pufferfish 8,9 . Previous studies of the yellowbelly pufferfish have mainly focused on behavioural 1 , morphological and growth characteristics, temperature and salinity effects on embryos and larval development 9 , and molecular marker development 10 . A genome of Takifugu flavidus was published in 2014 11 ; however, this genome was a fragmented draft with contig and scaffold N50 lengths of 2.7 kb and 305.7 kb, respectively. A high-quality reference genome of the yellowbelly pufferfish could facilitate and prompt conservation genetics research and investigation of the molecular mechanisms of important economic traits of the species.
The genomes of pufferfish have also played an important role in studies of vertebrate genome evolution due to the compactness of genus Takifugu genomes [12][13][14][15] . Previous studies have shown that the number of repetitive elements in pufferfish is significantly reduced 14,15 . The genome of Takifugu rubripes, the other pufferfish in genus Takifugu, exhibits conserved linkages with humans, implying the preservation of chromosomal segments from the common vertebrate ancestor [13][14][15] . Another pufferfish genome, that of Tetraodon nigroviridis, was the second pufferfish genome reported 14 . Although the yellowbelly pufferfish genome has been published, the genome was constructed using short reads from SOLiD next-generation sequencing, and the sequences have not been assembled at the chromosomal level 11 . The elucidation of genomic evolution among pufferfish species, and comparison with other vertebrates such as humans, will require further chromosomal genome assemblies for pufferfish species.
In this work, we applied a combined strategy involving Illumina, PacBio and Hi-C technologies to generate sequencing data for chromosomal genome construction for the yellowbelly pufferfish ( Fig. 1). More than 98% of BUSCO genes were detected, and the contig and scaffold N50 lengths reached 4.4 and 15.7 Mb, respectively, demonstrating the outstanding completeness and sequence continuity of the genome. A total of 29,416 protein-coding genes were predicted in the assembled genome, and more than 95% of those genes were successfully functionally annotated. We believe that the chromosomal genome assembly constructed in this work will not only be valuable for ecology, conservation and aquaculture studies of the yellowbelly pufferfish but will also be of general interest in the evolutionary investigation of teleosts and vertebrates.

Methods
Sample collection. A female yellowbelly pufferfish (Fig. 2), reared in the fish breeding centre of Fujian Normal University in Fuzhou City of Fujian Province was used for genome sequencing and assembly. Fresh white muscle, eye, skin, gonad, gut, liver, kidney, blood, gall bladder and air bladder tissues were collected and quickly frozen in liquid nitrogen for one hour. White muscle tissues were used for DNA sequencing for genome assembly, while all tissues were used for transcriptome sequencing.
DNA and rNA sequencing. Genomic DNA from white muscle tissue was extracted using the standard phenol/chloroform extraction method for DNA sequencing library construction. The integrity of the genomic DNA molecules was checked using agarose gel electrophoresis. Both the Illumina HiSeq X Ten platform and the PacBio SEQUEL platform were applied for genomic sequencing to generate short and long genomic reads, respectively. For the Illumina X Ten platform (San Diego, CA, USA), a paired-end library was constructed with an insert size of 250 base pairs (bp) according to the protocol provided by the manufacturer. As a result, 44 Gb (~120X coverage of the estimated genome size, Table 1) of accurate short reads were generated, which were further cleaned using the HTQC utility 16 . Adapter sequences and reads with more than 10% N bases or low-quality bases (≤5) were filtered from the sequencing data. After filtering, 41.8 Gb (~110X, Table 1) of cleaned data were retained for the following analysis. To obtain sufficient sequencing data for genome assembly, we constructed two 20 kb DNA libraries using the extracted DNA and the standard Pacific Biosciences (PacBio, Menlo Park, CA) protocol, and fragments were selected using the Blue Pippin Size-Selection System (Sage Science, MA, USA). The library was sequenced using the PacBio SEQUEL platform. After removing adaptors, we obtained 27.2 Gb subreads (~73X,   Table 1) for genome assembly. The genomic sequencing data used for subsequent genome assembly are summarized in Table 1.
We also performed RNA sequencing to generate transcriptome data for gene model prediction. To include as many tissue-specific transcripts as possible, multiple tissues were collected, as indicated in the Sample Collection section. TRIzol reagent (Invitrogen, USA) was used to separately extract RNA from all of the collected tissues, including white muscle, ocular, skin, gonad, intestine, liver, kidney, blood, gall bladder and air bladder tissues. RNA quality was checked with a NanoDrop ND-1000 spectrophotometer (Labtech, Ringmer, UK) and a 2100 Bioanalyzer (Agilent Technologies, CA, USA). Then, RNA molecules were equally mixed for transcriptome sequencing on the PacBio SEQUEL platform. First, cDNA was prepared using the SMARTer PCR cDNA Synthesis Kit (Clontech) from 1 μg of purified RNA. The Iso-Seq libraries were constructed from the BluePippin (Sage Science, MA, USA) size-selected cDNA with a size range of 0.6-3 kb according to the PacBio SEQUEL library construction protocol. Two SMRT flow cells were used for long-read transcriptome sequencing, and the resulting data used for gene prediction are summarized in Table 1.
De novo assembly of the yellowbelly pufferfish genome. For the Next Generation Sequencing (NGS) short reads, the Kmer-based method 17 was used to perform genome survey analysis to estimate the genome size, heterozygosity and repeat content of the yellowbelly pufferfish genome. We counted the number of each 17-mer with Jellyfish 18 , and the frequency distribution is plotted in Fig. 3. The yellowbelly pufferfish genome size was then estimated from the frequency distribution to be 377 Mb.
The Falcon package 19 was used to assemble the yellowbelly pufferfish genome with PacBio long reads, using the parameters length_cutoff = 10 kb and pr_length_cutoff = 8 kb. The procedures of long-read self-correction, corrected read alignment, sequence graph construction, and contig assembly were performed in Falcon. To correct random sequencing errors in the Falcon output, two steps of genome sequence polishing were applied: we first used arrow 20 to polish the genome using long sequencing data, and two rounds of polishing using NGS short reads were then applied with Pilon 21 . Finally, we obtained a final contig assembly of 366 Mb with a contig N50 length of 4.4 Mb ( Table 2).
To evaluate the quality of the assembled genome, the completeness and accuracy were assessed via BUSCO analysis and short-read mapping. The completeness of the assembled yellowbelly pufferfish genome was assessed by using BUSCO v3.0 22 with the vertebrata_odb9 database. We found that 95.7% and 2.7% of 2,586 BUSCO genes were completely and partially BUSCO genes were detected in the genome. We also aligned NGS short reads to the genome and found that more than 98.5% of the reads were reliably aligned, showing a high mapping ratio for the short-read sequencing data.
Chromosome construction using interaction information from Hi-C data. In this work, we applied the Hi-C technique for chromosome construction for the yellowbelly pufferfish. Although the Hi-C technique was first introduced to quantify genome-wide interactions 23 , the method exhibits suitability for chromosome assembly and has been successfully applied in many genomic projects 24 . In our study, we used 0.2 ml of blood from the same individual used for genome sequencing for Hi-C library construction and sequencing using the  www.nature.com/scientificdata www.nature.com/scientificdata/ same method as in a previous study 25,26 . From the Hi-C library sequencing, approximately 50.7 Gb of data were generated ( Table 1). The sequencing reads were mapped to the polished yellowbelly pufferfish genome with Bowtie 1.2.2 27 . We independently aligned the two read ends to the genome and only selected the read pairs for which both ends were uniquely aligned to the genome. The hiclib Python library 28 and a previously reported method 24 were applied to filter the Hi-C reads, and the interaction frequency was quantified and normalized among contigs. Lachesis 29 with default parameters was then applied to cluster contigs with the agglomerative hierarchical clustering method using the interaction matrix between sequences. Among the 169 million read pairs generated from Hi-C sequencing, 59 million read pairs (34.9%) provided valid interaction information for chromosome assembly. As a result, the contigs from the yellowbelly pufferfish were successfully clustered into 22 groups, which were further ordered and oriented into chromosomes. Finally, 271 contigs were reliably anchored on chromosomes, accounting for 95.9% of the total genome. The contig and scaffold N50 values reached 4.4 and 15.7 Mb (Tables 2 and 3), respectively, providing the first chromosomal genome assembly for the yellowbelly pufferfish.

Gene model prediction and functional annotations.
Repeat elements were annotated in the yellowbelly pufferfish before gene model annotation. We applied Tandem Repeat Finder (TRF) 30 , LTR_FINDER 31 , PILER 32 and RepeatScout 33 for the ab initio prediction of repeat elements in the genome. Thereafter, RepeatMasker and RepeatProteinMask (http://www.repeatmasker.org) were used to search the genome sequences for known repeat elements, with the genome sequences used as queries against the Repbase database 34 . The repetitive element annotations are listed in Table 4.
Gene annotation was performed on the repetitive-element-masked genome. A combined strategy of homology-based, ab initio and transcriptome-based gene prediction methods was used. Protein sequences of Astyanax mexicanus, Danio rerio, Gadus morhua, Ictalurus punctatus, Oryzias latipes, Takifugu rubripes, Tetraodon nigroviridis and Oreochromis niloticus were downloaded from Ensembl 35 . Proteins from the closely related fish species were mapped to the yellowbelly pufferfish genome using TBLASTN 36 . The alignments were joined with Solar, and GeneWise 37 was used to predict the exact gene structure of the corresponding genomic regions. Augustus 38 was also used for the ab initio prediction of genes in the repeat-masked genome. Finally, the full-length transcriptome sequences generated from PacBio sequencing were aligned to the genome using the Fig. 3 The 17-mer count distribution for the genome size estimation. Note that the peaks around the depths of 33, 66 and 132 represent the heterozygous, homozygous and repeated Kmers, respectively.  Table 2. Assembly statistics for the yellowbelly pufferfish. Note that the term contig here refers to the continuous sequences obtained after the Hi-C-data-based chromosome construction. Note that "new" represents the genome assembled in the present work and that "old" refers to the genome published in 2014. (2019) 6:267 | https://doi.org/10.1038/s41597-019-0279-z www.nature.com/scientificdata www.nature.com/scientificdata/ TopHat package 39 , and gene structure was predicted using Cufflinks 40 . All gene models were merged, and redundancy was removed by MAKER 41 , leading to a total of 29,416 protein-coding genes.

Content
The NCBI non-redundant protein (nr) database and the SwissProt database with an E-value threshold of 1e-5 were used for the functional annotation of the protein-coding genes using BLASTX and the BLASTN utility 42 . Functional ontology and pathway information from the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases was assigned to the genes using Blast2GO 43 . Ultimately, 28,017 genes (95.24% of the total) of the yellowbelly pufferfish were functionally annotated (Table 5).

Data records
The sequencing dataset and genome assembly were deposited in public repositories. The previous genome assembly can be accessed under accession number AOOT01000000 in NCBI 44 . The genomic Illumina sequencing data, genomic PacBio sequencing data, transcriptomic PacBio sequencing data, and the genomic Hi-C sequencing data were deposited as SRP162225 in NCBI 45 . The final chromosome assembly and genome annotation were submitted to NCBI Assembly with accession number RHFK00000000 46 . The functional annotation files are also available at figshare 47 .

Technical Validation
The quality of the DNA and RNA molecules and libraries used for genomic sequencing and transcriptome sequencing was validated before sequencing. The extracted DNA spectrophotometer ratios (SP) were 260/280 ≥1.6 for both Illumina and PacBio sequencing. DNA samples >2 μg and 20 μg were used for Illumina and PacBio sequencing, respectively. The concentration and quality of the total RNA were evaluated using a NanoVue Plus spectrophotometer (GE Healthcare, NJ, USA). RNA samples with a total RNA amount ≥10 μg, RNA integrity number ≥8, and rRNA ratio ≥1.5 were used to construct the sequencing library.
To validate our genome assembly, we compared the new genome to the previous genome. The new genome contained significantly fewer ambiguous bases (0.02 Mb) than the previous genome (67.9 Mb), but the size (366 Mb) of the new genome was approximately 20 Mb larger than that of the previous genome (346 Mb). Considering the estimated genome size of 377 Mb determined from the Kmer-based method, our new genome exhibited high completeness compared to the previous genome. The contig and scaffold N50 values of the newly assembled genome were almost 4,000 and 50 times higher than those of the previous genome, indicating a remarkable improvement in the sequence continuity of our assembly. We attributed the completeness and the continuity of the new genome to the application of PacBio long reads in the genome assembly. To further validate the improved continuity, we aligned genome fragments to our new genome with the NUCmer utility and found that more than 76% of the contigs were reliably mapped to the new genome with alignment ratios greater than 95%. Figure 4 provides two examples of alignments of genome sequences from the previous genome to our new assembly, showing that our new genome has significantly improved sequence continuity compared to the previous version.

Usage Notes
The contig sequences were assembled into chromosomes using interaction information from Hi-C sequencing data; therefore we used 100 bp to represent the unknown gap sizes among contigs in the chromosome sequences.  Table 5. The statistics of functional annotation of protein-coding genes. Note that "at least one database" here refers to genes with at least one hit in multiple databases.