lncRNA expression signatures in response to enterovirus 71 infection.

Outbreaks of hand, foot, and mouth disease caused by enterovirus 71 (EV71) have become considerable threats to the health of infants and young children. To identify the cellular long noncoding RNAs (lncRNAs) involved in the host response to EV71 infection, we performed comprehensive lncRNA and mRNA profiling in EV71-infected rhabdomyosarcoma cells through microarray. We observed the differential expression of more than 4800 lncRNAs during infection. Further analysis showed 160 regulated enhancer-like lncRNA and nearby mRNA pairs, as well as 313 regulated Rinn's lncRNA [M. Guttman I. Amit, M. Garber, C. French, M.F. Lin, D. Feldser, M. Huarte, O. Zuk, B.W. Carey, J.P. Cassady, M.N. Cabili, R. Jaenisch, T.S. Mikkelsen, T. Jacks, N. Hacohen, B.E. Bernstein, M. Kellis, A. Regev, J.L. Rinn, E.S. Lander. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature 458 (2009) 223-227, A.M. Khalil, M. Guttman, M. Huarte, M. Garber, A. Raj, D. Rivea Morales, K. Thomas, A. Presser, B.E. Bernstein, A. van Oudenaarden, A. Regev, E.S. Lander, J.L. Rinn. Many human large intergenic noncoding RNAs associate with chromatin-modifying complexes and affect gene expression. Proc. Natl. Acad. Sci. USA 106 (2009) 11667-11672] and nearby mRNA pairs. The results provided information for further research on the prevention and treatment of EV71 infection, as well as on distinguishing severe and mild EV71 cases.

Introduction

