1 INTRODUCTION

Half-smooth tongue sole (Cynoglossus semilaevis) is a commercial flatfish, ranking among the most valuable fishes in the coastal areas of China. Similar to other flatfish species, this animal possesses an asymmetric body shape with lateralization of the eyes to the same side during metamorphosis (Chen et al., 2014).

In contrast to other flatfish species, half-smooth tongue sole features some unique natural characteristics. In nature, the length and weight of mature females are 2–4 times of mature males, and the mature ovary reaches more than one thousand times the volume and weight of the testis (Sun et al., 2010; Li et al., 2018). This species employs a female heterogametic sex determination system (ZW female and ZZ male) (Chen et al., 2009, 2014; Shao et al., 2014). However, the molecular mechanism that regulates these special characteristics remains unknown.

Long noncoding RNAs (lncRNAs) are defined as noncoding RNAs longer than ~200 nt and lack protein-encoding capacity (Philippen et al., 2015; Santosh et al., 2015). To date, more than 10 000 lncRNAs have been identified in humans, and their number continually grows. However, information on fishes remains lacking. LncRNAs are classified into five subclasses based on their location and function, namely intergenic, intronic, sense overlapping, antisense, and bidirectional lncRNAs. In contrast to microRNAs, studies of lncRNAs are still in their infancy due to extremely complex and complicated gene regulation mechanisms (Jiang and Ning, 2015; Philippen et al., 2015; Santosh et al., 2015; Lorenzen and Thum, 2016). However, increasing evidence has shown that lncRNAs perform strong biological functions and may be involved in all biological and developmental processes (Sun and Kraus, 2015).

In this study, the lncRNA profile in adult halfsmooth tongue soles was investigated using highthroughput sequencing method. We also analyzed the expression patterns of some lncRNAs in different tissues by quantitative reverse transcriptionpolymerase chain reaction (qRT-PCR).

2 MATERIAL AND METHOD 2.1 Fish sampling and RNA isolation

Four mature female and four mature male halfsmooth tongue soles were obtained from a nearby hatchery (Laizhou City, Shandong Province, China). The fishes were anesthetized with MS-222 before being sampled. The gonad, liver, heart, kidney, gill, muscle, and brain tissues were collected individually and frozen at -80℃ (Table 1). In all experiments, total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). All tissues from mature female and male fishes (n=4) were pooled together to generate total RNA for high-throughput sequencing experiments. All animal experiments were approved by the Institutional Animal Care and Use Committee of Binzhou Medical University.

2.2 Library construction and high-throughput sequencing

The transcriptome of lncRNAs and mRNAs was determined via next-generation deep sequencing using Illumina platforms. Normal RNA-seq library was constructed. According to the manufacturer's transcriptome sequencing protocol, first strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase. DNA polymerase I and RNase H were used to synthesize the second strand cDNA. After adenylation of the 3ʹ ends of the DNA fragments, an adaptor with a hairpin loop structure was ligated to prepare for hybridization. To select cDNA fragments with the preferred length of 150–200 bp, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, CA, USA). Finally, quality was assessed on an Agilent Bioanalyzer 2100 system. The library was sequenced at the Annoroad Genome (Beijing, China) on an Illumina Hiseq 2000 platform.

2.3 Quality control and sequencing data analysis

