Identification of Potential Immune-Related circRNA–miRNA–mRNA Regulatory Network in Intestine of Paralichthys olivaceus During Edwardsiella tarda Infection

Olive flounder (Paralichthys olivaceus) is an important economical flatfish in Japan, Korea, and China, but its production has been greatly threatened by disease outbreaks. In this research, we aimed to explore the immune responsive mechanism of P. olivaceus against Edwardsiella tarda infection by profiling the expression of circRNA, miRNA, and mRNA by RNA-seq and constructing a regulatory circular circRNA–miRNA–mRNA network. Illumina sequencing of samples from normal control (H0), 2 h (H2), 8 h (H8), and 12 h (H12) post-challenge was conducted. Differentially expressed (DE) circRNA (DE–circRNAs), miRNAs (DE–miRNAs), and mRNAs [differential expression genes (DEGs)] between challenge and control groups were identified, resulting in a total of 62 DE–circRNAs, 39 DE–miRNAs, and 3,011 DEGs. Based on the differentially expressed gene results, miRNA target interactions (circRNA–miRNA pairs and miRNA–mRNA pairs) were predicted by MiRanda software. Once these paired were combined, a preliminary circRNA–miRNA–mRNA network was generated with 198 circRNA–miRNA edges and 3,873 miRNA–mRNA edges, including 44 DE–circRNAs, 32 DE–miRNAs, and 1,774 DEGs. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to evaluate the function of the DEGs in this network, and we focused and identified two important intestinal immune pathways (herpes simplex infection and intestinal immune network for IgA production) that showed statistical significance between the challenge and control groups. Furthermore, three critical DEGs (nectin2, MHC II α-chain, and MHC II β-chain) were identified, mapped into the preliminary circRNA–miRNA–mRNA network, and new circRNA–miRNA–mRNA regulatory networks were constructed. In conclusion, we, for the first time, identified circRNA–miRNA–mRNA network from P. olivaceus in the pathogenesis of E. tarda and provided valuable resources for further analyses of the molecular mechanisms and signaling networks.


