Identification and characterization of immune-related lncRNAs and lncRNA-miRNA-mRNA networks of Paralichthys olivaceus involved in Vibrio anguillarum infection

Abstract

Background

Long non-coding RNAs (lncRNAs) structurally resemble mRNAs and exert crucial effects on host immune defense against pathogen infection. Japanese flounder (Paralichthys olivaceus) is an economically important marine fish susceptible to Vibrio anguillarum infection. To date, study on lncRNAs in flounder is scarce.

Results

Here, we reported the first systematic identification and characterization of flounder lncRNAs induced by V. anguillarum infection at different time points. A total of 2,368 lncRNAs were identified, 414 of which were differentially expressed lncRNAs (DElncRNAs) that responded significantly to V. anguillarum infection. For these DElncRNAs, 3,990 target genes (named DETGs) and 42 target miRNAs (named DETmiRs) were identified based on integrated analyses of lncRNA-mRNA and lncRNA-miRNA expressions, respectively. The DETGs were enriched in a cohort of functional pathways associated with immunity. In addition to modulating mRNAs, 36 DElncRNAs were also found to act as competitive endogenous RNAs (ceRNAs) that regulate 37 DETGs through 16 DETmiRs. The DETmiRs, DElncRNAs, and DETGs formed ceRNA regulatory networks consisting of 114 interacting DElncRNAs-DETmiRs-DETGs trinities spanning 10 immune pathways.

Conclusions

This study provides a comprehensive picture of lncRNAs involved in V. anguillarum infection. The identified lncRNAs and ceRNA networks add new insights into the anti-bacterial immunity of flounder.

Background

Long non-coding RNAs (lncRNAs) of more than 200 nucleotides (nt) structurally resemble mRNAs, but exhibit poor sequence conservation and cannot be translated into functional proteins [1, 2]. LncRNAs exert vital effects on multiple biological processes, including development, reproduction, metabolism, and immunity [3,4,5,6]. Unlike microRNAs (miRNAs), whose post-transcriptional regulation mechanisms have been well characterized [7], the functional mechanisms of lncRNAs remain to be fully elucidated. Evidences show that lncRNAs modulate the expression of genes in close genomic proximity and distant transcriptional regulators via cis- and trans-acting, respectively [8]. Moreover, recent studies have revealed that lncRNAs can act as miRNA sponges to modulate the expressions of target mRNAs through common miRNA response elements (MREs) [9,10,11]. This regulatory mechanism is known as competitive endogenous RNA (ceRNA) activity, which generates a regulatory network across the transcriptome as a whole [12]. LncRNA-mediated ceRNA activity has been shown to be strongly relevant to cancer pathogenesis and provide important diagnostic biomarkers and therapeutic targets [13,14,15].

In teleost fish, infection-associated lncRNAs have been reported in a number of species. For example, it has been shown that tilapia lncRNAs were induced by Streptococcus agalactiae [16], rainbow trout lncRNAs were induced by Flavobacterium psychrophilum [17], orange-spotted grouper lncRNAs were induced by Pseudomonas plecoglossicida [18], and Atlantic salmon lncRNAs were induced by virus (infectious salmon anemia virus) [19], bacteria (Piscirickettsia salmonis) [19], and parasite (Caligus rogercresseyi) [19, 20]. The ceRNA activity controlled by circular RNA (circRNA)-miRNA-mRNA has also been reported to be implicated in fish immune regulation, for example, circRNA-mediated ceRNA was involved in anti-grass carp reovirus response and anti-Edwardsiella tarda response in grass carp and Japanese flounder, respectively [21,22,23]. To date, studies on lncRNA-mediated ceRNA in fish with pathogen infection are scare, except a recent report showing that a lncRNA regulated antiviral responses in miiuy croaker via ceRNA mechanism [24].

