1 INTRODUCTION

Large numbers of non-coding RNAs (ncRNAs) have been reported in mammalian and other eukaryotes genomes. Majority of these ncRNAs are by-products of RNA processing, performing various regulatory functions in cells (Mercer et al., 2009). NcRNAs were initially thought to be spurious transcriptional noise resulting from low RNA polymerase fidelity because their low protein-coding capacity compared with other RNAs (Struhl, 2007). However, previous research showed ncRNAs could take part in a variety of biological processes (McHugh et al., 2015). As for ncRNA types, they can be divided into sncRNA (length between 20–50 nt), mncRNA (length between 50–200 nt) and long non-coding RNAs (lncRNAs) (length over 200 nt) by different lengths. Numerous functional studies on ncRNAs have been conducted in recent years especially for lncRNA. For instance, tens of thousands lncRNAs have been identified in species including human (Homo sapiens), mouse (Mus musculus), zebrafish (Danio rerio), and Caenorhabditis elegans (Zhao et al., 2016). Depending on their relative positions to coding genes, lncRNAs can be classified into antisense lncRNAs, intronic lncRNAs, divergent lncRNAs, and intergenic lncRNAs (Ransohoff et al., 2018). lncRNAs play important roles in higher-order chromosome dynamics, chromatin modification, telomere biology, subcellular structural organization, transcriptional regulation, and post-transcriptional regulation. Mercer et al. (2009) found the down-regulation of a lncRNA called TROJAN impaired the proliferative potential of breast cancer cell lines in vitro, while over-expression of TROJAN promoted this ability (Jin et al., 2019). In addition, many lncRNA molecules related to human immune response were explored in early studies, such as Morrbid (Kotzin et al., 2016), lnc-Lsm3b (Jiang et al., 2018), NKILA (Huang et al., 2018). Moreover, innate immune responses of lncRNAs were widely investigated in other species. For example, mammals like sheep (Zhang et al., 2021a), arthropods like Drosophila (Zhang et al., 2020), several mollusks (Zhang et al., 2021b), even for plants (Gómez and Pallás, 2013). Therefore, the function of lncRNAs in the process of animal immune regulation is worth further exploration.

In teleost, the potential roles of lncRNAs have also been increasingly explored in recent years. For example, in common carp (Cyprinus carpio) undifferentiated gonads, 124 lncRNAs showed differential expression between experimental and control groups, which associated with gonadal development (Song et al., 2019). In the process of muscle growth, rainbow trout (Oncorhynchus mykiss) harbored 1 160 differentially expressed lncRNAs among different groups in the back muscle, involved in regulation of muscle growth (Paneru et al., 2018). Moreover, the lncRNAs in teleost involved in immune response against pathogenic infections was widely reported in previous studies. For example, 556 lncRNA showed differential expression when the rainbow trout infected by Flavobacterium psychrophilus. Several of them exhibited associations with immune related genes and pathways (Paneru et al., 2016). Furthermore, such immune related lncRNAs have also been identified in zebrafish (D. rerio) (Chen et al., 2018), olive flounder (Paralichthys olivaceus) (Xiu et al., 2021), and Nile tilapia (Oreochromis niloticus) (Li et al., 2018). However, knowledge of detailed immune regulation roles of lncRNAs in teleost is still limited.