INTRODUCTION
Circular RNAs (circRNAs), identified from RNA viruses in the 1970s (Sanger et al., 1976), were initially treated as viral genomes or by-products of rare mis-splicing, and thus, they have long been thought to be nonfunctional (Capel et al., 1993). CircRNAs are generated during the process of back-splicing and could be grouped into four categories: circular exonic RNAs (ecircRNAs), circular intronic RNAs (ciRNAs), exon-intron circRNAs (eiciRNAs), and intergenic circRNAs (Qu et al., 2015;Sonja and Sabine, 2015;Li et al., 2017a). During back-splicing, a downstream 5′ splice donor is joined with an upstream 3′ splice acceptor, and the resulting RNA circle is ligated by a 3′-5′ phosphodiester bond at the junction site (Lasda and Parker, 2014;Chen, 2016;Wilusz, 2018). Back-splicing is catalyzed by the canonical spliceosome machinery and modulated by both intronic complementary sequences and RNA binding proteins (Li et al., 2018b). Recent advancements to the high-throughput sequencing technology have benefited large amounts of circRNAs identified in succession from many organisms, such as plants, animals, human beings, fungi, and protists (Li et al., 2018b). Emerging researches suggest that some circRNAs are critical in many physiological and pathological conditions (Li et al., 2018b). For example, expression profiles and knockout experiments proved that circRNAs have been implicated in neuronal function (Rybak-Wolf et al., 2015;Piwecka and Glazar, 2017) and testes development (Capel et al., 1993;Hansen et al., 2013). Besides, more and more circRNAs have been found to be associated with human diseases, such as cancers (Conn et al., 2015;Guarnerio et al., 2016), Alzheimer's disease (Lukiw, 2013), neuronal diseases (Errichelli et al., 2017), and others. In addition, the most recent progresses reveal that some circRNAs are also involved in innate immune responses (Chen et al., 2017b;Li et al., 2017b;Wang et al., 2017). Collectively, considerable evidences proved that circRNAs are not simply accidental by-products but represent an essential part of noncoding RNA families.
Although relevant research on circRNAs is still in its infancy, it is becoming apparent that circRNAs play their regulatory roles through distinct mechanisms. Initially, circRNAs function as miRNA sponges through abundant binding sites for microRNAs and then modulate the activity of miRNAs on their target genes (Hyeon Ho et al., 2009). Remarkably, some circRNAs are strongly associated with cancer progression through competing with miRNAs to influence the expression of target genes that are involved in biological processes, for example, tumor cell proliferation, apoptosis, invasion, and migration (Zhong et al., 2018). Apart from acting as miRNA sponges, circRNAs play multiple functions through affecting splicing of their linear mRNA counterparts, regulating transcription of their parental genes, influencing splicing of their linear cognates, interacting with associated proteins, protein-coding genes, and generating pseudogenes (Li et al., 2018b).
Olive flounder (P. olivaceus) is an important economical flatfish that has been widely cultured in Japan, Korea, and China. The production of P. olivaceus has been greatly threatened by disease outbreaks, including bacteria, virus, and parasites (Kim et al., 2010). Edwardsiella tarda, associated with hemorrhagic septicemia of freshwater and marine fish, could also result in extensive economic losses to aquaculture industry of P. olivaceus (Mohanty and Sahoo, 2007;Xu and Zhang, 2014). It was reported that E. tarda is an important zoonotic and intestinal pathogen, and the intestine was likely the main route of entry to host Wang et al., 2012). Therefore, besides serving as the prime site for absorption of nutrients, intestine represents one of the first-line defense systems (Lauriano et al., 2019). It has been confirmed that intestinal hypo-immunity of fish favors E. tarda infection (Liu et al., 2014). Teleost fish possess a diffuse mucosa-associated immune system in the intestine where B cells act as one of the main responders (Parra et al., 2016). Moreover, IgT + B cells represent the predominant mucosal B-cell subset, and the accumulation of IgT + B cells has been detected in trout intestine after infection (Yu and Kong, 2018). Immunoglobulins produced by these B cells constitute a critical line of defense, which prevents the entrance of pathogens and commensal bacteria into the epithelium (Parra et al., 2016).
It is vitally necessary to understand and apply their immune mechanism against pathogen infection. Over the years, massive efforts have been conducted in exploring the immune mechanism of P. olivaceus at a molecular level (Hwang et al., 2018;Ma et al., 2018;Wu et al., 2018), among which a few researches have been conducted in non-coding RNA . However, there has been no report about the important roles of circRNAs during the immune process of P. olivaceus. In fish, few circRNA researches have been published on teleosts, including half-smooth tongue sole (Cynoglossus semilaevis) (Li et al., 2018a), large yellow croaker (Larimichthys crocea) (Xu et al., 2017), zebrafish (Danio rerio) (Shen et al., 2017), coelacanth (Anne et al., 2014) and grass carp (Ctenopharyngodon idella) (He et al., 2017). Interestingly, the most recent progresses reveal that circRNAs also take part in immune regulation and viral infection . It has been identified that the in vitro synthesized circRNAs would activate RIG-I-mediated innate immune responses, which will provide protection against viral infection (Chen et al., 2017b). Besides, the immune factors NF90/NF110 modulate circRNA biosynthesis and suppress viral infection by interacting with viral mRNAs (Li et al., 2017b). Moreover, it has been speculated that a circRNA-miRNA-mRNA network may be present in grass carp reovirus (GCRV)-infected grass carp, which provides new insights into the immune mechanism underlying grass carp against GCRV (He et al., 2017). However, there are still several intriguing questions that remained to be clarified: 1) Are circRNAs involved in antibacterial immune responses? 2) How do circRNAs contribute to antibacterial immune responses?
In this study, we examined interaction of circRNAs, miRNAs, and mRNAs of P. olivaceus in the pathogenesis of Edwardsiella tarda by high-throughput sequencing. We screened and identified differentially expressed circRNAs, miRNAs, and mRNAs; predicted the potential circRNA-miRNA-mRNA network; analyzed their significant enrichment pathways; and emphasized their implications in antibacterial immunity for the first time. August 2019 | Volume 10 | Article 731 Frontiers in Genetics | www.frontiersin.org

The Experimental Fish and Ethical Statement
Healthy olive flounders were obtained from Huanghai Aquaculture Company (Haiyang, Shandong, China). The fish were acclimatized in a recirculating water system (temperature 20 ± 1°C) for 1 week before processing, during which they were fed twice a day with commercial diet. In order to make sure that all the experimental fish were healthy, the olive flounders were monitored every day; after 1-week acclimation, they were randomly sampled for bacteriological examination. This study was carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health, Qingdao Agricultural University. The protocol was approved by the Committee on the Ethics of Animal Experiments of Qingdao Agricultural University IACUC (Institutional Animal Care and Use Committee).

Bacteria Challenge and Sample Collection
The bacterial of E. tarda were isolated from diseased olive flounders and kept by our laboratory. Before the challenge experiment, E. tarda were incubated in Luria broth (LB) medium at 28°C to mid-logarithmic stage. The concentration was determined by colony-forming unit (CFU) method. Overall, this experiment included four time points, and each time point contained three biological replicates (three fish for each biological replicate). Before E. tarda infection, nine fish were immersed in sterilized media, and their posterior intestines were sampled as the normal control group. The samples were designated as H0 (H0_1, H0_2, and H0_3). During challenge, the experimental groups were immersed in the bacteria solution with a final concentration of 6 × 10 7 CFU/ml for 2 h and then returned to the circulating water system. At 2, 8, and 12 h post-treatment, nine fish from each time point were collected, and their posterior intestines were collected for sequencing. The samples were designated as H2 (H2_1, H2_2, and H2_3), H8 (H8_1, H8_2, and H8_3), and H12 (H12_1, H12_2, and H12_3).