Japanese flounder (Paralichthys olivaceus) is an economically important marine fish in north Asia [25]. Flounder culture has been severely threatened by vibriosis, one of the most frequent aquaculture diseases caused by Vibrio spp, in particular Vibrio anguillarum [26]. Studies showed that some protein-coding genes of V. anguillarum, such as VAA [27], OmpK [28] and OmpR [29], are able to induce the immune responses of T and B lymphocytes in founder. Recently, transcriptome and micro-transcriptome analyses revealed that V. anguillarum induced the expression of a large amount immune related genes and miRNAs in Japanese flounder [30, 31]. However, no study on flounder lncRNA has been documented.

In this study, we systematically investigated the lncRNA expression profiles of flounder during V. anguillarum infection at 3 different time points. We identified the differentially expressed lncRNAs (DElncRNAs) induced by V. anguillarum, examined the integrative expressions of lncRNA-mRNA and lncRNA-miRNA, analyzed the target genes (termed DETGs) and target miRNAs (termed DETmiRs) of DElncRNAs, and characterized the immune-related ceRNA networks of DElncRNA-DETmiR-DETG. Our study provides a global profile of lncRNAs in flounder associated with bacterial infection, which adds new insights into the immune response of teleost during bacterial infection. The immune-related lncRNA-miRNA-mRNA networks identified in this study also can serve as potential targets for future investigations on the molecular mechanism of fish immune defense against bacterial pathogens.

Results

Identification and sequence characterization of lncRNAs

In a previous study, transcriptome was conducted to examine the mRNA profiles of Japanese flounder infected with V. anguillarum for 6, 12, and 24 h [30]. In the present study, the dataset was analyzed for lncRNA expression, and 2,368 lncRNAs were identified. Based on their physical locations in the genome, the lncRNAs were classified into intergenic and genic lncRNAs. Specifically, 1823 (76.98 %) lncRNAs are intergenic lncRNAs that overlap no protein-coding loci in the genome of flounder, and 545 (23.02 %) lncRNAs are from genic regions that overlap protein-coding genes in the sense or antisense orientation. Sequence conservation analysis showed that only 40 (1.69 %) lncRNAs had hits with known lncRNAs in other species, including 15 hits in humans (Homo sapiens), 9 hits in mouse (Mus musculus), 9 hits in zebrafish (Danio rerio), 4 hits in rat (Rattus norvegicus), 2 hits in cattle (Bos taurus), and one hit in opossum (Monodelphis domestica) (Additional file 1). Compared with the mRNAs detected in the same samples in a previous study [30], the lncRNAs were shorter in length (1861 bp on average) than the mRNAs (3367 bp on average) (Fig. 1a). The exon number contained in lncRNAs ranged from 2 to 13, with an average of 1.59, which was less than that contained in the mRNAs (average of 3.60) (Fig. 1b). The guanine-cytosine (GC) content of the lncRNAs ranged from 32.22 to 66.80, with an average of 46.44, which was lower than that of the mRNAs (average of 49.11) (Fig. 1c). LncRNAs exhibited a lower absolute value of the minimum free energy (MFE), an index evaluating the stability of the secondary structure of RNAs, indicating that the secondary structures of the lncRNAs were less stable than that of the mRNAs (Fig. 1d). The average expression level of the lncRNAs was 3.09, which was lower than that of the mRNAs (average FPKM of 9.31) (Fig. 1e).

Fig. 1

The features of lncRNAs versus mRNAs in Japanese flounder. a Sequence length. b The exon number contained in lncRNA or mRNA. c The guanine-cytosine (GC) content. d Minimum free energy (MFE). e Expression. FPKM, Fragments per kilobase of transcript per million mapped reads

Identification of V. anguillarum-induced lncRNAs