Black rockfish (Sebastes schlegelii) is one of the common economic fish species that widely distributed in the east coast of China, Korea, and Japan (Wang et al., 2020). The gram-negative pathogen Edwardsiella tarda in seawater that seriously endangers the aquaculture environment. The infection caused by E. tarda can result in systematic disease to fish, which can eventually lead to mass death (Xu and Zhang, 2013). To study the mechanism of black rockfish in response to pathogen invasion, several immune molecules were identified in recent studies including MyD88 (Shanaka et al., 2019), CTL (Du et al., 2018), CCL25 (Wang et al., 2020), TLRs (Cao et al., 2020), and NLRs (Cao et al., 2021a). Moreover, transcriptome analysis of non-coding RNAs in black rockfish, such as miRNAs (Im et al., 2020) and cricRNAs (Cao et al., 2021b) were conducted to explain their regulatory functions in black rockfish. However, the function of lncRNAs in black rockfish has not been studied so far. Accordingly, the differentially expressed lncRNAs and mRNAs between different infection groups and the control group were analyzed in the intestine of black rockfish. Potential lncRNA-mRNA interactions were thereafter analyzed based on the co-expression analysis of differentially expressed lncRNAs and mRNAs. Meanwhile, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed with the targeted mRNAs of lncRNAs. The results will provide novel knowledge about lncRNAs in the intestinal immune responses in black rockfish, and serve as important resources for further investigating roles of lncRNAs during pathogen infections in teleost.

2 MATERIAL AND METHOD 2.1 Bacterial challenge and sampling

Healthy black rockfish were purchased from a fish farm in Qingdao, Shandong Province, China with length of about 15±2 cm. The experimental protocols were approved by the Committee on the Ethics of Animal Experiments of Qingdao Agricultural University Institutional Animal Care and Use Committee (IACUC). E. tarda was isolated from diseased black rockfish. Subsequently, the purified E. tarda was inoculated in Luria-Bertani (LB) medium at 28 ℃ with 180 r/min. The concentration of E. tarda was 1×10⁷ CFU/mL in the experimental group tank. The survival rate for 1×10⁷ colony forming unit (CFU)/mL 4 h is about 62%, which was detected in our previous study (Cao et al., 2021a). In particular, individuals in the experimental group were cultured in tank with E. tarda for 4 h and then returned back to seawater. Meanwhile, the individuals in control group (CON) were cultured with seawater. The bath-challenge include four infection time points of 2 h, 6 h, 12 h, and 24 h. Therefore, the intestinal tissues from experimental groups were designated as EI2H, EI6H, EI12H, and EI24H, respectively. Furthermore, each time point had three replicates, and each replicate included 6 random individuals. All samples were flash-frozen in liquid nitrogen and then stored in a -80-℃ ultra-low freezer until preparation of RNA.

2.2 RNA isolation, cDNA library construction, and sequencing

RNAs were extracted and purified using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) and treated with RNase-Free DNase I reagent (TIANGEN, Beijing, China). Subsequently, 1% agarose gels were used to check the contamination and degradation. And the RNA concentration was examined in NanoPhotometer^® spectrophotometer. Meanwhile, the RNA Nano 6000 Assay Kit was used for measuring the integrity and quantity of RNA in Bioanalyzer 2100 system. According to manufacturer's suggestion, 5-μg total RNA per sample was used to construct sequence libraries with Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, USA). Then, the prepared libraries were sequenced on an Illumina Hiseq 4000 platform. When the raw data obtained, five sections of invalid data were removed including the reads with 5′ adapter, the reads without 3′ adapter, the reads with more than 10% N, the reads with poly DNA bases and the reads with low sequencing quality (bases with low quality values represented more than 50% of the entire reads). After above steps, the clean reads were obtained. And the parameters of clean data including Q20, Q30, and GC content were calculated.

2.3 Identification of lncRNAs

Firstly, the software HISAT2 was used to process the clean reads and mapped them onto the reference genome. For selecting lncRNAs in transcript libraries, filter conditions were set for all the clean-reads as follows: (1) selection of transcripts with exon ≥2, length > 200 bp; (2) filtering out the transcripts which overlapped with database annotated exon regions; (3) removal of the transcripts with protein-coding capability based on the information in CNCI, Pfam, and CPC2 database. The names of novel lncRNAs were identified according to the rules of HUGO Gene Nomenclature Committee (HGNC). Novel lncRNA sequences were also inspected and identified based on known lncRNA and mRNA sequences (Supplementary Dataset 1).