Histopathological Analysis
To observe histopathological changes of intestine in the E. tarda infected P. olivaceus, we took posterior intestines from nine fish at each time point to make pathological sections. Tissue samples were fixed in 4% paraformaldehyde in phosphate-buffered saline (PBS) and then further processed through the following steps: dehydrated in graded ethanol, cleared in xylene, embedded in paraffin, cut into 5-mm sections, and stained with hematoxylin and eosin (H&E) for examination by light microscopy (Licata et al., 2018). The histological measurements for the structures, height of mucosal folds, thickness of lamina propria, inner circular muscular layer, and outer longitudinal muscle were measured and analyzed. The mean ± standard error of mean (SEM) of each structure was compared among all of samples using the analysis of variance with Tukey LSD (SAS 9.4) at the significance level p < 0.05.

RNA Isolation, Library Construction, and Sequencing
Total RNA from samples was extracted by using the TRIzol reagent (Invitrogen, USA). The purity of total RNA was checked by using the NanoPhotometer ® spectrophotometer, its concentration was checked by using Qubit ® RNA Assay Kit in Qubit ® 2.0 Fluorometer, and its integrity was checked by Bioanalyzer 2100 system.
At different time points before (H0) and after (H2, H8, and H12) E. tarda infection, intestine tissues from nine olive flounders were respectively collected and used for circRNA sequencing. Three replicate samples were processed for each time point, and a total of 12 libraries were sequenced. For constructing the library of circRNAs or mRNAs, 5 μg of RNA for each sample was prepared. Then, Epicentre Ribo-Zero™ rRNA Removal Kit (Epicentre, USA) was used to remove rRNA, and ethanol precipitation was applied to clean up rRNA-free residue. Subsequently, the linear RNA was digested with 3 U of RNase R (Epicentre, USA) per microgram of RNA for mRNA library and without RNase R treatment for circRNAs, which was the only difference between library construction of circRNAs and mRNAs. The sequencing libraries were generated by using NEBNext ® Ultra™ Directional RNA Library Prep Kit for Illumina ® (NEB, USA) according to manufacturer's protocol. For circRNA, mRNA, and miRNAs, the library construction and sequencing were operated by Novogene Corporation (China) in the same way as previously reported (Lu et al., 2016). In consideration that the Ribo-Zero library contained both mRNA and lncRNA at the same time, but we were not interested in lncRNA in this research, we set up a series of strict screening conditions to identify and remove lncRNA according to their structural characteristics and functional characteristics. The screening process of lncRNA is shown as follows: 1) Select transcripts with exon number ≥2. 2) Select transcripts with length > 200 bp. 3) Screen transcripts that overlap with the annotated exon area. 4) Cuffquant was used to calculate the expression of each transcript, and the transcripts with fragments per kilobase of transcript per million mapped reads (FPKM) ≥0.5 were selected. (5) CNCI (coding-non-coding index) (v2), CPC (Coding Potential Calculator) (0.9-r2), Pfam Scan (v1.3), and PhyloCSF (phylogenetic codon substitution frequency) (v20121028) were used to predict the coding potential of transcripts, and the intersected transcripts without coding potential from these four software analysis were selected as the lncRNA.

