The sequence and de novo assembly of Oxygymnocypris stewartii genome

Animal genomes in the Qinghai-Tibetan Plateau provide valuable resources for scientists to understand the molecular mechanism of environmental adaptation. Tibetan fish species play essential roles in the local ecology; however, the genomic information for native fishes was still insufficient. Oxygymnocypris stewartii, belonging to Oxygymnocypris genus, Schizothoracinae subfamily, is a native fish in the Tibetan plateau living within the elevation from roughly 3,000 m to 4,200 m. In this report, PacBio and Illumina sequencing platform were used to generate ~385.3 Gb genomic sequencing data. A genome of about 1,849.2 Mb was obtained with a contig N50 length of 257.1 kb. More than 44.5% of the genome were identified as repetitive elements, and 46,400 protein-coding genes were annotated in the genome. The assembled genome can be used as a reference for future population genetic studies of O. stewartii and will improve our understanding of high altitude adaptation of fishes in the Qinghai-Tibetan Plateau.

Animal genomes in the Qinghai-Tibetan Plateau provide valuable resources for scientists to understand the molecular mechanism of environmental adaptation. Tibetan fish species play essential roles in the local ecology; however, the genomic information for native fishes was still insufficient. Oxygymnocypris stewartii, belonging to Oxygymnocypris genus, Schizothoracinae subfamily, is a native fish in the Tibetan plateau living within the elevation from roughly 3,000 m to 4,200 m. In this report, PacBio and Illumina sequencing platform were used to generate~385.3 Gb genomic sequencing data. A genome of about 1,849.2 Mb was obtained with a contig N50 length of 257.1 kb. More than 44.5% of the genome were identified as repetitive elements, and 46,400 protein-coding genes were annotated in the genome. The assembled genome can be used as a reference for future population genetic studies of O. stewartii and will improve our understanding of high altitude adaptation of fishes in the Qinghai-Tibetan Plateau.

Background & Summary
The Qinghai-Tibetan Plateau (QTP) is the largest and highest plateau in the world 1 . The upshift of QTP has formed complex mountain systems in Southwest China and greatly reshaped the drainage at this area 2 . The rapid alteration of topography in the QTP might act as significant barriers for gene flow of many species, leading to population isolations and initiating allopatric divergence and speciation 3 . Genomes of fish species in the QTP provide valuable resources for scientists to understand the molecular mechanism of environmental adaptation. Although we have successfully obtained the reference genome of Glyptosternon maculatum 4 , leading to the first high-quality fish genome in Tibet-plateau, the genome information of fish species in QTP is still lacking. The schizothoracine fishes (Schizothoracinae subfamily, Cyprinidae family, Cypriniformes order), also known as "mountain carps", which composed of approximately 100 species in 10-13 genera 5 . They can be diagnosed by two lines of enlarged scales along both sides of the urogenital opening and anus 6 . These fishes exhibit many unique traits that adapt to the extreme environment of the QTP 7 . Therefore, this taxon provides an excellent opportunity for investigating high altitude adaptation of teleost fishes.
Distributed in the QTP and its surrounding areas, they are the largest and most diverse taxon of the QTP icthyofauna 6 . Based on morphological traits, the schizothoracine fishes can be divided into three hierarchical groups that adapt to different environments of QTP: the primitive group (including Schizothorax, Schizocypris, and Aspiorhynchus), the specialized group (including Diptychus, Gymnodiptychus, and Ptychobarbus), and the highly specialized group (including Gymnocypris, Oxygymnocypris, Chuanchia, Herzensteinia, Platypharodon, and Schizopygopsis) 6 . The evolution of the three groups was proposed to be associated with the upshift history of the plateau 6,8 . Thus, schizothoracine fishes represent an excellent model for the study of speciation caused by geographical isolation, as well as a good model for the study of adaptive evolutions of fish species in the QTP.
Another prominent feature in the evolution of schizothoracine fishes is the complex chromosome compositions, and the majority of fishes in this taxon are considered to be polyploids 9 . Whole genome duplication (WGD) plays a vital role in the evolutionary history of plant and animals. There are at least three rounds of whole genome duplications early in teleost diversification 10,11 , and these events were suggested to be causally related to the evolutionary success of teleost 12,13 . The polyploid nature and rapid diversification of schizothoracine fishes make them a good model for the study of polyploidy driven speciation.
Oxygymnocypris stewartii (Lloyd, 1908) (NCBI Taxon ID: 361644, Fig. 1a), a highly specialized schzothoracine fish, is a one-time spawning fish species mainly distributed in the tributaries of the middle reaches in the YarlungZangbo River across an elevation ranging from roughly 3,000 m to 4,200 m 14 (Fig. 1b). O. stewartii is currently listed in the Red List by the World Conservation Union (IUCN) and identified as an endangered fish 15 . Therefore, it is imperative to protect and restore the population resources of the O. stewartii.
In this report, we provide the whole genome sequence of O. stewartii through the PacBio single molecule sequencing technique (SMRT). The availability of a fully sequenced and annotated genome is essential to support basic biological studies and will be helpful to the development of further protection strategies for this endangered species. Its whole genome sequence will also provide a foundation to explore the adaptive evolutionary processes of highland fishes, supplied as a starting point to study speciation mechanisms caused by the rapid rising of the QTP.