2.4 Differential expression analysis of lncRNAs

A direct embodiment of a gene's expression level was the abundance profile of its transcripts. And higher transcript abundance reflected higher levels of gene expression. The expression levels of genes were estimated by the count of sequencing sequences (reads) that mapped to a genomic regions or exonic regions. For quantification, we employed the StringTie's network stream algorithm to calculate FPKM values for each gene. Differential expression analysis was performed using EdgeR software. And the hierarchical clustering method was used to cluster the expression values of the samples. The padj-value were set less than 0.05 to filter differentially expressed genes (DEGs).

2.5 Co-expression analysis between lncRNA and their target mRNAs

The cis-acting target gene prediction and trans-acting target gene prediction were used for finding the target genes of lncRNAs. The specific method of cis-acting target gene prediction was used to select the protein-coding genes which located within 10–100 kb of lncRNAs. As for trans-acting target prediction, Pearson' correlations coefficients were used to calculate trans-acting regulatory elements between the coding genes and lncRNAs.

To better summarize the regulatory functions of lncRNAs in black rockfish, the interaction of lncRNA-mRNA was identified between the differentially expressed lncRNAs and mRNAs. When the number of samples was ≥5, the co-expression acting target genes were selected and merged with the differentially expressed mRNAs for analysis. The differentially expressed lncRNAs and mRNAs were connected by analyzing the association between the target genes of differentially expressed lncRNAs and the differentially expressed mRNAs through co-expression information.

2.6 The gene ontology and KEGG analysis

The gene ontology (GO) is a major bioinformatics database which defined one gene's function into three domains: cellular component (CC), molecular function (MF), and biological process (BP) (Young et al., 2010). Meanwhile, the KEGG (Kyoto Encyclopedia) database focuses on the pooling of pathway maps of genes, RNAs, chemical reactions, and the others (Kanehisa et al., 2008). In order to perform GO and KEGG enrichment analysis in black rockfish, the target genes of differentially expressed lncRNAs with length corrected were processed by the ClusterProfiler R package. And the enrichment results with corrected P-values less than 0.05 were considered significant.

2.7 Real-time quantitative PCR (qRT-PCR) validation of lncRNA and mRNA expression analysis

For validation of the expression analysis of lncRNAs and mRNAs in sequencing, total RNA was extracted from black rockfish. PrimeScript^TM RT reagent Kit (TaKaRa, Japan) was used to reversely transcribe RNA into cDNA according to the manufacturer's instructions. Elongation factor-1 alpha (EF1-α) was chosen as internal control. The validated qRT-PCR primers were designed by PrimerQuest Tool (https://www.idtdna.com/pages/tools/primerquest). Subsequently, CFX96 real-time PCR detection system (Bio-Rad Laboratories, Hercules, CA) was used to verify the expression levels of lncRNA and mRNA. The qRT-PCR reaction system is as follows: 10 μL of SYBR^® Premix Ex TaqTM II (TliRNaseH Plus), 0.6 μL of each primer, 2 μL of the 10-fold dilution cDNA, 6.8-μL nucleotide-free water. The qRT-PCR running program is set as follows: 5 min at 95 ℃, followed by 35 cycles of 95 ℃ for 5 s, 56 ℃ for 30 s and 72 ℃ for 30 s, then up to 95 ℃ with a rate of 0.1 ℃/s increment. After qRT-PCR, 2^-ΔΔCt method was used to calculate the expression fold change of lncRNA and mRNA (Livak and Schmittgen, 2001).

3 RESULT 3.1 Sequencing data