Data Analysis
Raw data (raw reads) were firstly processed through in-house perl scripts (for circRNAs and mRNAs) or customperl and python scripts (for miRNAs). More concretely, reads containing adapter, ploy-N, and low-quality reads from raw data were removed. Then, Q20, Q30, and GC contents of the clean data were calculated. All the downstream analyses were based on the clean data.
The circRNAs were detected and identified using find_circ (Memczak et al., 2013) and CIRI2 (Gao et al., 2018). The workflow of find_circ was as follows: 1) The reads that aligned contiguously to the genome were filtered out, and spliced reads were retained; 2) the terminal parts of each candidate read were mapped to the genome to find unique anchor positions; 3) candidate circRNAs were confirmed when their 3′ end of anchor sequence aligned to the upstream of 5′ end of anchor sequence, and the inferred breakpoint was flanked by GU/AG splice signals. The circos figures were constructed by using Circos software. Mapped small RNA tags were used to search for known miRNA. MiRBase20.0 was used as reference; modified software mirdeep2 (Friedländer et al., 2012) and srna-tools-cli were used to obtain the potential miRNA and draw the secondary structures. Custom scripts were used to obtain the miRNA counts as well as base bias on the first position of identified miRNA with certain length and on each position of all identified miRNA, respectively (Jian et al., 2018). The characteristics of hairpin structure of miRNA precursor can be used to predict novel miRNA. The available software miREvo (Wen, 2012) and mirdeep2 (Friedländer et al., 2012) were integrated to predict novel miRNA through exploring the secondary structure, the Dicer cleavage site, and the minimum free energy of the small RNA tags unannotated in the former steps. Meanwhile, custom scripts were used to obtain the identified miRNA counts as well as base bias on the first position with certain length and on each position of all identified miRNA, respectively . For transcriptome assembly, the mapped reads of each sample were assembled by StringTie (v1.3.1) (Pertea et al., 2016) in a reference-based approach. StringTie uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate fulllength transcripts representing multiple splice variants for each gene locus.
Differential Expression Analysis, Enrichment Analysis, and circRNA-miRNA-mRNA Network Analysis Differential expression analysis between two groups was performed using the DESeq R package (1.8.3). The p-value was adjusted using the Benjamini and Hochberg method. Corrected p-value of 0.05 was set as the threshold for significantly differential expression by default.
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used on significantly differential expressed genes, including host genes of differentially expressed circRNAs and the target gene candidates of differentially expressed miRNAs. Gene ontology (GO) enrichment analysis was implemented by the GOseq R package, in which gene length bias was corrected (Young et al., 2010). GO terms with corrected p-value of less than 0.05 were considered significantly enriched by differentially expressed genes. KEGG is a database resource for understanding high-level functions and utilities of the biological system (Kanehisa et al., 2008), such as the cell, the organism, and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http:// www.genome.jp/kegg/). We used KOBAS software to test their statistical enrichment in KEGG pathways (Mao et al., 2005).
The circRNA-miRNA-mRNA network was developed based on possible functional relationships between DE-circRNAs, DE-miRNAs, and differential expression genes (DEGs). Firstly, the target circRNAs of DE-miRNAs were predicted by scanning for conserved miRNA target sites with MiRanda (Enright et al., 2003); then the interactions between target circRNAs and DE-circRNAs were identified; and finally circRNA-miRNA regulation network was constructed. Secondly, the target mRNAs of DE-miRNAs were predicted by scanning for conserved miRNA target sites with MiRanda; then the interactions between target mRNAs and DEGs were identified; and finally miRNA-mRNA regulation network was constructed. At last, circRNA-miRNA-mRNA network was generated using a combination of circRNA-miRNA network and miRNA-mRNA network with Cytoscape 3.6.1 software (Su et al., 2014), and only the network follows the expression trend of "up-down-up" or "down-updown" was selected for further research. In conclusion, the construction of circRNA-miRNA-mRNA network followed the following principles: circRNAs served as bait, microRNAs served as core, and RNA served as target.

Confirmation of the Expression Level of circRNAs, miRNAs, and mRNAs
To validate the reliability of the data obtained from Illumina sequencing, real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR) and Sanger sequencing were conducted. To confirm the expression pattern of differentially expressed circRNAs, six circRNAs were randomly selected for qRT-PCR and Sanger sequencing. Primer Premier 5 software was used to design their divergent primers. In general, the divergent primers were designed to span the circRNA backsplice junction, and a fragment of 80-to 150-bp length was expected. Total RNA was extracted, digested using RNase-Free DNase (Promega), and then purified. A total of 1 μg of purified RNA was used to prepare first-strand cDNA by using a random 6 mers primers and the PrimeScript 1 st strand cDNA Synthesis Kit (Takara, Japan). qRT-PCR was carried out on a CFX 96 real-time PCR system (Bio-Rad, Hercules, CA, USA). Each RT-qPCR mixture contained 4.2 μl of ddH 2 O, 5 μl of 2× SYBR Green master mix (Aidlab), 0.4 μl of template, and 0.2 μl each of forward and reverse primers. The EF1α gene was used as an internal control. For each sample, three replicates were included. The program for qRT-PCR was as follows: 95°C for 3 min, followed by 40 cycles of 95°C for 10 s and 60°C for 34 s. Relative expression level was calculated using the Pfaffl method (Pfaffl, 2001). Data are shown as means ± SE of three replicates. To further prove that the amplified products were the template with circRNA rather than liner counterparts, products of qRT-PCR were subjected to pEASY-T1 cloning vector (TransGen, Beijing, China) and then sequenced by Tsingke Company (Tsingke, Qingdao, China).
Meanwhile, to further confirm the expression level of miRNAs and mRNAs, 8 miRNAs and 10 mRNAs were randomly selected for qRT-PCR experiments. For miRNAs analysis, small RNA (<200 nt) was harvested using the miRcute miRNA isolation kit (Tiangen Biotech). The amplification reactions were carried out using the miRcute miRNA qPCR detection kit (Tiangen Biotech) with the following conditions: 95°C for 15 min, 40 cycles of two steps (95°C for 5 s and 60°C for 30 s) within triplicate wells of each sample. The relative expression level of miRNA was normalized by 5S rRNA expression. For mRNAs analysis, reaction mixtures and the program for qRT-PCR were the same as circRNAs, and the EF1α gene was again used as an internal control for the normalization of gene expression.

