Embryonic transcriptome sequencing of the ocellate spot skate Okamejei kenojei

Chondrichthyans (cartilaginous fishes) exhibit highly variable reproductive styles, categorized as viviparity and oviparity. Among these, species with oviparity provide an enormous potential of molecular experimentation with stable sample supply which does not demand the sacrifices of live mothers. Cartilaginous fishes are divided into two subclasses, chimaeras (Holocephali) and elasmobranchs (Elasmobranchii), and the latter consists of two monophyletic groups, Batoidea (rays, skates and torpedoes) and Selachimorpha (sharks). Here we report transcriptome assemblies of the ocellate spot skate Okamejei kenojei, produced by strand-specific RNA-seq of its embryonic tissues. We obtained a total of 325 million illumina short reads from libraries prepared using four different tissue domains and assembled them all together. Our assembly result confirmed the species authenticity and high continuity of contig sequences. Also, assessment of its coverage of pre-selected one-to-one orthologs supported high diversity of transcripts in the assemblies. Our products are expected to provide a basis of comparative molecular studies encompassing other chondrichthyan species with emerging genomic and transcriptomic sequence information.


Background & Summary
Oviparous (egg-laying) chondrichthyans are distributed into three different selachimorphan (shark) orders, Carcharhiniformes, Heterodontiformes, and Orectolobiformes, and one batoid (skate) order Rajiformes (summarized in 1 ). The order Rajiformes consists of four families, Anacanthobatidae, Arhynchobatidae, Gurgesiellidae, and Rajidae. As of February 2018, molecular sequence data for rajiforms concentrate on species in Rajidae, such as the little skate Leucoraja erinacea, whose genome is being sequenced 2 , in the head of the list, and it is followed by Amblyraja radiata and Raja clavata. One of the genera in Rajidae is Okamejei, and among more than a dozen of species in the genus Okamejei, we focused on the ocellate spot skate Okamejei kenojei in the present study. Previously, this species was recognized as Raja porosa, Okamejei porosa, or Raja kenojei. The NCBI RefSeq entry for the whole mitochondrial genome sequence of this species Okamejei kenojei (NC_007173.1) was originally registered in GenBank under one of the former species names, Raja porosa (AY525783.1) 3 . The habitat of O. kenojei is the coasts of Japan, Korea, China, and Taiwan in the Northwest Pacific 4 . While the skate species with the most abundant sequence information, L. erinacea, inhabits a small area on the east coast of North America, O. kenojei is one of the most promising oviparous skate species available in East Asia for experimentation.
For modern life science studies, genomic and transcriptomic information serves as an indispensable fundamental resource. Especially for a species without abundant molecular sequences, even small-scale transcriptome sequencing with short-read data acquisition can provide a valuable start point, which can be achieved by simple tissue sampling and small financial investment. For efficient transcriptome sequencing data acquisition, a number of technical factors have been considered to optimize sample preparation, sequencing run design, and post-sequencing data analysis. Previously, the authors' group proposed some of those factors 5 , and in the present study, we extend those factors by incorporating latest reagents for sample preparation and tools for sequence data analysis.

Animal sampling and RNA extraction
All animal experiments were conducted in accordance with the Guideline of the Institutional Animal Care and Use Committee (IACUC) of RIKEN Kobe Branch (Approval ID: H16-11). An approximately 8 cm-long fertilized egg of the ocellate spot skate, O. kenojei (Fig. 1a), was purchased from a commercial marine organism supplier in Minami-ise town in Mie Prefecture, Japan, in February 2016. The embryo contained in the eggcase was 5 cm-long, corresponding morphologically to the stage 31 of typical shark embryos 6,7 and was dissected into six pieces labelled as head, gill, trunk, pectoral fins, cloaca, and caudal ( Fig. 1b). Total RNAs were extracted from the four parts of the embryo (head, trunk, cloaca, and caudal) using TRIzol reagent (Life Technologies) (Fig. 1c).