The sequencing libraries were named by their infection time points as: CON, EI2H, EI6H, EI12H, and EI24H. The sequenced reads were determined and classified into four data types include Raw_reads, Clean_reads, Raw_bases, and Clean_bases. A total of 2 935 704 384 raw reads were produced in 15 groups. Meanwhile, the 2 872 246 718 clean reads were generated after removing the low-quality, adapter-containing reads, and N-containing reads from the raw reads. All libraries had available sequencing quality with more than 15.06 Gb clean bases for each library. Meanwhile, in 15 sequencing libraries, the sequencing quality related parameters with the values of Q20 ≥97.93, Q30≥94.07, and error rate ≤0.03 (Supplementary Table S1).

A total of 9 311 lncRNAs were obtained from clean reads after quality control. The length of lncRNAs ranged from 201 to 71 349 bp, and the number of exons ranged from 2 to 37. Overall, the exon numbers of novel lncRNA were mostly in the range of 0–15, much smaller than that of mRNAs (Fig. 1c). The distribution tendency of open reading frames (ORF) is similar to that of exons (Fig. 1b), and the ORF lengths of lncRNA were all located within 1 000 bp, while those of mRNAs can reach up to 3 000 bp. In terms of length, novel lncRNAs had similar trends to annotated mRNAs, with the majority of the full length residing between 2 000–2 500 bp, and novel RNAs had a broader range and general longer full length (Fig. 1b). The number of included lncRNAs showed a decreasing trend as the length range increased stepwise (Fig. 1d). The proportion of exon, intergenic, and intronic in all samples were presented in Fig. 1e. Among the different types of lncRNAs, 65.7% of the discovered lncRNAs belong to long intergenic non-coding RNA (lincRNA), 16.2% of them belong to antisense, and 18.1% of them were overlapping. In addition, there were no sense intronic lncRNAs type existed (Fig. 1f).

3.2 The differential expression analysis of lncRNAs after E. tarda infection

A total of 102 lncRNAs underwent changes in expression trends, with 10 up-regulated and 12 down-regulated at 2 h, 34 up-regulated and 16 down-regulated at 6 h, 8 up-regulated and 16 down-regulated at 12 h, and 16 up-regulated and 12 down-regulated at 24 h (Fig. 2a). The most up-regulated lncRNAs were observed at 6 h with 34. Subsequently, 66 lncRNAs with significant expression changes at each time point were selected and subjected to expression heatmaps (Fig. 2b). Obviously, evm.TU.Chr8.478 was down-regulated at EI2H. Meanwhile, evm.TU.Chr5.1155 showed significant up-regulation at EI6H. Particularly, the lncRNA evm.TU.Chr22.794 showed obvious up-regulation trend at four infection time points, whereas XLOC_012887 showed obviously down-regulation trends at four infection time points. Overall, the number of up-regulated lncRNAs was slightly higher than down-regulated lncRNAs.

3.3 The enrichment analysis of lncRNA target mRNA

A total of 3 348 target mRNA were confirmed between groups to groups for further research. And the GO enrichment analysis was performed on all the target mRNAs (Fig. 3). It mainly contained three categories: cellular component, biological process and molecular function. The results of GO enrichment at 2 h, 12 h, and 24 h after infection with E. tarda showed higher number of genes involved in molecular function than the other two functions. Differently, the genes of the process function accounted for a higher proportion in EI6H group. Specifically, in EI2H group, the number of genes possessing catalytic activity was the greatest, followed by hydrolase activity. In EI6H group, the highest number of genes had the defense response function. Meanwhile, peptidase activity emerged as the function performed by the largest number of genes in EI12H group. Finally, in EI24H group, the number of genes with catalytic activity was significantly higher than the other enriched functions.

