Comprehensive analysis of circRNA expression pattern and circRNA–mRNA–miRNA network in Ctenopharyngodon idellus kidney (CIK) cells after grass carp reovirus (GCRV) infection

Aquaculture 258 (2019) 734349

Contents lists available at ScienceDirect

Aquaculture journal homepage: www.elsevier.com/locate/aquaculture

Comprehensive analysis of circRNA expression pattern and circRNA–mRNA–miRNA network in Ctenopharyngodon idellus kidney (CIK) cells after grass carp reovirus (GCRV) infection

T

Bo Liua,1, Rui Yuanb,1, Zi Lianga,1, Tingting Zhanga, Min Zhua, Xing Zhanga, Wei Genga, Ping Fangb, Mengsheng Jianga, Zhangyan Wanga, Yongjie Fenga,c, Xunmeng Liub, Yang Zhoud, ⁎ ⁎ Renyu Xuea,c, Guangli Caoa,c, Hui Chenb, Xiaolong Hua,c, , Chengliang Gonga,c, a

School of Biology and Basic Medical Sciences, Soochow University, Suzhou, Jiangsu 215123, China Jiangsu Center for Control and Prevention of Aquatic Animal Infectious Disease, Nanjing 210036, China Agricultural Biotechnology Research Institute, Agricultural Biotechnology and Ecological Research Institute, Soochow University, Suzhou 215123, China d Aquaculture Technical Extension Station of Dafeng City, Dafeng 224100, China b c

ARTICLE INFO

ABSTRACT

Keywords: Ctenopharyngodon idellus kidney cells circRNA Grass carp reovirus mRNA miRNA

Grass carp reovirus (GCRV) causes grass carp hemorrhage disease and is one of the most damaging factors in the culture of grass carp (Ctenopharyngodon idellus). Until now, the mechanisms underlying GCRV infection have been largely unknown. To explore the molecular responses to GCRV infection in C. idellus kidney (CIK) cells, we used high-throughput sequencing to investigate the circRNA, mRNA, and miRNA expression patterns and the circRNA–miRNA–mRNA network. A total of 76 circRNAs, 798 mRNAs, and 186 miRNAs were identified as differentially expressed (DE) RNAs in GCRV-infected cells compared with uninfected cells, of which 40 circRNAs, 658 mRNAs, and 146 miRNAs were upregulated, and 36 circRNAs, 140 mRNAs, and 40 miRNAs were downregulated (fold change > 2, p < .05). A Gene Ontology (GO) analysis of the DE circRNA parental genes, mRNAs, and miRNA target genes was performed. The parental genes of the DE circRNAs, the DE mRNAs, and the target genes of the DE miRNAs were mainly enriched in metal-ion-binding-, immune-, and DNA-replication-related items, respectively. A KEGG pathway analysis showed that the DE circRNA parental genes, DE mRNAs, and miRNA target genes were predominantly enriched in endocytosis-, immune-, and infection-related pathways, respectively. The miRNA–target interactions were also predicted and competing regulatory networks of endogenous RNAs were constructed based on the DE RNAs. Our prediction results suggest that the circRNA-2780/ cik–miRNA-7796/KAT6B, circRNA-1351/cik-miRNA–10648/B2L14 and circRNA-1591/cik-miR-8141/AMPD3 axes may be associated with GCRV infection. This study suggests that the crosstalk between circRNAs and their competing mRNAs play crucial roles in GCRV infection.

1. Introduction The grass carp Ctenopharyngodon idella is currently among the most widely cultured freshwater aquaculture fish species in China (He et al., 2017a; He et al., 2017b). Unfortunately, frequent disease has caused huge economic losses have occurred in the grass carp aquaculture industry (He et al., 2017a). The grass carp reovirus (GCRV) causes grass carp hemorrhage disease in cultured grass carp, with high mortality (Rao and Su, 2015). GCRV is classified in the genus Aquareovirus in the family Reoviridae, and has a genome of 11 double-stranded RNA

(dsRNA) segments (Pei et al., 2014). Twelve proteins are encoded by the 11 dsRNA segments, including seven structural proteins (VP1, VP2, VP3, VP4, VP5, VP6, and VP7) and five nonstructural proteins (NS80, NS38, NS31, NS26, and NS16) (Wang et al., 2012; Liu et al., 2016). In recent years, grass carp hemorrhagic disease, caused by GCRV, has received particular attention because this virus had widespread hosts. The rare minnow (Gobiocypris rarus), black carp (Mylopharyngodon piceus), and topmouth gudgeon (Pseudorasbora parva), all these cultured fish species, can be infected by GCRV, leading to hemorrhagic symptoms and death (He et al., 2017c). Because of the frequent outbreaks of

Corresponding authors at: School of Biology and Basic Medical Sciences, Soochow University, No.199 Ren'ai Road, Dushu Lake Higher Education Town, Suzhou Industrial Park, Suzhou 215123, PR China. E-mail addresses: [email protected] (X. Hu), [email protected] (C. Gong). 1 These authors contributed equally to this work. ⁎

https://doi.org/10.1016/j.aquaculture.2019.734349 Received 22 December 2018; Received in revised form 13 July 2019; Accepted 29 July 2019 Available online 29 July 2019 0044-8486/ © 2019 Elsevier B.V. All rights reserved.

Aquaculture 258 (2019) 734349

B. Liu, et al.

grass carp hemorrhagic disease in freshwater aquaculture, many studies of GCRV have been undertaken, but most have focused on the identification of immune genes (Chen et al., 2012; Su et al., 2012; Yang et al., 2012; Chen et al., 2013; Wang et al., 2013b; Jin et al., 2018), antiviral genes (Wang et al., 2013a; Wan and Su, 2015; Ma et al., 2018), and vaccine development (Xue et al., 2013; Wang et al., 2015; Chen et al., 2018a; Gao et al., 2018), and have sequenced the GCRV-infected grass carp transcriptome to identify immune response genes (Shi et al., 2014; Dai et al., 2017). A transcriptomic analysis provides a global perspective on the RNA expression patterns in the host in response to various pathogens, such as Hyphantria cunea larvae and Nuclear polyhedrosis virus (Sun et al., 2018), viral hemorrhagic septicemia virus in the olive flounder (Jeong et al., 2018), Aeromonas hydrophila in the black carp (Zhang et al., 2018a), and Aedes aegypti and Zika virus (Zhao et al., 2018). In the present study, whole-transcriptome and small RNA sequencing analyses of GCRV-infected CIK cells were used to determine the mechanisms underlying pathogenesis of GCRV and the hemorrhagic symptoms it induces. Circular RNAs (circRNAs) are a novel class of noncoding RNAs, which are abundant, widespread, and are tissue-specifically expressed in all kingdoms of organism (Memczak et al., 2013). This type of RNA characteristically has no distinct 5′ or 3′ terminus, but is covalently linked into a closed loop structure (Qu et al., 2015). Unlike mature mRNAs, the roles of most of circRNAs in various species are unknown. However, they play important roles in the regulation of gene expression at the posttranscriptional level, and have regulatory roles in the pathogenesis and progression of various diseases (Qu et al., 2015; Jin et al., 2016; Zhang et al., 2017a; Zhang et al., 2017b; Cheng et al., 2018; Lin and Chen, 2018; Zhang et al., 2018a). For example, the expression pattern of circRNAs in maize in response to maize Iranian mosaic virus infection was investigated, and 33 circRNAs were found to sponge 23 maize miRNAs that may affect maize metabolism and development (Ghorbani et al., 2018). Most circRNAs in tomato are derived from the exons of protein-coding genes, and the accumulation of Tomato yellow leaf curl virus (TYLCV) was reduced by silencing the parental genes of these circRNAs (Wang et al., 2018). The dysregulated circRNAs in the midgut of the silkworm after Bombyx mori cytoplasmic polyhedrosis virus (BmCPV) and B. mori nucleopolyhedrovirus (BmNPV) infection may be related to the pathogenesis of these viruses (Hu et al., 2018a; Hu et al., 2018b). Many of the parental genes of the circRNAs differentially expressed in GCRV-infected grass carp are involved in metal-ion binding, protein ubiquitination, enzyme activities, and nucleotide binding. Eight genes targeted by miRNAs that are sponged by the differentially expressed circRNAs are involved in immune responses, blood coagulation, hemostasis, and the complement and coagulation cascades (He et al., 2017a). However, most of the studies mentioned above have focused on other species, including the grass carp, but there has been no systematic study of the circRNAs, mRNAs, and miRNAs induced in response to GCRV infection or on the regulatory networks of competing endogenous RNAs (ceRNA) in Ctenopharyngodon idellus kidney (CIK) cells after GCRV infection. Ctenopharyngodon idellus kidney (CIK) cells derived from the grass carp were infected with GCRV, which characteristically induces a mild systemic immune response in the host, involving a variety of classical antiviral signal transduction pathways (He et al., 2017c; Yu et al., 2017; Ji et al., 2018). In this study, we used high-through sequencing to systematically investigate the virus-responsive DE circRNAs, mRNAs, and miRNAs in CIK cells after GCRV infection and predicted their putative roles during GCRV infection with a integration analysis of competing endogenous RNA moleclules (ceRNAs). CircRNA was reported as a class of ceRNAs that naturally sequester miRNAs to competitively suppress miRNA activity to regulate their mRNA targets (miRNA-dependent gene silencing involving Ago2-mediated cleavage of a circular antisense RNA; Circular RNAs are a large class of animal RNAs with regulatory potency). These results not only provide novel insights into the roles of circRNAs, mRNAs, and miRNAs in viral infection and