Histopathological Description and Cytokine Expression Analysis
As shown in Figure 1A, the healthy intestine sample contained tunica mucosa, submucosa, mucosal folds, circular muscular layer, and longitudinal muscular layer. The mucosal folds were lined with abundant of columnar epithelial cells and goblet cells. And all the cells appeared to be interrelated and uniformly arranged. For H2, in early infection, hyperplasia of intestinal mucosa was observed. Meanwhile, the blood vessel and the lymph vessel enlarged, which were associated with the increasing thickness of lamina propria (Supplementary Table 1). However, the intestine structure was still integrated, the height of the mucosal folds decreased (Supplementary Table 1), and most of the mucosal folds did not show significant lesions. The columnar epithelial cells and goblet cells were also lined tightly ( Figure 1B). For H8, cellular swelling and hydropic change were observed in intestinal epithelium, some intestinal epithelial cells were shedding, and some became disrupted. Inflammatory cells were infiltrated in connective tissue. Simultaneously, cells in the blood vessel and the lymph vessel became quantitatively more and more distinct. The intestine structure was injured severely ( Figure 1C). For H12, the infection led to destructive damage to the intestine structures. The mucosal folds suffered from further damage; fragmentation were observed. The blood vessel was also severely damaged, and only a few blood vessels could be observed. Necrosis was seen in the mucosa, submucosa, and muscle layers of intestinal wall ( Figure 1D). The lamina propria was significantly thicker at 12 h post-challenge (Supplementary Table 1; Figure 1D).
In order to clarify the biological importance for time course, the qRT-PCR check for expression of key cytokines for enteritis was conducted. The results showed that E. tarda could 1) up-regulate mRNA levels of intestinal pro-inflammatory cytokines interleukin 1β (IL-1β), IL-6, IL-8, IL-16, and IL-17D, tumor necrosis factor α (TNF-α), and granulocyte colonystimulating factor (G-CSF); and 2) down-regulate the mRNA levels of anti-inflammatory cytokines IL-10 (Figure 2). The primers of qRT-PCR are shown in Supplementary Table 2.

Overview of circRNA Sequencing Data
As shown in Table 1, raw reads, clean reads, clean bases, and Q20, Q30, and GC (guanine and cytosine) contents for each library were identified. All libraries gave a good quality base value ≥ 12.25 Gb, Q20 ≥ 96.99%, Q30 ≥ 92.32%, and an error rate ≤ 0.02. Therefore, all libraries proved to be suitable for further study. These data were deposited in NCBI database with the BioProject number of PRJNA511138.
Clean reads from the 12 libraries were used to identify circRNAs. After a series of selection, 5,478 novel circRNAs were obtained and termed from novel_circ_0000001 to novel_circ_0005478 (Supplementary Table 3). There was no single circRNA reported in olive flounders previously, so all the identified circRNAs were novel. A size distribution analysis revealed that the length of circRNAs ranged from 150 to 71,793 bp, but most (74.82%) were ≤5,000 bp ( Figure 3A). Most of the circRNAs were exonic and intronic ( Figure 3B).
To identify circRNAs that potentially participated in E. tarda infection, their expression profiles were examined at 0, 2, 8, and 12 h post-infection. As shown in Supplementary Table 4, a total of 34, 20, and 8 differentially expressed circRNAs (DE-circRNAs) were observed at 2, 8, and 12 h relative to 0-h control, respectively. A heatmap showing 62 DE-circRNAs ( Figure 3C) and a Venn diagram revealing one circRNA (novel_circ_0005065) were differentially expressed at all three comparisons ( Figure 3D).