KEGG enrichment of target gene also represents different functional enrichment at 2-h, 6-h, 12-h, and 24-h infection time points (Fig. 4). The EI2H group showed most genes played a role in metabolic pathways, especially for the inositol phosphate metabolism, glycerophospholipid metabolism, and carbon metabolism. There were also a number of immune-related pathways such as nucleotide oligomerization domain (NOD)-like receptor signaling pathway in enrichment analysis results. Meanwhile, the EI6H and EI24H groups showed similar results to the previous group. The difference was that the classic immune pathway, RIG-I-like receptor signaling pathway was appeared at these infection time points. And for EI12H group, in addition to metabolic and immunity pathways such as phagosome, calcium signaling pathway and focal adhesion, neuroactive ligand-receptor interaction and lysosome were also functionally clustered with a larger number of genes. Besides, few genes with high rich factor and low q-value displayed in steoid biosynthesis from EI2H, EI6H, and EI24H groups. Especially in EI24H group, the enrichment was significant in selenocompound metabolism.

3.4 Potential interaction between lncRNAs and DEmRNAs

Among the results of the GO and KEGG enrichment analysis, 38 differentially expressed mRNAs (DEmRNAs) in 4 groups (EI2H, EI6H, EI12H, and EI24H) which had functions in important pathway were further selected, such as metabolic or immune pathway. Finally, 8 DEmRNA (evm.TU.Chr9.463_ MMP13, evm.TU.Chr13.188_CARNMT2, and evm. TU.Chr16.796_BET1, evm.TU.Chr2.1072_GATM, evm.TU.Chr8.51_TNFSF13B, evm.TU.Chr16.518_ FCGBP, evm.TU.Chr2.1283_MUC5B, evm. TU.Chr7.943_DHCR3), which showed consistent trends of regulation with differentially expressed lncRNAs at two or more time points were selected. Meanwhile, 47 lncRNAs with strong co-expression for these DEmRNA were also selected (Fig. 5). In this analysis, most lncRNAs were co-expressed with one mRNA. Meanwhile, there was also a member of lncRNAs co-expressed with two and more DEmRNAs. Among them, the number of lncRNA co-expressed with two DEmRNA was thirteen, while the number of lncRNA co-expressed with three DEmRNAs was four (XLOC_067928, evm.TU.Chr1.1983, XLOC_030904, and XLOC_032649). In particular, only four DEmRNAs (Chr 13.188_CARNMT2, Chr 21.072_GATM, Chr 16.518_FCGBP, and Chr 2.1283_MUC5B) were the focal objects when one lncRNA was co-expressed with multiple DEmRNAs. The remaining DEmRNAs were all co-expressed with only a single lncRNA.

The identified eight DEmRNAs obtained from black rockfish were listed in Table 1. The function prediction of DEmRNAs were obtained according to the NR database. Among them, DEmRNA evm.TU.Chr8.51_TNFSF13B was identified as tumor necrosis factor like superfamily member 13b (TNF13B), a gene involved in the positive regulation of T-cell proliferation and associated with B-cell maintenance homeostasis in human and an important immune factor. Evm.TU.Chr9.463_ MMP13, evm.TU.Chr13.188_CARNMT2, and evm. TU.Chr7.943_DHCR3 were annotated as Collagenase 3 (MMP-13), Carnosine N-methyltransferase 2 and 3-beta-hydroxysteroid-Delta (8). Moreover, other DEmRNAs had different biological activities like mediating protein transport (evm.TU.Chr16.796_ BET1), involving in creatine biosynthesis (evm. TU. Chr2.1072_GATM), maintaining the mucosal architecture (evm.TU.Chr16.518_FCGBP), binding collagen and Ca²⁺ (evm.TU.Chr9.463_MMP13), and inhibited endopeptidases activity (evm.TU.Chr2.1283_ MUC5B) (Fig. 6).

3.5 qRT-PCR verification of selected lncRNA and mRNA

