Transcriptomic Analysis of Potential “lncRNA–mRNA” Interactions in Liver of the Marine Teleost Cynoglossus semilaevis Fed Diets With Different DHA/EPA Ratios

Long non-coding RNAs (lncRNA) have emerged as important regulators of lipid metabolism and have been shown to play multifaceted roles in controlling transcriptional gene regulation, but very little relevant information has been available in fish, especially in non-model fish species. With a feeding trial on a typical marine teleost tongue sole C. semilaevis followed by transcriptomic analysis, the present study investigated the possible involvement of lncRNA in hepatic mRNA expression in response to different levels of dietary DHA and EPA, which are two most important fatty acids for marine fish. An 80-day feeding trial was conducted in a flow-through seawater system, and in this trial three experimental diets differing basically in DHA/EPA ratio, i.e., 0.61 (D/E-0.61), 1.46 (D/E-1.46), and 2.75 (D/E-2.75), were randomly assigned to 9 tanks of experimental fish. A total of 124.04 G high quality genome-wide clean data about coding and non-coding transcripts was obtained in the analysis of hepatic transcriptome. Compared to diet D/E-0.61, D/E-1.46 up-regulated expression of 178 lncRNAs and 2629 mRNAs, and down-regulated that of 47 lncRNAs and 3059 mRNAs, while D/E-2.75 resulted in much less change in gene expression. The co-expression and co-localization analysis of differentially expressed (DE) lncRNA and mRNA among dietary groups were then conducted. The co-expressed DE lncRNA and mRNA were primarily enriched in GO terms such as Metabolic process, Intracellular organelle, Catalytic activity, and Oxidoreductase activity, as well as in KEGG pathways such as Ribosome and Oxidative phosphorylation. Overlap of co-expression and co-localization analysis, i.e., lncRNA–mRNA matches “XR_523541.1–solute carrier family 16, member 5 (slc16a5)” and “LNC_000285–bromodomain adjacent to zinc finger domain 2A (baz2a),” were observed in all inter-group comparisons, indicating that they might crucially mediate the effects of dietary DHA and EPA on hepatic gene expression in tongue sole. In conclusion, this was the first time in marine teleost to investigate the possible lncRNA–mRNA interactions in response to dietary fatty acids. The results provided novel knowledge of lncRNAs in non-model marine teleost, and will serve as important resources for future studies that further investigate the roles of lncRNAs in lipid metabolism of marine teleost.