After V. anguillarum challenge, 414 lncRNAs showed differential expressions at the 3 time points, and these lncRNAs were named DElncRNAs (Additional file 2). Specifically, at 6 h post-infection (hpi), 120 and 185 DElncRNAs were significantly up- and down-regulated, respectively; at 12 hpi, 78 and 100 DElncRNAs were significantly up- and down-regulated, respectively; at 24 hpi, 70 and 110 DElncRNAs were significantly up- and down-regulated, respectively (Fig. 2a, b). Seventy-nine (19.1 %) DElncRNAs were differentially expressed at all 3 time points (Fig. 2c). To validate the identified DElncRNAs, qRT-PCR was performed to determine the expressions of 8 DElncRNAs. The results showed high correlation coefficients (ranging from 0.86 to 0.99) with that of RNA-seq (Fig. S1), confirming the differential expression patterns of these DElncRNAs.

Fig. 2

Differentially expressed lncRNAs (DElncRNAs) detected at different time points. a Number of DElncRNAs at 6, 12 and 24 h post-infection (hpi). “Up” and “Down” indicate up- and down-regulated expression. b The expression profiles of DElncRNAs in different groups at different time points. For convenience, the control groups and V. anguillarum-infected groups were labeled with the capital letters “C” and “V”, respectively. The number after “C”/ “V” indicates hpi. c Venn diagram showing overlapping DElncRNAs at different time points

Identification of the target genes (mRNAs) of DElncRNAs and functional enrichment based on DElncRNA-target interactions

In order to identify interactive lncRNA-mRNA pairs, lncRNA and mRNA co-expression and co-localization analyses were performed to predict trans- and cis-interactions, respectively. A total of 7,140 mRNAs were found to exhibit significantly strong correlations (|r| > 0.9, and p < 0.05) with DElncRNAs in expression. Of these putative interacting mRNAs, 3,975 were differentially expressed after V. anguillarum infection, 3,922 of which were physically far away from DElncRNAs. These 3,922 genes were considered as trans-differentially expressed target genes of DElncRNAs (trans-DETGs). Genomic location analysis showed that 59 DElncRNAs were located near 278 mRNAs, 68 of these mRNAs were differentially expressed after V. anguillarum infection and were considered as cis-DETGs of DElncRNAs. Additionally, 53 of the 68 cis-DETGs were strongly correlated (|r| > 0.9, and p < 0.05) in expression with their respective DElncRNAs. In total, 3,990 DETGs (Additional file 2) were identified for the 414 DElncRNAs.

To gain insight into the biological processes in which the DElncRNAs were involved, functional enrichment network was constructed based on DElncRNA-DETG interactions and functional enrichment of the DETGs (Fig. 3). In the network, 25 pathways were highly enriched, at least 20 of which were associated with immunity, including TLR signaling pathway, TNF-α signaling pathway, NOD pathway, inflammatory response pathway, IL-6 signaling pathway, type II interferon signaling, complement activation pathway, complement and coagulation cascades, apoptosis, FasL and stress induction of HSP regulation, and BCR and TCR signaling pathways (Fig. 3).

Fig. 3

Functional enrichment network and the highly enriched pathways of DElncRNAs. DElncRNAs are indicated by green nodes. The pathways are indicated by squares colored from yellow to red according to their cumulative enrichment values (from low to high)

Identification of the target miRNAs of DElncRNAs

In a previous study of micro-transcriptome analysis, 1,218 miRNAs were identified in flounder [31]. In the present study, all the 1,218 miRNAs were predicted to be the targets of DElncRNAs. These miRNAs were subjected to correlation analysis with the expressions of DElncRNAs. The result showed that 74 miRNAs were significantly and negatively correlated with 164 DElncRNAs in expression. Further analysis was conducted with the 74 miRNAs based on their responses to V. anguillarum infection, and only the miRNAs with differential expression after V. anguillarum infection were selected. Finally, 42 miRNAs were identified as differentially expressed target miRNAs of DElncRNAs and were named DETmiRs (Additional file 2), whose expressions were both significantly regulated by V. anguillarum and significantly correlated (p < 0.05) with DElncRNAs expressions in a negative manner. The expression patterns of six pairs of DElncRNA-DETmiR, i.e., pol-lnc735-miR-21-y, pol-lnc735-pol-miR-21-3p, pol-lnc491-pol-miR-21-3p, pol-lnc131-pol-miR-n199-3p, pol-lnc163-pol-miR-n199-3p, and pol-lnc491-miR-221-x, were validated by qRT-PCR (Fig. S2). The results showed that in each pair, the expressions of DElncRNA and DETmiR were significantly (p < 0.05) and negatively correlated, with correlation coefficient r ranging from − 0.83 to − 0.96 (Fig. S2).