To verify the accuracy of the sequencing information, 4 DEmRNA and 8 lncRNAs were randomly selected for qRT-PCR verification. The results showed the gene expression trends of obtained by qRT-PCR were generally consistent with transcriptome data (Fig. 7). In detail, the expression patterns of most lncRNAs (XLOC_074289, XLOC_043541, XLOC_064520, XLOC_030548, XLOC_090070, evm.TU.Chr11.770, and XLOC_067928) calculated by qRT-PCR showed similar trends with those from sequencing results at two marked infection time points. For example, XLOC_074289 exhibited 6.22, 6.22, 6.14, and 6.20 fold down-regulation at 2 h, 6 h, 12 h, and 24 h infection time points based on the sequencing data. We found that the expression levels of XLOC_074289 from qRT-PCR downregulated with 3.03, 1.47, 2.00, and 2.88 fold, respectively. Similar trends were found in XLOC_043541 and XLOC_064520. As for mRNAs, the qRT-PCR fold change and transcriptome data for evm.TU.Chr8.51_TNFSF13B showed similar up-regulate trend from 2 h to 6 h. In contrast, evm. TU.Chr16.518_FCGBP showed down-regulate trend from 6 h time point to 24 h time point. In summary, the expression trends of lncRNAs and mRNAs obtained by qRT-PCR showed 71% consistency with transcriptome data. Therefore, the research on the co-expression of lncRNAs and DEmRNAs based on the transcriptome data were basically reliable.

4 DISCUSSION

Currently, E. tarda infection remains one of the pathogens that present a serious threat to black rockfish. In order to study the defense mechanism of black rockfish against this pathogen, we focused our research on long non-coding RNA which has been demonstrated with pivotal roles in various biological processes, especially the gene expression regulation, including transcriptional regulation, posttranscriptional control, and epigenetic processes. In recent years, the development of high-throughput sequencing technologies facilitated the discovery of lncRNAs. It has been found that lncRNAs play important roles in the regulation of immune responses in teleost (Geng and Tan, 2016). For example, the identified lncRNA-WAS and lncRNA-c8807 in grass carp (Ctenopharyngodon idellus) can interact with miRNA and activated NF-κB pathway to resist the invasion of pathogens (Fan et al., 2021). Moreover, the evolutionary conserved lncRNA were found in whitefish and could serve as biomarkers of liver injury (Florczyk et al., 2021). In the intestine of olive flounder, the identification of lncRNA from RNA-sequencing had also been carried out. And their target genes were enriched in numerous immune-related processes and exhibited a strong correlation with immune-related signaling pathways (Xiu et al., 2021). However, the relevant study of lncRNAs in black rockfish was still lacking. In this research, we obtained a total of 9 311 lncRNAs from one control group and four treatment groups infected with E. tarda. And the related sequencing quality parameters such as the value of Q20 ≥97.93, Q30 ≥94.07, and error rate ≤0.03 in all the 15 libraries can corroborate the sequencing data with high confidence.

The results of full length, ORF length, and number of exons for the novel lncRNAs presented different concentrated positions. Meanwhile, the ranges of the three parameters of lncRNAs were less than annotated mRNA and novel mRNA. It was related to the filter strategy of the transcriptome and the properties of the lncRNAs. From the identified lncRNAs classification charts, lincRNAs accounted for the highest proportion (65.7%). Based on the research recently, FANTOM5 consortium used cap analysis for the number of 27 919 human lncRNA, 13 105 (46.9%) of them were lincRNAs, which showed that lincRNAs also had a high proportion in human genome (Ransohoff et al., 2018).