INTRODUCTION
Long non-coding RNAs (lncRNAs) are a class of RNA molecules with more than 200 bases that function as RNAs with little or no protein-coding capacity (Spizzo et al., 2012). LncRNAs represent a new frontier in molecular biology. More and more studies have demonstrated that lncRNAs play critical roles in various biological processes, including chromatin modification, regulation of transcription, influence of nuclear architecture and regulation of gene expression at post-transcriptional and post-translational levels, and also interact with DNA, RNA and proteins during these processes Kornfeld and Brüning, 2014;Gardini and Shiekhattar, 2015;Lopez-Pajares, 2016;Smekalova et al., 2016;Delás and Hannon, 2017;Long et al., 2017).
In fish, however, little information has been available either about the identification and characterization of lncRNAs or about the roles of lncRNAs in lipid metabolism. Limited studies on lncRNAs in fish have been restricted to model fish species, and most of the studies focused on developmental biology (Kaushik et al., 2013;Liu et al., 2013;Haque et al., 2014;Dhiman et al., 2015;Al-Tobasei et al., 2016;Wang et al., 2016Wang et al., , 2017Sarangdhar et al., 2017;Hu et al., 2018). Therefore, following our previous studies on lipid/fatty acid nutrition of marine fish (Zuo et al., 2012a,b;Xu et al., 2016Xu et al., , 2017, the present study was aimed at investigating the potential involvement of lncRNAs in effects of dietary fatty acids on a non-model marine teleost tongue sole Cynoglossus semilaevis, which is also an important aquaculture species. A feeding trial was conducted in this study, followed by hepatic transcriptome analysis. Docosahexaenoic acid (DHA, 22:6n-3) and eicosapentaenoic acid (EPA, are the most important essential fatty acids for marine fish. Previous studies have widely demonstrated the critical roles of DHA and EPA in a series of physiological functions of marine fish such as functions of visual and neural systems (Bell et al., 1995;Furuita and Takeuchi, 1998;Ishizaki et al., 2000Ishizaki et al., , 2001Benítez-Santana et al., 2007;Noffs et al., 2009), bone development (Gapasin and Duray, 2001;Roo et al., 2009), pigmentation (Villalta et al., 2008;Vizcaíno-Ochoa et al., 2010), stress resistance (Kanazawa, 1997;Liu et al., 2002), immune response (Zuo et al., 2012b;Xu et al., 2016), and reproduction (Wilson, 2009;Xu et al., 2017). However, more fundamental mechanisms involved are still elusive. How DHA and EPA regulate physiological processes at transcriptional and post-transcriptional levels needs to be elucidated. Moreover, compared to total DHA+EPA contents, the effects of DHA/EPA ratio were interesting as well but relatively less studied (Sargent et al., 1999;Kim et al., 2002;Lee et al., 2003;Wu et al., 2003;Skalli and Robin, 2004;Hamre and Harboe, 2008;Wilson, 2009;Lund and Steenfeldt, 2011;ØStbye et al., 2011;Tocher, 2015). Our previous studies with marine species such as large yellow croaker Larmichthys crocea (Zuo et al., 2012a), Japanese seabass Lateolabrax japonicus (Xu et al., 2016) and tongue sole (Xu et al., 2017) have shown the significantly different roles of DHA and EPA in some physiological processes such as non-specific immune response and gonadal steroidogenesis. In the present study, taking advantage of completely sequenced genome data on tongue sole C. semilaevis, a genome-wide hepatic transcriptome analysis was conducted with tongue sole following an 80-day feeding trial with diets containing different DHA/EPA ratios. Potential lncRNA-mRNA interactions were thereafter analyzed based on the coexpression and co-localization analysis of differentially expressed lncRNAs and mRNAs among dietary groups. The results will provide novel knowledge about lncRNAs in non-model marine fish, and will serve as important resources for future studies that further investigate the roles of lncRNAs in lipid metabolism of marine teleost.

Experimental Diets, Experimental Fish, and Feeding Procedure
The experimental diets, experimental fish and feeding procedure has been described in a previous study of ours (Xu et al., 2018). Briefly, different levels of EPA enriched oil (containing 11.2% DHA and 52.0% EPA; in the form of triglyceride; Xi'an Renbang Biological Science and Technology Co., Ltd., Xi'an, China) and DHA enriched oil (containing 69.5% DHA and 6.6% EPA; in the form of triglyceride; Xi'an Renbang Biological Science and Technology Co., Ltd., Xi'an, China) were supplemented to the basal diet to obtain different DHA/EPA ratios, 0.61, 1.46, and 2.75, and the corresponding diets were designated as D/E-0.61, D/E-1.46, and D/E-2.75, respectively (Tables 1,2).   The feeding trial was conducted in a flow-through seawater system. Juvenile tongue sole C. semilaevis with an average initial weight of 23.40 ± 0.45 g were used in the feeding trial. Each diet was randomly assigned to triplicate tanks of 30 fish each. Fish were hand-fed to apparent satiation twice daily. The feeding trial lasted for 80 days. At the end of the feeding trial, after anesthetized with eugenol, liver samples from 6 fish each tank were collected, frozen with liquid nitrogen and stored at −80 • C prior to analysis. The feeding trial was carried out in accordance with the recommendations of Guidelines on Management of Experimental Animals, Animal Care and Use Committee of Yellow Sea Fisheries Research Institute. All sampling protocols, as well as fish rearing practices, were reviewed and approved by the Animal Care and Use Committee of Yellow Sea Fisheries Research Institute.

RNA Isolation, cDNA Library Construction, and Illumina Sequencing
Total RNA in liver samples was isolated using RNAiso Plus (TaKaRa Biotechnology (Dalian) Co., Ltd., Dalian, China). RNA degradation and contamination were monitored on 1% agarose gels. No genomic DNA contamination was observed. RNA purity was checked using the Nano Photometer R spectrophotometer (IMPLEN, Westlake Village, CA, USA). RNA concentration was measured using Qubit R RNA Assay Kit in Qubit R 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).
A total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations and all samples had RIN values > 8. Six individual samples from the same experimental tank were pooled in equal amounts to obtain a pooling sample for this replicate tank. Nine pooling samples were then used to prepare 9 separate Illumina sequencing libraries (three biological replicates for each dietary group).
After ribosomal RNA was removed with Epicentre Ribo-Zero TM rRNA Removal Kit (Epicentre, Madison, WI, USA), and rRNA free residues were cleaned up by ethanol precipitation, sequencing libraries were generated using the rRNA-depleted RNA by NEBNext R Ultra TM Directional RNA Library Prep Kit for Illumina R (NEB, USA) following the manufacturer's instructions. Briefly, fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H − ). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3 ′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of 150 ∼200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA). Then 3 µl USER Enzyme (NEB, Ipswich, MA, USA) was added to size-selected and adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min of 95 • C treatment before the PCR process. PCR was then carried out with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified with AMPure XP system (Beckman Coulter, Beverly, USA), and the library quality was assessed on Agilent Bioanalyzer 2100 system (Agilent, Santa Clara, USA).
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the libraries were sequenced on an Illumina Hiseq2500 platform (Illumina, Inc., San Diego, CA, USA) and 150 bp paired-end reads were generated.
Raw data (raw reads) in fastq format were firstly processed through in-house perl scripts. Clean data (clean reads) were obtained by removing reads containing adapter or poly-N, as well as low quality reads from raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. All the downstream analysis were based on clean data with high quality. The paired-end clean reads were then aligned to the reference genome (RefSeq assembly accession: GCA_000523025.1; https://www.ncbi.nlm.nih.gov/ genome/?term=Cynoglossus%20semilaevis) using TopHat v2.0.9. The mapped reads of each sample were assembled by both Scripture (beta2) (Guttman et al., 2010) and Cufflinks (v2.1.1) (Trapnell et al., 2010) in a reference-based approach. Cuffdiff (v2.1.1) was also used to calculate FPKMs (Fragments Per Kilo-base of exon per Million fragments mapped) of both lncRNAs and coding genes in each sample. FPKMs were calculated based on the length of the fragments and reads count mapped to this fragment. Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group. Cuffdiff provides statistical routines for determining differential expression in digital transcript or gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted with the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). Transcripts or genes with an adjusted P-value (P-adj) < 0.05 were assigned as differentially expressed (DE) between experimental groups (three biological replicates for each group).
Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) or lncRNA target genes was implemented by the GOseq R package, and GO terms with corrected P < 0.05 were considered significantly enriched by DEGs. KOBAS software was used to test the statistical enrichment of differentially expressed genes or lncRNA target genes in KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways (http://www.genome.jp/kegg/).

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) Validation of Illumina Sequencing Data
To validate the Illumina sequencing data, 10 mRNA and ten lncRNA, which were selected from the most potential "lncRNA-mRNA" interactions based on the co-expression and co-localization analysis of DElncRNAs and DEmRNA among dietary groups, were tested for qRT-PCR analysis, using the same RNA samples for the transcriptome profiling. Specific primers were designed based on the data from GenBank (Table 3), and beta-2-microglobulin (β-2-m) was used as the reference gene according to our previous screening. The real-time PCR was carried out with SYBR Green Real-time PCR Master Mix (TaKaRa Biotechnology (Dalian) Co., Ltd., Dalian, China) in a quantitative thermal cycler (Mastercycler eprealplex, Eppendorf, German). The amplification efficiency for all primers, which was estimated by standard curves based on a 6-step 10-fold dilution series of target template, was within 95 ∼105%, and the coefficients of linear regression (R 2 ) were more than 0.99. The detailed program was similar with Xu et al. (2014). The mRNA expression levels were studied by qRT-PCR method: 2 − CT (Livak and Schmittgen, 2001). The qRT-PCR data were subjected to one-way analysis of variance (ANOVA) in SPSS 16.0 (SPSS Inc., Chicago, USA) for Windows. Differences between means were tested by Tukey's multiple range test. The level of significance was chosen at P < 0.05, and the results were presented as means of triplicate tanks ± standard errors.

Availability of Materials
Detailed information about methods and materials used in the current study are available from the corresponding author on reasonable request.

Transcriptome Sequencing and Assembly
In the transcriptomic analysis, three pooled liver RNA samples were prepared for each dietary group (D/E-0.61, D/E-1.46, and D/E-2.75). Nine cDNA libraries were then constructed to perform Illumina sequencing. A total of 278,269,246, 269,325,362, and 279,307,386 clean reads were obtained for groups D/E-0.61, D/E-1.46, and D/E-2.75 respectively, giving rise to total clean bases of 41.74, 40.40 and 41.90 G, respectively (Supplementary Table 1). The average Q20 and Q30 (the sequencing error rate at 1 and 0.1% respectively) of the experimental samples was 96.17 and 90.17% respectively, indicating the high accuracy of the sequencing processes. Raw reads were deposited at the National Center for Biotechnology Information (NCBI)'s Sequence Read Archive under Accession No. SRP127310 (D1E2, D3E2, and D3E1 in the archived data represents groups D/E-0.61, D/E-1.46, and D/E-2.75 respectively). The reads were mapped on the genome of C. semilaevis and the average mapping rates of groups D/E-0.61, D/E-1.46, and D/E-2.75 was 77.62, 76.09, and 72.43%, respectively (Supplementary Table 2). The classification

Filtration and Characterization of lncRNA
Large-scale analysis of the mammalian transcriptome have shown that the number and types of lncRNAs far exceed those of protein-coding mRNAs, while a small proportion of lncRNAs have been reported to have biological functions (Carninci et al., 2005;Birney et al., 2007;Carninci and Hayashizaki, 2007;Kapranov et al., 2007). In fish, however, very little information has been available about the identification and function prediction of lncRNAs. Only limited studies have been published in model fish such as zebrafish Danio rerio and rainbow trout Oncorhynchus mykiss (Kaushik et al., 2013;Liu et al., 2013;Haque et al., 2014;Dhiman et al., 2015;Al-Tobasei et al., 2016;Wang et al., 2016). Moreover, previous studies in both mammals and zebrafish have shown that lncRNAs are less conserved than protein-coding genes (Guttman et al., 2010;Cabili et al., 2011;Ulitsky et al., 2011;Derrien et al., 2012). Therefore, specific investigation is needed to identify and characterize the lncRNA profile in a certain fish species.
In the present study with tongue sole, based on the transcript assembly results, lncRNA was filtered following five steps: (1) number of exon (≥2); (2) length of transcript (>200bp); (3) comparison with annotated transcript (Cuffcompare software): transcripts overlapped with annotated coding exon were eliminated, and transcripts overlapped with annotated lncRNAs were regarded as annotated lncRNA; (4) expression level (Cuffquant software, FPKM≥0.5); (5) prediction of coding potential: transcripts containing no coding potential predicted by any mainstream software (CPC, CNCI, and PFAM) was recognized as novel lncRNAs (see Supplementary Figure 1 for the results), and transcripts with coding potential predicted by at least one software were designated as Transcripts of Uncertain Coding Potential (TUCP). The filtration results following these steps were presented in Figure 1.
In the present study, the experimental diets significantly affected the expression levels of both lncRNAs and mRNAs  ( Table 4). The most significant difference existed between groups D/E-0.61 and D/E-1.46. Compared to group D/E-0.61, D/E-1.46 significantly (adjusted P < 0.05) up-regulated expression of 178 lncRNA and 2629 mRNA, and down-regulated that of 47 lncRNA and 3059 mRNA. However, compared to group D/E-0.61, D/E-2.75 only significantly up-regulated expression of 78 lncRNA and 1925 mRNA, and down-regulated that of 42 lncRNA and 2015 mRNA. Generally, this result indicated that dietary DHA/EPA ratio regulated gene transcription in tongue sole in a dose-dependent manner. This was in accordance with the effects of dietary DHA/EPA ratio on growth performances and other physiological status such as lipid accumulation (Xu et al., 2018). As reported previously, compared to group D/E-0.61, D/E-1.46 but not D/E-2.75 significantly increased the growth rate and whole-body lipid content of juvenile tongue sole. In addition, despite the large number of differentially expressed   (Figure 5). This indicated that each DHA/EPA ratio exerted a different regulation of the liver transcriptome. Moreover, the amount of DElncRNA was in proportion to that of DEmRNA among the inter-group comparisons. This suggested that the transcriptional activity of coding and non-coding RNAs in tongue sole liver was coexpressed in response to dietary DHA/EPA ratio, which was to be evidenced by the co-expression analysis described below.
To confirm the DEG results from the transcriptomic assay, 10 mRNA and 10 lncRNA, which were selected from most potential "lncRNA-mRNA" interactions based on the co-expression and co-localization analysis of DElncRNA and DEmRNA described below, were tested for quantitative RT-PCR analysis. The results showed that transcription of 18 out of 20 selected genes was in good accordance with the transcriptomic results (Table 5, Figure 6). The validation results confirmed the high accuracy of the transcriptomic results.

Co-expression and Co-localization of lncRNA and mRNA
In order to estimate the potential lncRNA-mRNA interactions in response to dietary DHA/EPA ratio and consequently to predict the potential roles of the obtained lncRNAs in mRNA expression, co-expression and co-localization relationship of DElncRNA and DEmRNA among dietary groups were analyzed ( Table 6).  (Supplementary Figures 4,5), indicating the possible wide involvement of lncRNA in the regulation of mRNA expression. On the other hand, these results also indicated the effects of dietary DHA/EPA on a wide range of    physiological processes in fish liver, which, however, was not the focus of the present study.
As some lncRNAs could regulate the expression of coding genes near genomic regions, the co-localization of DElncRNA and DEmRNA among dietary groups were also analyzed. In total, 684, 278, and 159 co-localization relationships were observed in comparisons D/E-1.46 vs. D/E-0.61, D/E-2.75 vs. D/E-0.61, and D/E-1.46 vs. D/E-2.75, respectively. The co-localization relationships were primarily enriched in GO terms such as RNA metabolic process, Membrane-bounded organelle, Nucleus, RNA polymerase, and Nucleoside-triphosphatase regulator activity, as well as in KEGG pathways such as Ribosome and Glycerolipid metabolism (Supplementary Figures 6-8). Compared to the co-expression relationships, the co-localization relationships were more related to gene expression regulation, indicating the distance from target mRNA probably influenced the transcription-regulating effects of lncRNAs.
Considering that concurrent existence of co-expression and co-localization relationships of a certain lncRNA-mRNA match must increase the possibility of lncRNA-mRNA interaction, the overlap of co-expression and co-localization relationships of DElncRNA and DEmRNA was analyzed ( Table 7). A total of 17, 13, and 4 overlapped lncRNA-mRNA matches was obtained for D/E-1.46 vs. D/E-0.61, D/E-2.75 vs. D/E-0.61, and D/E-1.46 vs. D/E-2.75, respectively. Two overlapped matches, "XR_523541.1solute carrier family 16, member 5 (slc16a5)" and "LNC_000285bromodomain adjacent to zinc finger domain 2A (baz2a)" were observed in all the three inter-group comparisons, indicating that they might play crucial roles in the effects of dietary DHA and EPA on hepatic gene transcription in tongue sole.
members (Halestrap and Meredith, 2004), and mainly catalyzes the rapid transport of many monocarboxylates across the plasma membrane. While function of some members of this family such as SLC16A1 (also known as MCT1), SLC16A7 (MCT2), SLC16A8 (MCT3), and SLC16A3 (MCT4) have been elucidated (Bröer et al., 1997(Bröer et al., , 1998Lin et al., 1998;Grollman et al., 2000;Manning Fox et al., 2000), functions of SLC16A5 (MCT6) are still unknown. Considering the possible roles of SLC16A5 in transmembrane transport of lactate, pyruvate, and ketone bodies, the significant difference in hepatic transcription of slc16a5 among groups with different DHA/EPA ratios indicated that DHA and EPA may have different efficiency and priority in energy supply in tongue sole. The co-expression and colocalization relationship of XR_523541.1 and slc16a5 suggested that the lncRNA XR_523541.1 may be involved in the regulation of slc16a5 transcription by dietary DHA/EPA. Precise functions of both XR_523541.1 and slc16a5 need to be elucidated by future studies.
Baz2a is an essential component of nucleolar remodeling complex, a complex that mediates silencing of a fraction of rDNA by recruiting histone-modifying enzymes and DNA methyltransferases, leading to heterochromatin formation and transcriptional silencing (Gu et al., 2015). The co-expression and co-localization relationship of "LNC_000285-baz2a" indicated that the lncRNA LNC_000285 probably mediated the regulatory effects of dietary DHA/EPA on gene expression through promoting heterochromatin formation. Inducing heterochromatin formation and suppressing mobile elements were major mechanisms of transcription-regulating effects of non-coding RNAs (Cusanelli and Chartrand, 2015;Meller et al., 2015;Acharya et al., 2017;da Rocha and Heard, 2017;Oliva-Rico and Herrera, 2017). In human pancreatic ductal adenocarcinoma, the lncRNA HOTAIR contributed in the silence of MicroRNA-34a by enchancer of zeste homolog 2 through induction of heterochromatin formation . Microsatellite repeat DXZ4-associated long non-coding RNAs also had developmental changes in expression coincident with heterochromatin formation at the human (Homo sapiens) macrosatellite repeat (Figueroa et al., 2015). Other evidence showed that lncRNA maturation to initiate heterochromatin formation in the nucleolus is required in exit from pluripotency in embryonic stem cells (Savić et al., 2014).
The overlapped match "XR_522182.1-lysine (K)-specific methyltransferase 2D (kmt2d)" was observed in inter-group Minus values in Pearson correlation denote negative correlation, i.e., "mRNA up-regulated but lncRNA down-regulated" or "mRNA down-regulated but lncRNA up-regulated." Frontiers in Physiology | www.frontiersin.org FIGURE 8 | Statistics of KEGG pathway enrichment for differentially expressed mRNA in lncRNA-mRNA co-expression analysis for D/E-1.46 vs. D/E-0.61. Rich factor is the ratio of number of differentially expressed genes in a certain pathway to number of all annotated genes in this pathway. qvalue is corrected P-value by multiple hypothesis test. q-value < 0.05 denotes significant differences.
comparisons D/E-2.75vs. D/E-0.61 and D/E-1.46vs. D/E-0.61. As a histone methyltransferase, Kmt2d methylates 'Lys-4' of histone H3 into H3K4me, which represents a specific tag for epigenetic transcriptional activation (Cho et al., 2007;Lan et al., 2007). This result indicated the potential roles of lncRNA in epigenetic transcriptional regulation. Other mechanisms of transcription-regulating effects of lncRNAs were also indicated from the present results. For example, several overlapped "lncRNA-zinc finger protein" matches, such as "XR_521587.1-zinc finger protein 319-like (zfp319), " "LNC_000360-zinc finger protein 574 (zfp574), " and "LNC_000360-zinc finger protein 585A-like (zfp585a), " were observed in the co-expression and co-localization analysis of DElncRNAs and DEmRNAs. Interacting with these zinc finger proteins may be a possible way of lncRNA functioning on gene expression.
Interactions of other lncRNAs with some mRNAs observed in overlapped "lncRNA-mRNA" matches of the current analysis have been reported in human and terrestrial animal studies. For example, the lncRNA-mRNA match "XR_523068.1-apolipoprotein A-IV (apoa4)" was overlapped in the co-expression and co-localization analysis in the present study. A study on mice has shown that an anti-sense lncRNA (APOA4-AS) concordantly and specifically regulated apoa4 expression both in vitro and in vivo with the involvement of the mRNA stabilizing protein HuR (Qin et al., 2016). The alignment analysis did not find homology between mice APOA4-AS and tongue sole XR_523068.1. Moreover, mice APOA4-AS was transcribed from the opposite strand of apoa4 locus (+ strand), and partially overlapped with apoa4 exon 3 and 3'-UTR, while tongue sole XR_523068.1 was far away from tongue sole apoa4 (67904 bp downstream). In mice, gene expression of another apolipoprotein APOC2 was also reported to be regulated by the lncRNA lncLSTR . These results indicated that lncRNAs may play certain roles in gene expression of apolipoproteins, which are key components of lipoproteins and play important roles in lipid transportation. In addition, the mRNA expression of apoa4 in response to dietary DHA/EPA ratio was negatively correlated with the hepatic lipid content of the experimental fish (Xu et al., 2018). The lower hepatic apoa4 expression in group D/E-1.46 may reduce the transportation of lipid out of liver, and thus contribute to higher hepatic lipid accumulation in this group.
Interactions between lncRNAs and suppressor of cytokine signaling 3-like (scs3) or methylenetetrahydrofolate dehydrogenase (NADP+ dependent) (mthfd) have also been observed in human studies. A recent study on postmenopausal osteoporosis observed a probable interaction between lncRNA LINC00963 and scs3 based on analysis of DElncRNA-DEmRNA co-expression network (Fei et al., 2018). However, different from the present study which showed that lncRNA LNC_000230 was 191 bp upstream relative to co-expressed scs3, in that study LINC00963 and scs3 were located in different chromosomes. Regarding mthfd, a study with MCF-7 breast cancer cells showed that a lncRNA, taurine-upregulated gene 1 (tug1), positively regulated the gene expression of mthfd2 (Zhao and Ren, 2016). Tug1 and mthfd2 were also located in different chromosomes, while in the present study, lncRNA XR_521789.1 was 54 bp downstream relative to co-expressed mthfd1.
In conclusion, this was the first time in marine teleost to investigate the possible lncRNA-mRNA interactions in response to dietary fatty acids. The results provided novel knowledge of lncRNAs in non-model marine teleost, and will serve as important resources for future studies that further investigate the roles of lncRNAs in lipid metabolism of marine teleost. Potential lncRNA-mRNA interactions "XR_523541.1-slc16a5" and "LNC_000285-baz2a" may play important roles in effects of dietary DHA/EPA on hepatic gene transcription in tongue sole.

DATA AVAILABILITY
The raw reads obtained from the transcriptomic analysis have been deposited at the National Center for Biotechnology Information (NCBI)'s Sequence Read Archive under Accession No. SRP127310. Other raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.