Construction of immune-related ceRNA networks of interactive DElncRNA-DETmiR-DETG

Integrated analyses of the interactions of DElncRNAs-DETmiRs and DEmiRs-DETGs, as well as the competitions of DElncRNAs-DETGs through MREs, were performed. As a result, 87 DElncRNAs, 28 DETmiRs, and 609 DETGs with interactive relationships were identified. Functional enrichment analysis revealed that these 609 competitive endogenous DETGs (ceDETGs) were involved in 10 immune-related pathways, including the signaling pathways of TLR, IL-6, IL-1, TNF-α, BCR, and TCR, as well as complement and coagulation cascades, apoptosis modulation and signaling, cytokines and inflammatory response, and lymphocyte TarBase (miR-targeted genes in leukocyte) (Fig. 4). The expression patterns of three DETGs (THBD, SERPINE1, and F10) involved in the pathway of complement and coagulation cascades were validated by qRT-PCR, which showed high correlations with that of RNA-seq (r ranging from 0.89 to 0.93) (Fig. S3). To gain insights into the role of the DElncRNA-DETmiR-DETG trinities in the immune response to V. anguillarum infection, an immune-related ceRNA network was constructed based on functional enrichment of the ceDETGs. The network consisted of 36 DElncRNAs, 16 DEmiRs, and 37 DETGs, which formed 114 interacting trinities that spanned 10 pathways associated with immunity (Fig. 5). The DETGs in the ceRNA trinities included SARM (sterile alpha and armadillo motif-containing protein), A2M (alpha-2-macroglobulin), F10 (coagulation factor Χ), CSF1 (colony stimulating factor 1), CREB1 (cAMP-responsive element-binding protein 1), CREBBP (CREB-binding protein), MAP3K2 (mitogen-activated protein kinase kinase kinase 2), MAP3K3 (mitogen-activated protein kinase kinase kinase 3), CASP2 (caspase-2), PTPN13 (protein tyrosine phosphatase-N13), GAB1 (GRB2-associated-binding protein 1), GAB2 (GRB2-associated-binding protein 2), PDCD4 (programmed cell death protein 4), and SLC25A1 (solute carrier family 25 member 1). In the ceRNA networks, six of the DEmiRs, i.e., pol-miR-n199-3p, pol-miR-n071-3p, miR-6240-x, miR-29-x, miR-11,987-x, and miR-194-y, were key immune-related DEmiRs identified in a previous study [31].

Fig. 4

Immune-related pathways enriched in the competitive endogenous DETGs of DElncRNAs

Fig. 5

Immune-related DElncRNA-DETmiR-DETG ceRNA networks. DElncRNA, differentially expressed lncRNA. DETmiR, differentially expressed target miRNA of DElncRNA. DETG, differentially expressed target gene of DElncRNA. CeRNA, competitive endogenous RNA. The round nodes indicate DETGs. The triangle nodes indicate DETmiRs, of which the key miRNAs identified in a miRNA transcriptome analysis of flounder infected with Vibrio anguillarum are labeled in red. The diamond nodes indicate DElncRNAs. “pol-lnc” indicates Japanese flounder lncRNAs

Discussion