GO and KEGG analyses were performed on 3 348 target mRNAs of lncRNAs, which were selected by cis-acting target gene prediction and trans-acting target gene prediction methods. Four kinds of most obvious enrichment items showed in GO analysis included catalytic activity, hydrolase activity, defense response, and peptidase activity. Actually, lncRNAs had also been identified with relevant activities in other species. For example, the catalytic activity of lncRNA was discovered in human endothelial cells which were adhered by prostate tumor cells PC-3M (Pan et al., 2021). In addition, we found that there were multiple lncRNAs co-expressed with the immune defense genes, such as lncRNA transcript-12631 and transcript-12631 in Spiroplasma eriocheiris (Ren et al., 2020). Meanwhile, a nuclear enriched lncRNA targeting ubiquitin carboxy terminal hydrolase L1 (Uchl1) could increase the transcription level of Uchl after neurodegenerative diseases in mouse (Carrieri et al., 2012). Not only in mammal, in Epinephelus coioides, lncRNAs can regulate genes related to serine hydrolase activity and mediate immune response (Tang et al., 2019). As for the peptidase activity, Liang et al. (2020) found lncRNA 2810403d21rik/ak007586/mirf can indirectly target usp15 (ubiquitin specific peptide 15) attenuated ischemic stress, which induced by cell death attenuated ischemic stress via microRNA competition pathway (Liang et al., 2020). The KEGG analysis showed the target mRNAs existed enrichment in several signaling pathways, especially in the metabolic pathway in black rockfish. For the differentially expressed lncRNAs selected from the gene panel acquired from hepatocellular carcinoma, the targeted mRNAs from them had obvious metabolic pathway enrichment (Wang et al., 2021b). And the metabolic-related pathway was further enriched by mRNA in olive flounder in previous research (Xiu et al., 2021). In addition, a partial of mRNA was also enriched in neuroactive ligand-receptor interaction and lysosome pathway with lower q-value.

In black rockfish, a total of 102 lncRNAs were identified after E. tarda infection. Among which, the lncRNA evm.TU.Chr5.1155 had particularly high expression at 6 h after E. tarda infection. Its target mRNA was G-protein coupled receptor activity. Recent studies had shown that lncRNA LINK-A in mouse mammary glands could facilitated crosstalk between the G-protein coupled receptor pathways and phosphatidylinositol-(3, 4, 5)-trisphosphate, finally causing metastatic mammary gland tumors (Hu et al., 2019). Thus, we can deduce that lncRNA evm. TU.Chr5.1155 can regulate the immune pathway of black rockfish with targeting mRNA. To search for functional DEmRNAs as well as lncRNAs with relevant regulatory roles in black rockfish. A total of 8 DEmRNAs with 47 related lncRNAs were selected from target mRNA and co-expression lncRNA at two or more E. tarda infection time points. As for selecting lncRNA with biological regulatory functions, the DEmRNAs involved in metabolism and immunity had been focused on. Especially, the DEmRNAs evm.TU.Chr8.51_TNFSF13B was identified as TNF13B which is also known as B-cell activating factor (BAFF) or B-lymphocyte stimulating factor (BLyS) (Ai et al., 2011). It was an important member of TNFSF family, the cytokines linked to diverse immunological and developmental pathway (Glenney and Wiens, 2007). The immune function of BAFF has also been characterized in various fish species, such as mefugu (Takifugu obscurus) (Ai et al., 2011), Japanese sea perch (Lateolabrax japonicus) (Cui et al., 2012), rainbow trout (Granja and Tafalla, 2019), and grass carp (Ctenopharyngodon idella) (Pandit et al., 2013).