cellular antiviral immunity, but also extend our understanding of the mechanisms used by GCRV-infected CIK cells. 2. Materials and methods 2.1. Cells, virus, and viral infection Ctenopharyngodon idellus kidney (CIK) cells are maintained in our laboratory and were grown in Medium 199 (Gibco, Gaithersburg, MD, USA) supplemented with 10% fetal bovine serum (FBS) at 28 °C. GCRV873 was kindly provided by Prof. Hui Chen (Jiangsu Center for Control and Prevention of Aquatic Animal Infectious Disease, Nanjing, China). The CIK cells were grown in 25 ml flasks and infected with GCRV-873 at a multiplicity of infection (MOI) of 5. Uninfected CIK cells were used as the control. 2.2. RNA isolation, library construction, and sequencing At 48 h postinfection (hpi), total RNAs was extracted from the GCRV-infected and control CIK cells with TRIzol Reagent (Invitrogen, CA, USA), according to the manufacturer's manual, and any contaminating DNA was digested with DNaseI (Epicentre, Charlotte, USA) at 37 °C for 15 min. The total RNAs were quantified with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Scotts Valley, CA, USA), and the integrity and quality of the total RNAs were evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The samples with an RNA integrity number (RIN) ≥ 7 were used in the subsequent analyses. Individual libraries of circRNA, mRNA, and miRNA were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold (Illumina, USA) and TruSeq Small RNA Sample Prep Kits (Illumina), according to the manufacturer's instructions. The libraries were then subjected to whole-transcriptome and small RNA sequencing on the Illumina sequencing platform (HiSeq™ 2500) at ShangHai Oebiotech Co. (Shanghai, China). Detailed information on the two sequencing protocols is given in previous reports (Zhang et al., 2017b, Zhang et al., 2018a,b). 2.3. Data analysis Clean data were obtained from the raw data by removing the adapter sequences. The Q20 and Q30 scores and GC contents of the clean data were calculated, and all downstream analyses were performed with high-quality clean data. The clean data were mapped to the grass carp reference genome (http://www.ncgr.ac.cn/grasscarp/ files/C_idella_female_scaffolds.fasta.v1.gz) with the TopHat2 software. The reads from the two samples were mapped to assembled transcripts (http://www.ncgr.ac.cn/grasscarp/files/C_idella_female_genemodels. v1.0.gz) with Bowtie2. As to unmapped reads, the junctions were picked out using back-splice algorithm (Gao et al., 2015). circRNAs were predicted with a software developed by OE which were considered as the reference sequence for further analysis(Zhang et al., 2017b). The expression level of gene expression was calculated as fragments per kilobase per million reads (Langmead and Salzberg, 2012). 2.4. Identification of differentially expressed (DE) circRNAs, mRNAs, and miRNAs The DE circRNAs, mRNAs, and miRNAs were identified with the DESeq software (http://bioconductor.org/packages/release/bioc/ html/DESeq.html). circRNAs, mRNAs, and miRNAs displaying changes of > 2-fold with adjusted p < .05 were considered DE RNAs. 2.5. Prediction of miRNA target genes and circRNA-binding miRNAs circRNAs affect miRNA activities by sponging natural miRNAs. The 2

Aquaculture 258 (2019) 734349

B. Liu, et al.

software miRanda (http://www.microrna.org/microrna/home.do) was used to predict the interactions between DE miRNAs and circRNAs, with a threshold of energy score ΔG ≤ −16 kcal/mol and a paring score S ≥ 160. The targets of the DE miRNAs were predicted with the miRanda (v3.3a) software in animals, with the parameters S ≥ 150, △G ≤ −30 kcal/mol, and strict 5′ seed pairing (John et al., 2004).

2.9. Response of circRNAs to GCRV challenge in CIK cells CIK cells (105) cultured in 3 cm dishes were infected with GCRV-873 (MOI = 5). The infected cells were harvested at 0, 12, 24, 36, and 48 hpi and their total RNAs extracted. cDNA synthesis and real-time PCR validation were performed as described above. The expression levels of eight circRNAs were elected to validate and all experiments were performed in triplicate. EF1a gene was used as the internal reference gene. The relative expression levels were calculated with the 2−ΔΔCt method (Livak and Schmittgen, 2001).

2.6. PCR amplification and Sanger sequencing To confirm the reliability of the circRNAs identified with sequencing, 16 circRNAs were randomly selected for PCR confirmation. Divergent primers for each selected circRNA were designed with the Primer Premier 5 software (Supplementary Table 1). The divergent primers for reverse transcription–PCR (RT–PCR) validation were designed near both side sequences at the circRNA junctions. The PCR products were confirmed with Sanger sequencing.

2.10. miRNA–mRNA interaction network analysis Based on the negative regulatory correlations between the expression levels of each miRNA and its mRNA targets, all the targets of the miRNAs were extracted and the miRNA–mRNA interaction network was constructed with the Cytoscape software (Shannon et al., 2003).

2.7. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses

2.11. ceRNA network analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the parental genes of the identified DE circRNAs, the mRNAs, and miRNA target genes were performed. The GO analysis (http://www.geneontology.org) was used to meaningfully annotate the RNAs, and identified the domains of biological processes (BPs), cellular components (CCs), and molecular functions (MFs). GO terms with a p value of < 0.05 were considered to be significantly enriched. The KEGG pathway analysis was performed to identify pathway clusters covering our knowledge of the molecular interactions and reaction networks during differentially regulated gene profiling. KEGG pathways with p < .05 were considered significantly enriched. The top 20 enriched GO terms and pathways of the DE RNAs were ranked according to the enrichment scores (−log10 p values).

The circRNAs, mRNAs, and miRNAs differentially expressed in the GCRV-infected and control CIKs cells that were meaningfully correlated were subjected to a ceRNA analysis. This describes the complex posttranscriptional regulatory network mediated by miRNAs: by sharing one or more miRNA response elements (MREs) (Kartha and Subramanian, 2014). We looked for potential MREs among the circRNA, mRNA, and miRNA sequences, and the occurrence of the same miRNA-seed-sequence-binding site in a circRNA and an mRNA sequence predicted an interaction between the circRNA and miRNA. Upand downregulated RNAs were selected to establish the ceRNA networks based on their coexpression patterns with negative analysis.

2.8. Real-time PCR

3.1. Fewer circRNAs identified in GCRV-infected cells than in uninfected cells

3. Results

To confirm the sequencing data, eight DE circRNAs, eight DE miRNAs, and nine DE mRNAs were randomly selected and validated with real-time PCR. Specific divergent primers for each circRNA were designed according to the sequence of the linear transcript. The primers for the mRNAs and miRNAs were designed according to the normal regulation. All the primers used in this study are listed in Supplementary Table 1. The total RNAs was then extracted from the GCRV-infected CIK cells, digested with RNase R (or not), and purified. The cDNA was synthesized with EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen Biotech Co. Ltd., China). Realtime PCR was performed on the CFX96 Touch™ Deep Well Real-Time PCR Detection System (Bio-Rad, Berkeley, USA) according to the manufacturer's instructions, with iTaq™ Universal SYBR Green Supermix (Bio-Rad). The expression levels of the circRNAs and mRNAs were normalized to that of EF1a mRNA (internal standard control) and the expression levels of the miRNAs were normalized to that of 18S rRNA (internal standard control). The relative expression levels were calculated with the 2−ΔΔCt method (Livak and Schmittgen, 2001). All experiments were performed in triplicate.

In total, the numbers of reads for circRNAs were 108,210,498 in the GCRV-infected CIK cells and 110,762,946 in the uninfected cells, which produced mapped read counts of 106,453,214 and 109,160,724, respectively (Table 1). The sizes of the candidate circRNAs in the CIK cells ranged from < 200 nucleotides (nt) to > 2000 nt. About 78.61% of the circRNAs had a predicted spliced length of < 2000 nt, and 43.95% had a length of 201–500 nt. The average length of the circRNAs was 3161.23 bp (Fig. 1A). The origins of the circRNAs were classified into five categories. Most circRNAs (78.53%, 2743 circRNAs) were spliced from a sense overlapping circRNA, and intronic circRNAs comprised only 0.46% (16 circRNAs) (Fig. 1B). In total, 3494 circRNAs were detected with sequencing, 1568 of which were identified in GCRV-infected CIK cells and 3013 in uninfected CIK cells (Fig. 1C). The splicing sites of the total circRNAs were analyzed and nearly all the circRNAs had the GT/AG splice site (Fig. 1D). To validate the authenticity of the circRNAs identified in the CIK cells, 16 circRNAs were randomly selected for validation with RT–PCR with divergent primers. The PCR products were further confirmed with Sanger sequencing. The results

Table 1 Statistics of circRNA data. Sample

Raw reads

Raw bases

Clean reads

Clean bases

Valid base

Q30

GC

Sample_CIK Sample_CIK_GCRV

110,762,946 108,210,498

13,845,368,250 13,526,312,250

109,160,724 106,453,214

13,642,692,949 13,303,890,697

98.53% 98.35%

96.94% 96.67%

44.00% 44.00%

Note: Raw reads: the total raw reads of each sample from sequencing. Raw bases: the total bases of each sample from sequencing. Clean reads: the reads from the filtering of the raw data. Clean bases: the total bases of clean reads. Q30(%): the percentage of Phred > 30 the bases in raw data to total number of bases. GC content (%): the percentage of the number of G and C to the total bases number. 3

Aquaculture 258 (2019) 734349

B. Liu, et al.

Fig. 1. Statistical analysis of the circRNAs identified in GCRV-infected CIK cells. (A) Length distributions of circRNAs. (B) Categories of circRNAs. (C) Total circRNAs identified in GCRV-infected and uninfected CIK cells. (D) Splicing signals in the identified circRNAs.

in response to GCRV infection, samples were analyzed with real-time PCR at multiple time points (0, 6, 12, 24, and 48 hpi). The expression of circRNA-3011, circRNA-1561, circRNA-1361, and circRNA-3305 was induced and tended to increase after GCRV infection. The expression of circRNA-1361, circRNA-3186, circRNA-0883, circRNA-0557, and circRNA-1351 was downregulated by GCRV infection (Fig. 3).

showed that the 5′ and 3′ termini of these circRNAs were covalently linked, forming closed loop structures (Table 2). 3.2. Expression of most (> 55%) circRNAs was inhibited by GCRV infection The numbers and distributions of circRNAs at the same level are shown in a volcano plot (Fig. 2A), and were required to satisfy both conditions mentioned above. The differences in expression of the circRNAs in the GCRV-infected and uninfected CIK cells were measured based on RPM (mapped back-splicing junction reads per million mapped reads). In total, 480 and 1925 circRNAs were specifically expressed in the GCRV-infected and uninfected CIK cells, respectively, and 1088 circRNAs were expressed in both groups (Fig. 2B). Seventy-six significantly DE circRNAs (fold change ≥2.0 and p < .05) were identified in the GCRV-infected CIKs relative to the uninfected CIK cells, including 40 upregulated and 36 downregulated circRNAs (Fig. 2C). The top 10 upregulated and downregulated circRNAs in the GCRV-infected cells relative to the uninfected CIK cells were identified based on the fold change and p values for their differential expression, and are shown in Table 3. To evaluate the accuracy of the expression profiling of the circRNAs, four circRNAs that were upregulated and four that were downregulated in the GCRV-infected CIK cells were randomly selected to confirm their expression levels with real-time PCR. The expression levels of these eight circRNAs determined with real-time PCR were consistent with the trends identified with sequencing (Fig. 2D). To determine the dynamic changes in the circRNA expression levels

3.3. More genes were upregulated than downregulated after GCRV infection of CIK cells In summary, 120,511,786 mRNA reads were detected in GCRV-infected CIK cells and 120,415,488 in uninfected cells, resulting in mapped read counts of 118,270,364 and 118,099,556, respectively (Table 4). In total, 20,200 mRNAs were detected in the GCRV-infected and uninfected CIK cells, and 798 DE mRNAs (fold change ≥2.0 and p < .05) were identified in the GCRV-infected CIKs relative to the uninfected cells, including 658 upregulated and 140 downregulated mRNAs (Fig. 4A). The GCRV-infected and uninfected CIK cells both expressed 708 mRNAs, whereas 64 and 26 mRNAs were specifically expressed in the GCRV-infected and uninfected CIK cells, respectively (Fig. 4B). The numbers and distributions of the mRNAs at the same level are shown in a volcano plot, and were required to satisfy both the above conditions (Fig. 4C). To validate the accuracy of the mRNA expression profiling, five mRNAs that were upregulated and four mRNAs that were downregulated in the GCRV-infected CIK cells were randomly selected to confirm their expression levels with real-time PCR. The realtime PCR data for the expression of the nine mRNAs were consistent with the trends determined with sequencing (Fig. 4D). 4

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 2 The sequencing results of circRNAs parental gene sequence and junction site.

(continued on next page) 5

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 2 (continued)

(continued on next page)

6

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 2 (continued)

(continued on next page)

7

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 2 (continued)

Fig. 2. Identification of circRNAs differentially expressed in GCRV-infected and uninfected CIK cells. (A) Volcano plot shows the numbers and distributions of circRNAs in the same plane. (B) Numbers of total and unique circRNAs in GCRV-infected and uninfected CIK cells. (C) Differentially regulated circRNAs detected with whole-transcriptome sequencing (fold change ≥2.0 and p < .05). (D) Validation of circRNA expression with real-time PCR. 8

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 3 The top 10 upregulated and downregulated circRNAs. id

log2FC

PValue

FDR

up_down

circRNA_chr

circRNA_strand

circRNA_start

circRNA_end

type

circRNA_3343 circRNA_2836 circRNA_0999 circRNA_1785 circRNA_2318 circRNA_1563 circRNA_0160 circRNA_0601 circRNA_3065 circRNA_0182 circRNA_1969 circRNA_0231 circRNA_0543 circRNA_2568 circRNA_0557 circRNA_1373 circRNA_1587 circRNA_0794 circRNA_3296 circRNA_1269

3.150665729 3.118210964 2.845111297 2.830824631 2.830824631 2.799403209 2.65535671 2.60898282 2.60898282 2.389551325 −1.947157613 −2.010827577 −2.130311197 −2.142992553 −2.186536039 −2.224695003 −2.224695003 −2.511850893 −2.511850893 −2.575520857

0.000145932 0.00330227 9.67E-08 0.0046627 0.0046627 0.016344234 0.001118588 0.027924617 0.027924617 0.047691132 0.030086015 0.030086015 0.014317568 0.032029795 0.009835004 0.020829363 0.020829363 0.014149991 0.014149991 0.009462925

0.026462 0.256634 0.000105 0.268724 0.268724 0.467961 0.152128 0.522921 0.522921 0.627771 0.522921 0.522921 0.432709 0.52488 0.365752 0.472132 0.472132 0.432709 0.432709 0.365752

Up Up Up Up Up Up Up Up Up Up Down Down Down Down Down Down Down Down Down Down

+ + − − − + + + − + + − − + − − − − − −

1,817,477 6,524,632 5,449,888 7,499,448 496,522 2,678,807 10,749,174 8,336,335 249,800 2,070,747 4,165,966 12,079,295 9,905,064 1,755,955 13,544,723 6,817,079 2,795,308 4,757,042 1,539,718 4,011,837

1,825,731 6,529,745 5,452,908 7,501,055 500,193 2,679,743 10,750,257 8,337,475 258,464 2,071,156 4,175,000 12,087,232 9,906,253 1,757,439 13,545,167 6,823,650 2,813,190 4,773,611 1,542,473 4,014,631

8255 5114 3021 1608 3672 937 1084 1141 8665 410 9035 7938 1190 1485 445 6572 17,883 16,570 2756 2795

Sense-overlapping Sense-overlapping Sense-overlapping Sense-overlapping Sense-overlapping Sense-overlapping Sense-overlapping Exonic Sense-overlapping Exonic Sense-overlapping Sense-overlapping Sense-overlapping Intergenic Sense-overlapping Sense-overlapping Sense-overlapping Sense-overlapping Sense-overlapping Sense-overlapping

Fig. 3. Expression levels of circRNAs in GCRV-infected CIK cells at different times after infection. 9

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 4 Statistics of mRNAs data. Sample

raw bases

clean reads

clean bases

valid bases

Q30

GC

GCRV infected CIK cells CIK cells

17,736,146,140 17,713,649,550

118,270,364 118,099,556

17,391,501,804 17,357,812,283

98.06% 97.99%

95.18% 95.06%

48% 49%

3.4. Upregulated miRNAs exceeded downregulated miRNAs in CIK cells infected by GCRV

Among these, 236 and 234 known miRNAs were identified in the GCRV-infected and uninfected CIK cells, respectively, and 222 known miRNAs were common to both groups. In the two groups, 1317 and 1318 miRNAs were predicted to be novel miRNAs, respectively, and 1198 novel miRNAs were common to both groups. Therefore, 1685 miRNAs were identified in the grass carp, 248 of which were already known and 1437 were novel miRNAs (Fig. 5B). According to the screening criteria for DE miRNAs, 186 DE miRNAs were identified when the GCRV-infected and uninfected CIK cells were compared, 146 of which were upregulated and 40 downregulated (p value < .05 and fold change > 2) (Fig. 5C). To determine the accuracy of the miRNA expression profiling, four upregulated and four downregulated miRNAs in the GCRV-infected CIK cells were randomly selected to validate their expression levels with real-time PCR. The real-time PCR data for the eight miRNAs were consistent with the trends determined with sequencing (Fig. 5D).

After the depletion of tRNA, small nuclear RNA (snRNA), rRNA, and low-quality sequences, 37,107,560 and 30,762,142 clean reads were acquired from the raw reads from the GCRV-infected and uninfected CIK cells (Table 5). When the length distributions of the clean reads from the two groups were analyzed, 29.7% and 17.8% of clean reads were 22 bp and 23 bp, respectively, in the GCRV-infected CIK cells, and 43.2% and 27.2% of clean reads were 22 bp and 23 bp, respectively, in the uninfected CIK cells (Fig. 5A), After quality control of the small RNA sequence data, the clean reads were analyzed with a BLAST search in miRBase to match them to known miRNAs, and the unmapped sequences (length > 18 bp) were analyzed with a BLAST search against the grass carp genome. The secondary structures of the miRNAs that were mapped to the grass carp genome were predicted with RNAfold and the sequences of the miRNA hairpin precursors formed from these sequences were considered to be novel miRNAs. In total, 1553 miRNAs were predicted in the GCRV-infected CIK cells and 1552 miRNAs in the uninfected cells, and 1420 miRNAs were shared by the two groups.

Fig. 4. Identification of differentially expressed mRNAs in GCRV-infected and uninfected CIK cells. (A) Differentially regulated mRNAs detected with whole-transcriptome sequencing (fold change ≥2.0 and p < .05). (B) Numbers of total and unique mRNAs in GCRV-infected and uninfected CIK cells. (C) Volcano plot shows the numbers and distributions of mRNAs in the same plane. (D) Validation of mRNA expression with real-time PCR. 10

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 5 Statistics of miRNAs data. Sample

Raw reads

Reads trimmed length

Reads rimmed Q20

Reads trimmed N

Clean reads

Unique reads

GCRV infected CIK cells CIK cells

37,769,394 31,428,967

37,218,029 30,857,002

37,189,080 30,833,865

37,107,560 30,762,142

37,107,560 30,762,142

879,575 638,992

Fig. 5. Identification of miRNAs in GCRV-infected and uninfected CIK cells. (A) Length distributions of miRNAs identified in GCRV-infected and uninfected CIK cells. (B) Numbers of total, unique, and shared circRNAs in GCRV-infected and uninfected CIK cells. (C) Differentially regulated miRNAs detected with small RNA sequencing (fold change ≥2.0 and p < .05). (D) Validation of miRNA expression with real-time PCR.

3.5. Parental genes of dysregulated circRNAs were enriched in GO items ‘calcitonin receptor activity’ and ‘metal ion binding’ and KEGG pathways ‘spliceosome’ and ‘endocytosis’

3.6. Dysregulated mRNAs enriched in GO items ‘innate immune response’ and ‘regulation of apoptotic process’ and KEGG pathways ‘TNF signaling’ and ‘steroid biosynthesis’

To determine the potential roles of the DE circRNAs induced in the CIK cells by GCRV infection, the parental genes of these DE circRNAs were analyzed with GO annotation and a KEGG pathway analysis. The top 20 GO processes of the parental genes of the dysregulated circRNAs were analyzed with routine GO classification algorithms. Many of these parental genes were enriched in the GO terms ‘response to toxic substance’, ‘calcitonin receptor activity’, ‘cAMP response element binding’, ‘protein binding’, and ‘lysosome’ (Table 6). A total of 21 KEGG pathways (p < .05) were identified for the upregulated and downregulated circRNAs. The top 20 KEGG pathways included ‘drug metabolism–cytochrome P450’, ‘spliceosome’, and ‘RNA degradation signaling pathway’ (Table 7).

To clarify the potential roles of the DE mRNAs induced in the CIK cells after GCRV infection, these DE mRNAs were analyzed with GO annotation and a KEGG pathway analysis. Based on the −log10 p values, 52 and 15 terms were related to upregulated and downregulated mRNAs, respectively. Most of the upregulated mRNAs were enriched in ‘transmembrane transport’, ‘inflammatory response’, ‘regulation of apoptotic process’, ‘innate immune response’, or ‘apoptotic process’ (Fig. 6A). Most of the downregulated mRNAs were enriched in ‘cholesterol biosynthetic process’, ‘cell division’, ‘DNA repair’, ‘chromatin’, ‘cytosol’, ‘endoplasmic reticulum membrane’, or ‘iron ion binding’ (Fig. 6B). Sixty-six and four KEGG pathways (p < .05) were identified among the upregulated and downregulated mRNAs, respectively. The upregulated mRNAs were enriched in immunity pathways and signal transduction pathways, such as the Toll-like receptor signaling pathway, Toll and IMD signaling pathways, T-cell receptor signaling pathway, TNF signaling pathway, HIF-1 signaling pathway, and JAK–STAT signaling pathway (Fig. 6C). The downregulated mRNAs 11

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 6 GO enrichment for differential expressed circRNAs with their parental genes. id

Term

Category

ListHits

ListTotal

PopHits

PopTotal

pval

GO:0003835 GO:0009311 GO:0004063 GO:0004064 GO:0009636 GO:0005635 GO:0006606 GO:0008140 GO:0032793 GO:0004385 GO:0007268 GO:0035255 GO:0043113 GO:0050808 GO:0097120 GO:0004948 GO:0007166 GO:0014069 GO:0051289 GO:0006683 GO:0005764 GO:0004336 GO:0045002 GO:0043142 GO:0043139 GO:0043137 GO:0033567 GO:0017108 GO:0016890 GO:0008409 GO:0006264 GO:0003678 GO:0000076 GO:0007589 GO:0006611 GO:0005049 GO:0010923 GO:0005975 GO:0006284 GO:0000729 GO:0046530 GO:0039023 GO:0035845 GO:0008594 GO:0043197 GO:0031527 GO:0016327 GO:0016323 GO:0042073 GO:0018279 GO:0045211 GO:0004518 GO:0030672 GO:0030992 GO:0032580 GO:0045740 GO:0051539 GO:0042384 GO:0006260 GO:0007010 GO:0016459 GO:0003774 GO:0007224 GO:0008360 GO:0030424 GO:0005929

Beta-galactoside alpha-2,6-sialyltransferase activity Oligosaccharide metabolic process Aryldialkylphosphatase activity Arylesterase activity Response to toxic substance Nuclear envelope Protein import into nucleus Camp response element binding protein binding Positive regulation of CREB transcription factor activity Guanylate kinase activity Synaptic transmission Ionotropic glutamate receptor binding Receptor clustering Synapse organization Receptor localization to synapse Calcitonin receptor activity Cell surface receptor signaling pathway Postsynaptic density Protein homotetramerization Galactosylceramide catabolic process Lysosome Galactosylceramidase activity Double-strand break repair via single-strand annealing Single-stranded DNA-dependent atpase activity 5′-3’ DNA helicase activity DNA replication, removal of RNA primer DNA replication, Okazaki fragment processing 5′-flap endonuclease activity Site-specific endodeoxyribonuclease activity, specific for altered base 5′-3′ exonuclease activity Mitochondrial DNA replication DNA helicase activity DNA replication checkpoint Body fluid secretion Protein export from nucleus Nuclear export signal receptor activity Negative regulation of phosphatase activity Carbohydrate metabolic process Base-excision repair DNA double-strand break processing Photoreceptor cell differentiation Pronephric duct morphogenesis Photoreceptor cell outer segment organization Photoreceptor cell morphogenesis Dendritic spine Filopodium membrane Apicolateral plasma membrane Basolateral plasma membrane Intraciliary transport Protein N-linked glycosylation via asparagine Postsynaptic membrane Nuclease activity Synaptic vesicle membrane Intraciliary transport particle B Golgi cisterna membrane Positive regulation of DNA replication 4 iron, 4 sulfur cluster binding Cilium assembly DNA replication Cytoskeleton organization Myosin complex Motor activity Smoothened signaling pathway Regulation of cell shape Axon Cilium

Molecular_function Biological_process Molecular_function Molecular_function Biological_process Cellular_component Biological_process Molecular_function Biological_process Molecular_function Biological_process Molecular_function Biological_process Biological_process Biological_process Molecular_function Biological_process Cellular_component Biological_process Biological_process Cellular_component Molecular_function Biological_process Molecular_function Molecular_function Biological_process Biological_process Molecular_function Molecular_function Molecular_function Biological_process Molecular_function Biological_process Biological_process Biological_process Molecular_function Biological_process Biological_process Biological_process Biological_process Biological_process Biological_process Biological_process Biological_process Cellular_component Cellular_component Cellular_component Cellular_component Biological_process Biological_process Cellular_component Molecular_function Cellular_component Cellular_component Cellular_component Biological_process Molecular_function Biological_process Biological_process Biological_process Cellular_component Molecular_function Biological_process Biological_process Cellular_component Cellular_component

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 6 7 8 9 9 13 13 13 14 15 16

701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701 701

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.37E-06 6.22E-05 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.000371 0.0011 0.0011 0.0011 0.0011 0.0011 0.0011 0.0011 0.0011 0.002175 0.002175 0.002175 0.002175 0.002175 0.002175 0.003583 0.003583 0.003583 0.003583 0.003583 0.005313 0.007354 0.009693 0.01232 0.01232 0.025496 0.025496 0.025496 0.029406 0.033544 0.0379

were only enriched in four pathways: ‘steroid biosynthesis’, ‘biosynthesis’, ‘cell cycle’, and ‘FOXO signaling pathway’ (Fig. 6D).

3.7. DE genes regulated by miRNA–mRNA axis were enriched in ‘hematopoietic cell lineage’, ‘AMPK signaling pathway’, and ‘PI3K–AKT signaling pathway’ To understand the functions of the DE miRNAs, the miRanda algorithm was used to predict the miRNA targets, and predicted 1972 target 12

Aquaculture 258 (2019) 734349

B. Liu, et al.

Table 7 KEGG enrichment for differential expressed circRNAs with their parental genes. id

term

class

ListHits

ListTotal

PopHits

PopTotal

pval

ko04142 ko04742 ko04974 ko04920 ko03018

Lysosome Taste transduction Protein digestion and absorption Adipocytokine signaling pathway RNA degradation

1 1 1 1 1

19 19 19 19 19

18 18 14 12 11

949 949 949 949 949

0.048025 0.048025 0.029958 0.022253 0.018767

ko00600 ko00500 ko03030 ko00140 ko00983 ko00040 ko00053 ko00830 ko00860 ko03040 ko04112 ko00604 ko00363 ko00627 ko00980

Sphingolipid metabolism Starch and sucrose metabolism DNA replication Steroid hormone biosynthesis Drug metabolism - other enzymes Pentose and glucuronate interconversions Ascorbate and aldarate metabolism Retinol metabolism Porphyrin and chlorophyll metabolism Spliceosome Cell cycle - Caulobacter Glycosphingolipid biosynthesis - ganglio series Bisphenol degradation Aminobenzoate degradation Metabolism of xenobiotics by cytochrome P450 Drug metabolism - cytochrome P450

Cellular Processes/Transport and catabolism Organismal Systems/Sensory system Organismal Systems/Digestive system Organismal Systems/Endocrine system Genetic Information Processing/Folding, sorting and degradation Metabolism/Lipid metabolism Metabolism/Carbohydrate metabolism Genetic Information Processing/Replication and repair Metabolism/Lipid metabolism Metabolism/Xenobiotics biodegradation and metabolism Metabolism/Carbohydrate metabolism Metabolism/Carbohydrate metabolism Metabolism/Metabolism of cofactors and vitamins Metabolism/Metabolism of cofactors and vitamins Genetic Information Processing/Transcription Cellular Processes/Cell growth and death Metabolism/Glycan biosynthesis and metabolism Metabolism/Xenobiotics biodegradation and metabolism Metabolism/Xenobiotics biodegradation and metabolism Metabolism/Xenobiotics biodegradation and metabolism

1 1 1 1 1 1 1 1 1 2 1 1 1 1 1

19 19 19 19 19 19 19 19 19 19 19 19 19 19 19

11 5 4 4 3 2 2 2 2 6 1 1 1 1 1

949 949 949 949 949 949 949 949 949 949 949 949 949 949 949

0.018767 0.003667 0.002227 0.002227 0.001127 0.00038 0.00038 0.00038 0.00038 0.000131 0 0 0 0 0

Metabolism/Xenobiotics biodegradation and metabolism

1

19

1

949

0

ko00982

genes for the 186 DE miRNAs. Most of the target genes were related to the immune and inflammation responses, such as interleukin 6 (IL-6), interferon regulatory factor 2 (IRF2), and complement receptor type 1 (CR1) (Supplementary Table 2). Furthermore, 1972 target genes were used in a functional analysis, and the 1136 GO terms obtained were classified into the three categories ‘biological process’ (BP), ‘cellular

component’ (CC) and ‘molecular function’ (MF), with 706, 201, and 229 of the 1136 GO terms enriched in each, respectively. The top 10 of terms in BP included ‘transcription’, ‘regulation of transcription’, and ‘apoptotic process’. The top 10 terms in CC included ‘nucleus’, ‘cytoplasm’, and ‘plasma membrane’. The top 10 terms in MF included ‘metal ion binding’, ‘ATP binding’, and ‘DNA binding’ (Fig. 7A). When the

Fig. 6. GO and KEGG pathway analyses of differentially expressed mRNAs. (A) GO terms for upregulated mRNAs. (B) GO terms for downregulated mRNAs. (C) KEGG pathways for upregulated mRNAs. (D) KEGG pathways for downregulated mRNAs. 13

Aquaculture 258 (2019) 734349

B. Liu, et al.

Fig. 7. GO and KEGG analyses of differentially expressed (DE) miRNAs. (A) GO terms for miRNA targets. (B) KEGG pathways for DE miRNAs. (C) Network of miRNAs and mRNAs. (D) KEGG enrichment analysis of DE mRNAs of miRNA targets.

KEGG pathways of the target genes of the miRNAs were analyzed, the top 20 enriched KEGG pathways included ‘HTLV-I infection’ and ‘microRNAs in cancer’ (Fig. 7B). To understand the miRNA regulatory mechanisms underlying the response to GCRV infection, the interactions between the miRNAs and their target mRNAs were investigated. Based on the regulatory mechanisms of miRNAs and their target genes (mRNAs), 43 mRNAs (39 upregulated and four downregulated) were identified as the targets of 23 miRNAs (15 upregulated and eight downregulated) (Fig. 7C). The miRNA–mRNA network constructed showed that one miRNA had several target mRNAs, or one miRNA had only one target mRNA. A KEGG pathway analysis showed that these 43 target mRNAs were enriched in the ‘glucagon signaling pathway’, ‘AMPK signaling pathway’, ‘PI3K–AKT signaling pathway’ etc. (Fig. 7D).

(Fig. 8). Thyroglobulin (THYG) gene expression was mediated by the interaction of cik-miR-5272* with circRNA_3415, and by the interaction of cik-miR-10719* with circRNA_3065 and circRNA_3011 together. The expression of the ephrin type-B receptor 2 (EPHB2), apoptosis facilitator Bcl-2-like protein 14 (B2L14), membrane frizzled-related protein (MFRP), serine/threonine-protein kinase SIK2 (SIK2), POU domain, class 2, transcription factor 2 (PO2F2), and transient receptor potential cation channel subfamily M member 5 (TRPM5) genes was regulated by the sponging of circRNA_2528/circRNA_1315 to cik-miRNA-10648. The expression of AMP deaminase 3 (AMPD3) was mediated by cik-miR8141 through circRNA-1591/circRNA-1602. The circRNA-2780/cikmiRNA-7796/KAT6B and circRNA-1351/cik-miRNA-10648/B2L14 axes may be key regulatory pathways in viral infection. These results indicate that there are complex regulatory relationships among circRNAs, miRNAs, and mRNAs (Table 8), and suggest that circRNAs act in ceRNAs to regulate mRNAs by sponging their shared miRNAs during viral infection.

3.8. Construction of ceRNA network To specify the key circRNA–miRNA–mRNA regulatory modules involved in the process of GCRV infection, 76 DE circRNAs, 786 DE mRNAs, and 186 DE miRNAs were integrated into a ceRNA (circRNA–miRNA–mRNA) regulatory network. The network was composed of nine circRNA–miRNA pairs and 18 miRNA–mRNA pairs, which included nine DE circRNAs, seven DE miRNAs, and 17 DE mRNAs

4. Discussion In this study, whole-transcriptome and small RNA sequencing were used to systematically investigate the responses of circRNA, mRNA, and miRNA expression profiles to GCRV infection in CIK cells. After a series 14

Aquaculture 258 (2019) 734349

B. Liu, et al.

Fig. 8. Interaction network of differentially expressed (DE) circRNAs–DE miRNAs–DE mRNAs.

of strict selection processes, 3494 circRNAs, 20,200 mRNAs, and 1553 miRNAs were identified in CIK cells infected with GCRV and uninfected CIK cells. Among these, 1568 and 3013 circRNAs were identified in the GCRV-infected and uninfected CIK cells. The huge number of differences in the circRNAs expressed in these two groups may be attributable to the nuclear export of NF90/NF110 to the cytoplasm during the process of viral infection and the inhibition of viral replication by the NF90/NF110 released from circRNP complexes, to bind viral mRNAs (Li et al., 2017). Interestingly, fewer circRNAs were identified in the CIK cells than in the spleen cells of the grass carp (5052 circRNAs), suggesting that the expression of circRNAs has inherent features, such as temporal and tissue-specific expression patterns (Jeck et al., 2013). Previous studies have demonstrated that many dysregulated circRNAs are involved in many biological events. For example, the circRNA expression patterns in the silkworm were altered after

infection with BmNPV or BmCPV. These dysregulated circRNAs may be associated with viral infection and many immunity-related genes are regulated by these circRNAs (Hu et al., 2018a; Hu et al., 2018b). After GCRV infection, the circRNA expression profiles in the grass carp spleen were dysregulated and the functions of the target miRNAs of these dysregulated circRNAs may involve endothelial and blood cell damage and hemorrhagic symptoms (He et al., 2017a). However, the regulatory mechanisms involving circRNAs, mRNAs, and miRNA during GCRV infection and their roles in the hemorrhagic symptoms of GCRV-infected fish are still unclear. The functions of many circRNAs in diseases of human and other species have been comprehensively analyzed. Not only can these circRNAs be used as biomarkers for the diagnosis of disease (Cui et al., 2016; Ouyang et al., 2018; Yuan et al., 2018), but also play important roles in the pathogenesis of diseases (Jin et al., 2016; Zhang et al., 2017b; Cheng et al., 2018; Lin and Chen, 2018;

Table 8 The interaction network of circRNA-host miRNA-mRNA. circRNA

miRNA

Target gene

Target gene Description

Target gene ID

circRNA_0988

cik-miR-3526*

circRNA_1351 circRNA_2528

cik-miR-10,648

circRNA_1591

cik-miR-6846 cik-miR-8141

GFOD1 NLRP3 B2L14 EPHB2 MFRP PO2F2 SIK2 TRPM5 TGBR3 C8AP2 AMPD3 NLRP3 KAT6B THYG

Glucose-fructose oxidoreductase domain-containing protein 1 NACHT, LRR and PYD domains-containing protein 3 Apoptosis facilitator Bcl-2-like protein 14 Ephrin type-B receptor 2 Membrane frizzled-related protein POU domain, class 2, transcription factor 2 Serine/threonine-protein kinase SIK2 Transient receptor potential cation channel subfamily M member 5 Transforming growth factor beta receptor type 3 CASP8-associated protein 2 AMP deaminase 3 NACHT, LRR and PYD domains-containing protein 3 Histone acetyltransferase KAT6B Thyroglobulin

CI01000265_00039287_00064707 CI01000167_00131588_00209846 CI01000013_05905378_05906438 CI01000001_07137510_07187079 CI01000092_04413965_04426684 CI01000319_05269495_05286542 CI01000037_03946645_03954021 CI01000026_16965731_16979755 CI01000340_00927319_01001912 CI01000325_03637765_03645294 CI01000013_06038814_06043689 CI01000152_01042484_01077436 CI01000012_12381049_12417489 CI01000027_09002993_09040903

EDRF1 PPTC7 WC11

Erythroid differentiation-related factor 1 Protein phosphatase PTC7 homolog Antigen WC1.1

CI01000310_03054431_03066827 CI01000010_09902329_09911570 CI01073364_00000079_00000892

circRNA_1602 circRNA_2780 circRNA_3011 circRNA_3065 circRNA_3415

cik-miR-7796 cik-miR-10,719 cik-miR-5272*

15

Aquaculture 258 (2019) 734349

B. Liu, et al.

Zhang et al., 2018a). circRNAs are transcribed from their parental genes, and analyzing the roles of these parental genes and their related circRNAs should allow a comprehensive understanding of their functions. In this study, the parental genes of the dysregulated circRNAs were determined for GO annotation and a KEGG signaling pathway analysis. Interestingly, many parental genes of the DE circRNAs were related to ‘ATP binding’, ‘metal ion binding’, ‘DNA binding’, ‘regulation of transcription’, ‘regulation of transcription’, and ‘apoptotic process’. Metal ions are important immune factors in the defense against pathogen invasion (Jesse et al., 2014). In the present study, a GO enrichment analysis of the linear transcripts of the dysregulated circRNAs suggested that the metal ion balance is disrupted in cells after GCRV infection. A KEGG signaling pathways analysis of the linear transcripts showed that most of the linear transcripts related to dysregulated circRNAs were enriched in ‘spliceosome’, ‘endocytosis’, and ‘cell cycle’. The roles of circRNAs that contribute to viral pathogenesis or the infection process may be facilitated by miRNA-mediated mRNA interactions. In this study, to fully understand the roles of circRNAs in the CIK cell response to GCRV infection, 798 DE mRNAs and 186 DE miRNAs (fold change ≥2.0 and p < .05) were identified in the GCRVinfected CIKs relative to uninfected control cells. Interesting, 47.5% and 70.4% of the miRNAs had lengths of 22 bp and 23 bp, respectively, in the GCRV-infected CIK cells and uninfected CIK cells. This suggests that miRNA biosynthesis was inhibited by GCRV infection, which may contribute to the replication of GCRV in susceptible cells. miRNA is a kind of noncoding RNA that has been highly conserved during evolution, and the targeting of an mRNA by an miRNA can be blocked by a single-nucleotide changes in the RNA target. Therefore, it is possible that viruses have evolved to be largely refractory to inhibition by endogenous miRNAs (Cullen, 2013). The DE mRNAs and miRNAs identified in the CIK response to GCRV infection were analyzed for their potential functions. Most of the upregulated mRNAs were enriched in immunity-related GO terms and most of the downregulated mRNAs were enriched in metabolism-related GO terms. The upregulated mRNAs were also enriched in immunity-related signaling pathways and the downregulated mRNAs were mainly enriched in biosynthesis-related signaling pathway. These data suggest that the cellular immunity and metabolism are greatly altered by GCRV infection. It is possible that these changes in cells after GCRV infection are involved in defending the cells against viral infection. These important immunity pathways, including the Toll signaling pathway and the T-cell receptor signaling pathway, are also enriched in the large yellow croaker during Aeromonas hydrophila infection (Mu et al., 2010). GCRV infection is a gradual process, involving the adsorption of the virus to the cell surface, followed by its endocytosis into the cell, and its transport by lysosomes, ultimately resulting in cell necrosis and/or apoptosis (Chen et al., 2018b). It is known that fish have complete innate and adaptive immune systems. Nine significantly upregulated genes were enriched in the Toll-like receptor signaling pathway, and myeloid differentiation factor 88 (MYD88) is also reportedly upregulated as an antiviral factor in the spleen of the grass carp (Su et al., 2011). Four DE genes were enriched in the T-cell receptor signaling pathway, which contain E3 ubiquitin protein ligase (CBL), phosphatidylinositol 3-kinase (PI3K), mitogen-activated protein kinase kinase 7 (MKK7), and activator protein 1 (AP1). These results suggest that CBL, PI3K, MKK7, and AP1 perform specific functions in CIK cells after GCRV infection. Other important enriched signaling pathways included the TNF signaling pathway, MAPK signaling pathway, JAK–STAT signaling pathway, and HIPPO signaling pathway, all of which are associated with growth, development, proliferation, differentiation, apoptosis, and disease. The target genes of miRNAs were analyzed in a KEGG signaling pathway analysis, and the results showed that the HTLV-I infection and miRNAs in cancer were enriched. Furthermore, 43 target mRNAs that interacted with 23 miRNAs were enriched in the glucagon signaling

pathway, AMPK signaling pathway, PI3K–AKT signaling pathway, etc. These findings indicate that the roles of these immunity pathways are conserved among different species and play important roles in the process of GCRV infection. It has also been reported that when the grass carp spleen is infected with GCRV, the targets of DE miRNAs are predominantly mRNAs associated with the immune response, blood coagulation, hemostasis, and the complement and coagulation cascades (He et al., 2017b). One of the most important peculiarities of circRNAs is its role as an “miRNA sponge”. It can recruit, bind, and regulate miRNAs, further influencing the expression of their target genes, by competing for miRNA sites (Hansen et al., 2013). ceRNA is strongly implicated in many important biological processes (Mitra et al., 2018). To further investigate the functions of the dysregulated circRNAs identified in this study, the DE mRNAs and miRNAs were identified with sequencing and the targets of the miRNAs and circRNAs were predicted with the MiRanda software. In total, 76 DE circRNAs, 798 mRNAs, and 186 miRNAs were differentially expressed in the GCRV-infected CIK cells relative to the uninfected cells. With strict selection and a correlation analysis, a ceRNA network was constructed with these dysregulated RNAs. The network was composed of nine circRNA–miRNA pairs and 18 miRNA–mRNA pairs, which included nine DE circRNAs, seven DE miRNAs, and 17 DE mRNAs. According to the regulatory mechanisms of the circRNA/miRNA and miRNA/mRNA pairs, individual circRNAs sponged to several miRNAs and miRNAs had more than one mRNA target. Here, THYG was mediated by cik-miR-5272* and cik-miR10719* together. cik-miRNA-10648 had six targets (EPHB2, B2L14, MFRP, SIK2, PO2F2, and TRPM5 mRNAs). cik-miR-10719* was sponged by two miRNAs (circRNA-3011 and circRNA-3065). The expression of AMP deaminase 3 (AMPD3) was mediated by cik-miR-8141 through circRNA-1591 and circRNA-1602. The circRNA-2780/cikmiRNA-7796/KAT6B and circRNA-1351/cik-miRNA-10648/B2L14 axes may be key regulatory pathways in viral infection. AMP deaminase (AMPD) plays an important role in purine/urate metabolism, catalyzing AMP to inosine monophosphate. Skeletal muscle and red blood cell deficiencies are caused by AMPD1 and AMPD3 deficiencies, respectively, and AMPD1 deficiency had been detected in patients with metabolic myopathy (Morisaki and Morisaki, 2008). Therefore, the hemorrhagic symptoms and reduced movement capacity of GCRVinfected fish may be associated with the circRNA-1591/circRNA-1602/ cik-miR-8141/AMPD3 regulatory axis. circRNA-0988 is predicted to interact with cik-miR-3526, and the expression of cik-miR-3526 was negative with its targets glucose–fructose oxidoreductase domain-containing protein 1 (GFOD1), NACHT, LRR, and PYD domain-containing protein 3 (NLRP3). NLRP3 is the main component of the inflammasome, which acts as the activation platform for caspase 1 and promotes the maturation of the inflammatory cytokines IL-1β, IL-18, and IL-33. The activation of the NLRP3 inflammasome can be protective during infection, but an unintended pathology can occur in response to nonpathogenic endogenous or exogenous stimuli when the activation of the NLRP3 inflammasome is unregulated (Elliott and Sutterwala, 2015). In the process of GCRV infection, upregulated circRNA-0988 and cik-miR3526 promote the expression of NLRP3, activate the inflammasome, and initiate the apoptosis to inhibit GCRV infection. GFOD1 plays important roles in energy metabolism, and the circRNA-0988/cik-miR3526/GFOD1 axis is associated with its expression. These findings are consistent with the GO terms predicted in this study. The ceRNA analysis identified several dysregulated circRNAs associated with cell apoptosis. This regulatory network includes circRNA-1351/cik-miR1064/B2L4 (apoptosis facilitator Bcl-2-like protein 14), circRNA-1591/ cik-miR-6846/transforming growth factor beta receptor type 3 (TGBR3), and circRNA-1591/cik-miR-6846/CASP8-associated protein 2 (C8AP2), which may play important roles in the defense against the GCRV infection through circRNA-1351-mediated cellular apoptosis. The expression patterns of genes in GCRV-infected fish were altered relative to those in the control fish. PO2F2 is a transcription factor, and 16

Aquaculture 258 (2019) 734349

B. Liu, et al.

it expression was significantly upregulated by GCRV infection. In this study, the expression of PO2F2 was regulated by circRNA-1351 and circRNA-2528, and we predict that many genes are regulated by these two circRNAs after GCRV infection. The expression of erythroid differentiation-related factor 1 (EDRF1), which is involved in the differentiation of red blood cells, was upregulated in the GCRV-infected grass carp, so we predict that the induced expression of EDRF1 is associated with the hemorrhagic symptoms of infected fish. From our analysis of the potential roles of these dysregulated circRNAs, mRNAs, and miRNAs, we predict that the hemorrhagic symptoms of fish infected with GCRV are related with these dysregulated RNAs, and that knowledge of miRNA-mediated circRNAs–mRNAs will allow us to comprehensively understand the mechanisms of GCRV infection and its pathogenesis.

major depressive disorder. Biomark. Med 10, 943–952. Cullen, B.R., 2013. How do viruses avoid inhibition by endogenous cellular microRNAs? PLoS Pathog. 9, e1003694. Dai, Z., Li, J.C., Hu, C.Y., Wang, F., Wang, B.H., Shi, X., Hou, Q.H., Huang, W.G., Lin, G., 2017. Transcriptome data analysis of grass carp (Ctenopharyngodon idella) infected by reovirus provides insights into two immune-related genes. Fish & shellfish immunology 64, 68–77. Elliott, E.I., Sutterwala, F.S., 2015. Initiation and perpetuation of NLRP3 inflammasome activation and assembly. Immunol. Rev. 265, 35–52. Gao, Y., Wang, J.F., Zhao, F.Q., 2015. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 16. Gao, Y., Pei, C., Sun, X., Zhang, C., Li, L., Kong, X., 2018. Novel subunit vaccine based on grass carp reovirus VP35 protein provides protective immunity against grass carp hemorrhagic disease. Fish & shellfish immunology 75, 91–98. Ghorbani, A., Izadpanah, K., Peters, J.R., Dietzgen, R.G., Mitter, N., 2018. Detection and profiling of circular RNAs in uninfected and maize Iranian mosaic virus-infected maize. Plant science : an international journal of experimental plant biology 274, 402–409. Hansen, T.B., Jensen, T.I., Clausen, B.H., Bramsen, J.B., Finsen, B., Damgaard, C.K., Kjems, J., 2013. Natural RNA circles function as efficient microRNA sponges. Nature 495, 384–388. He, L., Zhang, A., Xiong, L., Li, Y., Huang, R., Liao, L., Zhu, Z., Wang, A.Y., 2017a. Deep circular RNA sequencing provides insights into the mechanism underlying grass carp reovirus infection. Int. J. Mol. Sci. 18. He, L.B., Zhang, A.D., Chu, P.F., Li, Y.M., Huang, R., Liao, L.J., Zhu, Z.Y., Wang, Y.P., 2017b. Deep Illumina sequencing reveals conserved and novel microRNAs in grass carp in response to grass carp reovirus infection. BMC Genomics 18. He, L.B., Zhang, A.D., Pei, Y.Y., Chu, P.F., Li, Y.M., Huang, R., Liao, L.J., Zhu, Z.Y., Wang, Y.P., 2017c. Differences in responses of grass carp to different types of grass carp reovirus (GCRV) and the mechanism of hemorrhage revealed by transcriptome sequencing. BMC Genomics 18. Hu, X., Zhu, M., Liu, B., Liang, Z., Huang, L., Xu, J., Yu, L., Li, K., Jiang, M., Xue, R., Cao, G., Gong, C., 2018a. Circular RNA alterations in the Bombyx mori midgut following B. mori nucleopolyhedrovirus infection. Mol. Immunol. 101, 461–470. Hu, X., Zhu, M., Zhang, X., Liu, B., Liang, Z., Huang, L., Xu, J., Yu, L., Li, K., Zar, M.S., Xue, R., Cao, G., Gong, C., 2018b. Identification and characterization of circular RNAs in the silkworm midgut following Bombyx mori cytoplasmic polyhedrosis virus infection. RNA Biol. 15, 292–301. Jeck, W.R., Sorrentino, J.A., Wang, K., Slevin, M.K., Burd, C.E., Liu, J., Marzluff, W.F., Sharpless, N.E., 2013. Circular RNAs are abundant, conserved, and associated with ALU repeats. Rna 19, 141–157. Jeong, J.M., Jeswin, J., Bae, J.S., Park, C.I., 2018. Differential expression of olive flounder (Paralichthys olivaceus) transcriptome during viral hemorrhagic septicemia virus (VHSV) infection at warmer and colder temperature. Data in brief 19, 2023–2028. Jesse, H.E., Roberts, I.S., Cavet, J.S., 2014. Metal ion homeostasis in Listeria monocytogenes and importance in host-pathogen interactions. Adv. Microb. Physiol. 65, 83–123. Ji, J., Rao, Y., Wan, Q., Liao, Z., Su, J., 2018. Teleost-specific TLR19 localizes to endosome, recognizes dsRNA, recruits TRIF, triggers both IFN and NF-kappaB pathways, and protects cells from grass carp Reovirus infection. J. Immunol. 200, 573–585. Jin, X., Feng, C.Y., Xiang, Z., Chen, Y.P., Li, Y.M., 2016. CircRNA expression pattern and circRNA-miRNA-mRNA network in the pathogenesis of nonalcoholic steatohepatitis. Oncotarget 7, 66455–66467. Jin, S., Zhao, X., Wang, H., Su, J., Wang, J., Ding, C., Li, Y., Xiao, T., 2018. Molecular characterization and expression of TLR7 and TLR8 in barbel chub (Squaliobarbus curriculus): responses to stimulation of grass carp reovirus and lipopolysaccharide. Fish & shellfish immunology 83, 292–307. John, B., Enright, A.J., Aravin, A., Tuschl, T., Sander, C., Marks, D.S., 2004. Human MicroRNA targets. PLoS Biol. 2, e363. Kartha, R.V., Subramanian, S., 2014. Competing endogenous RNAs (ceRNAs): new entrants to the intricacies of gene regulation. Front. Genet. 5, 8. Langmead, B., Salzberg, S.L., 2012. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. Li, X., Liu, C.X., Xue, W., Zhang, Y., Jiang, S., Yin, Q.F., Wei, J., Yao, R.W., Yang, L., Chen, L.L., 2017. Coordinated circRNA biogenesis and function with NF90/NF110 in viral infection. Mol. Cell 67, 214–227 e217. Lin, X., Chen, Y., 2018. Identification of potentially functional CircRNA-miRNA-mRNA regulatory network in hepatocellular carcinoma by integrated microarray analysis. Med. Sci. Monit. Basic Res. 24, 70–78. Liu, B., Gong, Y.C., Li, Z., Hu, X.L., Cao, G.L., Xue, R.Y., Gong, C.L., 2016. Baculovirusmediated GCRV vp7 and vp6 genes expression in silkworm and grass carp. Mol. Biol. Rep. 43, 509–515. Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2(−Delta Delta C(T)) method. Methods 25, 402–408. Ma, J., Fan, Y.D., Zhou, Y., Liu, W.Z., Jiang, N., Zhang, J.M., Zeng, L.B., 2018. Efficient resistance to grass carp reovirus infection in JAM-A knockout cells using CRISPR/ Cas9. Fish & shellfish immunology 76, 206–215. Memczak, S., Jens, M., Elefsinioti, A., Torti, F., Krueger, J., Rybak, A., Maier, L., Mackowiak, S.D., Gregersen, L.H., Munschauer, M., Loewer, A., Ziebold, U., Landthaler, M., Kocks, C., le Noble, F., Rajewsky, N., 2013. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 495, 333–338. Mitra, A., Pfeifer, K., Park, K.S., 2018. Circular RNAs and competing endogenous RNA (ceRNA) networks. Transl. Cancer Res. 7, S624–S628. Morisaki, H., Morisaki, T., 2008. AMPD genes and urate metabolism. Nihon rinsho Japanese journal of clinical medicine 66, 771–777. Mu, Y.N., Ding, F., Cui, P., Ao, J.Q., Hu, S.N., Chen, X.H., 2010. Transcriptome and

5. Conclusions 3493 circRNAs were identified in CIK cells after GCRV infection, 76 of which were differentially expressed. Many parental genes of the DE circRNAs identified in CIK cells infected with GCRV are involved in ATP binding, metal ion binding, DNA binding, the regulation of transcription, and the apoptotic process. A total of 20,200 mRNAs were identified in the GCRV-infected CIK cells and uninfected CIK cells, 798 of which were differentially expressed, and these genes are involved in the immune responses, blood coagulation, hemostasis, and the complement and coagulation cascades. In total, 1685 miRNAs were detected in our study, and 186 DE miRNAs were identified in the CIK cells infected with GCRV. The targets of these miRNAs were associated with immunity. The complex regulatory relationships identified here indicate that circRNAs act in ceRNAs to regulate mRNAs by sponging their shared miRNAs during the process of infection. These results provided novel insights into the mechanisms underlying GCRV infection. Declaration of Competing Interest All authors declare that they have no conflicts of interest. Acknowledgments This work was supported by the Key Research and Development Program of Jiangsu Province (Modern Agriculture) (BE2016322), the Triple-New Project of Aquaculture of Jiangsu Province of China (D2017-3), and a project funded by the Priority Academic Program for the Development of Jiangsu Higher Education. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.aquaculture.2019.734349. References Chen, J., Li, C., Huang, R., Du, F., Liao, L., Zhu, Z., Wang, Y., 2012. Transcriptome analysis of head kidney in grass carp and discovery of immune-related genes. BMC Vet. Res. 8, 108. Chen, X.H., Wang, Q., Yang, C.R., Rao, Y.L., Li, Q.M., Wan, Q.Y., Peng, L.M., Wu, S.Q., Su, J.G., 2013. Identification, expression profiling of a grass carp TLR8 and its inhibition leading to the resistance to reovirus in CIK cells. Dev. Comp. Immunol. 41, 82–93. Chen, D.D., Yao, Y.Y., Cui, Z.W., Zhang, X.Y., Peng, K.S., Guo, X., Wang, B., Zhou, Y.Y., Li, S., Wu, N., Zhang, Y.A., 2018a. Comparative study of the immunoprotective effect of two DNA vaccines against grass carp reovirus. Fish & shellfish immunology 75, 66–73. Chen, G., He, L.B., Luo, L.F., Huang, R., Liao, L.J., Li, Y.M., Zhu, Z.Y., Wang, Y.P., 2018b. Transcriptomics sequencing provides insights into understanding the mechanism of grass carp reovirus infection. Int. J. Mol. Sci. 19. Cheng, J., Zhuo, H., Xu, M., Wang, L., Xu, H., Peng, J., Hou, J., Lin, L., Cai, J., 2018. Regulatory network of circRNA-miRNA-mRNA contributes to the histological classification and disease progression in gastric cancer. J. Transl. Med. 16, 216. Cui, X., Niu, W., Kong, L., He, M., Jiang, K., Chen, S., Zhong, A., Li, W., Lu, J., Zhang, L., 2016. hsa:circRNA_103636: potential novel diagnostic and therapeutic biomarker in

17

Aquaculture 258 (2019) 734349

B. Liu, et al. expression profiling analysis revealed changes of multiple signaling pathways involved in immunity in the large yellow croaker during Aeromonas hydrophila infection. BMC Genomics 11. Ouyang, Q., Huang, Q., Jiang, Z., Zhao, J., Shi, G.P., Yang, M., 2018. Using plasma circRNA_002453 as a novel biomarker in the diagnosis of lupus nephritis. Mol. Immunol. 101, 531–538. Pei, C., Ke, F., Chen, Z.Y., Zhang, Q.Y., 2014. Complete genome sequence and comparative analysis of grass carp reovirus strain 109 (GCReV-109) with other grass carp reovirus strains reveals no significant correlation with regional distribution. Arch. Virol. 159, 2435–2440. Qu, S., Yang, X., Li, X., Wang, J., Gao, Y., Shang, R., Sun, W., Dou, K., Li, H., 2015. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 365, 141–148. Rao, Y.L., Su, J.G., 2015. Insights into the antiviral immunity against grass carp (Ctenopharyngodon idella) Reovirus (GCRV) in grass carp. J Immunol Res, 670437. Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., Ramage, D., Amin, N., Schwikowski, B., Ideker, T., 2003. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. Shi, M., Huang, R., Du, F., Pei, Y., Liao, L., Zhu, Z., Wang, Y., 2014. RNA-seq profiles from grass carp tissues after reovirus (GCRV) infection based on singular and modular enrichment analyses. Mol. Immunol. 61, 44–53. Su, J., Dong, J., Huang, T., Zhang, R., Yang, C., Heng, J., 2011. Myeloid differentiation factor 88 gene is involved in antiviral immunity in grass carp Ctenopharyngodon idella. J. Fish Biol. 78, 973–979. Su, J., Heng, J., Huang, T., Peng, L., Yang, C., Li, Q., 2012. Identification, mRNA expression and genomic structure of TLR22 and its association with GCRV susceptibility/resistance in grass carp (Ctenopharyngodon idella). Dev. Comp. Immunol. 36, 450–462. Sun, L., Liu, P., Sun, S., Yan, S., Cao, C., 2018. Transcriptomic analysis of interactions between Hyphantria cunea larvae and nucleopolyhedrovirus. Pest management science 75, 1024–1033. Wan, Q.Y., Su, J.G., 2015. Transcriptome analysis provides insights into the regulatory function of alternative splicing in antiviral immunity in grass carp (Ctenopharyngodon idella). Sci Rep-Uk 5. Wang, Q., Zeng, W.W., Liu, C., Zhang, C., Wang, Y.Y., Shi, C.B., Wu, S.Q., 2012. Complete genome sequence of a reovirus isolated from grass carp, indicating different genotypes of GCRV in China. J. Virol. 86, 12466. Wang, H., Shen, X., Xu, D., Lu, L., 2013a. Lipopolysaccharide-induced TNF-alpha factor in grass carp (Ctenopharyngodon idella): evidence for its involvement in antiviral innate immunity. Fish & shellfish immunology 34, 538–545. Wang, L., Shang, N., Feng, H., Guo, Q., Dai, H., 2013b. Molecular cloning of grass carp

(Ctenopharyngodon idellus) T-bet and GATA-3, and their expression profiles with IFN-gamma in response to grass carp reovirus (GCRV) infection. Fish Physiol. Biochem. 39, 793–805. Wang, Y., Liu, G.L., Li, D.L., Ling, F., Zhu, B., Wang, G.X., 2015. The protective immunity against grass carp reovirus in grass carp induced by a DNA vaccination using singlewalled carbon nanotubes as delivery vehicles. Fish & shellfish immunology 47, 732–742. Wang, J.Y., Yang, Y.W., Jin, L.M., Ling, X.T., Liu, T.L., Chen, T.Z., Ji, Y.H., Yu, W.G., Zhang, B.L., 2018. Re-analysis of long non-coding RNAs and prediction of circRNAs reveal their novel roles in susceptible tomato following TYLCV infection. BMC Plant Biol. 18. Xue, R., Liu, L., Cao, G., Xu, S., Li, J., Zou, Y., Chen, H., Gong, C., 2013. Oral vaccination of BacFish-vp6 against grass carp reovirus evoking antibody response in grass carp. Fish & shellfish immunology 34, 348–355. Yang, C., Su, J., Li, Q., Zhang, R., Rao, Y., 2012. Identification and expression profiles of ADAR1 gene, responsible for RNA editing, in responses to dsRNA and GCRV challenge in grass carp (Ctenopharyngodon idella). Fish & shellfish immunology 33, 1042–1049. Yu, F., Wang, L., Wang, H., Sheng, J., Lu, L., 2017. Repression of SUMOylation pathway by grass carp reovirus contributes to the upregulation of PKR in an IFN-independent manner. Oncotarget 8, 71500–71511. Yuan, Y., Liu, W., Zhang, Y., Zhang, Y., Sun, S., 2018. CircRNA circ_0026344 as a prognostic biomarker suppresses colorectal cancer progression via microRNA-21 and microRNA-31. Biochem. Biophys. Res. Commun. 503, 870–875. Zhang, X., Yan, Y., Lei, X., Li, A., Zhang, H., Dai, Z., Li, X., Chen, W., Lin, W., Chen, F., Ma, J., Xie, Q., 2017a. Circular RNA alterations are involved in resistance to avian leukosis virus subgroup-J-induced tumor formation in chickens. Oncotarget 8, 34961–34970. Zhang, X., Zhu, M., Yang, R., Zhao, W., Hu, X., Gan, J., 2017b. Identification and comparison of novel circular RNAs with associated co-expression and competing endogenous RNA networks in pulmonary tuberculosis. Oncotarget 8, 113571–113582. Zhang, F., Zhang, R., Zhang, X., Wu, Y., Li, X., Zhang, S., Hou, W., Ding, Y., Tian, J., Sun, L., Kong, X., 2018a. Comprehensive analysis of circRNA expression pattern and circRNA-miRNA-mRNA network in the pathogenesis of atherosclerosis in rabbits. Aging 10, 2266–2283. Zhang, X., Zhu, M., Hu, X., 2018b. Integrated miRNA and mRNA expression profiling to identify mRNA targets of dysregulated miRNAs in pulmonary tuberculosis. Epigenomics 10, 1051–1069. Zhao, L., Alto, B.W., Shin, D., Yu, F., 2018. The effect of permethrin resistance on Aedes aegypti transcriptome following ingestion of Zika virus infected blood. Viruses 10.

18