Overview of miRNA Sequencing Data
Meanwhile, to study the miRNA profile of P. olivaceus after E. tarda infection, four sRNA libraries (i.e., H0, H2, H8, and H12) were also constructed and sequenced. Altogether, 33,957,978, FIGURE 2 | Relative expression of interleukin 1β (IL-1β), IL-6, IL-8, IL-10, IL-16, IL-17D, tumor necrosis factor α (TNF-α), and granulocyte colony-stimulating factor (G-CSF) in posterior intestines of E. tarda-infected P. olivaceus. The data present the mean ± standard error (SE) derived from triplicate experiments. * represents the values that are significantly different (p < 0.05) from H0 time point.  ,067,562, 34,149,754, and 33,590,933 raw reads were acquired from H0, H2, H8, and H12, respectively. These data were deposited in NCBI database with the BioProject number of PRJNA510916. After the low-quality reads, adaptor sequences, and reads with sequences < 1 or >35 nt were filtered, all sRNAs were obtained ( Table 2). As shown in Figure 4A, the majority of the sRNAs from the four libraries ranged from 20 to 23 nt, and the peak distribution was for sequences that were 21 nt long. Furthermore, the sRNAs of all samples were compared with their reference genomes. The results showed that most of sRNAs (>89%) were able to be mapped onto the olive flounder genome ( Table 2). After Rfam and genome databases were searched, other non-coding RNAs (rRNA, tRNA, snRNA, and snoRNA) (Supplementary Table 5) and repeat sequences (Supplementary Table 6) were annotated. After a series of selections, a total of 303 miRNAs were identified, including 33 known miRNAs and 270 putative novel miRNAs ( Table 3). The expression levels of miRNAs were calculated based on the read count and were subsequently normalized to TPM (number of transcripts per million clean tags). Among the 303 miRNAs, approximately 22% miRNAs were in a higher level (TPM interval > 60). However, approximately 25% miRNAs showed a much lower expression level (0 < TPM interval < 0.1) (Supplementary Table 7). Specifically, the TPM of novel miRNAs was mostly in a lower level, and the known miRNAs were mostly in a higher level (Supplementary Table 8). While comparison was performed between different groups, the overall expression pattern of miRNAs among the four groups was highly consistent (Figure 4B). Compared with the miRNA expression levels of uninfected group, 39 miRNAs showed significantly differential expression (p < 0.05), including 26, 7, and 6 DE-miRNAs in H2 vs H0, H8 vs H0, and H12 vs H0 comparisons, respectively ( Figure 4C). Unsupervised hierarchical clustering revealed that all samples were clustered according to their respective groups, proving that the miRNA expression signatures were able to differentiate challenge groups from normal group (Figure 4D). Venn diagram revealed that some DE-miRNAs were differentially expressed at two or three comparisons ( Figure 4E).
In comparison with the H0 library, 2,100, 400, and 511 genes were identified as DEGs in H2, H8, and H12 libraries, respectively ( Figure 5A, Supplementary Table 9). As shown in the Venn diagram, some DEGs were differentially expressed at two or three comparisons ( Figure 5B). Considering that circRNAs could regulate transcription of their parental genes, we took the intersection of DE-circRNAs parental genes and DEGs. As shown in Supplementary  Table 9, 13 DEGs served as parental genes of 13 DE-circRNAs. To further explore the functions of the DEGs in response to E. tarda infection, GO and KEGG enrichment analyses were conducted. For GO analysis, the dominant functions in each of the three main categories were metabolic process (GO:0008152) in the biological process (BP) category, cell (GO:0005623) in the cellular component (CC) category, and structural constituent of ribosome (GO:0003735) in the molecular function (MF) category ( Figure 5C). In addition, the DEGs were aligned against the KEGG pathways database to identify pathways that were responsive to E. tarda infection. As shown in Figure 5D, the pathways of ribosome, proteasome, oxidative phosphorylation, and spliceosome were mostly activated.

Construction of the Potential circRNA-miRNA Network
As we all know, circRNAs serving as competing endogenous RNA (ceRNA) of miRNA could regulate the expression of corresponding genes. In order to construct circRNA-miRNA regulation network, the MiRanda software was used to predict the relationship between DE-circRNAs and DE-miRNAs. The circRNA-miRNA network contained 325 circRNA-miRNA pairs, including 51 circRNAs and 32 miRNAs (Supplementary Table 10). As shown in Figure 6A, some circRNAs were predicted to combine several miRNAs; for example, novel_ circ_0003296 could combine 20 miRNAs. Meanwhile, some miRNAs could also bind to several circRNAs; for example, novel_51 could link 30 circRNAs ( Figure 6B). Considering the importance of hub genes in a network, we employed an MCODE approach (Bader and Hogue, 2003) to screen hub genes from the protein-protein interaction (PPI) network. With the k-core = 2, one subnetwork with 11 nodes and 18 edges was identified, including four circRNAs (novel_ circ_0002248, novel_circ_0002267, novel_circ_0001856, and novel_circ_0000799) and seven miRNAs (novel_149,  novel_154, novel_171, novel_204, novel_272, pol-miR-144-5p, and pol-miR-182-5p) ( Figure 6C).

