Next Article in Journal
Morphological and Micromorphological Description of the Larvae of Two Endemic Species of Duvalius (Coleoptera, Carabidae, Trechini)
Next Article in Special Issue
Comprehensive CircRNA Profiling and Selection of Key CircRNAs Reveal the Potential Regulatory Roles of CircRNAs throughout Ovarian Development and Maturation in Cynoglossus semilaevis
Previous Article in Journal
Ozone Gel in Chronic Periodontal Disease: A Randomized Clinical Trial on the Anti-Inflammatory Effects of Ozone Application
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Biology
Volume 10
Issue 7
10.3390/biology10070626
biology-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Revealing New Landscape of Turbot (Scophthalmus maximus) Spleen Infected with Aeromonas salmonicida through Immune Related circRNA-miRNA-mRNA Axis

by
Ting Xue
,
Yiping Liu
,
Min Cao
,
Mengyu Tian
,
Lu Zhang
,
Beibei Wang
,
Xiaoli Liu
and
Chao Li
*
School of Marine Science and Engineering, Qingdao Agricultural University, Qingdao 266109, China
*
Author to whom correspondence should be addressed.
Biology 2021, 10(7), 626; https://doi.org/10.3390/biology10070626
Submission received: 31 May 2021 / Revised: 30 June 2021 / Accepted: 2 July 2021 / Published: 6 July 2021
(This article belongs to the Special Issue The Application of Genetic and Genomic Biotechnology in Aquaculture)
Browse Figures
Review Reports Versions Notes

Abstract

:

Simple Summary

In this study, the expression of circRNAs, miRNAs, and mRNA in the immune organs spleen of turbot (Scophthalmus maximus) infected with Aeromonas salmonicida was analyzed by high-throughput sequencing, and circRNA-miRNA-mRNA network was constructed, so as to explore the function of non-coding RNA in the immune system of teleost. A total of 119, 140, and 510 differential expressed circRNAs, miRNAs, and mRNAs were identified in the infected groups compared with the uninfected group. The qRT-PCR verified the reliability and accuracy of the Illumina sequencing data. Fifteen triple networks of circRNA-miRNA-mRNA were presented in the form of “up (circRNA)-down (miRNA)-up (mRNA)” or “down-up-down”. Immune-related genes were also found in these networks. These results indicate that circRNAs and miRNAs may regulate the expression of immune-related genes through the circRNA-miRNA-mRNA regulatory network and thus participate in the immune response of turbot spleen after pathogen infection.

Abstract

Increasing evidence suggests that non-coding RNAs (ncRNA) play an important role in a variety of biological life processes by regulating gene expression at the transcriptional and post-transcriptional levels. Turbot (Scophthalmus maximus) has been threatened by various pathogens. In this study, the expression of circular RNAs (circRNAs), microRNAs (miRNAs), and mRNA in the immune organs spleen of turbot infected with Aeromonas salmonicida was analyzed by high-throughput sequencing, and a circRNA-miRNA-mRNA network was constructed, so as to explore the function of non-coding RNA in the immune system of teleost. Illumina sequencing was performed on the uninfected group and infected group. A total of 119 differential expressed circRNAs (DE-circRNAs), 140 DE-miRNAs, and 510 DE-mRNAs were identified in the four infected groups compared with the uninfected group. Most DE-mRNAs and the target genes of DE-ncRNAs were involved in immune-related pathways. The quantitative real-time PCR (qRT-PCR) results verified the reliability and accuracy of the high-throughput sequencing data. Ninety-six differentially expressed circRNA-miRNA-mRNA regulatory networks were finally constructed. Among them, 15 circRNA-miRNA-mRNA were presented in the form of “up (circRNA)-down (miRNA)-up (mRNA)” or “down-up-down”. Immune-related genes gap junction CX32.2, cell adhesion molecule 3, and CC chemokine were also found in these networks. These results indicate that ncRNA may regulate the expression of immune-related genes through the circRNA-miRNA-mRNA regulatory network and thus participate in the immune response of turbot spleen after pathogen infection.

1. Introduction

Turbot (Scophthalmus maximus, S. maximus) is a kind of cold-water benthic fish that originated from Europe, which is an important economic species in many countries. However, in recent years, the aquaculture of turbot was frequently interfered with by bacterial diseases, especially Aeromonas salmonicida (A. salmonicida), Edwardsiella tarda (E. tarda) and Vibrio anguillarum (V. anguillarum), which seriously restricted the development of turbot aquaculture [1]. Therefore, it is very urgent to study the immune system and disease mechanism of turbot.
Gene regulatory network is a complex dynamic network system, which is the overall performance of the interaction between various factors. In recent years, with the development of next-generation sequencing technology and the integration analysis and application of multi-omics data, a variety of non-coding RNAs (ncRNAs), such as long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and circular RNAs (circRNAs), have been confirmed to be involved in the complex biological process of fish resistance to bacterial immunity.
CircRNAs are a class of single-stranded closed RNA with 3′ and 5′ ends linked to each other [2]. The circRNA was first discovered in 1979 using electron microscopy in eukaryotic cells [3]. CircRNAs were initially considered to be the abnormally spliced transcripts [4,5]. However, thousands of circRNAs have been found in various organisms [6,7,8,9]. With abundant miRNA binding sites, circRNAs can act as miRNA sponges by binding and adsorbing with miRNA [10]. MicroRNA is a short (~22 nt), endogenous, evolutionarily conserved non-coding RNA that can down-regulate the expression of a target gene by binding to the 3′ UTR of mature mRNA [11]. Therefore, the adsorption of circRNAs to miRNAs prevents the binding of miRNA to its target mRNAs, thus indirectly regulating the expression of miRNA target genes [12].
Recent studies showed that circRNAs are involved in a variety of biological processes, such as cell regulation [13], signaling pathway regulation [14], tissue development [15], and so on. Nevertheless, other numerous researches revealed that circRNAs also plays an important role in the occurrence and treatment of diseases such as tumors, tissue injury, and depression [16,17,18]. Furthermore, some of these studies revealed that circRNAs are also involved in innate immune responses [19,20]. CircRNA (102911) is involved in the proliferation and apoptosis of CD3(+)/CD4(+) T lymphocytes through adsorbing miR-129-5p [21]. A high level of circ_0020710 could upregulate the CXCL12 expression via sponging miR-370-3p and promotes melanoma growth in vivo [22]. Circ-RelL1 could regulate inflammatory response in ox-LDL-induced endothelial cells through miR-6873-3p/MyD88/NF-κB axis [23].
The identification of circRNAs in fish has also gradually increased, including those participating in the process of teleost resistance to pathogen infection [24,25,26]. In grass carp, 41 differential expressed circRNAs (DE-circRNAs) were identified after GCRV (Grass Carp Reovirus) infection, and these circRNAs may be related to the hemorrhagic symptoms [24]. circPIKfyve activated the NF-κB/IRF3 pathway by sponging miR-21-3p, thereby enhancing the innate antiviral response of teleost [25]. As a competitive endogenous RNA (ceRNA) of miR-21-3p, circrasgF1b alleviated the inhibitory effect of miR-21-3p on its target, MITA, and thus enhanced the innate anti- SCRV (Sininiperca Chuatsi Rhabdovirus) response [26].
In addition to the functions of circRNAs through translation and pseudogene derivation influencing the splicing of its linear cognates, regulating transcription and splicing, circRNAs mainly play their role through binding with miRNA [27]. CircRNAs affect the expression of mRNAs by combining with miRNAs, forming circRNA-miRNA-mRNA regulatory network, and then affecting downstream genes or pathways, thus playing a role in a variety of biological processes, such as autoimmune disease [28], stress response [29], cancer [30], and cardiovascular disease [31]. Emerging studies have shown that circRNA-miRNA-mRNA regulatory networks have also been identified in pathogen-infected teleost [9,32,33].
In this study, we examined the expression of circRNA, miRNA, and mRNA in the spleen of turbot infected with A. salmonicida at different time points by high-throughput sequencing. The differentially expressed circRNAs, miRNAs, and mRNAs were screened and the regulatory networks of circRNA-miRNA, miRNA-mRNA, and circRNA-miRNA-mRNA were constructed. Our aim is not only to broaden our understanding of turbot immune regulatory network, but also to provide new ideas to enhance turbot resistance to pathogenic bacteria invasion.

2. Materials and Methods

2.1. Turbot Treatment and Sample Collection

