Identification of microRNAs involved in the host response to enterovirus 71 infection by a deep sequencing approach.

Role of microRNA (miRNA) has been highlighted in pathogen-host interactions recently. To identify cellular miRNAs involved in the host response to enterovirus 71 (EV71) infection, we performed a comprehensive miRNA profiling in EV71-infected Hep2 cells through deep sequencing. 64 miRNAs were found whose expression levels changed for more than 2-fold in response to EV71 infection. Gene ontology analysis revealed that many of these mRNAs play roles in neurological process, immune response, and cell death pathways, which are known to be associated with the extreme virulence of EV71. To our knowledge, this is the first study on host miRNAs expression alteration response to EV71 infection. Our findings supported the hypothesis that certain miRNAs might be essential in the host-pathogen interactions.

1. Introduction

Hand, foot, and mouth disease (HFMD), a common febrile illness in children, is usually caused by human enteroviruses. Enterovirus 71 (EV71) and coxsackievirus A16 (CA16) are two major causative agents of HFMD. EV71 and CA16 infections manifesting as lesions on the skin and oral mucosa are clinically similar, but EV71 infection is more frequently associated with serious neurological diseases such as aseptic meningitis, encephalitis, and acute flaccid paralysis and fatalities [1-4] while the CA16-associated HFMD has a milder outcome [5]. More than 500,000 HFMD cases caused by EV71 were reported nationwide, including 176 fatal cases in China since March 2008 [6]. Elucidating the cellular events following EV71 infection will facilitate the development of strategies to prevent and treat this virus. However, the molecular mechanisms of the host response to EV71 infection are not completely understood.

microRNAs (miRNAs) have emerged as key regulators in many biological processes, from development to defense, at almost all organismal levels through mRNA degradation or translational repression of their targets [7, 8]. Recently, their role has been highlighted in pathogen-host interactions. Not only the miRNA encoded by viral genomes but host encoded miRNAs have been found participating in host-virus interactions. They can stimulate as well as suppress viral infections. For example, liver-specifc miR-122 is an indispensable factor in supporting hepatitis C virus (HCV) replication [9], whereas, miR-125b and miR-223 directly target human immunodeficiency virus-1 (HIV-1) mRNA, thereby attenuating viral gene expression in resting CD4+ T cells [10]. It has been showed that Epstein-Barr virus (EBV) encodes five miRNAs in its large DNA genome. These miRNAs can potentially regulate several genes encoding proteins involved in apoptosis, cell proliferation, signal transduction, transcription regulation, and immune response [11]. In contrast, the global changes, miRNAs expression during EV71 infection have not yet to be extensively elucidated. To determine which cellular microRNAs play a role in the host response to enterovirus infection, in this study, we performed a comprehensive miRNA profling in EV71-infected Hep2 cells through deep sequencing.

2. Materials and Methods

2.1. Cell Culture and Virus Infection

Human epidermoid carcinoma (Hep2) cells were grown in RPMI 1640 medium (Invitrogen) supplemented with 10% fetal calf serum (FCS). When the Hep2 cells were grown to 70% confluence in 25 cm2 flasks, they were infected with EV71 at a multiplicity of infection (m.o.i.) of 0.03 50% tissue culture infectious doses (TCID50) and maintained after infection at 37°C in RPMI 1640 medium with 2% FCS.

2.2. Isolation of RNA and miRNA

Hep2 cells were infected with EV71 as described above. At 6 h, 24 h, 48 h, 72 h, and 96 h post infection, total RNA or miRNA was extracted from cells using the mirVana kit according to the manufacturer's protocol (Ambion). RNA was also extracted from noninfected control cells at the same times. The quality and the concentration of the RNA samples were monitored by gel electrophoresis and absorbance at A260/280 ratio.

2.3. Library Construction, SOLiD Sequencing and Analysis

The miRNA sequencing library construction followed standard procedure of SOLiD small RNA expression kit (Applied Biosystem). All SOLiD run parameters followed standard Applied Biosystems protocols. Different barcodes were introduced to two samples in the polymerase chain reaction of library construction, and all the samples were sequenced in a single sequencing run.