Library preparation and sequencing
Using 1 μg of each of the extracted total RNAs, four strand-specific RNA-seq libraries were prepared with KAPA Stranded mRNA-Seq Kit (Kapa Biosystems, cat. No. KK8420) according to its standard protocol unless stated otherwise below. Before the total volume PCR amplification was performed, we performed a preliminary PCR using a 1.5 μl aliquot of 10 μl DNA from the previous step, with KAPA Real-Time Library Amplification Kit (Kapa Biosystems, cat. No. KK2702). This demonstrated that the amplification of the products reached Standard 1 accompanying this kit between three and four PCR cycles, which instructed us to perform the full-volume PCR with three PCR cycles, introducing the minimal amplification (Fig. 1d)

Read trimming and assembly
The obtained sequence reads in the fastq files (Data Citation 1) were processed with the program Trim Galore! v0.3.3 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) with the options '--phred33 --stringency 2 --quality 30 --length 25 --paired'. The reads after adaptor trimming were assembled with the program Trinity v2.5.1 8 with the options '--trimmomatic --min_kmer_cov 2 --SS_lib_type RF'. Among the resultant contig sequences, those matching PhiX, mitochondrial DNA, and genome and transcriptome assemblies of the species sequenced in the same HiSeq run in the BLASTN v2.2.30+ 9 results executed with the options '-perc_identity 95' were removed. We designated the resultant sequence set as the redundant nuclear transcriptome assembly (Data Citation 2). In this assembly, protein-coding regions were predicted with the program Transdecoder v5.0.2 10 following its user documentation, which employed the protein-level similarity to SWISSPROT 11  (Data Citation 4). Annotation of the putative peptide sequences were performed with similarity searches using BLASTP v2.2.30+ towards RefSeq Protein sequences for the human (GRCh38, 113,620 sequencesas of July 12, 2018) and the Callorhinchus milii (NCBI RefSeq, 28,600 sequences -as of July 12, 2018). The whole post-sequencing procedure is outlined in Fig. 2.

Code availability
No custom computer code was employed in this study.

Data Records
The approximately 325 million raw sequence reads from four different portions of a O. kenojei embryo were released in the NCBI Sequence Read Archive (Table 1 and Data Citation 1). The nuclear transcriptome assembly using all the obtained reads consisted of 1,081,614 sequences (Data Citation 2), which could include spurious intergenic transcripts and alternative splicing variants. Mapping of the reads employed in the assembly to the assembled transcript contig sequences (Data Citation 2) yielded the mapping rates of as high as 90.2 to 93.1%. The nuclear protein-coding transcriptome assembly consisted of 167,783 sequences available both in the nucleotide and amino acid sequences (Data Citation 3). Out of those putative protein-coding sequences, 88,376 (52.7%) were predicted to have complete ORFs    AB195842.1, AB201248.1, and AB201247.1), and for each of those sequences, a transcript contig in our assembly showed no less than 98% similarity (no more than 1% gaps) ( Table 2). These high similarities authentificate that the species used for this study was O. kenojei.

Transcript diversity measured by one-to-one ortholog coverage
We also employed a completeness assessment on the web server gVolante 15 in which the ortholog search pipeline BUSCO 16 is implemented. This method evaluates the coverage of one-to-one protein-coding orthologs selected in advance. For this purpose, we used the ortholog set CVG 5 introduced for more accurate assessment for vertebrate sequence sets than performed with other ortholog sets, as well as the Vertebrata ortholog set released together with BUSCO 16 . Although the assessment results showed a slight fluctuation depending on the type of input data (nucleotide or amino acid sequences), out of the 233 components of CVG, our resultant data were shown to always contain at least 218 full-length ortholog sequences ('complete') and 232 partial sequences ('fragmented') ( Table 3). Our assessment with the BUSCO's Vertebrata ortholog set also produced comparable scores, mostly >92% in percentages. These high scores ascertain the high coverage of protein-coding genes in our resultant transcript sequence data set.