In this study, we examined the lncRNAs of Japanese flounder induced by V. anguillarum. We found that more than 98 % of the flounder lncRNAs had no orthologs in other species, which concurred with the notion that most lncRNAs lack primary sequence conservation [32]. Compared with mRNAs, lncRNAs exhibited lower MFE values, suggesting a flexibility in their secondary structures, which may facilitate their binding to miRNAs and mRNAs [33]. Other characteristics of flounder lncRNAs, including exon number, GC content, and expression level, were similar to that observed in the lncRNAs of yellow croaker, rainbow trout, Atlantic salmon, tilapia, and tongue sole [16, 17, 19, 34, 35]. However, the average length of founder lncRNAs (1861 bp) was longer than that of the lncRNAs of zebrafish (1113 bp), tilapia (764 bp), rainbow trout (400 bp), and Atlantic salmon (400 bp) [19, 34, 36, 37]. Since we used only three individual samples in each group at each time point, there is a possibility that some of the characteristics observed in our study may be different if the sample size is increased.

In this study, a total of 414 DElncRNAs were identified to be affected by V. anguillarum infection in a time-dependent fashion, with more DElncRNAs detected in the early infection stage. LncRNAs play an important role in host immune defense against pathogen infection [38,39,40]. Congruously, we found that flounder DElncRNAs, by targeting DETGs, were highly involved in the pathways associated with immunity. The enrichment of the pathways of TLR signaling and complement activation indicated that V. anguillarum stimulated the pathogen recognition process and initiated host immune response. The pathogen recognition process was also found to play an important role in our previous study on V. anguillarum-induced core immune genes of founder [30]. Pathways associated with inflammation and apoptosis were also strongly enriched, suggesting an involvement of these responses in the clearance of the invading pathogen as reported previously [41, 42]. In rainbow trout and tilapia, pathogen-induced host lncRNAs were shown to be engaged in adaptive immunity, such as TCR signaling and MHC protein complex [16, 17]. In our study, we found that the DETGs of DElncRNAs were enriched in BCR and TCR signaling pathways. These results suggest that fish lncRNAs are able to stimulate both innate and adaptive immune responses, which likely provide more efficient protections against pathogen infection. In mammals, lncRNAs are known to regulate host immunity by serving as ceRNAs to modulate mRNA expression [43,44,45]. In our study, immune-related ceRNA networks were found to be formed by 36 DElncRNAs and their corresponding DETmiRs and DETGs involved in 10 immune pathways, including pathogen recognition, inflammation, apoptosis, and adaptive immunity response, which are discussed below.

Pathogen recognition and killing

In this study, pol-lnc78 was localized in the immune-related ceRNA network and targeted pol-miR-n199-3p, which regulated SARM, a newly identified TLR adaptor. In mammals, SARM is known to be a negative regulator of immune response through the TLR signaling pathway, and suppress LPS- and poly (I:C)-mediated AP-1 activation during pathogen infection [46,47,48]. In fish, SARM expression was down-regulated by LPS, and over-expression of SARM inhibited GCRV (grass carp reovirus) triggered IFN-I response [49]. In accordance with the previous reports, we found that SARM was enriched in the TLR signaling pathway and significantly reduced during V. anguillarum infection at all 3 time points, suggesting a persisted activation of TLR signaling and the pathogen recognition process. In addition to SARM, A2M and F10 were also identified in our study as the targets of pol-lnc126 (via pol-miR-n199-3p/miR-6240-x) and pol-lnc131 (via miR-6240-x), respectively. A2M has a common evolutionary origin with complement components C3 and C4, and was shown to contribute to bacteriostatic activity in amphioxus [50,51,52]. F10 is involved in antibacterial infection through initiating the coagulation cascade [53], and induced by both bacterial (Aeromonas hydrophila) and fungal (Aphanomyces invadans) infections in fish [54]. Taken together, these observations indicate that in flounder, V. anguillarum infection induces ceRNA trinities that likely promote pathogen recognition and bacterial killing.

NF-κB-regulated inflammation response