2.4. Data Analysis

SOLiD data were first analyzed by SOLiD System Small RNA Analysis Pipeline Tool (RNA2MAP). The miRBase sequences (Sanger) of human being were downloaded from miRBase (http://www.mirbase.org/).The number of bases to use when generating initial seeds locations was 18 with a tolerance of 3 mismatches. After extension step, at most 6 mismatches were allowed in full length mapping.

Potential conserved target genes of differentially expressed miRNAs were firstly predicted by targetscan (http://www.targetscan.org/) [12-14]. Target genes of some miRNAs, such as mir-1972, mir-1974, mir-1975, mir-1979, and mir-764, could not be predicted in targetscan database, we further predicted those miRNAs targets using DIANA-microT v3.0 [15, 16]. In brief, each differentially expressed miRNA was submitted to targetscan individually and all of its targets predicted in targetscan or microT v3.0 were used for the following gene ontology (GO) analysis (http://www.babelomics.bioinfo.cipf.es/). All targets of induced and repressed miRNAs were submitted to FatiGO program [17]. Functional category enrichment based on the GO terms was evaluated on the targets of these differentially expressed miRNAs.

2.5. Confirmation of Differentially Expressed miRNAs by Real-time Quantitative RT-PCR

To confirm the expression of miRNAs by deep sequencing approach, stem-loop quantitative RT-PCR (qRT-PCR) was performed. In brief, cDNA was synthesized from total RNA by using AMV reverse transcriptase (TaKaRa). The 20 μl reactions were incubated for 15 min at 16°C, 30 min at 42°C and 5 min at 85°C, and then held at 4°C. Subsequently, real-time quantification was performed using an Applied Biosystems 7500 Sequence Detection system (Applied Biosystems). The 20 μl PCR reactions included 1 μl RT-PCR product, 10 μl Premix Ex Taq (TaKaRa), and 1 μl SYBR green (Invitrogen). The reactions were incubated in a 96-well optical plate at 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. All reactions were run in triplicate. After reaction, the threshold cycle value (CT) data were determined using default threshold settings, and the mean CT was determined from the duplicate PCRs. Human small nuclear RNA U6 was used for normalization. The expression levels of miRNAs were measured in terms of CT and normalized to U6 using 2−ΔΔCT [18].

2.6. Statistical Analysis

We firstly used the Z test to determine the statistical significance of the differences between the two libraries [19]. This approach is to look at the number of copies of a specifc miRNA per cell as a fraction or proportion of the total number of miRNA molecules in that cell. The same proportion of specifc tags should be present in the miRNA library of all sequenced tags. In this test a false discovery rate less than 5% was selected. miRNAs were considered significantly altered only when they fulflled three criteria: (1) mean fold change  >2  or <0.5, (2) having at least 10 copies by SOLiD sequencing, and (3) Z score>1.96 or <−1.96.

3. Results

3.1. Replication of EV71 in Hep2 Cells

Cytopathic effect (CPE) was first observed at 24 h post infection and progressed to moderate and severe CPE at 72 and 96 h, respectively. miRNA extracted from cells after 72 h infection was used for SOLiD deep sequencing.

3.2. Deep Sequencing

A total of 411,419 and 479,414 filtered high quality reads were obtained from infected and control cells by deep sequencing, respectively. After filtrating reads contaminated by rRNA, tRNA, snRNA, and snoRNA, 35,272 and 78,143 reads were obtained from infected and non-infected cells, respectively. Out of these reads, 11,188 of these high quality reads were exact matches while other 68,041 reads were loose matches to known human miRNAs (Table 1). Loose matches were defined by sequence reads that aligned with human miRNA consensus sequence with 1–4 mismatches. These may represent sequencing errors (when occurring in low copy numbers), mutations, and/or RNA editing events.

Table 1

Number of reads of miRNAs from EV71 infected and non-infected Hep2 cells.

	Infected cells	Control cells
High quality/both adapter	35,272	78,143
Exact match to known human miRNAs	2,731	8,457
Loose match to known human miRNAs	20,389	47,652

The size distribution of sequence reads showed that the majority of miRNAs was 18–25 for both libraries (>90%), with 22 nt small RNA being the most abundant (Figure 1), which is within the typical size range of human miRNAs.

Figure 1

Size distribution of sequenced short RNAs.

After reads were compared with an miRBase database (release 14.0), 569 miRNAs were detected in EV71-infected cells while 540 miRNAs were detected in non-infected control cells.

3.3. Aberrant miRNAs Expression in EV71-Infected Cells

On the basis of differentially expressed miRNA, we found 64 that miRNAs were differentially expressed between infected and non-infected cells. More miRNAs (42 out of 64 miRNAs) were upregulated than down-regulated during EV71 infection in Hep2 cells (Tables 2 and 3).

Table 2

Upexpressed miRNAs in infected cells compared with non-infected cells.

miRNA	Reads in infected	Reads in non-infected	Z score	Ratio infected/non-infected (Normalized)
hsa-mir-636	11	0	5.05958439	#
hsa-mir-619	305	41	23.59443534	17.31099765
hsa-mir-1302-8	387	54	26.46428291	16.67720702
hsa-mir-1302-2	1560	228	53.05155216	15.92193571
hsa-mir-1302-6,7	387	57	26.28151611	15.79945928
hsa-mir-1302-3	390	58	26.34977656	15.64741958
hsa-mir-1302-1	366	59	25.24116716	14.4356116
hsa-mir-1302-5	305	50	22.98061978	14.19501807
hsa-mir-1302-4	300	50	22.73564953	13.96231285
hsa-mir-1273	162	37	15.80907527	10.18871479
hsa-mir-1290	34	8	7.194848151	9.889971606
hsa-mir-1268	17	4	5.08718048	9.889971606
hsa-mir-1178	11	3	3.959263353	8.532524522
hsa-mir-1246	382	121	22.48247616	7.346561309
hsa-mir-1972	34	12	6.488332484	6.593314404
hsa-mir-1285-1	96	35	10.79491649	6.382771591
hsa-mir-566	13	5	3.901900612	6.05033557
hsa-mir-1226	10	4	3.376339014	5.817630356
hsa-mir-635	12	5	3.645008994	5.584925142
hsa-mir-627	16	7	4.132798028	5.318976326
hsa-mir-518c	13	6	3.647443806	5.041946309
hsa-mir-1185-2	20	10	4.373804587	4.654104285
hsa-mir-1289-1	14	7	3.659280977	4.654104285
hsa-mir-539	11	6	3.116008646	4.266262261
hsa-mir-1247	18	10	3.95066374	4.188693856
hsa-mir-764	12	7	3.147056335	3.989232244
hsa-mir-1979	330	201	16.15562506	3.820533368
hsa-mir-1272	11	7	2.873705427	3.656796224
hsa-mir-324	28	18	4.558616136	3.619858888
hsa-mir-1975	829	548	24.5078521	3.520303332
hsa-mir-744	28	19	4.414239693	3.429339999
hsa-mir-541	19	16	3.126630669	2.763374419
hsa-mir-1254	13	11	2.576283222	2.750152532
hsa-mir-421	42	38	4.379270931	2.572005
hsa-mir-1237	11	10	2.231381474	2.559757357
hsa-mir-320c-1	74	77	5.081196578	2.236387773
hsa-mir-615	16	17	2.309101549	2.190166722
hsa-mir-720	14	15	2.140061026	2.171915333
hsa-mir-140	184	208	7.299247353	2.058546126
hsa-mir-451	188	214	7.31642594	2.044326181
hsa-mir-21	3445	3992	31.35678948	2.008190038

Table 3

Down-expressed miRNAs in infected cells compared with non-infected cells.

miRNA	Reads in infected	Reads in non-infected	Z score	Ratio infected/non-infected (Normalized)
hsa-mir-584	0	10	−2.07306	0
hsa-mir-221	4	70	−4.6261	0.132974
hsa-let-7e	23	199	−6.40515	0.268956
hsa-mir-1180	5	43	−2.96815	0.270587
hsa-mir-1259	6	46	−2.91293	0.303529
hsa-mir-338	128	905	−12.4248	0.32913
hsa-mir-30a	304	2083	−18.6008	0.339618
hsa-mir-19b-1	487	2952	−20.5604	0.383901
hsa-mir-19b-2	491	2946	−20.3931	0.387842
hsa-mir-545	7	41	−2.33845	0.397302
hsa-mir-433	6	35	−2.15407	0.398923
hsa-mir-582	93	531	−8.27217	0.407563
hsa-mir-26a-1	129	735	−9.72425	0.408421
hsa-mir-26a-2	132	746	−9.73602	0.411757
hsa-mir-27b	123	677	−9.08034	0.422788
hsa-mir-452	31	158	−4.09666	0.456574
hsa-mir-1974	38	188	−4.34508	0.470362
hsa-mir-27a	280	1382	−11.8087	0.471472
hsa-mir-26b	35	172	−4.12863	0.473528
hsa-mir-30e	86	419	−6.39521	0.477629
hsa-mir-222	50	236	−4.64215	0.49302
hsa-mir-660	48	225	−4.4992	0.496438

3.4. Confirmation of Differentially Expressed miRNAs

Quantitative RT-PCR assays were used to confirm the expression pattern of differentially expressed miRNAs in Hep2 cells. There was general consistency between quantitative RT-PCR assay and deep sequence analysis in four miRNAs (miR-1246, miR-1237, miR-30a, and miR-222) in terms of directions of regulation and significance. Specifically, there was a 1.92-fold upregulation (7.35-fold in deep sequencing analysis) in miR-1246, 1.58-fold upregulation (2.58-fold in deep sequencing analysis) in mir-1237, (P < .05), 2.45-fold down-regulation (2.94-fold in deep sequencing analysis) in miR-30a, and 1.42-fold down-regulation (2.02-fold in deep sequencing analysis) in miR-222.

3.5. Gene Ontology Analysis

Targets were predicted for all identified differentially expressed miRNA families. In total 5765 unique target genes were predicted for 64 of the differentially expressed miRNAs ( see Table S1 in Supplementray Material available online at doi:115/2010/425939). On the basis of the biological functions described by FatiGO program (http://www.babelomics.bioinfo.cipf.es/), these target genes can be grouped into 72 categories (S2). The top 30 Gene Ontology terms are shown in Figure 2. The majority of targets fall into the category of metabolic process, regulation of biological process, and cell communication indicating intense biological change in Hep2 cells after EV71 infection. Several other groups contain genes regulating death (255 target genes including 239 apoptosis-related target genes), neurological process (211 target genes), and immune response (120 target genes).

Figure 2

Predicted targets gene ontology terms in the biological process category.

Hundreds of altered transcripts in response to EV71 infection were found in two transcriptomic studies [20, 21]. To determine whether miRNAs might be modulators of mRNAs that were differentially expressed, we investigated whether those differentially expressed mRNAs are enriched for predicted targets and interrogated their inversely correlated targets for functional associations. To our surprise, targets for differentially expressed miRNAs or miRNA families in the present study account for 12.1% (19 out of 157) of the differentially expressed mRNAs in one study [20] whereas they account for 22.4% (13 out of 58) of the transcripts that were significantly altered in another study (Tables 4 and 5) [21].

Table 4

Reverse regulatory association of up-regulated miRNAs and their predicted mRNA targets.

Up-regulated miRNA	Target genes	Functions	Reference
hsa-mir-1289-1	Interferon regulatory factor 2	IFN response	[20]
hsa-mir-140,hsa-mir-320c	v-yes-1 Yamaguchi sarcoma viral oncogene homolog 1	Oncogene	[20]
hsa-mir-421	v-Ki-ras2 Kirsten rat sarcoma viral oncogene homolog		[21]
hsa-mir-1226	Interleukin 10 receptor, beta	Immune response	[20]
hsa-mir-1290	Nuclear receptor coactivator 2	Transcriptional regulators	[20]
hsa-mir-320c,hsa-mir-1290	Transcription factor AP-2 beta		[20]
hsa-mir-320c	pre-B-cell leukemia homeobox 3		[21]
hsa-mir-636,hsa-miR-1979,hsa-mir-1302,hsa-mir-518c,hsa-mir-1226,hsa-mir-1290,hsa-mir-421,hsa-mir-21,hsa-mir-140	ribosomal protein S6 kinase	Kinases and phosphatases	[20]
hsa-mir-1289-1,hsa-mir-1272,hsa-mir-421	Splicing factor, arginine/serine-rich 1	RNA synthesis and modification	[20]
hsa-mir-1272,hsa-mir-140,hsa-mir-636	Neuro-oncological ventral antigen 1	Neuron specific	[20]
hsa-mir-539	DEAD (Asp-Glu-Ala-Asp) box polypeptide 5		[21]
hsa-mir-1254,hsa-mir-627,hsa-mir-320c,hsa-mir-421,hsa-mir-21,hsa-mir-140	cyclin-dependent kinase 6	cell cycle	[21]
hsa-mir-619	protocadherin 9	cell adhesion	[21]
hsa-mir-627	RER1 retention in endoplasmic reticulum 1 homolog (S. cerevisiae)	trafficking and targeting proteins	[21]
hsa-mir-140	sorting nexin 2		[21]
hsa-mir-1289-1	xenotropic and polytropic retrovirus receptor	Cellular receptors	[21]
hsa-mir-21	interleukin 12A (natural killer cell stimulatory factor 1, cytotoxic lymphocyte maturation factor 1, p35)	cell signalling, extracellular communication	[21]
hsa-mir-1178,hsa-mir-140	transforming growth factor, alpha		[21]
hsa-mir-1178	dual specificity phosphatase 6	intracellular trasducers/effector/modulators	[21]

Table 5

Reverse regulatory association of down-regulated miRNAs and their predicted mRNA targets.

Down-regulated miRNA	Target genes	Functions	Reference
hsa-mir-27	Cytochrome P450,superfamily I polypeptide 1	Mitochondrial function	[20]
hsa-mir-26a,hsa-mir-27	BCL2-antagonist/killer 1	Apoptosis	[20]
hsa-mir-26a	Programmed cell death 10		[20]
hsa-mir-582	Amyloid beta (A4) precursor protein		[20]
hsa-mir-19	BCL2-associated athanogene 5		[20]
hsa-mir-1180,hsa-mir-19,hsa-mir-582,hsa-mir-26a,hsa-mir-30a,hsa-let-7e,hsa-mir-545	calcium channel, voltage-dependent, L type, alpha 1C subunit	Membrane transporters	[20]
hsa-mir-433,hsa-let-7e	Cathepsin C	Protein degradation	[20]
hsa-mir-26a,hsa-mir-545	Muscle RAS oncogene homolog	Oncogene	[20]
hsa-mir-19	RAB5B, member RAS oncogene family		[20]
hsa-mir-19	neurotrophic tyrosine kinase, receptor, type 2	Kinases and phosphatases	[20]
hsa-mir-19	Fas-activated serine/threonine kinase		[20]
hsa-mir-582,hsa-let-7e	adaptor-related protein complex 1, sigma 1 subunit	trafficking and targeting proteins	[21]
hsa-mir-27	actin, alpha 2, smooth muscle, aorta	cytoskeleton/motility	[21]

4. Discussion

Many technologies have been developed for miRNA profiling, including real-time quantitative RT-PCR [22], northern blotting [23, 24], and microarray analyses based on either direct hybridization or hybridization coupled with enzymatic extension [25, 26]. These methods have been used successfully in a variety of studies. But they still have some technical limitations. For example, first two methods were not high-throughput while microarray method needs large amounts of starting materials. Recent progress in high-throughput sequencing technologies allows deep sequencing of large libraries of short RNAs [27-30]. The longest reads are obtained by the 454 technology, which currently gives reads of about 400 base pairs (bp). However, this technology yields much less reads than other techniques (about 400 000 per sample) [27]. The Solexa (Illumina) and SOLiD platform (Applied Biosystems) generates shorter reads (up to 35 bp) but yields 1–3 million reads per sample [28, 29]. The other high-throughput technique, such as massively parallel sequencing (MPSS), gives more reads than Solexa but the reads are even shorter, only 17 bp [30]. In the present study, we used SOLiD platform to explore differential miRNA profling of Hep2 cells response to EV71 infection. We obtained about 0.4 million reads in both samples, well below the capacity of SOLiD. One possibility is that we used cell line as infection model but not tissue. Numbers of small RNAs were not expressed in Hep2 cells. Even now, we found that there was general consistency between quantitative RT-PCR assay and deep sequence analysis. So, SOLiD deep sequencing successfully revealed miRNA profling in EV71-infected and control Hep2 cells.

Based on a comprehensive examination of miRNA expression from EV71-infected and non-infected control Hep2 cells, we identified 64 miRNAs that were differentially expressed, with most of them (65.6%) upregulated in EV71-infected cells. The altered patterns of cellular miRNAs we observed for EV71-infected cells are similar in some respects to changes seen for cells infected with Hepatitis C virus (HCV) [31] but contrast with those with Epstein-Barr virus (EBV) [32] and human cytomegalovirus (HCMV) in which more miRNAs were dow-regulated in response to virus infection [33]. Most of the miRNAs affected by EV71 are different from those affected by HCV, EBV, and HCMV. We observed miR-636 and miR-584 expressed in only infected or non-infected cells. Though their functions during virus infection have not been explored, their predicted target genes have been identified to be involved in virus entry, replication and propagation. For example, reticulon 3 (RTN3), one predicted target gene of miR-636, can bind the 2C protein of enterovirus 71 and is required for viral replication, [34]. Poliovirus receptor-related 1 (herpesvirus entry mediator C) (PVRL1, also known as nectin-1), another target gene of miR-636, can serve as receptor for herpes simplex virus and pseudorabies virus entry [35]. Abhydrolase domain containing 2 (ABHD2), one predicted target gene of miR-584, is essential for Hepatitis B virus propagation [36]. Thus, it will be important to investigate the mechanisms of regulation of miRNA levels during virus infection, which could be at the stages of transcription, maturation, and/or degradation.

The host response to viral infection represents complex orchestration of divergent pathways deigned to eliminate the virus and protect the host. Viruses impact on many aspects of the host cell's biology and function [20]. As seen from gene ontology analysis, the top 3 gene ontology terms related to metabolic process which indicated EV71 infection have enormous effect on Hep2 cells metabolism. Intriguingly, neurological process, apoptosis, and immune response related GO terms in biological process were also enriched from the predicted targets. Considering the fact that EV71 can cause lethal encephalitis or myocarditis, both apoptosis and immune response contribute to EV71 pathogenesis. Myocarditis, for example, represents an intricate interplay between virus and patient responses, in which both direct viral injury and immunopathologic damage caused by innocent-bystander phenomena affect the disease course [37]. Further study of the functions of those underlying miRNAs related to neurological process, apoptosis, and immune response will help to elucidate the molecular mechanisms of EV71 pathogenesis.

Our study has some limitations that will need to be addressed in future studies. We did not assess the roles in infection of miRNAs whose expression were altered after infection. As microRNAs predominately function as repressors of target gene expression, we indeed found numbers of targets in an another transcriptomic research whose expression was inversely correlated with the expression of dysregulated miRNAs in the present study. Although all assays were executed when cells demonstrated similar CPE, the virus may interact with miRNA regulatory pathways differently in different cell types in which human neural SF268 cells and rhabdomyosarcoma cells were used in those studies, respectively. Even then, these findings suggest that expression modification of host miRNAs during EV71 infection could be related to a number of cellular physiological processes that eventually control the cell fate.

In summary, in this study, we identified the miRNAs involved in the host response to EV71 infection using deep sequencing technology. Our findings provide a deeper understanding of the mechanisms underlying EV71 infection. Once the role of these miRNAs in the regulation of host-EV71 interaction has been determined, it will improve the protection and treatment strategies in enterovirus infection.

Supplementary materials contains “Gene Ontology Analysis”. Targets were predicted for all identified differentially expressed miRNA families. In total 5765 unique target genes were predicted for 64 of the differentially expressed miRNAs.

Click here for additional data file.

Click here for additional data file.