Healthy turbot (100 ± 10 g) were purchased from Huanghai Aquaculture Company (Haiyang, Shandong, China). The experimental protocols were approved by the Committee on the Ethics of Animal Experiments of Qingdao Agricultural University IACUC (Institutional Animal Care and Use Committee). The experiment was carried out after 2 weeks of temporary breeding in the laboratory. The feeding was stopped 3 days before the experiment. According to Coscelli et al., A. salmonicida was injected into turbot by intraperitoneal injection [34]. The A. salmonicida was first reactivated from the infected turbot and then cultured by LB culture medium. The bacteria were resuspended in PBS at a concentration of 1 × 107 CFU and inoculated intraperitoneally at 0.1 mL per tail. The turbot spleen tissue samples from 15 fish injected with bacteria for 6 h (ASS6), 12 h (ASS12), 24 h (ASS24), and 96 h (ASS96) (five fish per pool) were sampled. All fish were anesthetized with 0.01% ethyl-m-aminobenzoate until breathing ceased before the samples were collected. The PBS-injected turbot was taken as control (ASS0). The extracted tissue samples were frozen in liquid nitrogen and subsequently analyzed by high-throughput sequencing.

2.2. Total RNA Isolation, Library Construction, and Sequencing

Total RNA was extracted from the collected samples with TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s manual. The purity, concentration, and quality of the total RNA were measured using NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). The integrity of RNA was accurately detected by Agilent 2100 bioanalyzer system (Agilent Technologies, CA, USA).
The preparation steps of mRNA and circRNA libraries followed manufacturer’s recommendations of NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA). First, after the RNA was tested and qualified, 3 μg RNA was prepared for each sample and ribosomal RNA (rRNA) was removed from the total RNA using the Epicenter Ribo-Zero rRNA Removal Kit (Illumina, WI, USA). mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Different from mRNA library, the construction of circRNA library required the addition of 40U RNase R to the rRNA removal system and incubation at 37 °C for 3 h to remove the linear RNA. Next, the treated RNA samples were fragmented to 250–300 bp.
The first strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H−). The DNA Polymerase I was then used for the synthesis of second strand cDNA. The library fragments were purified and selected using AMPure XP system (Beckman Coulter, CA, USA). The size-selected, adaptor-ligated cDNA was degraded by USER enzyme (NEB, MA, USA). Then PCR was performed to obtain the chain specific cDNA library and the quality of the purified library was assessed on the Agilent Bioanalyzer 2100 system.
The NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, MA, USA) was used to construct the miRNA library. In brief, approximately 3 μg of total RNA was used as the initial sample, based on the special structure of sRNA with hydroxyl group at the 3′ end and phosphate group at the 5′ end; sRNA was specificity ligated with adaptors at both ends, and then the cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H−). Sequencing adaptors were ligated to double-stranded cDNA using PCR and the final miRNA library contains 140–160 bp fragments that were separated and recovered from an 8% polyacrylamide gel. The library quality was detected on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. At last, the mRNA, circRNA, and sRNA library were clustered and then sequenced on an Illumina Hiseq 2500/2000 platform.

2.3. Data Analysis

2.3.1. Quality Control

Clean data was generated by removing reads with an adapter, reads with ploy-N, and low quality reads from raw data. Simultaneously, Q20, Q30, GC-content, and sequence duplication level of the clean data were calculated. All the downstream analyses were based on the high-quality clean data. Then, the clean data was aligned to the assembled genome of S. maximus using STAR (v2.5.1b) or Bowtie [35].

2.3.2. RNA Identification and Annotation