Sample collection and sequencing
A healthy female fish captured from Gongga Country, Lhasa, Tibet (Fig. 1a,b) was used for genome sequencing. Genomic DNA was isolated using Qiagen DNA purification kit (Qiagen, Valencia, CA, USA) from the white muscular tissue as in our previous studies 4 .
To generate enough read data for the genome assembly, both the PacBio SEQUEL and the Illumina HiSeq 4000 platform were used for the sequencing. Long reads generated from the PacBio platform were used for genome assembly, and the short but accurate reads from the Illumina platform were analyzed for genome survey and base level correction after the assembly. For the PacBio platform, genomic sequencing libraries were constructed according to the PacBio suggested protocol and 141.1 Gb long sequencing reads were obtained from 27 SMRT cells. A total of 140.7 Gb (coverage of 74.3×) subreads were obtained after removing adaptors in polymerase reads ( Table 1). The subreads N50 and average lengths were 14.2 and 9,0 kb, respectively. For the Illumina HiSeq 4000 sequencing platform, one ug genomic DNA molecules were used for sequencing library construction. DNA molecules were fragmented, end-paired and ligated to the adaptor, which was further fractionated on agarose gels and purified by PCR amplification. To improve the representativeness of reads for the O. stewartii genome, 11 paired-end sequencing libraries were constructed with insert length of 250 bp according to Illumina's protocol (Illumina, San Diego, CA, USA). Finally, a total of 145.4 Gb (coverage of 70.8×) short sequencing reads were generated. Reads with the adaptors and a quality value lower than 20 (corresponding to a 1% error rate) were filtered out. As a result, we obtained 144.3 Gb cleaned reads for the k-mer analysis and base correction of the genome ( Table 1).
The individual used for the genomic sequencing was also used for the transcriptome sequencing, providing necessary gene expression data for the genome sequence annotation. Given that gene www.nature.com/sdata/ SCIENTIFIC DATA | 6:190009 | https://doi.org/10.1038/sdata.2019.9 expression exhibited clear tissue-specificity, 12 tissues, including skin, eye, swim bladder, muscle, brain, gill, heart, liver, gut, ovary, fat tissue and kidney were collected for the following transcriptome sequencing. As per the similar method in our previous study 4 , RNA molecules were extracted using RNAiso Pure RNA Isolation Kit (Takara, Japan) for all samples, and DNase I treatment was performed to eliminate DNA contamination. After the quality assessment of the extracted RNAs using NanoVue Plus spectrophotometer (GE Healthcare, NJ, USA), RNA-seq libraries were constructed according to the protocol 4 and were sequenced by Illumina HiSeq 4000 in paired-end 150 bp mode, resulting in a total of 50 Gb transcriptome data. All genome and transcriptome sequencing data were summarized in Table 1.
De novo assembly of Oxygymnocypris stewartii genome Genome size was estimated using Illunima sequencing data with the Kmer-based method 16 . As per our previous study 4 Table 3).
The completeness and the accuracy of the genome were evaluated by CEGMA, BUSCO and read mapping. The completeness of the genome assembly was assessed by the single copy orthologs (BUSCO, version 3.0) 21 and CEGMA 22 software. 94.2% complete and 3.6% partial of the 2,586 vertebrate BUSCO genes were identified in the final assembly. Using CEGMA 22 , we revealed that 95.56% of the 248 core genes were evolutionarily conserved genes identified in the genome. Both BUSCO and CEGMA confirmed the completeness of the genome assembly. The accuracy of the genome was evaluated by the Illumina short read mapping with BWA 23 and the transcript alignment with BLAT 24 . More than 98.6% of the reads were aligned to the genome, and the insert length distribution exhibited a single peak that was     consistent with the experimental design. Meanwhile, the transcriptome was de novo assembled by Trinity 25 , and the transcripts were mapped to the genome assembly using BLAT 24 with default parameters. We found that the alignment coverage (alignment length to transcript length) of expressed genes ranged from 96.44 to 99.95% in the genome assembly.      Repetitive element and non-coding gene annotation in the O. stewartii genome To annotate repeat elements in the O. stewartii genome, both homologous comparison and ab initio prediction were applied. The similar annotation process in our previous work 4 was employed. For ab initio repeat annotation, LTR_FINDER 26 , RepeatScout 27 , and RepeatModeler (http://repeatmasker.org/ RepeatModeler/) were used to construct a de novo repetitive element database, and the RepeatMasker 28 (http://repeatmasker.org/RMDownload.html) were used to annotate repeat elements with the database. Then, RepeatMasker and RepeatProteinMask 28 were used for known repeat element types by searching against Repbase database 29 . Tandem repeats were also ab initio predicted using TRF tool 30 . A total of 822.84 Mb repetitive elements were identified in the O. stewartii genome by those repeat annotation processes, accounting for 44.50% of the whole genome (Tables 4 and 5 and Fig. 3). For non-coding genes, 24,208 tRNAs were predicted using tRNAscan-SE 31 , and 1,363 rRNA genes were annotated using BLASTN tool with an E-value of 1E-10 32 against human rRNA sequence. Small   nuclear and nucleolar RNAs in the O. stewartii genome were also annotated by the infernal tool 33 using Rfam database 34 (Table 6).

Protein-coding gene prediction and functional annotation
The gene model prediction method in our previous study 4 was applied to the protein-coding gene annotation in the O. stewartii genome. We merged the evidence of the gene prediction from multiple methods, including homolog based, ab initio and RNA-seq based annotations. The protein and coding sequences were obtained from the Ensembl database 35 for the following species, including human (Homo sapiens, GCF_000001405.37), mouse (Mus musculus, GCF_000001635.26), zebrafish (Barchydanio rerio var, GCF_000002035.5), common carp (Cyprinus carpio, GCF_000951615.1), tiger puffer (Takifugu rubripes, GCF_000180615.1), channel catfish (Ictalurus punctatus, GCF_001660625.1), Sinocyclocheilus graham (GCF_001515645.1) and grass carp 36 (Ctenopharyngodon idellus). The protein sequences were aligned against the O. stewartii genome using TBLASTN 37 search with parameters of e-value 1e-5. After filtering low-quality records, the gene structure was predicted by GeneWise 38 (referred to "Homology" in Table 7). Secondly, transcripts assembled from twelve tissues RNA-Seq data were aligned against the O. stewartii genome using Program to Assemble Spliced Alignment (PASA) 39 (referred to "PASA" in Table 7). Augustus 40 , GeneID 41 , GeneScan 42 , GlimmerHMM 43 , and SNAP 44 were used for ab initio prediction with the optimized parameters that trained using high-quality proteins that derived from the PASA gene models. RNA-seq reads were also aligned to the O. stewartii genome directly using TopHat 45 v2.0.9, and the gene models were constructed by Cufflinks 46 v2.2.1 (referred to Cufflinks in Table 7). Finally, EvidenceModeler 39 was applied to combine all gene models that were predicted by various methods with the identical weights with our previous work 4 . Untranslated regions (UTRs) and alternative splicing variations were annotated using PASA2 39 (referred to "PASA-update" in Table 7). Finally, 46,400 protein-coding genes with a mean of 8.41 exons per gene (Table 7) were annotated in the O. stewartii genome. The statistics of gene models, including lengths of a gene, CDS, intron, and exon in O. stewartii were comparable to those for close-related species (Table 8 and Fig. 4).
Public biological function databases of SwissProt 47 , InterPro 48 , NR from NCBI and Kyoto Encyclopedia of Genes and Genomes (KEGG) 49 were used for the functional annotation of the predicted genes. BLASTX utility 32 were used for the homolog search with an E-value threshold of 1E-5. InterPro database 48 was used to predict protein function based on the conserved protein domains by InterproScan tool 50 . A total of 45,991 genes (99.1%) were successfully annotated by at least one public database. (Table 9 and Fig. 5).

Code Availability
The sequence data were generated using the software provided by the sequencing platform manufacturer and the sequencing data were processed with commands with the guidance from the public software that is cited in the manuscript. No custom computer codes were generated in this work.

Data Records
All PacBio long-read sequencing data and Illumina short-read sequencing data have been deposited to NCBI Sequence Read Archive (SRA) (Data Citation 1).
The transcriptome data are available through the NCBI SRA (Data Citation 2). The assembled genome version is available at GenBank (Data Citation 3). The annotation gff3 file of the assembled genome is available at Figshare (Data Citation 4).

Species Gene Number
Average  Technical Validation

RNA integrity
The transcriptomes for twelve tissues from three fish individuals were sequenced. Before constructing RNA-Seq libraries, the concentration and quality of total RNA were evaluated using NanoVue Plus spectrophotometer (GE Healthcare, NJ, USA). The total amount of RNA, RNA integrity and rRNA ratio were used to estimate the quality, content and degradation level of RNA samples. In the present study, RNAs samples with a total RNA amount ≥10 μg, RNA integrity number ≥8, and rRNA ratio ≥1.5 were finally subjected to construct the sequencing library.    Quality filtering of Illumina sequencing raw reads The raw sequencing reads generated from the Illumina platform were rigorously cleaned by the following procedures as in the previous study 4 . Firstly, adaptors in the reads were filtered out; secondly, reads with more than 10% of N bases were filtered out; thirdly, reads with more than 50% of the low-quality bases (phred quality score o= 5) were filtered out. If any end pair was classified as low quality, both pairs were discarded. The initially generated raw sequencing reads were also evaluated for quality distribution, GC content distribution, base composition, average quality score at each position and other metrics.