Particularly, the DEmRNA evm.TU.Chr2.1283_ MUC5B and evm.TU.Chr16.518_FCGBP were previously considered to be important molecules in intestinal immunity in teleost. The evm. TU.Chr2.1283_MUC5B was annotated as Mucin-5B. In previous research, to block the invasion of pathogenic organisms, the fish bodies were always covered with a mucus layer which especially existing in the tissues like gut, epidermis, and gill (Subramanian et al., 2008). And the composition of mucus was very complex, containing mucins, immunoglobulins, complement, lectins, lysozyme, and others (Voynow and Rubin, 2009). Among them, mucin was expressed by epithelial cells and played a role in maintaining homeostasis in various organs. And the MUC2, MUC5AC, and MUC5B have received attention as three protein molecules that have been implicated in a variety of infection as cancer (Conze et al., 2010; Dong et al., 2020), gastric tumors (Mejías-Luque et al., 2010), and idiopathic interstitial pneumonia (Lou et al., 2020). In teleost, several mucin related gene segments were also explored in sequencing information of multiple fish species like Scophthalmus maximus (Gao et al., 2021), Salmo salar (Micallef et al., 2012), and Sparus aurata (Pérez-Sánchez et al., 2013). Especially for common carp, the MUC5B showed down-regulation after CyHV-3 associated disease in the skin which can cause secondary infection among other fish (Adamek et al., 2013). Moreover, the carp MUC5B was identified to have high similarity to its counterpart in mammals and birds (van der Marel et al., 2012). It means the MUC5B had a certain degree of conservation in the evolutionary process. The DEmRNA evm.TU.Chr16.518_FCGBP was identified as IgGFc-binding protein, which also can be named as FCGBP, an important intestinal mucus protein interacted with the Fc portion of IgG and with MUC2 (Johansson et al., 2009; Ehrencrona et al., 2021). Moreover, the FCGBP was ubiquitous in vertebrates, has a conserved N-terminal domain. Although its specific mechanism of action was not yet clear in teleost, its differential expression in numerous tumors made researchers interested in its immune function (Wang et al., 2021a). As for aquatic, in hepatopancreas of Sinonovacula constricta, the IgGFc-binding protein showed significant expression after Vibrio parahaemolyticus infection through the analysis of transcriptomes (Zhao et al., 2017). There is no doubt that the immune function of two DEmRNAs evm.TU.Chr2.1283_MUC5B and evm. TU.Chr16.518_FCGBP in black rockfish need further study.

For metabolism function, the DEmRNA evm. TU.Chr9.463_MMP13, evm.TU.Chr13.188_ CARNMT2, and evm.TU.Chr7.943_DHCR3 were identified as Collagenase 3, Carnosine N-methyltransferase 2, and 3-beta-hydroxysteroid-Delta(8), Delta(7)-isomerase, which performed collagen catabolism, carnosine metabolism, and sterol metabolism respectively in biological processes. Although there are limited studies on these metabolism-related molecules in teleost, their biological functions were also proved in other species. For example, Collagenase 3 played an important role in cartilage erosion independent of proteoglycan loss in mutant mice. And Carnosine N-methyltransferase 2 can transfer a methyl group to carnosine (β-alanyl-l-histidine) and closely related to the production of anserine (Cecilia Berin and Chehade, 2010), which act to retard and reverse non-enzymatic glycation (Szwergold, 2005). For 3-beta-hydroxysteroid-delta (8), delta (7)-isomerase, it is an enzyme that mainly played a role in cholesterol synthesis. The mutations in this gene can cause mosaic patchy epidermal hyperkeratinization (Man et al., 2016). These metabolic molecules and the co-expressed lncRNAs provided material for in-depth study of fish metabolic system.

5 CONCLUSION

The expression patterns of lncRNAs in the intestine of black rockfish were investigated by high-throughput sequencing technique, and the functional enrichment analysis such as GO and KEGG were performed. The results show that DEGs are enriched in immune and metabolic-related pathways. Totally, 8 DEmRNAs and 47 lncRNAs were obtained from 3 348 target mRNAs and 1 998 associated lncRNAs. Meanwhile, the expression patterns of DEmRNAs and lncRNAs from RNA-seq were verified by qRT-PCR. Finally, the functions of immune and metabolic-related DEmRNAs and lncRNAs in black rockfish were discussed. This study laid the foundation for further investigation of the regulatory roles of lncRNAs in the intestinal immune and metabolic response of black rockfish.

6 DATA AVAILABILITY STATEMENT

The mRNA analyzed during this study is included in the published article available at https://doi.org/10.1007/s00343-021-1230-7. The lncRNA dataset analyzed during the current study is available in the Supplementary Dataset 1.

Electronic supplementary material

Supplementary material (Supplementary Dataset 1 & Table S1) is available in the online version of this article at https://doi.org/10.1016/B978-1-4377-0271-2.00044-4.