In consideration of the high false positives in circRNA identification, two software (find_circ and CIRI2) [36,37] were used for the identification of circRNA and the results from these two software were intersected.
The mapped small RNAs were first identified in the miRBase 20.0 database to find known miRNAs. Repeatmasker and RFAM software were used to remove the possible rRNA, tRNA, small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA). Novel miRNA was analyzed using miREvo and mirdeep2 software [38,39], which depended on the hairpin structure of miRNA precursors. A total rRNA ratio of less than 40% was considered as a sample quality indicator. The base bias of seed region was calculated using Perl scripts. miFam.dat (http://www.mirbase.org/ftp.shtml, accessed on 3 September 2019) and Rfam (http://rfam.sanger.ac.uk/search/, accessed on 5 September 2019) database were used to look for families of known miRNA and novel miRNA precursors, respectively. miRanda [40] was used to predict the miRNA binding sites to target genes.

2.3.3. Quantification of Gene Expression Levels

The expression levels of mRNA, circRNA, and miRNA were estimated using TPM (number of transcripts per million clean tags) [41]. DESeq2 R package (1.10.1) [42] was used for the different expression analysis of mRNA, circRNA, and miRNA under the condition of threshold p-value < 0.05.

2.3.4. GO and KEGG Analyses

GOseq R package [43] and KOBAS (http://www.genome.jp/kegg/, accessed on 15 September 2019) were used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to further explore the function of differential expressed mRNA (DE-mRNA), DE-circRNA, or DE-miRNA in S. maximus.

2.3.5. Co-Expression Network Analysis

A plurality of miRNA-mRNA and circRNA-miRNA interaction pairs were generated by target prediction of that miRNA to the circRNA and the mRNA. According to the miRNA that targeted circRNA and mRNA simultaneously, circRNA-miRNA-mRNA networks were constructed by Cytoscape 3.4.0 software [44]. In the established circRNA-miRNA-mRNA network, the expression trend of “up (circRNA)-down (miRNA)-up (mRNA)” or “down-up-down” was retained. The immune-related circRNA-miRNA-mRNA networks will be selected as candidates for future study.

2.3.6. Validation the Expression of RNAs with Quantitative Real-Time PCR (qRT-PCR)

To confirm the reliability of sequencing data, four DE-circRNAs, four DE-miRNAs, and eight DE-mRNAs were randomly selected to verify their expressions by qRT-PCR. Oligo7 software was used for the specific primers design. The design of circRNA divergent primers spanned the circRNA backsplice junction. The forward primers of miRNA were designed according to the instructions of miRcute miRNA isolation kit (Tiangen Biotech, China). All the primers used in this study are listed in Supplementary Table S1. The total RNA extracted forms of A. salmonicida-infected and control S. maximus spleen were reverse-transcribed into cDNA by using the PrimeScript™ RT reagent Kit and 6 mers primers (Takara, Japan) for circRNAs and mRNAs detection. The total RNA extracted forms of A. salmonicida-infected and control S. maximus spleen were reverse-transcribed into cDNA by using the miRcute miRNA First-Stand cDFNA Synthesis Kit (Tiangen Biotech) for miRNA detection. qRT-PCR was performed on a CFX 96 real-time PCR system (Bio-Rad, Hercules, CA, USA) with SYBR® Premix Ex TaqTM II (TliRNaseH Plus) (Takara, Japan).
The 18S rRNA was used as an internal control of circRNAs and mRNAs and U6 was used as an internal control of miRNA. Relative expression level was calculated using the 2−ΔΔCt method [45]. All experiments were performed in triplicate and the data were shown as means ± SE. Correlation analysis of differential expression of RNAs in qRT-PCR and high-throughput sequencing was performed using Spearman’s nonparametric correlation (SPSS20.0).

3. Results

3.1. Statistical Analysis of circRNAs Data

In rRNA-depleted non-infected (ASS0) and infected libraries from different time points (6, 12, 24, and 96 h), average 87.30, 86.59, 96.44, 94.79, and 94.98 million clean reads were identified, respectively (Table 1). The results showed that the quality of the clean reads of all libraries was high, with Q20 > 95% and Q30 > 90%, and most (>85%) clean reads were mapped to the reference genome (Table 1).
CircRNAs were identified in all samples, and they were widely distributed on different chromosomes (Figure 1A and Supplementary Table S2). Of all the identified circRNAs, 78.56% were located in exons, 5.92% in introns, and 15.52% in intergene regions (Figure 1B). The size of most circRNAs ranged from 200 to 400 nt, and 82.27% of them had a predicted spliced length of <400 nt (Figure 1C). As shown in Figure 1D and Supplementary Table S3, a total of 28, 34, 42, and 47 DE-circRNAs were identified at 6 h, 12 h, 24 h, and 96 h relative to control, respectively. Of these, 59 circRNAs were up-regulated and 60 were down-regulated. Further, we performed a hierarchical cluster analysis and Venn diagram of all DE-circRNAs (Figure 1E,F). The results showed that the expression patterns of uninfected and A. salmonicida-infected samples were classified into different clusters (Figure 1E). Nevertheless, the types of DE-circRNA were changed with the increase of infection time. One circRNA (novel_circ_0005673) was shown differentially expressed at all four comparisons in the Venn diagram (Figure 1F).
Go and KEGG were performed to elucidate the biological function of circRNAs in S. maximus after A. salmonicida infection. The GO terms of all DE-circRNAs from four time points post-infection were classified into three types, biological processes, molecular functions, and cellular components (Figure 1G). The biological processes include cellular process, metabolic process, response to stimulus, signal transduction, cell communication, and cytokinesis. Protein binding, molecular function regulator, identical protein binding, and transmembrane receptor protein were the important functional molecules in A. salmonicida infection. Extracellular region, cell periphery, cell cortex, envelope, and outer membrane were functionally enriched in the cellular component. KEGG pathways analysis showed that the DE-circRNAs were involved in multiple biological processes, such as TGF-β signaling pathway, melanogenesis, tight junction, and adherens junction (Figure 1H).

3.2. Statistical Analysis of miRNAs Data

After removing the adaptor-contained and low-quality raw reads, approximately 44.15, 48.26, 40.76, 46.50, and 51.29 million clean reads were identified in 15 sRNA libraries (ASS0, ASS6, ASS12, ASS24, ASS96), respectively (Table 2). In order to analyze the distribution of small RNA on the genomes, filtered sRNA was compared with reference genomes. The results showed that 93.20%, 93.08%, 92.76%, 85.98%, and 94.72% small RNA reads of the five samples were successfully mapped onto the S. maximus genome (Table 2). Clean reads of each sample with a length range of 18~35 nt were filtered for subsequent analysis. From the statistical diagram of the length distribution of sRNA, the length of sRNA from the five groups varieties ranged from 21 to 23 nt (Figure 2A). rRNA, tRNA, snRNA, and snoRNA were annotated with Rfam and removed from sRNA.
The known miRNAs were analyzed by the miRBase database and novel miRNAs were predicted through the signature hairpin structure of miRNA precursors, respectively. A total of 910 miRNAs were identified, including 252 known miRNAs and 668 novel miRNAs (Table 3). The expression levels of known and novel miRNAs in each sample were counted and normalized by TPM. The overall expression pattern of miRNAs among the five groups was exhibited in Figure 2B and was highly consistent. Among the 910 miRNAs, 140 miRNAs showed significant changes in expression compared with the control group. There were 44, 63, 39, and 96 miRNAs differentially expressed in ASS6, ASS12, ASS24, and ASS96 compared with ASS0, respectively (Figure 2C). Hierarchical cluster analysis revealed that the DE-miRNAs formed obvious clustering, and there were significant differences between the uninfected and A. salmonicida-infected group (Figure 2D). As can be seen from the Venn diagram, one part of the miRNAs was differentially expressed at different time points of infection, and the other parts of the miRNAs were differentially expressed only at certain time points after infection (Figure 2E). More miRNAs were induced at 96 h after infection (Figure 2E).
The function of these DE-miRNAs was predicted through GO and KEGG analysis of their target genes. The GO terms of all target genes were also classified into three types. Most of the target genes involved biological processes, such as response to stimulus, signal transduction, cell communication, response to stimulus, and metabolic process. The majority of genes were important components of nucleoplasm, organelle inner membrane, mitochondrial inner membrane, and microtubule organizing center. These target genes mainly had a protein binding function, receptor activity, tumor necrosis factor receptor binding, and MHC protein binding function (Figure 2F). KEGG analysis of differentially expressed miRNAs showed that DE-miRNAs were annotated to wnt signaling pathway, TGF-β signaling pathway, MAPK signaling pathway, and adherens junction (Figure 2G).

3.3. Statistical Analysis of mRNAs Data

Through comparison, splicing, and screening, the expression levels of all transcripts were quantified, and 243, 187, 126, and 232 differential expressed genes were identified in four infected groups (ASS6, ASS12, ASS24, ASS96) compared with control group, respectively (Figure 3A). With the increase of infection time, the types of differentially expressed genes also changed (Figure 3B). Among them, 510 genes were differentially expressed at two or more infected time points, and 20 genes were differentially expressed at all four infected time points (Figure 3C).
Go analysis was performed to analyze the biological function of DE-mRNAs. The results showed that as the components of the cell membrane, extracellular region, or extrinsic component of membrane, these DE-mRNAs were mainly involved in a variety of biological processes, such as biological regulation, response to stimulus, cell communication, signal transduction, and transport. In these biological processes, they exhibited important functions in transporter activity, transmembrane transporter activity, signal transducer activity, receptor activity, and so on (Figure 3D). KEGG analysis showed that DE-mRNAs were annotated to various metabolic pathways, including glycolysis/gluconeogenesis, purine metabolism, carbon metabolism, RNA degradation, and biosynthesis of amino acids (Figure 3E).

3.4. Construction of the circRNA-miRNA Regulatory Network

CircRNAs regulate the expression of mRNAs by acting as miRNA sponges to antagonize the inhibition of miRNA on their target genes. The close relation of miRNA to the occurrence of disease by regulating mRNA expression suggests that circRNA also plays an important regulatory role in the occurrence of disease. Therefore, circRNA-miRNA regulatory networks were constructed to analyze the role of circRNAs as miRNA sponges in the process of A. salmonicida infection. A total of 1051 circRNA-miRNA interactions including 116 DE-circRNAs and 290 miRNAs were identified in the four infected groups (Supplementary Table S4).
Analysis showed that a single circRNA could adsorb more than one miRNA, while a single miRNA could bind to more than one circRNA. For example, novel_circ_0002735 could adsorb 28 miRNAs, and dre-miR-1306 could bind to 18 circRNAs. These miRNAs targeted a variety of immune-related genes, such as dre-miR-15b-5p targeted interleukin-20, dre-miR-138-2-3p targeted interleukin-8, dre-miR-138-5p targeted interferon regulatory factor-4, dre-miR-133a-3p and dre-miR-133a-3p targeted C-X-C motif chemokine 14, and dre-miR-133b-3p targeted tumor necrosis factor β (Supplementary Table S4). This finding suggested that these circRNA-miRNA pairs act as important regulatory mediators during the infection of A. salmonicida.

3.5. Construction of the miRNA-mRNA Regulatory Network

MicroRNAs usually downregulate gene expression in a variety of manners through partially complementary to one or more mRNA molecules. The miRNA-mRNA regulation networks were constructed to analyze the main functions of miRNA-mRNA pairs in S. maximus infected with A. salmonicida. In four infected groups, a total of 9605 mRNAs that were identified could be targeted by 140 DE-miRNAs (Supplementary Table S5). The absolute majority of miRNAs could target more than one mRNA, while only 10 miRNAs specifically targeted a single mRNA, such as dre-miR-126a-3p, dre-miR-203a-3p, dre-miR-142a-5p, dre-miR-1, and novel_1136. Similarly, the same one mRNA could also be bound by multiple or single miRNAs. The results showed that the mRNAs targeted by these miRNAs contained a variety of immune-related genes, including tight junction protein, interferon regulatory factor, interleukin, complement component, tumor necrosis factor receptor, integrin, C-X-C motif chemokine, gap junction, NF-κB inhibitor, and so on. C-X-C motif chemokine 14 could be targeted by dre-miR-133a-3p, dre-miR-133b-3p, and dre-miR-210-3p; gap junction α-5 could be targeted by novel_513; NF-κB inhibitor δ could be targeted by dre-miR-107b; complement C1q could be targeted by novel_753 and novel_880; interleukin-17 receptor B, tight junction protein ZO-1, and NF-κB inhibitor ζ all could be targeted by novel_324. The above results suggested that miRNA-mRNA pairs also play an important role in the S. maximus spleen infected with A. salmonicida.

3.6. Construction of the circRNA-miRNA-mRNA Regulatory Network

Gene regulatory network is a complex dynamic network system, which is the overall performance of the interaction between various factors, including circRNA and miRNA. A large number of circRNA-miRNA pairs and miRNA-mRNA pairs were identified in the spleen of S. maximus infected with A. salmonicida. The circRNA-miRNA-mRNA regulatory network has been verified to play an important role in various bacterial infected fish tissues. Therefore, the triple networks of circRNA-miRNA-mRNA were constructed for exploring underlying pathogenesis of A. salmonicida infection by combining the data from circRNA-miRNA interactions and miRNA-mRNA interactions. After a lot of in-depth analysis, 96 circRNA-miRNA-mRNA regulatory networks (name as DE-circRNA-miRNA-mRNA) consisting of 24 DE-circRNAs, 15 overlap DE-miRNAs, and 44 target DE-mRNAs were finally obtained (Figure 4A and Supplementary Table S6). Among all the DE-circRNA-miRNA-mRNA regulatory networks, 15 of them showed significant differences in the expression of all three RNAs in the form of “up (circRNA)-down (miRNA)-up (mRNA)” or “down-up-down” after being infected with A. salmonicida, six were expressed in the form of “up-down-up”, and nine were expressed in the form of “down-up-down” (Figure 4B and Supplementary Table S6).
Six hours post-infection, the downregulated circRNAs novel_circ_0006815 interacted with dre-miR-138-5p and affected the expression of two target genes of dre-miR-138-5p, ras-related protein rab-6A, and desmin. The same situation also occurred at 12 h after infection. The novel_circ_0001612 and novel_circ_0000208 that interacted with upregulated dre-miR-138-5p were downregulated, and the downregulation of dre-miR-138-5p targeted keratin type I cytoskeletal 13 (krt13) and phosphatase-related protein type 1 (ppr1) was also observed. Acts as novel_324 adsorbate, the low expression of novel_circ_0002055 can hardly weaken the inhibition of novel_324 on transcription factors mafB, ataxin-2, SH3 and cysteine-rich domain-containing protein 2 (stac2), and cell adhesion molecule 3 (cam3) at 24 h post-infection. Furthermore, the novel_circ_0001830, novel_circ_0000627, novel_circ_0001881, and novel_circ_0003494 downregulated the expression of dre-miR-193a-3p, novel_187, dre-miR-210-5p, novel_685, and then the miRNA-targeted genes, including gap junction Cx32.2, glucagon-1, and CC chemokine (ccl19), were upregulated in the 96 h post-infection. The detailed information of circRNA-miRNA-mRNA triple networks of different time points is shown in Supplementary Table S6.

3.7. qRT-PCR Verification of Selected circRNAs, miRNAs, and mRNAs

qRT-PCR was used to validate the expression profiles of circRNAs, miRNAs, and mRNAs obtained by Illumina sequencing. The expression profiles of four detected DE-circRNAs (novel_circ_0000498, novel_circ_0004184, novel_circ_0002683, and novel_circ_0003361), four detected DE-miRNAs (dre-miR-233, dre-miR-200a-3p, dre-miR-200b-3p, and dre-miR-122), and eight detected DE-mRNAs (gene14346, gene658, gene9492, gene1381, gene16286, gene16422, gene19654, and gene9752) in A. salmonicida-infected S. maximus spleen were randomly selected and measured (Figure 5, Figure 6 and Figure 7).
As shown in Figure 5, most of the qRT-PCR results of four selected circRNAs were consistent with those of Illumina sequencing (correlation coefficient rs ≥ 0.9, except for novel_circ_0004184 = 0.3). As shown in Figure 6, there was similarity between the quantitative assay and high-throughput sequencing analysis of the four miRNAs in terms of fold change and significance of differential expression (correlation coefficient rs ≥ 0.6). Although there were few differences in fold change of expression, the variation trend was identical. As shown in Figure 7, all eight randomly selected DE-mRNAs showed a similar expression pattern between qRT-PCR and Illumina sequencing (correlation coefficient rs ≥ 0.66), although there were slight differences in the fold change. The results confirmed the reliability and accuracy of the high-throughput sequencing data.
For example, both quantitative assay and high-throughput sequencing analysis showed that three immune-related genes (gene16286 (CD83), gene16422 (tumor necrosis factor α -induced protein 2), and gene9752 (CC chemokine)) were upregulated in S. maximus spleen when infected with A. salmonicida (Figure 7). Novel_circ_0004184 showed downregulated expression and the other three detected circRNAs showed upregulated expression in all time points post-infection in response to A. salmonicida infection (Figure 5). The similar expression trends of miRNAs were also found by quantitative assay and high-throughput sequencing analysis (Figure 6). The results showed that the expression trends of most genes in qRT-PCR were in agreement with the Illumina sequencing data.

4. Discussion

The spleen is the primary immune tissue of teleost except kidney, thymus, and mucosa-associated lymphoid tissues. The spleen not only contains immune T cells and immune B cells that can destroy or neutralize antigens to prevent further manifestations during the disease process, but also plays an important role in antigen presentation and adaptive immune response [46,47]. Therefore, expression profiling of the spleen of S. maximus infected with A. salmonicida will be helpful to understand the mechanism of disease resistance and disease defense of S. maximus.
The analysis of the whole transcriptome results showed that we have obtained high quality and credible expression profile data of circRNA, miRNA, and mRNA. qRT-PCR results of randomly selected circRNAs, miRNAs, and mRNAs also confirmed that.
Non-coding RNAs are composed of a variety of RNAs that do not encode proteins and have long been known as “junk DNA”. However, many studies have shown that ncRNAs are functional RNA molecules that regulate gene expression at the transcriptional and post-transcriptional levels [48]. According to their length, ncRNAs can be divided into two categories according to the length: >200 nt and <200 nt. The ncRNAs of >200 nt are mainly long non-coding RNA, and those of <200 nt include microRNA, short interfering RNA (siRNA), and piwi-interacting RNA (piRNA) [49]. Different from the linear ncRNAs mentioned above, circRNA is a covalently closed loop of ncRNAs with large length variations [2]. These ncRNAs perform their function in various physiological and pathological conditions, including embryonic development, tissue differentiation, tissue homeostasis, disease development, and disease suppression [15,50,51,52].
These ncRNAs participate in life processes by forming regulatory networks, such as circRNA-miRNA-mRNA regulatory networks. MiR-195a bound to circRNA_28313 and CSF1 to form circRNA-miRNA-mRNA network. As a result, circRNA_28313 released the inhibitory effect of mir-195a on CSF1 and regulated osteoclast differentiation of BMM (bone marrow monocyte/macrophage) cells [53]. CircRNA_09505 could promote AKT1 expression in macrophages via IκBα/NF-κB signaling pathway by acting as a miR-6089 sponge, thereby affecting the inflammatory response involved in macrophages [54]. circRNAs showed significant regulation in Type 1 Diabetes Mellitus. Hsa_circ_0060450 could release the target gene (SHP2 and PTPN11) of miR-199a-5p and further suppressed the JAK-STAT signaling pathway triggered by type I interferon (IFN-I) to inhibit macrophage-mediated inflammation [55]. In a word, circRNA-miRNA-mRNA regulatory networks existed and played a pivotal role in a variety of physiological and pathological processes. However, there were relatively few studies on the regulatory role of circRNA and miRNA in teleost fish after pathogen infection. Therefore, we systematically analyzed the circRNA, miRNA, and mRNA expression profiles of S. maximus spleen infected with A. salmonicida.
A total of 119 DE-circRNAs, 140 DE-miRNAs, and 510 DE-mRNAs were identified at four infection time points in S. maximus exposed to A. salmonicida infection. It has been found that these differential RNAs constituted multiple circRNA-miRNA-mRNA regulatory networks.
Similarly, 87 DE-circRNAs were identified by expression profile analysis in the intestine of Sebastes schlegelii after pathogen infection. The circRNA-miRNA-mRNA regulation network was found involved in the regulation of the NF-kB signaling pathway and the expression of interleukin 1 receptor [32]. Novel_circ_0004195, serving as ceRNA of novel_530, could activate the downstream NF-kB signaling pathway. It was demonstrated that circRNAs play an important antibacterial role in the intestinal tissues of Sebastes schlegelii challenged with E. tarda by forming the circRNA-miRNA-mRNA network.
In our study, circRNA was also found to regulate the expression of immune genes acting as miRNA sponges. Just as found in the constructed circRNA-miRNA and miRNA-mRNA regulatory networks, immune-related genes (gap junction Cx32.2, cell adhesion molecule 3, and CC chemokine) were also found in DE-circRNA-miRNA-mRNA regulatory networks (Figure 4B). These results indicated that circRNAs and miRNAs might play an important role in pathogenic bacteria-infected S. maximus in the form of circRNA-miRNA-mRNA network.
Connexin32 (Cx32) is a member of the gap junction protein family, which serves a useful function in the mammalian immune system. There is research suggesting that the expression of LPS-induced pro-inflammatory cytokines IL-8 and TNF-α was significantly reduced by Cx32 in Japanese flounder (Paralichthys olivaceus), indicating the important immunological functions of Cx32 [56]. The expression of dre-miR-193a-3p was decreased in the spleen of S. maximus infected with A. salmonicida for 96 h, and novel_circ_0001830 and novel_circ_0000627, which interacted with dre-miR-193a-3p, were highly expressed. Meanwhile, the expression of connexin32, which was predicted and targeted by dre-miR-193p, was also increased. Therefore, we speculated that novel_circ_0001830 and novel_circ_0000627 could indirectly regulate connexin32 to resist the invasion of pathogenic bacteria.
Chemokines are a family of cytokines that regulate the migration and immune response of leukocytes [57,58]. Chemokine (CC motif) ligand 19 (CCL19) is a cytokine that belongs to the CC chemokine subfamily. It is a key regulator of T cell activation, immune tolerance, and inflammatory response [59]. CCL19 could be targeted and inhibited by miR-325-3p to attenuate renal inflammation in diabetic nephropathy mice [60]. The expression level of ccl19 in the spleen of S. maximus infected with A. salmonicida was increased. ccl19 was reported to promote inflammation by activating Mo/Mφ and lymphocytes in pathogen-infected teleost [61]. Recombinant ccl19 could stimulate the expression of cytokines in the mucosal system and the proliferation of leukocytes in head kidney [62]. In S. maximus, ccl19 was also reported to induce leukocyte trafficking and promote anti-viral and anti-bacterial defense [63].
In our study, the expression of ccl19 was up-regulated in the spleen of S. maximus infected with A. salmonicida for 96 h, and the novel_685 that targeted ccl19 was down-regulated, while the predicted sponge of novel_685 and novel_circ_0001881 was up-regulated at the same time. ccl19, novel_685, and novel_circ_0001881 formed a circRNA-miRNA-mRNA regulatory network. This indicated that circRNAs and miRNAs could participate in the anti-pathogen immune response in S. maximus infected with A. salmonicida through the formation of circRNA-miRNA-mRNA regulatory network.

5. Conclusions

In conclusion, we analyzed the expression profile of Scophthalmus maximus spleen infected with A. salmonicida and constructed the circRNA-miRNA-mRNA regulatory networks. The establishment of circRNA-miRNA and miRNA-mRNA regulatory networks provided evidence that circRNA could participate in the immune regulation of fish through the circRNA-miRNA-mRNA regulatory network. It also provided a new insight into the immune response pathway of fish after pathogen infection. These immune-related circRNA-miRNA-mRNA regulatory networks still need to be further validated in order to understand the immune mechanism of S. maximus against A. salmonicida from a broader perspective.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/biology10070626/s1, Table S1: Primers used in this study, Table S2: circRNAs distributed on chromosomes, Table S3: DE-circRNA, Table S4: circRNA-miRNA, Table S5: miRNA-mRNA, Table S6: circRNA-miRNA-mRNA.

Author Contributions

T.X. has contributed to data curation, formal analysis, project administration, visualization, writing—original draft, and writing—review & editing; C.L. has contributed to funding acquisition, methodology, project administration, validation, visualization, and writing—review & editing; M.C. has contributed to investigation, resources, and visualization; Y.L. has contributed to formal analysis, investigation, and visualization; M.T., L.Z., B.W. and X.L. have contributed to investigation, resources, and visualization. All authors have read and agreed to the published version of the manuscript.

Funding

This study was supported by the Joint Funds of the National Natural Science Foundation of China and Shandong Province (No. U1706205), National Science Foundation of China (Grant No.: 31902403), Young Experts of Taishan Scholars (No.: tsqn201909130); Science and Technology Support Plan for Youth Innovation of Colleges and Universities in Shandong Province (No.: 2019KJF003); the advanced Talents Foundation of QAU grant (Grant No.: 663-1120029), the “First class fishery discipline” program in Shandong Province, and Breeding Plan of Shandong Provincial Qingchuang Research Team (2019), China.

Institutional Review Board Statement

The animal study was reviewed and approved by Qingdao Agricultural University.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

The data presented in this study are available in this article or Supplementary Materials.

Conflicts of Interest

We also affirm that there are no financial or other relations that could lead to a conflict of interest.

References

  1. Yang, N.; Wang, B.; Yu, Z.; Liu, X.; Fu, Q.; Cao, M.; Xue, T.; Ren, Y.; Tan, F.; Li, C. Characterization of a novel lncRNA (SETD3-OT) in turbot (Scophthalmus maximus L.). Fish Shellfish Immunol. 2020, 102, 145–151. [Google Scholar] [CrossRef]
  2. Jeck, W.R.; Sorrentino, J.A.; Wang, K.; Slevin, M.K.; Burd, C.E.; Liu, J.; Marzluff, W.F.; Sharpless, E.N. Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA 2013, 19, 141–157. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Hsu, M.-T.; Coca-Prados, M. Electron microscopic evidence for the circular form of RNA in the cytoplasm of eukaryotic cells. Nature 1979, 280, 339–340. [Google Scholar] [CrossRef] [PubMed]
  4. Cocquerelle, C.; Mascrez, B.; Hétuin, D.; Bailleul, B. Mis-splicing yields circular RNA molecules. FASEB J. 1993, 7, 155–160. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Nigro, J.M.; Cho, K.R.; Fearon, E.R.; Kern, S.E.; Ruppert, J.; Oliner, J.D.; Kinzler, K.W.; Vogelstein, B. Scrambled exons. Cell 1991, 64, 607–613. [Google Scholar] [CrossRef]
  6. Li, H.; Li, K.; Lai, W.; Li, X.; Wang, H.; Yang, J.; Chu, S.; Wang, H.; Kang, C.; Qiu, Y. Comprehensive circular RNA profiles in plasma reveals that circular RNAs can be used as novel biomarkers for systemic lupus erythematosus. Clin. Chim. Acta 2018, 480, 17–25. [Google Scholar] [CrossRef] [PubMed]
  7. Rybak-Wolf, A.; Stottmeister, C.; Glažar, P.; Jens, M.; Pino, N.; Giusti, S.; Hanan, M.; Behm, M.; Bartok, O.; Ashwal-Fluss, R.; et al. Circular RNAs in the Mammalian Brain Are Highly Abundant, Conserved, and Dynamically Expressed. Mol. Cell 2015, 58, 870–885. [Google Scholar] [CrossRef] [Green Version]
  8. Sai, L.; Li, L.; Hu, C.; Qu, B.; Guo, Q.; Jia, Q.; Zhang, Y.; Bo, C.; Li, X.; Shao, H.; et al. Identification of circular RNAs and their alterations involved in developing male Xenopus laevis chronically exposed to atrazine. Chemosphere 2018, 200, 295–301. [Google Scholar] [CrossRef]
  9. Xiu, Y.; Jiang, G.; Zhou, S.; Diao, J.; Liu, H.; Su, B.; Li, C. Identification of Potential Immune-Related circRNA-miRNA-mRNA Regulatory Network in Intestine of Paralichthys olivaceus During Edwardsiella tarda Infection. Front. Genet. 2019, 10, 731. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Kulcheski, F.R.; Christoff, A.P.; Margis, R. Circular RNAs are miRNA sponges and can be used as a new class of biomarker. J. Biotechnol. 2016, 238, 42–51. [Google Scholar] [CrossRef]
  11. Berezikov, E. Evolution of microRNA diversity and regulation in animals. Nat. Rev. Genet. 2011, 12, 846–860. [Google Scholar] [CrossRef] [PubMed]
  12. Panda, A.C. Circular RNAs Act as miRNA Sponges. Adv. Exp. Med. Biol. 2018, 1087, 67–79. [Google Scholar] [PubMed]
  13. Chen, L.L. The expanding regulatory mechanisms and cellular functions of circular RNAs. Nat. Rev. Mol. Cell Biol. 2020, 21, 475–490. [Google Scholar] [CrossRef] [PubMed]
  14. Yang, B.; Zhao, J.; Huo, T.; Zhang, M.; Wu, X. Effects of CircRNA-ITCH on proliferation and apoptosis of hepatocellular carcinoma cells through inhibiting Wnt/β-catenin signaling pathway. J. BUON 2020, 25, 1368–1374. [Google Scholar] [PubMed]
  15. Zhu, Z.; Ma, Y.; Li, Y.; Li, P.; Cheng, Z.; Li, H.; Zhang, L.; Tang, Z. The comprehensive detection of miRNA, lncRNA, and circRNA in regulation of mouse melanocyte and skin development. Biol. Res. 2020, 53, 4. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Zhang, Q.; Wang, W.; Zhou, Q.; Chen, C.; Yuan, W.; Liu, J.; Li, X.; Sun, Z. Roles of circRNAs in the tumour microenvironment. Mol. Cancer 2020, 19, 14. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Peng, P.; Zhang, B.; Huang, J.; Xing, C.; Liu, W.; Sun, C.; Guo, W.; Yao, S.; Ruan, W.; Ning, G.; et al. Identification of a circRNA-miRNA-mRNA network to explore the effects of circRNAs on pathogenesis and treatment of spinal cord injury. Life Sci. 2020, 257, 118039. [Google Scholar] [CrossRef] [PubMed]
  18. Liu, N.; Wang, Z.-Z.; Zhao, M.; Zhang, Y.; Chen, N.-H. Role of non-coding RNA in the pathogenesis of depression. Gene 2020, 735, 144276. [Google Scholar] [CrossRef] [PubMed]
  19. Yan, L.; Chen, Y.G. Circular RNAs in Immune Response and Viral Infection. Trends Biochem. Sci. 2020, 45, 1022–1034. [Google Scholar] [CrossRef] [PubMed]
  20. Zhou, Z.; Sun, B.; Huang, S.; Zhao, L. Roles of circular RNAs in immune regulation and autoimmune diseases. Cell Death Dis. 2019, 10, 503. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  21. Zhao, H.; Shao, X.; Liu, H.; Liu, Q.; Lu, J.; Li, W. The circRNA_102911/miR-129-5p/SOX6 axis is involved with T lymphocyte immune function in elderly patients with laparoscopic left hepatectomy for hepatolithiasis. Exp. Ther. Med. 2020, 21. [Google Scholar] [CrossRef] [PubMed]
  22. Wei, C.-Y.; Zhu, M.-X.; Lu, N.-H.; Liu, J.-Q.; Yang, Y.-W.; Zhang, Y.; Shi, Y.-D.; Feng, Z.-H.; Li, J.-X.; Qi, F.-Z.; et al. Circular RNA circ_0020710 drives tumor progression and immune evasion by regulating the miR-370-3p/CXCL12 axis in melanoma. Mol. Cancer 2020, 19, 84. [Google Scholar] [CrossRef] [PubMed]
  23. Huang, H.-S.; Huang, X.Y.; Yu, H.Z.; Xue, Y.; Zhu, P.L. Corrigendum to first author affiliations of circular RNA circ-RELL1 regulates inflammatory response by miR-6873-3p/MyD88/NF-kB axis in endothelial cells. Biochem. Biophys. Res. Commun. 2020, 532, 682. [Google Scholar] [CrossRef]
  24. He, L.; Zhang, A.; Xiong, L.; Li, Y.; Huang, R.; Liao, L.; Zhu, Z. Deep Circular RNA Sequencing Provides Insights into the Mechanism Underlying Grass Carp Reovirus Infection. Int. J. Mol. Sci. 2017, 18, 1977. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Su, H.; Chu, Q.; Zheng, W.; Chang, R.; Gao, W.; Zhang, L.; Xu, T. Circular RNA circPIKfyve acts as a sponge of miR-21-3p to enhance antiviral immunity through regulating MAVS in teleost fish. J. Virol. 2021, 95. [Google Scholar] [CrossRef]
  26. Chu, Q.; Zheng, W.; Su, H.; Zhang, L.; Chang, R.; Gao, W.; Xu, T. A Highly Conserved Circular RNA circRasGEF1B Enhances Antiviral Immunity by Regulating miR-21-3p/MITA Pathway in Lower Vertebrates. J. Virol. 2021, 95. [Google Scholar] [CrossRef] [PubMed]
  27. Li, X.; Yang, L.; Chen, L.-L. The Biogenesis, Functions, and Challenges of Circular RNAs. Mol. Cell 2018, 71, 428–442. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Zhang, J.; Liu, Y.; Shi, G. The circRNA-miRNA-mRNA regulatory network in systemic lupus erythematosus. Clin. Rheumatol. 2021, 40, 331–339. [Google Scholar] [CrossRef]
  29. Wu, T.; Li, Y.; Liang, X.; Liu, X.; Tang, M. Identification of potential circRNA-miRNA-mRNA regulatory networks in response to graphene quantum dots in microglia by microarray analysis. Ecotoxicol. Environ. Saf. 2021, 208, 111672. [Google Scholar] [CrossRef] [PubMed]
  30. Zou, F.W.; Cao, D.; Tang, Y.F.; Shu, L.; Zuo, Z.; Zhang, L.Y. Identification of CircRNA-miRNA-mRNA Regulatory Network in Gastrointestinal Stromal Tumor. Front. Genet. 2020, 11, 403. [Google Scholar] [CrossRef] [PubMed]
  31. Su, Q.; Lv, X. Revealing new landscape of cardiovascular disease through circular RNA-miRNA-mRNA axis. Genomics 2020, 112, 1680–1685. [Google Scholar] [CrossRef]
  32. Cao, M.; Yan, X.; Su, B.; Yang, N.; Fu, Q.; Xue, T.; Song, L.; Li, Q.; Li, C. Integrated Analysis of circRNA-miRNA-mRNA Regulatory Networks in the Intestine of Sebastes schlegelii Following Edwardsiella tarda Challenge. Front. Immunol. 2021, 11, 618687. [Google Scholar] [CrossRef]
  33. Shen, Y.; Liang, W.; Lin, Y.; Yang, H.; Chen, X.; Feng, P.; Zhang, B.; Zhu, J.; Zhang, Y.; Luo, H. Single molecule real-time sequencing and RNA-seq unravel the role of long non-coding and circular RNA in the regulatory network during Nile tilapia (Oreochromis niloticus) infection with Streptococcus agalactiae. Fish Shellfish Immunol. 2020, 104, 640–653. [Google Scholar] [CrossRef] [PubMed]
  34. Coscelli, G.; Bermúdez, R.; Ronza, P.; Losada, A.P.; Quiroga, M.I.; García, A.P.L. Immunohistochemical study of inducible nitric oxide synthase and tumour necrosis factor alpha response in turbot (Scophthalmus maximus) experimentally infected with Aeromonas salmonicida subsp. salmonicida. Fish Shellfish Immunol. 2016, 56, 294–302. [Google Scholar] [CrossRef]
  35. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10, R25. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Memczak, S.; Jens, M.; Elefsinioti, A.; Torti, F.; Krueger, J.; Rybak, A.; Maier, L.; Mackowiak, S.D.; Gregersen, L.H.; Munschauer, M.; et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 2013, 495, 333–338. [Google Scholar] [CrossRef] [PubMed]
  37. Gao, Y.; Zhang, J.; Zhao, F. Circular RNA identification based on multiple seed matching. Brief. Bioinform. 2018, 19, 803–810. [Google Scholar] [CrossRef]
  38. Wen, M.; Shen, Y.; Shi, S.; Tang, T. MiREvo: An integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinform. 2012, 13, 140. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Friedländer, M.R.; Mackowiak, S.; Li, N.; Chen, W.; Rajewsky, N. MiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012, 40, 37–52. [Google Scholar] [CrossRef]
  40. Enright, A.J.; John, B.; Gaul, U.; Tuschl, T.; Sander, C.; Marks, D.S. MicroRNA targets in Drosophila. Genome Biol. 2003, 5, R1. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Zhou, L.; Chen, J.; Li, Z.; Li, X.; Hu, X.; Huang, Y.; Zhao, X.; Liang, C.; Wang, Y.; Sun, L.; et al. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS ONE 2010, 5, e15224. [Google Scholar] [CrossRef] [PubMed]
  42. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Young, M.D.; Wakefield, M.J.; Smyth, G.K.; Oshlack, A. Gene ontology analysis for RNA-seq: Accounting for selection bias. Genome Biol. 2010, 11, R14. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Su, G.; Morris, J.H.; Demchak, B.; Bader, G.D. Biological network exploration with Cytoscape 3. Curr. Protoc. Bioinform. 2014, 47, 8.13.1–8.13.24. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Livak, K.J.; Schmittgen, T.D. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  46. Ali, A.; Rexroad, C.E.; Thorgaard, G.H.; Yao, J.; Salem, M. Characterization of the rainbow trout spleen transcriptome and identification of immune-related genes. Front. Genet. 2014, 5, 348. [Google Scholar] [CrossRef] [Green Version]
  47. Muñoz-Atienza, E.; Díaz-Rosales, P.; Tafalla, C. Systemic and Mucosal B and T Cell Responses Upon Mucosal Vaccination of Teleost Fish. Front. Immunol. 2021, 11, 622377. [Google Scholar] [CrossRef]
  48. Chakrabarti, J.; Mitra, S. (Eds.) Chapter 1—Introduction. In Cancer and Noncoding RNAs; Academic Press: Cambridge, MA, USA, 2017; Volume 1, pp. 1–23. [Google Scholar]
  49. Hombach, S.; Kretz, M. Non-coding RNAs: Classification, Biology and Functioning. Adv. Exp. Med. Biol. 2016, 937, 3–17. [Google Scholar]
  50. Fu, Q.; Liu, C.-J.; Zhai, Z.-S.; Zhang, X.; Qin, T.; Zhang, H.-W. Single-Cell Non-coding RNA in Embryonic Development. Adv. Exp. Med. Biol. 2018, 1068, 19–32. [Google Scholar]
  51. Schmitz, S.U.; Grote, P.; Herrmann, B.G. Mechanisms of long noncoding RNA function in development and disease. Cell. Mol. Life Sci. 2016, 73, 2491–2509. [Google Scholar] [CrossRef] [Green Version]
  52. Chistiakov, D.A.; Orekhov, A.N.; Bobryshev, Y.V. The role of miR-126 in embryonic angiogenesis, adult vascular homeostasis, and vascular repair and its alterations in atherosclerotic disease. J. Mol. Cell. Cardiol. 2016, 97, 47–55. [Google Scholar] [CrossRef]
  53. Chen, X.; Ouyang, Z.; Shen, Y.; Liu, B.; Zhang, Q.; Wan, L.; Yin, Z.; Zhu, W.; Li, S.; Peng, D. CircRNA_28313/miR-195a/CSF1 axis modulates osteoclast differentiation to affect OVX-induced bone absorption in mice. RNA Biol. 2019, 16, 1249–1262. [Google Scholar] [CrossRef]
  54. Yang, J.; Cheng, M.; Gu, B.; Wang, J.; Yan, S.; Xu, D. CircRNA_09505 aggravates inflammation and joint damage in collagen-induced arthritis mice via miR-6089/AKT1/NF-κB axis. Cell Death Dis. 2020, 11, 833. [Google Scholar] [CrossRef] [PubMed]
  55. Yang, L.; Han, X.; Zhang, C.; Sun, C.; Huang, S.; Xiao, W.; Gao, Y.; Liang, Q.; Luo, F.; Lu, W.; et al. Hsa_circ_0060450 Negatively Regulates Type I Interferon-Induced Inflammation by Serving as miR-199a-5p Sponge in Type 1 Diabetes Mellitus. Front. Immunol. 2020, 11, 576903. [Google Scholar] [CrossRef]
  56. Li, S.; Wang, N.; Zhang, T.; Feng, Y.; Wang, L.; Sun, J. Characterization of three connexin32 genes and their role in inflammation-induced ATP release in the Japanese flounder Paralichthys olivaceus. Fish Shellfish Immunol. 2020, 106, 181–189. [Google Scholar] [CrossRef] [PubMed]
  57. Moser, B.; Wolf, M.; Walz, A.; Loetscher, P. Chemokines: Multiple levels of leukocyte migration control. Trends Immunol. 2004, 25, 75–84. [Google Scholar] [CrossRef] [PubMed]
  58. Laing, K.J.; Secombes, C.J. Chemokines. Dev. Comp. Immunol. 2004, 28, 443–460. [Google Scholar] [CrossRef]
  59. Yan, Y.; Chen, R.; Wang, X.; Hu, K.; Huang, L.; Lu, M.; Hu, Q. CCL19 and CCR7 Expression, Signaling Pathways, and Adjuvant Functions in Viral Infection and Prevention. Front. Cell Dev. Biol. 2019, 7, 212. [Google Scholar] [CrossRef]
  60. Sun, J.; Wang, J.; Lu, W.; Xie, L.; Lv, J.; Li, H.; Yang, S. MiR-325-3p inhibits renal inflammation and fibrosis by targeting CCL19 in diabetic nephropathy. Clin. Exp. Pharmacol. Physiol. 2020, 47, 1850–1860. [Google Scholar] [CrossRef]
  61. Chen, F.; Lu, X.-J.; Nie, L.; Ning, Y.-J.; Chen, J. Molecular characterization of a CC motif chemokine 19-like gene in ayu (Plecoglossus altivelis) and its role in leukocyte trafficking. Fish Shellfish Immunol. 2018, 72, 301–308. [Google Scholar] [CrossRef] [PubMed]
  62. Sepahi, A.; Tacchi, L.; Casadei, E.; Takizawa, F.; LaPatra, S.E.; Salinas, I. CK12a, a CCL19-like Chemokine That Orchestrates both Nasal and Systemic Antiviral Immune Responses in Rainbow Trout. J. Immunol. 2017, 199, 3900–3913. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Chen, C.; Hu, Y.-H.; Xiao, Z.-Z.; Sun, L. SmCCL19, a CC chemokine of turbot Scophthalmus maximus, induces leukocyte trafficking and promotes anti-viral and anti-bacterial defense. Fish Shellfish Immunol. 2013, 35, 1677–1682. [Google Scholar] [CrossRef] [PubMed]
Biology 10 00626 g001 550
Figure 1. Analysis of circRNA sequencing data. (A) The distribution of circRNAs on turbot different chromosomes; the box represents the length of the chromosomes. (B) The number of circRNAs arising from different genomic locus (exon, intronic, and intergenic). (C) Size distribution of unique circRNA sequences from ASS0, ASS6, ASS12, ASS24, and ASS96 groups. (D) Volcano Plots of DE-circRNAs among ASS0 and infected groups. Red blocks represent up-regulated circRNAs and green blocks represent down-regulated circRNAs. (E) Heatmap was used to assess the expression of DE-circRNAs. Red and blue denote high and low expression, respectively. Each DE-circRNA is represented by a single row of colored boxes, and each group is represented by a single column. (F) Pairwise comparisons of differentially expressed circRNAs between two groups. The Venn diagrams display the distribution of the 119 unique DE-circRNAs between ASS6, ASS12, ASS24, ASS96, and ASS0, respectively. (G) Go term analysis of DE-circRNAs. (H) KEGG analysis of DE-circRNAs.
Figure 1. Analysis of circRNA sequencing data. (A) The distribution of circRNAs on turbot different chromosomes; the box represents the length of the chromosomes. (B) The number of circRNAs arising from different genomic locus (exon, intronic, and intergenic). (C) Size distribution of unique circRNA sequences from ASS0, ASS6, ASS12, ASS24, and ASS96 groups. (D) Volcano Plots of DE-circRNAs among ASS0 and infected groups. Red blocks represent up-regulated circRNAs and green blocks represent down-regulated circRNAs. (E) Heatmap was used to assess the expression of DE-circRNAs. Red and blue denote high and low expression, respectively. Each DE-circRNA is represented by a single row of colored boxes, and each group is represented by a single column. (F) Pairwise comparisons of differentially expressed circRNAs between two groups. The Venn diagrams display the distribution of the 119 unique DE-circRNAs between ASS6, ASS12, ASS24, ASS96, and ASS0, respectively. (G) Go term analysis of DE-circRNAs. (H) KEGG analysis of DE-circRNAs.
Biology 10 00626 g001
Biology 10 00626 g002 550
Figure 2. Analysis of miRNA sequencing data. (A) The length of miRNAs from ASS0, ASS6, ASS12, ASS24, and ASS96 groups. (B) TPM density distribution of sRNAs from 15 samples. (C) Volcano Plots of DE-miRNAs among ASS0 and infected groups. Red blocks represent up-regulated miRNAs and green blocks represent down-regulated miRNAs. (D) Heatmap was used to assess the expression of DE-miRNAs. Red and blue denote high and low expression, respectively. Each DE-miRNA is represented by a single row of colored boxes, and each group is represented by a single column. (E) Pairwise comparisons of differentially expressed miRNAs between two groups. The Venn diagrams display the distribution of the 140 unique DE-miRNAs between ASS6, ASS12, ASS24, ASS96, and ASS0, respectively. (F) Go term analysis of DE-miRNAs. (G) KEGG analysis of DE-miRNAs.
Figure 2. Analysis of miRNA sequencing data. (A) The length of miRNAs from ASS0, ASS6, ASS12, ASS24, and ASS96 groups. (B) TPM density distribution of sRNAs from 15 samples. (C) Volcano Plots of DE-miRNAs among ASS0 and infected groups. Red blocks represent up-regulated miRNAs and green blocks represent down-regulated miRNAs. (D) Heatmap was used to assess the expression of DE-miRNAs. Red and blue denote high and low expression, respectively. Each DE-miRNA is represented by a single row of colored boxes, and each group is represented by a single column. (E) Pairwise comparisons of differentially expressed miRNAs between two groups. The Venn diagrams display the distribution of the 140 unique DE-miRNAs between ASS6, ASS12, ASS24, ASS96, and ASS0, respectively. (F) Go term analysis of DE-miRNAs. (G) KEGG analysis of DE-miRNAs.
Biology 10 00626 g002
Biology 10 00626 g003 550
Figure 3. Analysis of mRNA sequencing data. (A) Volcano Plots of DE-mRNAs among ASS0 and infected groups. Red blocks represent up-regulated mRNAs and green blocks represent down-regulated mRNAs. (B) Heatmap was used to assess the expression of DE-mRNAs. Red and blue denote high and low expression, respectively. Each DE-mRNA is represented by a single row of colored boxes, and each group is represented by a single column. (C) Pairwise comparisons of differentially expressed mRNAs between two groups. The Venn diagrams display the distribution of the 530 unique DE-mRNAs between ASS6, ASS12, ASS24, ASS96, and ASS0, respectively. (D) Go term analysis of DE-mRNAs. (E) KEGG analysis of DE-mRNAs.
Figure 3. Analysis of mRNA sequencing data. (A) Volcano Plots of DE-mRNAs among ASS0 and infected groups. Red blocks represent up-regulated mRNAs and green blocks represent down-regulated mRNAs. (B) Heatmap was used to assess the expression of DE-mRNAs. Red and blue denote high and low expression, respectively. Each DE-mRNA is represented by a single row of colored boxes, and each group is represented by a single column. (C) Pairwise comparisons of differentially expressed mRNAs between two groups. The Venn diagrams display the distribution of the 530 unique DE-mRNAs between ASS6, ASS12, ASS24, ASS96, and ASS0, respectively. (D) Go term analysis of DE-mRNAs. (E) KEGG analysis of DE-mRNAs.
Biology 10 00626 g003
Biology 10 00626 g004 550
Figure 4. CircRNA-miRNA-mRNA regulatory network. (A) The DE-circRNA-miRNA-mRNA regulatory networks found in ASS6, ASS12, ASS24, and ASS96 groups compared with ASS0. (B) The DE-circRNA-miRNA-mRNA regulatory network in the form of “up-down-up” (circRNA-miRNA-mRNA) or “down-up-down” after infection. CircRNAs, miRNAs, and mRNA are indicated by square, circle, and triangle, respectively; red represents up regulation and blue represents down regulation.
Figure 4. CircRNA-miRNA-mRNA regulatory network. (A) The DE-circRNA-miRNA-mRNA regulatory networks found in ASS6, ASS12, ASS24, and ASS96 groups compared with ASS0. (B) The DE-circRNA-miRNA-mRNA regulatory network in the form of “up-down-up” (circRNA-miRNA-mRNA) or “down-up-down” after infection. CircRNAs, miRNAs, and mRNA are indicated by square, circle, and triangle, respectively; red represents up regulation and blue represents down regulation.
Biology 10 00626 g004
Biology 10 00626 g005 550
Figure 5. Verification of circRNA expression level by qRT-PCR. Four circRNAs were randomly selected for qRT-qPCR analysis and compared with data obtained using RNA-sequencing (RNA-seq). The relative expression levels of circRNAs in ASS6, ASS12, ASS24, and ASS96 groups were calculated as the ratio of the gene expression level (qRT-PCR) or normalized transcripts per million (TPM) (RNA-seq) relative to ASS0 group. The qRT-PCR data are given as means ± standard deviation of three replicates.
Figure 5. Verification of circRNA expression level by qRT-PCR. Four circRNAs were randomly selected for qRT-qPCR analysis and compared with data obtained using RNA-sequencing (RNA-seq). The relative expression levels of circRNAs in ASS6, ASS12, ASS24, and ASS96 groups were calculated as the ratio of the gene expression level (qRT-PCR) or normalized transcripts per million (TPM) (RNA-seq) relative to ASS0 group. The qRT-PCR data are given as means ± standard deviation of three replicates.
Biology 10 00626 g005
Biology 10 00626 g006 550
Figure 6. Verification of miRNA expression level by qRT-PCR. Four miRNAs were randomly selected for qRT-qPCR analysis and compared with data obtained using RNA-sequencing (RNA-seq). The relative expression level of miRNAs in ASS6, ASS12, ASS24, and ASS96 groups was calculated as the ratio of the gene expression level (qRT-PCR) or normalized transcripts per million (TPM) (RNA-seq) relative to the ASS0 group. The qRT-PCR data are given as means ± standard deviation of three replicates.
Figure 6. Verification of miRNA expression level by qRT-PCR. Four miRNAs were randomly selected for qRT-qPCR analysis and compared with data obtained using RNA-sequencing (RNA-seq). The relative expression level of miRNAs in ASS6, ASS12, ASS24, and ASS96 groups was calculated as the ratio of the gene expression level (qRT-PCR) or normalized transcripts per million (TPM) (RNA-seq) relative to the ASS0 group. The qRT-PCR data are given as means ± standard deviation of three replicates.
Biology 10 00626 g006
Biology 10 00626 g007 550
Figure 7. Verification of mRNA expression level by qRT-PCR. Eight mRNAs were randomly selected for qRT-qPCR analysis and compared with data obtained using RNA-sequencing (RNA-seq). The relative expression level of mRNAs in ASS6, ASS12, ASS24, and ASS96 groups was calculated as the ratio of the gene expression level (qRT-PCR) or normalized transcripts per million (TPM) (RNA-seq) relative to the ASS0 group. The qRT-PCR data are given as means ± standard deviation of three replicates. tnfaip2: tumor necrosis factor alpha-induced protein 2.
Figure 7. Verification of mRNA expression level by qRT-PCR. Eight mRNAs were randomly selected for qRT-qPCR analysis and compared with data obtained using RNA-sequencing (RNA-seq). The relative expression level of mRNAs in ASS6, ASS12, ASS24, and ASS96 groups was calculated as the ratio of the gene expression level (qRT-PCR) or normalized transcripts per million (TPM) (RNA-seq) relative to the ASS0 group. The qRT-PCR data are given as means ± standard deviation of three replicates. tnfaip2: tumor necrosis factor alpha-induced protein 2.
Biology 10 00626 g007
Table
Table 1. Information list of circRNA sequencing data.
Table 1. Information list of circRNA sequencing data.
Sample_NameRaw_ReadsClean_ReadsClean_Bases
(Gb)		Error Rate
(%)		Q20 (%)Q30 (%)GC_Content
(%)		Mapped to Genome (%)
ASS0189,427,87087,674,31813.150.0397.2292.5854.1189.81
ASS0288,835,29686,925,03413.040.0397.3292.8451.3390.62
ASS0391,281,97687,308,90813.100.0397.0192.2353.2286.94
ASS6184,909,27478,465,24611.770.0396.8192.1255.2890.44
ASS6295,743,38693,531,87014.030.0397.3792.9352.1787.85
ASS6389,110,60287,785,72813.170.0397.5193.2152.9490.53
ASS12197,847,96896,072,44614.410.0397.4292.9353.2185.35
ASS122107,925,458106,515,64815.980.0397.2692.6053.0487.26
ASS12388,551,05686,735,52813.010.0397.2492.6152.5490.15
ASS241108,968,470107,506,05016.130.0397.2692.6755.7480.36
ASS24293,831,13692,152,78013.820.0397.4192.9651.6985.96
ASS24386,686,37484,702,79812.710.0397.2992.7653.0888.55
ASS961100,039,32698,296,58614.740.0397.6893.6154.4788.51
ASS96287,993,17083,798,12812.570.0396.9892.2453.9084.57
ASS963104,529,398102,835,52615.430.0397.6293.4353.8184.48
Table
Table 2. Information list of miRNA sequencing data.
Table 2. Information list of miRNA sequencing data.
SampleTotal_ReadsClean ReadsError Rate (%)GC Content
(%)		Mapped sRNA (%)Uniq Reads
ASS0116,992,25716,776,8910.0148.9194.56225,970
ASS0213,585,72913,445,1930.0149.7291.22256,756
ASS0314,090,92613,924,0280.0149.1793.82233,122
ASS6116,759,63116,480,9230.0148.4394.80292,697
ASS6214,970,40814,487,8620.0149.2790.79297,090
ASS6317,691,02717,287,8380.0148.9793.65380,538
ASS12115,944,20715,759,0310.0148.9095.73280,956
ASS12214,993,05514,650,2900.0149.4392.89331,440
ASS12314,495,40810,351,4180.0151.0889.65172,990
ASS24119,726,92919,222,4540.0148.9595.82276,881
ASS24213,696,25312,062,2660.0154.2866.84301,900
ASS24315,432,97415,212,0950.0148.6095.28305,646
ASS96116,295,41316,051,4910.0148.7095.38326,589
ASS96218,209,62817,873,0270.0148.3594.96376,083
ASS96317,714,90817,364,7450.0149.4493.82319,525
Table
Table 3. Information of identified miRNAs.
Table 3. Information of identified miRNAs.
TypesKnown miRNANew miRNA
Mapped MatureMapped HairpinMapped Uniq sRNAMapped Total sRNAMapped MatureMapped StarMapped HairpinMapped Uniq sRNAMapped Total sRNA
ASS01214258253910,109,0193541594221721849,828
ASS0220625523947,201,0713301434001504486,725
ASS0321225725598,142,4783551734321714577,499
ASS121218263264310,991,7263891804641910721,131
ASS12221025724688,730,4533441674041690641,866
ASS12320223718343,157,3602341123051036234,731
ASS241216245290612,630,4123931854792055974,386
ASS24217922312721,934,747207117282810154,078
ASS243211258257310,726,0493971744741935776,575
ASS61208271244610,849,3913961684621829779,494
ASS6220324924017,794,5813181513921613356,696
ASS63211256255911,483,8323991674701867693,347
ASS961211241269010,345,7873951864642005880,402
ASS962212259260212,439,7174151984932073983,840
ASS963214256272010,828,0353941784721941845,700
Total25228665,488245,032,27666844272645,78416,824,392
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Xue, T.; Liu, Y.; Cao, M.; Tian, M.; Zhang, L.; Wang, B.; Liu, X.; Li, C. Revealing New Landscape of Turbot (Scophthalmus maximus) Spleen Infected with Aeromonas salmonicida through Immune Related circRNA-miRNA-mRNA Axis. Biology 2021, 10, 626. https://doi.org/10.3390/biology10070626

AMA Style

Xue T, Liu Y, Cao M, Tian M, Zhang L, Wang B, Liu X, Li C. Revealing New Landscape of Turbot (Scophthalmus maximus) Spleen Infected with Aeromonas salmonicida through Immune Related circRNA-miRNA-mRNA Axis. Biology. 2021; 10(7):626. https://doi.org/10.3390/biology10070626

Chicago/Turabian Style

Xue, Ting, Yiping Liu, Min Cao, Mengyu Tian, Lu Zhang, Beibei Wang, Xiaoli Liu, and Chao Li. 2021. "Revealing New Landscape of Turbot (Scophthalmus maximus) Spleen Infected with Aeromonas salmonicida through Immune Related circRNA-miRNA-mRNA Axis" Biology 10, no. 7: 626. https://doi.org/10.3390/biology10070626

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop