circRNA expression patterns and circRNA-miRNA-mRNA networks during CV-A16 infection of SH-SY5Y cells.

Coxsackievirus A16 (CV-A16) has caused worldwide epidemics of hand, foot, and mouth disease (HFMD) in infants and preschool children. Circular RNAs (circRNAs), a class of noncoding RNA molecules, participate in the progression of viral infectious diseases. Although the function of circRNAs has been a heavily researched topic, their role in CV-A16 infection is still unclear. In this study, the viral effects of CV-A16 on the cellular circRNA transcriptome were investigated using next-generation sequencing technology. The results showed that a total of 8726, 8611, and 6826 circRNAs were identified at 0, 12, and 24 h postinfection, respectively. Moreover, it was found that 1769 and 1192 circRNAs were differentially expressed in at 12 and 24 h postinfection, respectively. The common differentially expressed circRNAs were used for functional annotation analysis, and it was found that the parent genes of differentially expressed circRNAs might be associated with the viral infection process, especially the "Immune system process" in GO analysis and the "Inflammation mediated by chemokine and cytokine signaling pathway" in KEGG analysis. Subsequently, circRNA-miRNA-mRNA regulatory networks were constructed, and the hsa_circ_0004447/hsa-miR-942-5p/MMP2, hsa_circ_0078617/hsa-miR-6780b-5p/MMP2 and hsa_circ_0078617/hsa-miR-5196-5p/MMP2 regulatory axes were identified by enrichment analysis as important networks during the progression of CV-A16 infection. Finally, six dysregulated circRNAs were selected for validation and were verified to be consistent with the sequencing results. Considering all of these results, to the best of our knowledge, this study is the first to present a comprehensive overview of circRNAs induced by CV-A16 infection, and this research demonstrated that a network of enriched circRNAs and circRNA-associated competitive endogenous RNAs (ceRNAs) is involved in the regulation of CV-A16 infection, thereby helping to elucidate the mechanisms underlying CV-A16-host interactions.

Introduction

Hand, foot, and mouth disease (HFMD) is a common infectious illness in infants and young children under 5 years of age that includes such symptoms as 3–4 days of fever, reduced appetite, and the appearance of maculopapular or vesicular rash on the mouth, hands, feet, and/or buttocks [7, 28]. The two principal pathogenic agents of HFMD are enterovirus 71 (EV-A71) and coxsackievirus A16 (CV-A16) [2, 11]. In general, EV-A71 is more frequently associated with severe cases and fatalities accompanied by serious complications, especially those involving the central nervous system (CNS) [4], while CV-A16 infection usually results in a mild outcome [11, 22]. Hence, most studies have focused on EV-A71, and prophylactic EV-A71 vaccines developed by three companies (Beijing Vigoo Biological, Sinovac Biotech Co. Ltd, and the Institute of Medical Biology) have all shown high efficacy and satisfactory safety in providing protection against EV-A71-associated diseases [2]. Unfortunately, these vaccines do not have cross-strain protective activity against HFMD caused by other enteroviruses, including CV-A16 [2, 4]. However, in the past several decades, the number of clinical cases caused by CV-A16 has been observed to increase steadily in the Asia–Pacific region [25]. Moreover, CV-A16 infections have also been reported to have more severe CNS complications, such as aseptic meningitis, cerebellar ataxia, poliomyelitis-like paralysis, acute brainstem encephalitis, cardiopulmonary failure, and fulminant neurogenic pulmonary edema [25, 36]. Importantly, no antiviral therapies or vaccines for CV-A16 are currently available; therefore, several vaccine companies and academic institutions are launching projects to develop monovalent or bivalent CV-A16 vaccines.

Circular RNAs (circRNAs) are a novel class of noncoding RNAs characterized by the formation of covalently closed continuous loops without 5′–3′ polarities or poly(A) tails [21, 29]. Numerous studies have clearly indicated that circRNAs are widely expressed in plants, animals and humans with tissue- and developmental-stage-specific patterns [5, 19]. More importantly, circRNAs are involved in the regulatory process of gene expression by serving as microRNA (miRNA) sponges, regulators of splicing and transcription, platforms for the assembly of protein complexes, and templates for translation [6]. Therefore, circRNAs play an important role in many cellular processes, such as embryonic development, control of the cell cycle, cellular senescence, cell signaling, and response to cellular stress [15]. The perturbation of circRNA expression might be correlated with various diseases, such as cancer, heart disease, neurological disorders, diabetes, and atherosclerosis, and even viral diseases [13, 34]. Hence, in this study, deep circRNA sequencing technology was used to explore the potential functions of aberrantly expressed circRNAs in CV-A16-infected SH-SY5Y cells, and these results might help to elucidate the mechanisms governing the interactions between CV-A16 and the host and to identify new therapeutic targets against CV-A16 infection.

Materials and methods

Cell culture and virus inoculation

The human neuroblastoma cell line SH-SY5Y, obtained from Jennino Biological Technology, was cultivated in Roswell Park Memorial Institute-1640 medium (RPMI-1640; Corning, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, USA), streptomycin (100 mg/ml), and penicillin (100 U/ml) in a humidified atmosphere with 5% CO2 and 90% humidity at 37 °C. Cell passage was performed once the cells had reached 90% confluence by rinsing twice with phosphate-buffered saline (PBS) and digesting with 1 ml of 0.25% trypsin (Beyotime Biotechnology Co., China). Next, the single-cell suspension was placed into culture flasks or 6-well plates for further experiments.

The virus CV-A16-G20 strain (subgenotype B, GenBank no. JN590244.1), isolated from an HFMD patient in Guangxi, China, in 2010, was used in this study. Six-well tissue culture plates were seeded at ∼ 1 × 106 SH-SY5Y cells/well. When the cells reached approximately 90% confluence, they were washed twice with PBS, and CV-A16 was incubated with the cells at a multiplicity of infection (MOI) of 1. Following infection for 0 h, 12 h, or 24 h, the cells were collected, and the cells treated with CV-A16 for 0 h were defined as the control group.

Analysis of viral replication and viral protein synthesis by quantitative real-time polymerase chain reaction (qRT-PCR) and western blot (WB)

Total viral RNA was extracted from the samples using an RNA extraction kit (TIANGEN, China), following the manufacturer’s instructions. The regions of standards generated from CV-A16-VP1 RNA were amplified using PrimeSTAR HS reagent (TAKARA, Japan) with the primers shown in Supplementary Table S1. The target 682-bp region (CV-A16-VP1) was cloned into the pGM-T vector (TIANGEN, China) following the manufacturer’s guidelines. A single recombinant clone containing the insert under the control of a T7 promoter was selected, and the insert was transcribed in vitro using a TranscriptAid T7 High Yield Transcription kit (Thermo Scientific, USA) according to the manufacturer’s recommendations. The in vitro-transcribed RNA was purified using a GeneJET RNA Purification Kit (Thermo Scientific, USA) and quantified using a NanoDrop 2000 Spectrophotometer (Thermo Scientific, USA). qRT-PCR was performed using a One Step PrimeScriptTM RT-PCR Kit (TAKARA, Japan). To generate a standard curve for cycle threshold (Ct) vs. virus copy number, the purified RNA was serially diluted to known concentrations in a range of 101–108 copies/µl and assayed concurrently with the test samples. The PCR conditions were as follows: 5 min at 42 °C for reverse transcription and 10 s of denaturation at 95 °C, followed by 40 cycles of 95 °C for 5 s and 60 °C for 34 s. The reactions were run on a 7500 Fast Real-Time PCR System (Applied Biosystems, USA). The qRT-PCR primers used for determination of CV-A16 copy numbers are listed in Table 1. Each sample was tested in triplicate.

Table 1

Primers used for qRT-PCR analysis of circRNA and reference genes

CircBase ID or gene name	Primer sequences
hsa_circ_0000205	Forward primer: 5’-AGTTGGCTCTCACTGCTTCT-3’ Reverse primer: 5’-GGTCACTCCTGCAATAAGACT-3’
hsa_circ_0002485	Forward primer: 5’-CGTTTTCAGCGTGACAAGGA-3’ Reverse primer: 5’-CCAGCTGGTACCTCTTACTGT-3’
hsa_circ_0007059	Forward primer: 5’-GAGACAGTAGCCATCCAGC-3’ Reverse primer: 5’-TGATCTGAGTCCAGGTGTT-3’
hsa_circ_0002483	Forward primer: 5’-TGCCAAAAGGATTTCTAAACCAGT-3’ Reverse primer: 5’-TTGGGGTCAAGGTAAGCAGC-3’
hsa_circ_0002301	Forward primer: 5’-TATATGGTCAACTGCAACTTGGC-3’ Reverse primer: 5’-TCACATGTCTCCACCCTTGT-3’
hsa_circ_0002141	Forward primer: 5’-TCATTGTCAAGAGAGCAGATACT-3’ Reverse primer: 5’-TCCTTGTCTTTTCTGCATCTTGA-3’
GAPDH	Forward primer: 5’-ATGAGGTCCACCACCCTGTT-3’ Reverse primer: 5’-CTCAAGGGCATCCTGGGCTA-3’
CV-A16-VP1-1 (used for the construction of standard RNA)	5’-AACACTGAGGCTAGTAGTCAC-3’ (sense) 5’-CGTGTTTGATTCTCATGTACACC-3’ (anti-sense)
CV-A16-VP1-2 (used for qRT-PCR examination)	5’-GTTTGTGAAAATGACGGACCC-3’ (sense) 5’-GTCATTTGCTTGAAGGTGCTC-3’ (anti-sense) Probe: FAM-CAGCTCAAGTGTCAGTCCCCT-TAMRA

WB analysis was performed according to a standard method. Briefly, total proteins from cell samples (including control, CV-A16-12 h, and CV-A16-24 h) were extracted using radioimmunoprecipitation assay (RIPA) buffer containing protease inhibitors, and the protein concentration was determined using a BCA protein assay kit (Beyotime, China). Thirty μg of protein lysates from each time point was resolved by 8–10% sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) at 80 V for 2 h and transferred to polyvinylidene difluoride (PVDF) membranes (Millipore, USA) at a constant current of 200 mA for 60 min. After blocking with 5% skimmed milk at room temperature for 1 h, the PVDF membranes were incubated with primary antibodies against VP1 (1:1000 dilution; Millipore, USA) and GAPDH (1:500 dilution; Boster, China) overnight at 4 ℃. After washing three times with 1× TBST, the membranes were incubated with the corresponding secondary antibodies (1:10,000 dilution; Boster, China) for 1 hour at 37 ℃. The membranes were again washed three times with TBST and finally visualized using an ECL Western Blot Detection Kit (Amersham, USA) and exposed to X-ray films (Kodak, Japan). Relative expression levels of each protein were normalized to the endogenous control, GAPDH, using ImageJ 1.8.0 software.

RNA extraction and circRNA sequencing

Total RNA was extracted at each time point from three independent experiments and pooled into one group, following the instructions of TRIzol Reagent (Invitrogen, USA). RNA integrity was analyzed using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA), while the RNA concentration was determined using a Qubit RNA Assay Kit in a Qubit fluorometer (Invitrogen, USA). Total RNA samples that met the following requirements were used in the subsequent experiments: total RNA level > 5 μg, RNA integrity number (RIN) ≥ 7.0, and a 28S:18S ratio ≥ 1.5. Finally, a Ribo-Zero™ Magnetic Kit (Epicentre Technologies, USA) was used to remove ribosomal RNA (rRNA), and RNase R (Epicentre Technologies, USA) was used to remove linear RNA and enrich circular RNA.