In this study, we found that pol-lnc79, pol-lnc163, and pol-lnc449 competitively targeted pol-miR-n199-3p against CSF1, and pol-lnc291, pol-lnc449, and pol-lnc1387 acted as ceRNAs that regulated CREBBP by sponging miR-29-x. Pol-miR-n199-3p and miR-29-x were identified as the key miRNAs induced by V. anguillarum infection in our previous micro-transcriptome analysis of the same samples [31]. CSF1 is a cytokine highly implicated in inflammation [55, 56]; CREBBP is a crucial cofactor that binds numerous transcription factors, e.g., human CREBBP cooperates with NF-κB to regulate IL-6 promoter [57, 58]. MAP3K1-3 also exerts a crucial effect on NF-κB activation [59, 60]. In mice, MAP3K3 affects both TNF- and IL-1-induced NF-κB activation [61]. Moreover, MAP3K3 was reported to activate NF-κB through the lncRNA-MALAT1-miR-424-MAP3K3 axis [62]. In our study, pol-lnc198, pol-lnc765, and pol-lnc2111 sponged miR-194-y to regulate MAP3K2 associated IL-1 signaling pathway, and pol-lnc126 sponged miR-146-x to regulate MAP3K3 associated IL-1 and TNF-α signaling pathways. Together, these results suggest that lncRNA-miRNA-mRNA trinities are involved in the regulation of NF-κB-mediated inflammatory response.

Apoptosis

Four lncRNAs, i.e., pol-lnc73, pol-lnc233, pol-lnc491, and pol-lnc2061, were found to regulate CASP2 by competitively binding miR-221-x, and four other lncRNAs, i.e., pol-lnc491, pol-lnc536, pol-lnc735, and pol-lnc1234, regulated PTPN13 by competitively binding pol-miR-21-3p or miR-21-y. CASP2 is an initiator for apoptosis execution and involved in bacteria-induced immune defense in striped murrel [63]. PTPN13 is known to inhibit Fas-induced apoptosis, and its down-regulation significantly increased the survival of human papillomavirus infected patients with squamous cell carcinoma [64, 65]. In this study, we found that CASP2 and PTPN13 were the targets of lncRNAs and down-regulated during infection. The suppression of these target genes suggests that they may be part of the host defensive response against pathogen, or the consequence of bacterial manipulation of the host immune system.

Adaptive immunity

BCR signaling is pivotal for the activation of B cells, and TCR recognition of pathogen-derived peptides is a hallmark of the adaptive immunity [66, 67]. In this study, we found that pol-lnc78, pol-lnc125, and pol-lnc291/pol-lnc491 regulated GAB2, GAB1, and CREB1, respectively, by sponging the target miRNAs, and GAB1/2 and CREB1 were involved in the signaling pathways of BCR and TCR. GAB1/2 are key adaptor/scaffolding factors mediating signal transductions down-stream of BCR and TCR [68], and CREB1 is a CREB/ATF family protein controlling the transcription of various viruses (human T-lymphotropic virus, herpes simplex virus, Epstein-Barr virus, and cytomegalovirus) in humans [69]. These observations indicate that in flounder, V. anguillarum probably induces T and B cell-mediated adaptive immunity via lncRNA-regulated ceRNA axes.

Conclusions

This study identified systematically the lncRNAs of Japanese flounder induced by V. anguillarum. A total of 414 DElncRNAs were detected, which target 3,990 DETGs and 42 DETmiRs. Thirty-six DElncRNAs compete with DETmiRs to regulate DETGs and form immune related ceRNA networks. These observations indicate that in teleost, bacteria-induced immune response is controlled by complicated regulatory networks involving non-coding RNAs and ceRNAs. The results of our study add new insights into the anti-microbial defense of flounder and provide a valuable dataset for future studies on the immune mechanism and disease control of flounder based on non-coding RNAs and ceRNAs.

Methods

LncRNA identification