Enterovirus 71 (EV71) is a typical positive-strand RNA virus belonging to the Picornaviridae family [1]. EV71 infection is a major cause of hand, foot, and mouth disease (HFMD) in infants and young children [2], [3] This infection can also cause neurological diseases such as aseptic meningitis, encephalitis, and acute flaccid paralysis, which can lead to permanent paralysis and even death [4], [5]. Outbreaks of EV71 infection have been reported in many countries [6]. In recent years, EV71 infection has become a considerable threat to public health in China [7], where the government has reported 1,619,706 cases of HFMD (with 509 deaths) in 2011 and 1,587,849 cases (with 463 deaths) from January to August 2012 (http://www.moh.gov.cn/publicfiles//business/htmlfiles/mohjbyfkzj/s3578/list.htm). However, to date, no effective vaccine or therapy is available to prevent or treat this infection. In humans, cellular immunity is important in preventing the development of serious complications after EV71 infection [8], [9]. Thus, understanding the cellular events after EV71 infection can facilitate the development of new strategies for preventing and treating this infection.

In 2004, Shih [10] reported that EV71 infection leads to increased levels of mRNAs encoding chemokines, proteins involved in protein degradation, complement proteins, and proapoptotis proteins. The infection also results in the decreased expression of genes encoding proteins involved in host RNA synthesis. In 2006, Leong et al. [11] identified 152 down-regulated genes and 39 up-regulated genes in rhabdomyosarcoma (RD) cells infected with EV71. In 2010, 64 microRNAs were up- or down-regulated more than twofold in response to EV71 infection [12].

In the last decade, long noncoding RNAs (lncRNAs) have been shown to play important roles in gene expression regulation, dosage compensation, genomic imprinting, nuclear organization and compartmentalization, as well as nuclear-cytoplasmic trafficking [13], [14], [15], [16], [17]. Recent studies have demonstrated the changes in host lncRNA expression in response to virus infection. After infectious bursal disease virus and Marek’s disease virus infection in chicken, eight and two lncRNAs are differently expressed, respectively [18]. In 2010, Peng et al. [19] reported the differential expression of more than 500 lncRNAs in mice after severe acute respiratory syndrome coronavirus (SARS-CoV) infection. Other studies have shown that most lncRNAs are similarly regulated in response to influenza virus infection, and that they have distinctive kinetic expression profiles in type I interferon receptor and STAT1 knockout mice during SARS-CoV infection [19]. These findings suggest the widespread differential expression of lncRNAs in response to virus infection and their involvement in regulating the host response, including innate immunity [19]. To determine which cellular lncRNAs play roles in the host response to EV71 infection, we performed lncRNA and mRNA microarray analyses in mock- and EV71-infected RD cells.

Materials and methods

Cell culture and virus infection

Human RD cells were grown in Eagle’s minimum essential medium (MEM; Gibco) supplemented with 10% fetal bovine serum (FBS; Gibco). When the cells had grown to 90% confluence in 25 cm2 flasks, they were infected with EV71 (GDFS-3; isolated and identified from Guangdong Province, China in 2008) [20] at a multiplicity of infection of 100 50% tissue culture infectious doses (TCID50). After adsorption for 1 h at 37 °C, the inoculum was removed and Eagle’s MEM with 2% FBS was added. The culture was maintained at 37 °C.

Isolation of RNA

At 24 h post-EV71 infection, the medium was removed, an appropriate volume of buffer RLT was added, and the total RNA was extracted using an RNeasy Mini Kit according to the manufacturer’s protocol (Qiagen). RNA was also extracted from mock-infected cells. The quality and the concentration of the RNA samples were monitored at absorbance ratios of A260/A280 and A260/230 using a NanoDrop ND-1000 spectrophotometer and standard denaturing agarose gel electrophoresis.

RNA labeling and array hybridization

Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology) with minor modifications. Briefly, mRNA was purified from total RNA after the removal of rRNA (mRNA-ONLY™ Eukaryotic mRNA Isolation Kit, Epicentre). Each sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without a 3′ bias using the random priming method. The labeled cRNAs were purified by using an RNeasy Mini Kit (Qiagen). The concentration and specific activity of the labeled cRNAs (pmol Cy3/μg cRNA) were measured using NanoDrop ND-1000. About 1 μg of each labeled cRNA was fragmented by adding 5 μl of 10× blocking agent and 1 μl of 25× fragmentation buffer, followed by heating the mixture at 60 °C for 30 min. Finally, 25 μl of 2× GE hybridization buffer was added to dilute the labeled cRNA. About 50 μl of hybridization solution was dispensed into the gasket slide and assembled onto an Arraystar Human lncRNA Array v2.0 slide (Arraystar, USA). The slides were incubated for 17 h at 65 °C in an Agilent hybridization oven. The hybridized arrays were washed, fixed, and scanned using an Agilent DNA Microarray Scanner (part number G2505C).

Data analysis

Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.0 software package (Agilent Technologies). After quantile normalization of the raw data, lncRNAs and mRNAs with “Present” or “Marginal” (“All Targets Value”) flags in mock- and EV71-infected samples were subjected to further data analysis. Differentially expressed lncRNAs and mRNAs between the two samples were identified through fold-change (greater than twofold) filtering.

LncRNA classification and subgroup analysis

lncRNAs with enhancer-like function were identified using the GENCODE annotation [21] of human genes [22]. The selection of lncRNAs with enhancer-like function involved the exclusion of transcripts mapped to the exons and introns of annotated protein coding genes, as well as the natural antisense transcripts overlapping with the protein coding genes and all known transcripts. Rinn’s lncRNAs were identified based on the studies of Rinn [23], [24]. Human homeobox transcription factors (HOX) cluster lncRNAs were also identified based on the study of Rinn [25].

Confirmation of differentially expressed lncRNAs by real-time quantitative RT-PCR (qPCR)

qPCR was performed to confirm the expression of lncRNAs by microarray analysis. Briefly, cDNA was synthesized from total RNA using a PrimeScript RT reagent kit with a gDNA Eraser (TaKaRa). Primers for four lncRNAs were designed and synthesized (Table 1 ). Then, qPCR was performed using a LightCycler 480 (Roche Applied Science). The 10 μl PCR reactions included 1 μl of cDNA product and 5 μl of SYBR Premix Ex Taq II (TaKaRa). The reactions were incubated at 95 °C for 1 min, followed by 40 cycles at 95 °C for 5 s, 60 °C for 10 s, and 72 °C for 15 s. 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 18S rRNA was used for normalization. The expression levels of lncRNAs were measured in terms of CT, and then normalized to 18S using 2−ΔΔCT [26].

Table 1

Primers used for qPCR.

lncRNAs	Sense primer (5′-3′)	Anti-sense primer(5′-3′)	Product (bp)
AP000688.29	ttgtgctgctaaccactgagact	aggataatggtgcttgcttgac	87
AC002511.1	tttacaagcatcagccaccac	agatccagcatgaggaaccc	93
RP5-843L14.1	ttgtggagacaggatttggacc	gcaggaaccgggacttgaa	84
RP4-620F22.3	cccacctgtgaagtgaagcc	ttcctcatcaagtgagaagggtt	102
GAPDH	acccactcctccacctttgac	accaccctgttgctgtagcc	107

Results

lncRNA and mRNA microarray data

Cytopathic effects were observed 24 h post-infection. Total RNA was extracted from mock- and EV71-infected cells 24 h post-infection. The OD260/OD280 ratios were approximately 2.1, and the OD260/OD230 ratios were more than 1.9, which suggested that the total RNAs were sufficiently pure for the succeeding experiments. Subsequently, mRNA was purified, cRNA was prepared, and array hybridization was performed using Arraystar Human LncRNA Array v2.0. After quantile normalization of the raw data, the expression profiles of 22971 lncRNAs and 18194 mRNAs were obtained from mock- and EV71-infected cells (Tables S1 and S2). The distributions of the log2 ratios of lncRNAs and mRNAs between EV71- and mock-infected samples were almost the same. Fig. 1 shows the heat maps of the expression ratios (log2 scale) of lncRNAs and mRNAs.

Fig. 1

Heat maps of the EV71-/mock-infected expression ratios (log2 scale) of lncRNAs (A) and mRNAs (B) in RD cells. “Red” denotes high relative expression and “blue” denotes low relative expression. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Aberrant lncRNA and mRNA expression in EV71-infected cells

To identify differentially expressed lncRNAs and mRNAs, we performed fold-change filtering between mock- and EV71-infected cells (fold change > 2.0). We found that 2990 lncRNAs and 1718 mRNAs were up-regulated, whereas 1876 lncRNAs and 2552 mRNAs were down-regulated in EV71-infected cells (Tables S3 and S4).

Confirmation of some differentially expressed lncRNAs

We performed qPCR assays to confirm the expression pattern of four differentially expressed lncRNAs in RD cells. A general consistency between the qPCR and microarray analysis results was confirmed in four lncRNAs (AP000688.29, AC002511.1, RP5-843L14.1, and RP4-620F22.3) in terms of regulation direction and significance. Specifically, a 3.31-fold down-regulation (2.25-fold in microarray analysis) was observed in AP000688.29, 3.33-fold up-regulation (2.77-fold in microarray analysis) in AC002511.1, 2.29-fold up-regulation (2.40-fold in microarray analysis) in RP5-843L14.1, and 2.99-fold up-regulation (2.22-fold in microarray analysis) in RP4-620F22.3 (Fig. 2 ).

Fig. 2

Comparison between microarray data and qPCR results. AP000688.29, AC002511.1, RP5-843L14.1, and RP4-620F22.3 differentially expressed in EV-71-infected cells compared with mock-infected cells by microarray were validated by qPCR. The heights of the columns in the chart represent the log-transformed median fold changes (T/N) in the expression between EV71- and mock-infected cells, and the bars represent standard errors. The validation results of the four lncRNAs indicated that the microarray data well correlated with the qPCR results.

lncRNA classification and subgroup analysis