The sequencing data discussed in this publication have been deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession number GSE98425. To assess the quality of RNA-seq data, each base in the reads was assigned a quality score (Q) by Fast QC with a Phred-like algorithm. Clean data were obtained by removing reads containing adapters, reads containing over 10% poly-N, rRNA mapping reads and low-quality reads (> 50% of bases whose Q scores were < 20%) from the raw data. The clean data were mapped to the C. semilaevis reference genome (https://www.ncbi.nlm.nih.gov/genome/?term=Cynoglossus+semilaevis, GCA_000523025.1 Cse_v1.0) by TopHat2 (v2.0.13) with "-library-type fr-firststrand" as a parameter (Trapnell et al., 2012). The transcripts were assembled using Cufflinks (v2.2.0) according to the instructions provided (Trapnell et al., 2012).

2.4 LncRNA prediction

To identify lncRNA genes from the assembly results, the following criteria were considered. Firstly, transcripts shorter than 200 nt and single exons of transcripts were removed. Secondly, transcripts with an open reading frame that was longer than 300 nt were removed. Next, the transcripts with low expression levels (fragments per kilobase of transcript per million fragments mapped, FPKM < 1) and minimal read coverage threshold below 3 were discarded. Transcripts that overlapped with known protein-coding genes on the same strand were discarded. Then, Coding-Potential Assessment Tool (CPAT) (Wang et al., 2013), CodingNon-Coding Index (CNCI) (Sun et al., 2013) and Coding Potential Calculator (CPC) (Kong et al., 2007), were used to distinguish mRNA from lncRNAs. CPAT was used to assess the long noncoding transcripts by a logistic regression model trained by CPAT software based on the coding and noncoding transcripts of C. semilaevis from NCBI. The coding probability score from CPAT of C. semilaevis was optimized according to the instructions provided. Finally, the optimum cutoff for C. semilaevis gene annotation was determined as 0.356. Scores < 0.356 were retained as putative noncoding RNAs. CNCI software was downloaded from the web (https://github.com/www-bioinfo-org/CNCI) and used with the default parameters to predict noncoding RNAs. CPC (http://cpc.cbi.pku.edu.cn/) was used to predict the protein-coding potential for transcripts. Only transcripts with a CPC score ≤-1 were retained. The transcripts that remained were considered putative lncRNAs. The transcripts were submitted to identify known and novel lncRNAs by comparison against the NCBI-NR database.

2.5 qRT-qPCR of selected lncRNAs

Using specific primers designed for selected novel lncRNAs isolated from the library (Table 2), RT-PCR was performed to verify differential expression in seven tissues (gill, heart, brain, liver, muscles, kidneys and ovary) of female fish. The specific primers were designed by Primer Premier 5.0 software based on the sequencing results. Total RNAs from different tissues from female fishes (n=6) were reverse-transcribed to cDNA using SuperScript Ⅳ Reverse Transcriptase (ThermoFisher, USA). LncRNAs was performed on cDNAs using TaqMan^TM Fast Advanced Master Mix (ThermoFisher, USA). β-actin was used for internal control. The relative gene expression of selected lncRNAs was calculated by comparing the threshold cycle (C_t).

3 RESULT 3.1 Sequencing results and quality control

A total of 109 680 782 raw reads were produced from the cDNA libraries. After quality control, 96 726 222 clean reads were obtained. The proportion of reads with a Phred quality value higher than 30 among the clean reads totaled 96.03% (Table 3). Overall, 88.00% of the clean reads were aligned against the Cynoglossus semilaevis reference genome. Among these mapped reads, 80.7% of the reads were mapped to exon regions, 7.4% to intron regions, and 11.9% to intergenic regions (Fig. 1). According to their genome positions, these reads were widely distributed in all chromosomes of the half-smooth tongue sole (Fig. 2).

3.2 Identification of lncRNAs in half-smooth tongue soles

Transcripts were reconstructed using Cufflinks. To obtain a unique set of isoforms, all the reconstructed transcripts were assembled using Cuffcompare against the C. semilaevis reference genome, Cse_ v1.0. Potential lncRNAs were identified based on their sequences, amino acid peptide lengths, and protein-coding potential using CPAT, CNCI, and CPC methods. In this report, 1 694 known lncRNAs were identified from the sequencing library by BLAST in the NCBI-NR database. Table 4 lists the top 20 most abundant known lncRNAs. In addition to known lncRNAs, novel lncRNAs were predicted from the total assembled transcripts.

After removing low-abundance transcripts, known transcripts, CPAT-removed transcripts, CNCIremoved transcripts and CPC-removed transcripts, 1 412 potential lncRNA transcripts were obtained, and the results are shown in Table 5. Finally, 803 novel lncRNAs were identified; these novel lncRNAs included 287 novel intergenic, 1 novel intron, 488 novel senses, and 27 novel antisense types (Table 6). We analyzed the sequence length and exon number of novel lncRNAs. Results show that sequence length was mainly distributed in the range of 200 nt to 3 000 nt, and most lncRNAs contained 2–6 exons (Fig. 3). A total of 22 104 protein-coding transcripts were identified, and analysis of sequence length revealed that lncRNAs were shorter in sequence length and featured fewer exons than mRNAs (Fig. 4).

3.3 Validation of differentially expressed lncRNAs in tissues of female fish

To validate our deep sequencing results, six novel lncRNAs identified by deep sequencing analysis were selected. Their expression levels in different tissues of female fish were further determined by qRT-PCR. Table 2 lists the primers used in the study. Expression patterns of these lncRNAs showed dissimilarity in the tissues of female fish (Fig. 5). Lnc_673 and lnc_40 showed broad expressions and were both expressed in all tissue samples. Lnc_190, lnc_493, lnc_770 and lnc_150 showed characteristic tissue-specific expressions. Lnc_190 was strongly expressed in the brain and comparatively weaker in the liver of females. Lnc_493 was mainly expressed in brain and ovary compared with other female tissues. Lnc_770 and lnc_150 were mainly expressed in the ovary and slightly expressed in the liver and muscles of the female fish. These results imply that the two novel lncRNAs may be involved in the development of ovary in the fish.

4 DISCUSSION

In the past decades, significant genetic improvements have been achieved in half-smooth tongue soles. However, all aspects were focused on mRNAs and microRNA research (Liu et al., 2016a, b; Xu et al., 2016; Yan et al., 2016). Information on lncRNAs in half-smooth tongue soles remains lacking. The mechanisms and functions of lncRNAs have also remained unclear in marine fish. LncRNAs are noncoding RNAs longer than 200 nt; significant progress has been achieved in the study of these molecules in mice, humans and other model organisms in the last few years. Studies have shown that lncRNAs regulate metabolic tissue development and function, including muscle differentiation, microvascular dysfunction, colorectal cancer progression and vascular development (Gong et al., 2015; Viereck et al., 2015; Yan et al., 2015; Ma et al., 2016). Although numerous studies have indicated the importance of lncRNAs in gene regulation, little is known about their biological function in marine fishes. Our study is the first to screen for lncRNAs in half-smooth tongue sole (C. semilaevis) using largescale deep sequencing. After quality control, 96 726 222 clean reads were obtained, and 88.00% of the clean reads were successfully mapped to the C. semilaevis reference genome. As the first study of lncRNAs in C. semilaevis, we identified 803 novel lncRNAs. Six novel lncRNAs were selected for expression analysis by qRT-PCR. Some lncRNAs showed tissue-specific expression, indicating that lncRNAs may be involved in tissue development.

According to expression profile results by qRTPCR, the novel lncRNA lnc_190 was strongly expressed in the brain, lnc_493 was mainly expressed in the brain and ovary, whereas lnc_770 and lnc_150 were mainly expressed in the ovary. Numerous studies showed that lncRNA genes are expressed in tissuespecific patterns (Roberts et al., 2014). In studies of human cell lines, results show that 29% of lncRNAs were expressed specifically in a single-cell type, and only 10% were expressed in all cell types (Roberts et al., 2014). For example, in mice, the linc-Brn1b expression is primarily restricted to specific brain regions and predominantly restricted in the developing cerebral cortex. linc-Brn1b shows spatiotemporally regulated patterns of expression during cortical development (Sauvageau et al., 2013). These results support the attribution of potential function based on tissue-specific expression. The expression of lncRNAs is regulated both at the transcriptional and posttranscriptional levels. However, how tissue-specific expression is achieved is remains unclear, especially in fish. Transcription factors and miRNAs targeting each lncRNA might have some contribution to tissuespecific expression. LncRNA transcription could affect the expression of adjacent genes in cis, hybridize to the overlapping sense transcript, act as ceRNA, or function in more complex ways to regulate gene expression. In cotton (Gossypium hirsutum L.), expression patterns analysis showed that most intronic lncRNAs (75.95%) exhibit a consistent expression pattern with their adjacent protein-coding genes (Lu et al., 2016). C2 calcium-dependent domaincontaining protein 4C-like (C2cd4c) is an adjacent gene to lnc_190 on chromosome 6 in C. semilaevis. C2cd4c is overexpressed in the whole brain in other species (Yue et al., 2014). This implies lnc_190 and C2cd4c might have a certain regulation relationship. The hypothalamic-pituitary-gonadal axis plays an important role in reproduction regulation in fish (Shao et al., 2014). However, the regulatory mechanism remains unclear, especially in C. semilaevis. These novel lncRNAs showed tissue-specific expression patterns in the brain or ovary. Thus, lncRNAs may be involved in reproduction regulation in C. semilaevis. The detailed function and mechanisms require further studies.

Compared with those in mouse and humans, the function of lncRNAs in marine fish poses considerable research difficulties. However, with some newly emerging technologies, the function of lncRNAs in marine fish will be easily studied. Recently the technology of genome editing has been successfully used in C. semilaevis to study gene functions (Cui et al., 2017).

5 CONCLUSION

Non-coding RNAs comprise the majority of the transcriptome, whose important functions have been discovered in model organisms. To date, few lncRNAs have been characterized in marine fish, including half-smooth tongue sole. In this study, hundreds of lncRNAs were identified by highthroughput sequencing and bioinformatics-analyses. A total of 1 694 known lncRNAs and 803 novel lncRNAs were identified. Some lncRNAs showed differential tissue distribution. These results suggest that lncRNAs might be important regulatory factors in the development of the fish.

6 DATA AVAILABILITY STATEMENT

The raw reads of the present study were also uploaded to the Sequence Read Archive databases of NCBI with accession number SRR5494302 (NCBI BioProject accession number, PRJNA385074). This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession number GFXK00000000. The version described in this paper is the first version, GFXK01000000.

7 COMPETING INTEREST

The authors declare there are no competing interests.