LncRNA analysis was performed with the previously reported raw RNA sequencing data of 18 spleen libraries of flounder infected with and without (control) V. anguillarum [30]. In this previous study [30], Japanese flounder (purchased from a local commercial fish farm in Shandong Province, China) were injected with V. anguillarum C312 (named group V) or PBS (control, named group C); at 6, 12, and 24 h post infection (hpi), three fish from each group were collected and euthanized with tricaine methanesulfonate (Sigma, St. Louis, USA); spleen was then collected from the fish and used for RNA sequencing [30]. In the present study, the software fastp (v0.12.4) was used to filter the raw reads to obtain high quality reads, which were then mapped to the reference genome of Japanese flounder (GenBank project accession PRJNA369269) using TopHat 2 (v2.0.3.12) [70] as previously reported [30]. Cufflinks (v2.2.1) [71] was used to reconstruct transcripts from the mapped reads as previously reported [30]. The transcripts with the length of more than 200 bp were analyzed using the software coding-non-coding-index (CNCI, v2), a powerful tool to effectively distinguish protein-coding and non-coding sequences independent of known annotations, to filter out the protein-coding mRNAs [72]. The software coding potential calculator (CPC) was used to evaluate the protein-coding potential of the remaining non-coding sequences [73]. Moreover, the non-coding sequences were also searched against the UniProtKB/SwissProt database, and the sequences with blast hits of E-values < 1e-6 were removed. Finally, the transcripts assessed to have no protein-coding potential and have no hits in the database were identified as lncRNAs of Japanese flounder and named pol-lnc1 to pol-lnc2368. The conservation of the lncRNAs was evaluated by aligning with the lncRNAs of other species in database NONCODE v5.0 [74].

Differential expression analysis