Rinn et al. [25] identified numerous transcripts from the four HOX loci, which included 101 mRNA, 75 introns, and 231 intergenic transcripts. These lncRNAs, which were expressed in temporal and site-specific manners, possibly used the same enhancers as HOX genes and may have the same global regulating functions as HOX. In this paper, the profiling data of all probes targeting these 407 discrete transcripts are presented in Table S5. The data showed that 43 mRNAs can be detected in RD cells, with seven of them being differentially expressed. Then, 300 transcribed noncoding RNAs (including introns and intergenic transcripts) were detected, with 64 of them being differentially expressed.

Using chromatin-state maps, Rinn’s studies [23], [24] have identified 3019 lncRNAs with clear evolutionary conservation and association with distinct and diverse biological processes, such as cell proliferation, RNA binding complexes, immune surveillance, embryonic stem cell pluripotency, neuronal processes, morphogenesis, gametogenesis, and muscle development. The profiling data of all probes for these lncRNAs are provided in Table S6, which indicated that 477 from the detected 2200 lncRNAs were differentially expressed. Among them, 190 were down-regulated and 287 were up-regulated. Further analysis resulted in 313 differentially expressed lncRNAs and nearby coding gene pairs (distance < 300 kb) for each comparison between mock- and EV71-infected cells (Table S7). Among the 163 pairs, lncRNAs and nearby coding genes were regulated in the same direction (up or down), whereas 150 pairs were regulated in the opposite direction.

In 2010, using the GENCODE annotation [21] of human genes, Orom et al. [22] identified a set of lncRNAs with enhancer-like function from human cell lines. The depletion of these lncRNAs resulted in a concomitant decrease in the expression of neighboring genes. The profiling data of all probes for lncRNAs with enhancer-like function are shown in Table S8; 1025 enhancer-like lncRNAs were detected and 252 of them were differentially expressed (fold change > 2). Among these 252 lncRNAs, 160 had differentially expressed nearby coding genes (distance < 300 kb) for each comparison, as shown in Table S9. In 70 of the 160 pairs, lncRNAs and nearby coding genes were regulated in the same direction (up or down); in 90 of the 160 pairs, they were regulated in the opposite direction.

Nearby coding gene function analysis

Inferring the putative functions of protein-coding genes located near lncRNAs is an important approach to lncRNA research [23], [27]. In this paper, we combined differentially expressed nearby mRNA pairs of 313 differentially expressed Rinn’s lncRNAs and 160 enhancer-like lncRNAs, and abandoned the duplicated mRNAs. Then, we used the DAVID Functional Annotation Chart [28], [29] for functional enrichment analysis of the differentially regulated protein-coding gene and lncRNA pairs. The most significant functional groups consisted of annotation terms of alternative splicing, splice variant, phosphoprotein, nucleus, cytoplasm, and acetylation (Fig. 3 ). We hypothesized that the lncRNAs can modulate host responses through nearby protein-coding genes.

Fig. 3

Functional enrichment analysis on differently regulated mRNAs which were with differently expressed nearby lncRNAs. The mRNAs were from the 160 regulated enhancer-like lncRNA and nearby mRNA pairs, and 313 regulated Rinn’s lncRNA and nearby mRNA pairs. The functional enrichment analysis was performed by utilizing the DAVID Functional Annotation Chart [28], [29].

Discussion

Continuous interactions between viruses and hosts during their co-evolution have shaped their immune system. Consequently, viruses have manipulated host immune-control mechanisms to facilitate their propagation. Previous studies on virus–host interactions and viral pathogenesis have largely focused on protein-coding genes. In the last decade, evidence of host-cellular microRNA modulation of the expression of various viral genes has been reported. This modulation plays a pivotal role in the host–pathogen interaction network [30]. Recent studies have shown that virus infection alters the expression profiles of host lncRNAs. For example, eight lncRNAs are differentially expressed in virus-infected birds [18]. In 2010, comprehensive deep sequencing showed that more than 500 lncRNAs were differentially expressed in mice after SARS-CoV infection [19]. However, lncRNAs in other virus infections are not well documented. In the present study, using Arraystar microarray analysis, we identified the differentially expressed lncRNAs in RD cells after EV71 infection, together with nearby differentially expressed mRNA pairs.