Construction of the Potential miRNA-mRNA Network
In order to construct miRNA-mRNA regulation network, the MiRanda software was used to predict the relationship between DE-miRNAs and DEGs. We predicted the potential target genes of the 39 DE-miRNAs using MiRanda software. Then, we took the intersection of potential target genes and DEGs. We obtained in total 3,873 possible miRNA-mRNA target pairs under three comparisons. As shown in Supplementary Table 11, all miRNAs had more than one intersected DEGs. Five significant miRNAs, novel_51 (degree = 510), novel_144 (degree = 417), novel_171 (degree = 294), pol-miR-144-3p (degree = 218), and novel_318 (degree = 207), had the most target genes. Moreover, many mRNAs were associated with more than one miRNA, such as FIGURE 5 | Analysis of mRNA sequencing data. Volcano plots for differential expression genes (DEGs) between H2, H8, H12, and H0, respectively. (B) Venn diagrams display the distribution of the 2,531 unique DEGs between H2, H8, H12, and H0, respectively. (C) Gene ontology (GO) analysis for all of DEGs between H2, H8, H12, and H0. X-axis represents the number of genes. The Y-axis on the left represents the GO term, and the Y-axis on the right represents GO type. The green column indicates the biological process, the red column the cellular component, and the gray column the molecular function. (D) Statistics of pathways enrichment of all DEGs between H2, H8, H12, and H0. Colors of the points refer to the q-value of the respective signaling pathway. Size of the point refers to the number of genes within each pathway. serpin H1 (gene_id, 109634651; transcript_id, XM_020095324.1), which was targeted by novel_265, novel_495, pol-miR-144-3p, novel_204, novel_17, novel_59, pol-miR-206-3p, novel_3, pollet-7a-5p, novel_127, novel_154, and pol-miR-144-3p. Construction of the Potential circRNA-miRNA-mRNA Network According to the differentially expressed results, circRNA-miRNA pairs and miRNA-mRNA pairs were predicted by MiRanda software. Then, circRNA-miRNA-mRNA network was generated using a combination of data obtained from circRNA-miRNA pairs and miRNA-mRNA pairs (Figure 7). This network contained 198 circRNA-miRNA pairs and 3,873 miRNA-mRNA pairs, including 44 circRNAs, 32 miRNAs, and 1,774 mRNAs (Supplementary Table  12). Among 198 circRNA-miRNA pairs, four pairs of circRNA-miRNA existed in multiple comparison groups; for example, novel_circ_0005065-novel_171 existed in all three comparisons, novel_circ_0003068-novel_51 and novel_circ_0003068-novel_144 both existed in H2 vs H0 and H12 vs H0 comparisons, and novel_ circ_0005065-pol-miR-144-5p existed in H2 vs H0 and H8 vs H0 comparisons. Among 3,873 miRNA-mRNA pairs, 178 miRNA-mRNA repeats existed in multiple comparison groups; for example, novel_171-109646742, novel_171-109646311, and novel_171-109644261 pairs existed in all of the three comparisons.
Next, GO and KEGG analyses were performed to evaluate the function of the DEGs in the network. GO analysis revealed that there were 107, 53, and 96 enriched GO terms with statistical significance (p < 0.05) in the biological process, cellular component, and molecular function categories, respectively.
GO analysis results suggested that some of the DEGs might play important biological roles during olive flounder against the E. tarda infection. As shown in Figure 8A, specific GO items were mainly involved in biological processes (e.g., metabolic process, organic substance metabolic process, and primary metabolic process), cell components (e.g., cell, cell part, and intracellular), and molecular function (e.g., catalytic activity, organic cyclic compound binding, and heterocyclic compound binding).
qRT-PCR Verification of Selected circRNAs, miRNAs, and mRNAs To further confirm the expression level of circRNAs obtained by Illumina sequencing, six DE-circRNAs (novel_circ_0001462, novel_circ_0002610, novel_circ_0002746, novel_circ_0003643, novel_circ_0003068, and novel_circ_0002248) were selected for qRT-PCR. Divergent primers were designed for each selected circRNA, and EF1α was used as internal control ( Supplementary  Table 13). Specifically, divergent primers could only amplify circular RNA forms, but not genomic DNA or linearized mRNAs, and Sanger sequencing further confirmed the amplified products to be circRNAs ( Figure 9A). Their relative expression level at different time points (H2, H8, and H12) was compared with that of H0. As shown in Figure 9B, most of the qRT-PCR results were consistent with those of Illumina sequencing. Therefore, the Sanger sequencing and qRT-PCR results confirmed the reliability and accuracy of the circRNA sequencing data.
To further confirm the expression level of miRNAs obtained by Illumina sequencing, eight novel_318,novel_171,novel_561,novel_154,novel_272,and novel_54) were selected for qRT-PCR. Different primers were designed for each selected miRNA, and 5S rRNA was used as internal control (Supplementary Table 14). As shown in Figure 10, there was similarity between the quantitative assay and high-throughput sequencing analysis of the eight miRNAs in terms of fold change and significance of differential expression. Although there were few differences in fold change of expression, the variation trend was identical.
To validate the expression patterns of the mRNAs, qRT-PCR was utilized to detect the expression of 10 randomly selected DEGs (XM_020102825, XM_020094518, XM_020094521, XM_ 020094517, XM_020091231, XM_020092535, XM_020113897, XM_020112506, XM_020093123, and XM_020106476). Different primers were designed for each selected mRNA, and EF1α was again used as internal control (Supplementary Table 15). As shown in Figure 11, all of 10 randomly selected DEGs showed a similar expression pattern between qRT-PCR and Illumina sequencing, although there were slight differences in the fold change.