The abundance of each lncRNA was counted using RSEM [75], and the expression level was normalized with FPKM (Fragments Per Kilobase of transcript per Million mapped reads) [76]. The R package edgeR (v3.12.1) (http://www.r-project.org/) was applied to conduct differential expression analysis. False discovery rate (FDR) < 0.05 and fold change (FC) in expression > 2 (log₂|FC| > 1) were set as the threshold to identify the differentially expressed lncRNAs (DElncRNAs).

Experimental validation

To verify the DElncRNAs obtained by RNA-seq, quantitative real-time reverse transcription-PCR (qRT-PCR) was employed to test the expression of eight randomly selected DElncRNAs at 6, 12, and 24 hpi (three fish/time point) with the collected spleen tissues that had been used for RNA-seq in the previous study [30]. The expression profiles of DETmiRs were evaluated by qRT-PCR with the collected spleen tissues that had been used for micro-transcriptome in the previous study [31]. Similarly, three DETGs (THBD, SERPINE1, and F10) of DElncRNAs were validated for expression by qRT-PCR. The PCR reactions were performed with SYBR Premix Ex TaqII (TaKaRa, Dalian, China) in QuantStudio 3 Real-Time PCR Systems (Thermo Fisher Scientific, CA, USA) according to the manufacturer’s protocol. At the end of each reaction, the melting curve analysis was conducted to verify the specific of PCR product. The expression levels of lncRNAs, mRNAs and miRNAs was evaluated with 2^−ΔΔCt comparative Ct method [77]. Correlations between qRT-PCR and RNA-seq results, as well as between DElncRNAs and their corresponding DETmiRs, were analyzed using cor.test in R (v3.5.2). The primer sequences were listed in Table 1.

Table 1 Summary of primers used for qRT-PCR in this study

Target genes identification

Based on the trans- and cis-acting of lncRNAs, the target genes for DElncRNAs were predicted by co-expression and co-location analysis [78]. The co-expression analysis was performed using cor.test in R (v3.5.2) based on the lncRNA expression data in this study and the mRNA expression data in a previous study with the same samples [30]. All of the 2,368 lncRNAs identified in this study and the 43,494 mRNAs identified previously [30] were used in the analysis. The expression of mRNAs was counted using RSEM (v1.2.19) and normalized across libraries using FPKM (Fragments Per Kilobase of transcript per Million mapped reads) method as described previously [30]. The genes correlated (correlation coefficient r > 0.9 or r < − 0.9, and p < 0.05) with DElncRNAs in expression were identified as the candidate trans-regulated target genes [79], which were then subjected to analysis of differential expression induced by V. anguillarum using R (v3.5.2). As reported previously [30], the exact negative binomial test in the R package edgeR (v3.12.1) (http://www.r-project.org/) was used to conduct pairwise comparison between group C and group V at each time point, and the threshold for significant difference was set as false discovery rate (FDR) < 0.05 and log2|FC| >1. The genes that both correlated with DElncRNAs expressions and differentially induced by V. anguillarum were subjected to distance evaluation. Finally, the genes that were differentially expressed, correlated with DElncRNAs, and located far from DElncRNAs (> 10 kb) were considered as trans-regulated target genes of DElncRNAs (named trans-DETGs). The co-location analysis was conducted according to the physical distance between the lncRNAs and potential target genes in the genome of flounder. The genes that physically overlapped or were near (< 10 kb) to DElncRNAs were identified as candidate cis-regulated target genes, which were then analyzed for differential expression induced by V. anguillarum using R (v3.5.2). The genes that were both co-located with DElncRNAs and differentially expressed after V. anguillarum infection were considered as cis-regulated target genes (named cis-DETGs). All the trans- and cis-regulated target genes were termed differentially expressed target genes of DElncRNAs (named DETGs).

Functional enrichment network analysis

Enrichment analysis and visualization (ENViz), a Cytoscape plug-in, is generally used to predict biological functions engaged by ncRNAs [80]. In this study, ENViz of Cytoscape (v3.7.1) was applied to the analysis of DElncRNAs and their DETGs to construct the lncRNA-pathway enrichment networks. The pathways with a cumulative enrichment cutoff value > 3 were displayed and colored from yellow to red according to their cumulative enrichment values ranging from low to high, and the pathway with a cumulative enrichment cutoff value > 250 was considered a highly enriched pathway.

Identification of the target miRNAs of DElncRNAs

Three softwares, i.e., RNAhybrid (v2.1.2) + svm_light (v6.01), Miranda (v3.3a), and TargetScan (v7.0), were employed to predict the interacting miRNAs of DElncRNAs using 1,218 miRNAs of flounder identified in a previous study with the same samples [31] and the DElncRNAs identified in this study. The common results predicted by the three algorithms were identified as candidate target miRNAs. For these candidate target miRNAs, R (v3.5.2) was employed to detect the differentially expressed miRNAs induced by V. anguillarum, and cor.test in R (v3.5.2) was used to detect the differentially expressed miRNAs that negatively correlated (r < − 0.7, p < 0.05) with their DElncRNAs based on the lncRNA data in this study and the miRNA data in a previous study with the same samples [31]. The candidate target miRNAs that were both differentially expressed during V. anguillarum infection and negatively correlated with DElncRNAs expressions were identified as the differentially expressed target miRNAs of DElncRNAs and were named DETmiRs.

ceRNA network construction

Integrated analysis was conducted with several datasets, including the lncRNA expression data in this study and the mRNA expression data and miRNA expression data in previous studies with the same samples [30, 31]. Based on the ceRNA hypothesis, i.e., lncRNAs act as miRNA sponges to relieve the expression inhibition of miRNAs on the target mRNAs [12], the interacting DElncRNA-miRNA-mRNA networks were identified as follows: (1) detection of the negatively correlated DElncRNA-DETmiR pairs; (2) detection of the negatively correlated DETG-DEmiR pairs; (3) detection of the positively correlated DElncRNA-DETG pairs, in which, the DElncRNA and DETG of each pair share the same miRNA response elements (MREs). To be stringent, the hypergeometric cumulative distribution function test was applied to obtain the significantly enriched DElncRNA-DETG pairs that targeted the common miRNAs, with the threshold of p < 0.05. After above procedures, the remaining were the intertwined DElncRNA-DETmiR-DETG ceRNA trinities. The immune-related ceRNA network was constructed with Cytoscape (v3.7.1) [81] using the competitive endogenous DETGs enriched in the immune-related pathways and their corresponding DElncRNA-DETmiR in the intertwined DElncRNA-DETmiR-DETG trinities.