In a study on the molecular mechanisms of host response to EV71 infection, Hsih et al. [11] found the up-regulation of mRNAs encoding chemokines, complement proteins, proteins involved in protein degradation, and proapoptosis proteins. They also observed the down-regulation of several genes encoding proteins involved in host RNA synthesis in EV71-infected SF268 cells. Further investigations have shown that EV71 infection alters the transcription of genes encoding components of cytoskeleton, protein translation, and protein modification; cellular transport proteins; protein degradation mediators; cell death mediators; mitochondrial-related and metabolism proteins; as well as cellular receptors and signal transducers in RD cells. Recent microRNA profiling analysis in Hep2 cells has identified 64 microRNAs whose expression levels changed more than twofold in response to EV71 infection. Their potential conserved target genes encode proteins with functions in neurological processes, immune responses, and cell death pathways. These proteins are known to be associated with the extreme virulence of EV71 [12]. In this study, lncRNA and mRNA expression analyses identified differentially expressed lncRNA and mRNA pairs. Functional enrichment analysis further indicated that the mRNAs were associated with alternative splicing, splice variant, phosphoprotein, nucleus, cytoplasm, and acetylation. These results and those of future works can expand the molecular mechanisms of the host response to EV71 infection.

Some lncRNAs reportedly serve as enhancers and have positive effects on gene expression. Evf-2 ncRNA forms a complex with the homeodomain-containing protein Dlx2 and causes transcriptional enhancement [31]. The binding of heat-shock RNA-1 ncRNA with heat-shock transcription factor 1 leads to the induction of heat-shock proteins [32]. An isoform of ncRNA steroid receptor RNA activator is also co-activated with steroid receptor responsiveness [33]. Recently, Orom et al. [22] showed that noncoding RNA-activating 1–7 enhance the expression of nearby genes. In our study, we identified 160 differentially expressed enhancer-like lncRNA and mRNA pairs, and 4l.5% (70/160) of these pairs were regulated in the same direction. We speculated that some of these lncRNAs function as enhancers that activate nearby genes; however, further research is needed to prove this hypothesis.

Identifying the putative functions of nearby genes of lncRNAs may aid in understanding the functional roles of lncRNA [23], [27]. Peng et al. [19] performed functional enrichment analysis on the nearby protein-coding genes of differentially expressed lncRNAs in SARS-CoV infected mouse. They found that the most significant functional group consisted of annotation terms related to gene expression, including transcription regulation, nuclear and DNA-binding transcription factor activity, as well as regulation of RNA metabolic process. In this study, functional enrichment analysis of differentially expressed mRNAs with differentially expressed nearby lncRNA partners showed that they were functionally related with alternative splicing, splice variant, phosphoprotein, nucleus, cytoplasm, and acetylation. Although the functions of lncRNAs during virus infection have not been explored, we speculate that EV71 infection may change some lncRNA expression levels that further regulate the expression of proteins related with alternative splicing and signal transduction.

EV71 infection may be asymptomatic or may cause diarrhea, rashes, and HFMD. EV71 can also cause severe neurological disease [4], [5]. Distinguishing mild and severe cases of EV71 infection in the early infection phase is very important in determining the appropriate treatment process. Recent reports have shown that the expression levels of five microRNAs significantly increase in Coxsackievirus A16 (CAV16)-infected patients compared with EV71-infected ones. The combination of three microRNAs also shows a moderate ability to differentiate between CVA16 and EV71. These data indicate that microRNA expression profiles can serve as supplementary biomarkers for diagnosing and classifying enteroviral HFMD infections [34]. We suggest that further studies on the putative differential expression profiles of lncRNAs in different EV71-infected cases can help distinguish mild and severe cases.

In summary, we identified the lncRNAs putatively involved in the host response to EV71 infection, which provided deeper insight into the mechanisms underlying EV71 infection. After determining the role of lncRNAs in the regulation of host-EV71 interactions, the protection and treatment methods for EV71 infection can be improved, and severe cases can be distinguished from mild ones in the earlier phase.