DISCUSSION
Histopathological evaluation and cytokine expression analysis on the posterior intestine of P. olivaceus helped in understanding and determining the pathogenesis of E. tarda. This research described the morphological changes and cytokine expression changes of the posterior intestine from olive flounder after E. tarda infection at FIGURE 10 | Confirmation of miRNAs by qRT-PCR analysis. qRT-PCR analysis result (orange) was compared with data obtained from Illumina sequencing (crimson). The expression rates of qRT-PCR between H2, H8, H12, and H0 samples are shown by fold change. The data present the mean ± standard error (SE) derived from triplicate experiments. The relative expression levels by high-throughput sequencing analysis are represented by 2 log2(treatment/control) .
FIGURE 11 | Confirmation of mRNAs by qRT-PCR analysis. qRT-PCR analysis result (orange) was compared with data obtained from Illumina sequencing (crimson). The relative expression level between H2, H8, H12, and H0 samples are shown by fold change. The data present the mean ± standard error (SE) derived from triplicate experiments. The relative expression levels by high-throughput sequencing analysis are represented by 2 log2(treatment/control) . different time points. As the infection developed, the integrity of the intestinal mucosal structures showed pathological changes, such as cellular swelling, thickness of the lamina propria, shedding of the epithelial cells, and fragmentation of mucosal folds. Furthermore, structures of the posterior intestine, which were for immunity, especially the adaptive immunity, had seen an increased number of inflammatory cells (predominantly lymphocytes) and goblet cells at late infection (i.e., 8 and 12 h post E. tarda infection), which was consistent with qRT-PCR results of cytokine expression. Significantly thicker lamina propria, which was the structure for immunity, was observed at 8 and 12 h post-infection, an indication that the immune system actively participated at that particular time period. Pathogenesis study integrating morphological histology approach and next-generation sequencing approach would benefit a better understanding and elucidation of the intestinal immune response in the course of host-bacterial interaction.
More and more researches have proved that circRNAs, acting as miRNA sponges, counteract miRNA and eventually mediate expression of mRNAs (Hansen et al., 2013). Recently, the regulatory network of circRNA-miRNA-mRNA has been increasingly demonstrated in different kinds of diseases (Zheng et al., 2016;Chen et al., 2017a). In this research, we constructed a global circRNA-miRNA-mRNA network based on predicted circRNA-miRNA and miRNA-mRNA pairs in the pathogenesis of E. tarda in olive flounder. The integrated circRNA-miRNA-mRNA network consisted of 44 circRNAs, 32 miRNAs, and 1,774 mRNAs. GO and KEGG analyses were performed to evaluate the function of the DEGs in the network. KEGG analysis revealed that two important intestinal immune pathways (herpes simplex infection pathway and intestinal immune network for IgA production pathway) showed statistical significance between challenge and control groups.
Humans are the natural host of herpes simplex virus (HSV), of which HSV-1 resides in greater than 60% of the world's population and causes peripheral disease (Farooq and Deepak, 2012). HSV-1 and HSV-2 initiate the infection process in epithelial cells of mucosal surfaces. Binding glycoproteins gB and gC of HSV-1 with heparan sulfate proteoglycans from host cell surface allows attachment of the viral glycoproteins gB, gD, and gL to host cellular receptors, such as nectin1, herpes virus entry mediator, or 3-O-sulfated HS for membrane fusion and viral entry (Agelidis and Shukla, 2015;Menendez and Carr, 2017). In addition, nectin2 also functions as a receptor for HSV-2, although the binding to the gD is weak. Consistent with these conclusions, CHO cell line expressing hNectin2 was susceptible to HSV-2 infection (Fujimoto et al., 2016). Coincidentally, we also found nectin2 (mRNA id, XM_020093024.1; gene id, 109633274) in the current circRNA-miRNA-mRNA network, and it was down-regulated between H2 and H0 comparison. Similar to HSV, the challenged fish were immersed in the bacteria solution so that their mucosal surfaces (gill, skin, and gastrointestinal tract) became important sites of bacterial exposure and colonization. Therefore, we concluded that nectin2 of P. olivaceus may be also a receptor of E. tarda, and the down-regulation of nectin2 is an active regulation of the host for self-protection. Then, we found several circRNA-miRNA-RNA networks where nectin2 is located, including 1) novel_circ_0002248/novel_circ_0002610/ novel_circ_0003068/novel_circ_0004671-novel_144-nectin2 and 2) novel_circ_0000686/novel_circ_0002248/novel_ circ_0005065-novel_149-nectin2. In the following study, we focused on the function of nectin2 and its circRNA-miRNA-mRNA regulatory network.
In conclusion, by employing Illumina sequencing, bioinformatics, and qRT-PCR technologies, we constructed circRNA-miRNA-mRNA networks and found two important intestinal immune pathways (herpes simplex infection pathway and intestinal immune network for IgA production pathway). In addition, three critical DEGs (nectin2, MHC II α-chain, and MHC II β-chain) were identified, and their circRNA-miRNA-mRNA networks were also constructed and discussed. Our study provides a novel insight into the immune response for P. olivaceus to E. tarda infection from the circRNA-miRNA-mRNA view. Future research on the specific mechanism of action should be investigated in the pathology of E. tarda.

AUTHOR CONTRIBUTIONS
YX: constructed the circRNA-miRNA-mRNA network and wrote this paper; GJ: raised Paralichthys olivaceus and conducted