Sequencing libraries were generated using an NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer’s protocol. The libraries were controlled for quality and quantified using a BioAnalyzer 2100 system. Fragments were sequenced on an Illumina HiSeq2000 platform using 150-bp paired-end reads. We removed adapters from the raw FASTQ files using Trim Galore. After 3’ adapters were trimmed and low-quality reads were removed using Cutadapt software, the high-quality trimmed reads were used for analysis of circRNAs. We aligned the filtered sequence data to a human reference genome sequence (GRCh37/hg19) obtained from the UCSC genome database (http://genome.ucsc.edu/), using TopHat2. The sequencing data that could not be aligned directly to the reference genome sequence were subjected to subsequent circRNA analysis by recognition of the reverse splicing events using Find_circ (v1.0) and CIRCexplorer2 software. We normalized the circRNA content as the number of uniquely mapped fragments per kilobase of exon per million fragments mapped (FPKM).

Bioinformatic analysis

Detection and characterization of overall circRNA data

Raw junction reads for all samples were normalized by total read number and were log2 transformed. Next, the distribution of circRNAs on human chromosomes was explored, and the positive and negative strands of circRNAs located on different chromosomes were also analyzed. The log2 ratio was used for quantile normalization. circRNA sequencing data were deposited in the NCBI Gene Expression Omnibus (GEO) under accession number GSE151890.

Identification of differently expressed circRNAs

The Limma (v3.32.10) package of R software was used to analyze the differential expression of circRNAs. circRNAs with fold changes > 2 and P-values < 0.05 were selected as significantly differential expression. According to the source of circRNA formation, the classification of circRNAs was divided into four main categories, namely, exonic, intronic, antisense, and intragenic, which are summarized in a pie chart.

Screening of common circRNAs that are differentially expressed during CV-A16 infection

To identify the majority of CV-A16-dependent circRNAs that were induced in both the CV-A16-12 h and CV-A16-24 h groups, Venn diagram analysis was performed. Hierarchical clustering was performed to obtain an overview of the expression profiles of circRNAs among the samples using these selected common circRNAs.

Functional analysis of the differentially expressed circRNAs and their host genes

To investigate the potential functions of the parent genes of the differentially expressed circRNA, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway based on the Database for Annotation, Visualization and Integrated Discovery (DAVID) were used in this study. GO terms provided information about the biological process (BP), cellular component (CC), and molecular function (MF), whereas KEGG analysis was performed to explore the pathways related to circRNA-hosting genes.

Construction of a circRNA-miRNA-mRNA coexpression network

To establish a circRNA-miRNA-mRNA network, we first searched for miRNA response elements to circRNAs using miRanda and RNAhybrid software. A hit between circRNAs and a target miRNA with a perfect seed match was considered when the miRanda score was 150 or higher and the RNAhybrid energy was less than -25. The putative miRNAs were then further applied to the prediction of their target genes based on the TargetScan and miRDB algorithms. In order to classify the target genes corresponding to the miRNAs according to their function, we screened the miRNA target genes according to the functions regulated by the circRNAs, but the screening of the miRNA target gene had to provide a perfect seed match between the miRNA and the target gene in TargetScan, with a target score ≥ 80 in miRDB.

Confirmation of differential expression of circRNAs by qRT-PCR

We used qRT-PCR to validate our RNA sequencing results, focusing on circRNAs that differed significantly in their expression levels. Total RNA from cells and tissue samples was isolated using TRIzol Reagent (Invitrogen, USA). Specific divergent primers spanning the back-splice junction site of circRNAs were designed by Sangon Biotech (Shanghai, China) and are listed in Table 1. GAPDH, a housekeeping gene, was employed as an internal control. To quantify the amount of circRNA, complementary DNA (cDNA) was synthesized using a PrimeScript™ RT Reagent Kit (TaKaRa, Japan). qRT-PCR analysis was performed using TB Green™ Premix Ex Taq™ II (TaKaRa, Japan) in an ABI 7500 real-time PCR system (Applied Biosystems, USA). The formula 2−ΔΔCt was used to analyze the relative expression levels of different circRNAs.

Statistical analysis

All of the qRT-PCR experiments were conducted in triplicate, and all data are reported as the mean ± SD. Differences between two groups were analyzed by Student’s t-test by using SPSS 22.0 (IBM-SPSS, USA), and a P-value < 0.05 indicated a statistically significant difference.

Results

Growth curves of CV-A16 in SH-SY5Y cells

To study CV-A16 replication kinetics in SH-SY5Y cells, viral RNA and protein synthesis was monitored by qRT-PCR and WB, respectively, at each time point. The growth curves showed that intracellular viral RNA gradually increased over time (Supplementary Fig. S1A). It was also observed that the intracellular viral protein VP1 was detected at 12 hpi and increased at 24 hpi (Fig. S1B). These data show that CV-A16 could actively replicate in SH-SY5Y cells.

Overview of circRNA expression

We identified circRNAs in CV-A16-infected SH-SY5Y cells and uninfected SH-SY5Y cells using RNA-Seq data generated from rRNA-depleted samples. We identified 8726, 8611, and 6826 circRNAs in the control, CV-A16-12 h, and CV-A16-24 h group, respectively. Moreover, the genomic loci from which the circRNAs were derived were observed to be widely distributed across the chromosomes, except the Y chromosome, and the positive and negative strands of circRNAs located on different chromosomes are indicated in Fig. 1. These differences at different time points indicate that CV-A16 infection triggered changes in circRNA expression in SH-SY5Y cells.

Fig. 1

View of circRNAs expression profiles showing the chromosomal distribution of the circRNAs in the control group and CV-A16-infected groups

Analysis of differentially expressed circRNAs

Differential expression analysis showed that compared with the controls, 1769 differentially expressed circRNAs were found in the CV-A16-12 h group, 821 of which were upregulated and 948 of which were downregulated, while 1192 differentially expressed circRNAs were found in the CV-A16-24 h group, 653 of which were upregulated and 539 of which were downregulated (Fig. 2A). Mapping of the circRNAs to the source showed that most of the significantly differentially expressed circRNAs were transcribed from the exons of protein coding regions (Fig. 2B). To investigate the roles of these distinct circRNAs in the progression of CV-A16 infection, circRNAs that were differentially expressed in both the CV-A16-12 h group and the CV-A16-24 h group were identified. A Venn diagram revealed that the cross-section represented 554 common circRNAs in each sample, and the independent section represents the number of sample-specific circRNAs during CV-A16 infection (Fig. 2C). Furthermore, these common circRNAs were utilized to perform hierarchical clustering analysis, which illustrated the distinguishable circRNA expression pattern among the samples (Fig. 2D). These data suggest that changes in numerous circRNAs are induced by CV-A16 infection, and the common differentially expressed circRNAs at different time points during CV-A16 infection might be regulatory factors in the pathogenesis of CV-A16.

Fig. 2

Differences and characteristics of circRNA expression profiles between the CV-A16-infected groups and the control group. A Number of dysregulated differentially expressed circRNAs in the CV-A16-infected groups. B The percentage of significantly differentially expressed circRNAs arising from different genomic loci, presented as a pie chart. C Venn diagram of differentially expressed circRNAs in the CV-A16-infected groups showing similarities and differences. D Hierarchical clustering results of circRNA expression profiles of the CV-A16-infected groups and the control group. “Red” indicates high relative expression, and “green” indicates low relative expression

Biomathematically predicted biological function of host linear transcripts

To investigate the potential functional implications of the circRNAs, we performed functional enrichment analysis for the parent genes of these common differentially expressed circRNAs. The GO terms and KEGG pathway analysis were presented and ranked by P-value. The results showed that for GO-BP terms, cellular process, metabolic process, biological regulation, localization, response to stimulus, signaling, developmental process, multicellular organismal process, biological adhesion, locomotion, immune system process, reproductive process, interspecies interaction between organisms, reproduction, biological phase, multi-organism process, and biomineralization were enriched (Fig. 3A). For GO-MF terms, the majority of the enriched terms were related to binding, catalytic activity, molecular function regulator, transporter activity, molecular transducer activity, molecular adaptor activity, structural molecule activity, and translation regulator activity (Fig. 3B). For GO-CC enrichment terms, host linear transcripts of differentially expressed circRNAs involved in cellular anatomical entity, intracellular, and protein-containing complex were enriched (Fig. 3C). Additionally, KEGG pathway analysis demonstrated that the top 10 regulated pathways were the gonadotropin-releasing hormone receptor pathway, EGF receptor signaling pathway, Wnt signaling pathway, FGF signaling pathway, angiogenesis, integrin signaling pathway, inflammation mediated by chemokine and cytokine signaling pathway, PDGF signaling pathway, CCKR signaling map, and p53 pathway (Fig. 4). Hence, the results implied that vital biological functions regulated by the host linear transcripts of differentially expressed circRNAs might be closely associated with the pathogenic mechanism of CV-A16 infection.

Fig. 3

Enrichment map of GO annotation for the host genes of the common differentially expressed circRNAs. A Biological process. B Molecular function. C Cellular component. The number of genes and GO terms are shown on the x- and y- axes, respectively

Fig. 4

KEGG enrichment analysis performed on the host genes of common differentially expressed circRNAs in the CV-A16-infected groups versus the control groups. The top 10 significantly enriched pathways are displayed. The vertical lines correspond to the names of the pathways, and the horizontal line represents the number of genes

Prediction of circRNA-miRNA-target gene associations

A large number of reports have shown that circular RNAs affect miRNA-mediated regulation of gene expression through miRNA sequestration [20]. In order to identify circRNAs involved in the pathogenic mechanism of CV-A16 infection, we selected circRNAs related to “Immune system process” in GO-BPs analysis and circRNAs related to “Inflammation mediated by chemokine and cytokine signaling pathway” in KEGG analysis. Detailed information on these circRNAs is listed in Tables 2 and 3, respectively. Potential circRNA-binding miRNAs were predicted based on the miRanda and RNAhybrid algorithm. Moreover, to further elucidate circRNA functions, the target genes of the putative miRNAs were predicted using TargetScan and miRDB. The final result of the prediction was used to draw a circRNA-miRNA-mRNA regulatory network map, which showed that that five circRNAs, 69 miRNAs, and 17 mRNAs were in the “Immune system process”-related network (Fig. 5A), while nine circRNAs, 134 miRNAs and 23 mRNAs were in the “Inflammation mediated by chemokine and cytokine signaling pathway”-related network (Fig. 5B). Taken together, the regulatory relationships among circRNAs, miRNAs and mRNAs might suggest that dysregulated circRNAs play a role in CV-A16 infection.

Table 2

Differentially expressed circRNAs associated with “Immune system process” in GO-BPs analysis

Acc ID	CV-A16-12 h		CV-A16-24 h		CHROM	Gene name	Circ type	Seq length
	Log2FC	Style	Log2FC	Style
chr13_46728984_46723221_+5763-LRCH1	20	Up	1	Up	chr13	LRCH1	Exonic	248
chr4_73125312_73118688_-6624-ANKRD17	2	Up	− 2.584962501	Down	chr4	ANKRD17	Exonic	954
chr2_195681050_195680032_+1018-SLC39A10	20	Up	− 20	Down	chr2	SLC39A10	Exonic	1019
chr6_47286595_47283938_-2657-TNFRSF21	1.378511623	Up	− 20	Down	chr6	TNFRSF21	Exonic	1147
chr5_108798389_108768093_+30296-FER	20	Up	− 20	Down	chr5	FER	Exonic	412

Table 3

Differentially expressed circRNAs associated with “Inflammation mediated by chemokine and cytokine signaling pathway” in KEGG analysis

Acc ID	CV-A16-12 h		CV-A16-24 h		CHROM	Gene name	Circ type	Seq length
	Log2FC	Style	Log2FC	Style
chr10_94032252_94030683_+1569-PLCE1	20	Up	− 20	Down	chr10	PLCE1	Exonic	1570
chr16_69677335_69626403_+50932-NFAT5	20	Up	− 20	Down	chr16	NFAT5	Exonic	1563
chr16_68174514_68121987_+52527-NFATC3	1.736965594	Up	− 3.321928095	Down	chr16	NFATC3	Exonic	1812
chr21_33432871_33414888_+17983-IFNGR2	20	Up	− 20	Down	chr21	IFNGR2	Exonic	806
chr1_243843282_243637611_-205671-AKT3	− 20	Down	20	Up	chr1	AKT3	Exonic	673
chr21_33349543_33345246_+4297-IFNAR1	20	Up	− 20	Down	chr21	IFNAR1	Exonic	470
chr2_46010517_46001404_+9113-PRKCE	20	Up	− 20	Down	chr2	PRKCE	Exonic	614
chr2_39058804_39051144_-7660-SOS1	20	Up	− 20	Down	chr2	SOS1	Exonic	651
chr6_161049979_161034259_+15720-MAP3K4	1.94753258	Up	− 3.169925001	Down	chr6	MAP3K4	Exonic	1555

Fig. 5

circRNA-miRNA-mRNA interaction network diagram. The network consists of circRNAs (red round rectangles), miRNAs (yellow ellipses), and their target genes (blue triangles). Based on the prediction and bioinformatics analysis, the circRNA-miRNA-mRNA triple network was found to be associated with “Immune system process” in GO-BPs (A) and “Inflammation mediated by chemokine and cytokine signaling pathway” in KEGG analysis (B). Moreover, the overlapping target genes were further screened and the key ceRNA interactions were found (C)

Innate antiviral and inflammatory responses can immediately be elicited upon viral infection [1, 14]. These responses must be finely regulated to prevent viral dissemination and present the appropriate effective immune responses [18]. Therefore, we searched for targeted genes that were present in both network diagrams, and 16 genes were identified (Fig. 5C). One of these 16 genes was MMP2. Therefore, the network axis of MMP2 regulated by circRNA was screened, including hsa_circ_0004447/hsa-miR-942-5p/MMP2, hsa_circ_0078617/hsa-miR-6780b-5p/MMP2, and hsa_circ_0078617/hsa-miR-5196-5p/MMP2 (Fig. 5C), which simultaneously participated in “Immune system process” and “Inflammation mediated by chemokine and cytokine signaling pathway” during CV-A16 infection. Thus, the three key circRNA-associated ceRNA networks pointed out a new direction for our future research.

Validation of the differential expression levels of circRNAs

In order to verify the reliability of the circRNA-seq results, six circRNAs that showed differential expression, hsa_circ_0000205, hsa_circ_0002485, hsa_circ_0007059, hsa_circ_0002483, hsa_circ_0002301, and hsa_circ_0002141, were selected for validation by qRT-PCR. The results verified that the data obtained by qRT-PCR were highly consistent with the high-throughput data for all six of these circRNAs (Fig. 6).

Fig. 6

Validation of selected circRNAs. Six circRNAs were significantly amplified by qRT-PCR, and the results were basically consistent with the microarray results

Discussion

circRNAs were originally thought to be by-products of aberrant splicing with little functional potential [10]. Novel bioinformatic approaches coupled with biochemical enrichment strategies and deep sequencing have enabled comprehensive studies of circRNAs to be performed [9], and several studies have addressed the identification, characterization, and function of circRNAs during viral infection [33]. For example, human cytomegalovirus (HCMV) infection causes significant changes in host circRNA expression, and it was found that circSP100 potentially bound to 10 HCMV-encoding proteins, which further influenced the process of HCMV infection [8]. The host circRNA-miRNA-mRNA network was found to be significantly changed in human lung adenocarcinoma epithelial (Calu-3) cells infected with the highly pathogenic Middle East respiratory syndrome coronavirus (MERS-CoV), and circFNDC3B and circCNOT were also shown to reduce the MERS-CoV load by modulating various biological pathways, including the mitogen-activated protein kinase (MAPK) and ubiquitination pathways [37]. Nevertheless, the properties and potential roles of circRNAs during CV-A16 infection have not been thoroughly elucidated to date. In this study, we applied a next-generation sequencing technique to systematically analyze circRNA expression profiles in SH-SY5Y cells infected with CV-A16. Here, we chose SH-SY5Y cells, mainly because they are often used as susceptible cells to investigate the neuropathic mechanism of neurotropic viruses in vitro. For example, SH-SY5Y cells have been used as an in vitro adult human neuronal cell-based model, for studying the biology of Zika virus (ZIKV) [24]. Likewise, CV-A16 is also a neurotropic virus that has been reported to induce neurological complications.

In this study, thousands of circRNAs from diverse genomic locations at different times postinfection were identified. Among these circRNAs, it was found that 1769 and 1192 were differentially expressed at 12 and 24 h postinfection, respectively. While a small fraction of these circRNAs were intronic in origin, more than 90% were exonic. The large number of different circRNAs present at different times indicated that the transcription of circRNAs was very active after CV-A16 infection. Functional circRNAs can affect host gene expression directly by regulating transcription or interfering with splicing, mainly due to the fact that circRNAs are usually generated from host genes through backsplicing. Thus, the function of these dysregulated circRNAs was assessed by GO and KEGG enrichment analysis with their host genes. Furthermore, to identify circRNAs that are involved in the process of CV-A16 infection, we screened out only the circRNAs that were differentially expressed at both 12 and 24 h postinfection. The host genes that were related to these 554 overlapping dysregulated circRNAs participated in various biological processes, including cellular process, metabolic process, biological regulation, localization, response to stimulus, signaling, developmental process, multicellular organismal process, biological adhesion, locomotion, immune system process, reproductive process, interspecies interaction between organisms, reproduction, biological phase, multi-organism process, and biomineralization, most of which are associated with the initiation and progression of viral infections. Additionally, the KEGG annotation also showed that these circRNAs might primarily regulate the gonadotropin-releasing hormone receptor pathway, EGF receptor signaling pathway, Wnt signaling pathway, FGF signaling pathway, angiogenesis, integrin signaling pathway, inflammation mediated by chemokine and cytokine signaling pathway, PDGF signaling pathway, CCKR signaling map, and p53 pathway, which were involved in host regulatory mechanisms during virus infection. For instance, EGF receptor inhibitors have been found to inhibit hepatitis B virus (HBV) and hepatitis C virus (HCV) replication via downregulation of signal transducers and activator of transcription 3 (STAT3) phosphorylation [12]. Activation of the Wnt/β-catenin pathway has been shown to increase influenza virus mRNA and virus production in vitro in mouse lung epithelial E10 cells and mRNA expression of influenza virus genes in vivo in mouse lungs [26]. Thence, the above results suggested that the alterations in circRNAs induced by CV-A16 infection in SH-SY5Y cells might play important regulatory roles by directly influencing host biological functions and pathways. We used qRT-PCR methods to test and verify the reliability of the RNA-seq data and found that the expression of randomly selected circRNAs was consistent with the sequencing results, although the fold change was slightly different from that deduced from the sequencing results.

As highly conserved endogenous RNAs, many circRNAs harbor abundant miRNA binding sites, indicating that they can sequester the corresponding miRNAs and thus function as competing endogenous RNAs (ceRNAs) to regulate gene expression [31]. For instance, circRNA Cdr1as was shown to function as a ceRNA to promote hepatocellular carcinoma progression by acting as a sponge for miR-1270 and enhancing the expression of AFP [30]. In another study, the circRNA_15698/miR-185/TGF-beta1 axis was found to aggravate the extracellular matrix of diabetic nephropathy mesangial cells, which promoted the pathogenesis of diabetic nephropathy [16]. Thus, in order to investigate the potential functions of key circRNAs, we focused on the major significantly changed GO terms (i.e., “Immune system process”) and most correlated pathways (i.e., “Inflammation mediated by chemokine and cytokine signaling pathway”), and a circRNA-associated-ceRNA network was constructed. Most viruses are controlled sufficiently by the immune system that limited damage is done to host tissues [32]. However, during the evolution of the viruses, they also can enhance their survival in the host cell through an immune escape mechanism [27]. Thus, exploring the “Immune system process” induced by CV-A16 infection might provide more clues for understanding its pathogenic mechanism. In addition, there is accumulating evidence that inflammation is a major driving force in virus infection [3, 17], and the development of neurological symptoms caused by CV-A16 is closely associated with inflammation of the nervous system. Hence, the investigation of “Inflammation mediated by chemokine and cytokine signaling pathway” could help us to discover potential mechanisms and therapeutic agents in CV-A16 infection progression. In the circRNA-miRNA-mRNA triple network map, it was seen that miRNAs that potentially bind to circRNAs, and the most likely target genes for each miRNA were identified. Among these complicated networks, we paid particularly close attention to the circRNAs that may act as ceRNAs to regulate the expression of MMP2, which has previously been reported to participate in the regulation of viral infections. For example, MMP2 can mediate viral clearance during HBV infection by cleaving membrane-bound CD100 into soluble CD100 from T cells [35]. Dengue virus can infect macrovascular endothelial cells, resulting in overproduction of MMP-2, which might contribute to the pathogenesis of severe dengue infection [23]. In the current study, it was discovered that there were three circRNAs-regulated axes on MMP2, namely hsa_circ_0004447/hsa-miR-942-5p/MMP2, hsa_circ_0078617/hsa-miR-6780b-5p/MMP2, and hsa_circ_0078617/hsa-miR-5196-5p/MMP2. This gave us a clear direction to study the specific mechanism of these circRNAs in the pathogenesis of CV-A16 infection.

In conclusion, we have identified a population of circRNAs that are differentially expressed in CV-A16-treated and CV-A16-untreated SH-SY5Y cells by high-throughput sequencing analysis and verified the expression of six dysregulated circRNAs by qRT-PCR. Further analysis suggested that these common differentially expressed circRNAs might play a role in viral pathogenicity by participating in the regulation of viral-infection-associated biological processes and signaling pathways. Moreover, the establishment of a ceRNA network further outlined the regulatory function of circRNAs that could potentially restrict or facilitate EV-A71 infection through regulation of gene expression, especially three circRNAs regulatory axes involved in MMP2 regulation. Thus, the results of this study may be helpful for future studies investigating the molecular functions of circRNAs in viral pathogenesis and virus-host interactions involved in CV-A16 infection.

Below is the link to the electronic supplementary material.

Supplementary Fig. S1 Growth curves of CV-A16 in SH-SY5Y cells at different time points after CV-A16 infection at an MOI of 1. (A) The kinetics of CV-A16 replication. (B) The kinetics of virus VP1 protein synthesis (TIF 866 KB)