Characterization of Long Non-coding RNA Profile in RD and SH-SY5Y Cells Infected by CV-B5 using RNA Sequencing


 Hand, foot, and mouth disease caused by Coxsackievirus B5 (CV-B5) poses considerable threats to the health of infants especially in neurological damage. And the long noncoding RNAs (lncRNAs) act pivotal factors in regulating and participating in virus-host interactions. However, the role of lncRNAs in CV-B5-host interactions has not yet been elucidated. In this study, we used the RNA sequencing to determine the expression profiles of lncRNAs in CV-B5 infected human rhabdomyosarcoma (RD) and SH-SY5Y cells. Our results identified that in the differentially expressed lncRNAs a total of 508 up-regulated and 760 down-regulated lncRNAs in RD cell, with 46.2% were lincRNAs, 28.6% were anti-sense lncRNAs, 24.1% were sense overlapping lncRNAs, and 1.0% were sense intronic lncRNAs. Moreover, 792 lncRNAs were significantly increased and 811 lncRNAs were greatly decreased in SH-SY5Y cell including 48.6% were lincRNAs, 34.7% were anti-sense lncRNAs, 16.0% were sense overlapping lncRNAs, and 0.8% were sense intronic lncRNAs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway showed that the differentially expressed lncRNAs participated in the occurrence of disease in RD cell and associated with signaling pathway in SH-SY5Y cell after CVB5 infection respectively. In addition, similar results were obtained when seven lncRNAs were selected for validation using RT-qPCR assays. Moreover, we conducted the candidate lncRNA-IL12A secondary structures and found that it inhibits viral replication through Wnt signaling pathway. Our results reveal that lncRNAs can become a possible novel molecular target for the prevention and treatment of CV-B5 infection, and provided information for distinguishing neurogenic CV-B5 disease.


Introduction
Hand, foot, and mouth disease (HFMD) is a common contagious disease among young children especially under the age of ve years. HFMD usually occurs in the spring to fall and follows a mild and self-limiting course, such as fever and papulovesicular rash on the hands and feet, mouth [1]. Severe complications such as the central nervous system and signi cant mortality have also been observed in a small proportion of children suffering from HFMD [2,3]. The main pathogens of HFMD are enterovirus and coxsackievirus. Studies conducted in the past ten years had reported that enterovirus A71 (EV-A71) and coxsackievirus A16 (CV-A16) were the two major pathogens worldwide [4,5]. However, in recent years the incidence of HFMD caused by coxsackievirus B5 (CV-B5) has increased. CV-B5 had been widely associated with both sporadic cases and outbreaks of HFMD worldwide, particularly in China, South East Asia and Europe [6][7][8]. CV-B5 infection not only caused severe HFMD, but also the neurological damages, including aseptic meningitis, viral encephalitis, and acute accid paralysis [9][10][11]. CV-B5 belongs to the Enterovirus genus of the Picornaviridae family and is a positive-sense single-stranded RNA with approximately 7500 bp. However, there is no effective vaccine for the diseases caused by CV-B5 nor the mechanism of CV-B5 association with the nervous system has yet been elucidated.
Long non-coding RNAs (lncRNAs) are non-coding RNAs (ncRNAs) that are greater than 200 nucleotides in length and do not have protein coding ability. They are classi ed into several categories including sense, antisense, bidirectional, intronic and intergenic lncRNA based on the location of lncRNAs and target gene on chromosome [12,13]. LncRNAs are versatile molecules that can interact physically and functionally with DNA, other RNAs, and proteins, either through nucleotide base pairing or via formation of structural domains generated by RNA folding [14]. Previous studies have reported that lncRNAs regulate various aspects of an organism's life including the cell cycle and differentiation, epigenetic regulation of genome expression, transcription and post-transcriptional mRNA regulation, and protein translation [15,16].
High throughput RNA sequencing has proved that lncRNAs also regulated and participated in virus-host interactions. Some studies have also reported that lncRNAs are involved in innate immune responses by regulating innate immune cells and transcription of gene expression programs [17,18]. In addition, they regulated viral life cycle and modulated host transcription, thereby promoting viral replication. Indeed, some of the earliest characterized lncRNAs were those transcribed from viral genomes [19,20]. LncRNAs can also affect downstream signaling pathways and change the interaction between host and virus [21,22]. Therefore, many researches have showed that lncRNAs were associated with virus replication and infection. However, the function of lncRNAs in CV-B5 infection is still unknown. Therefore, understanding the relationship between the CV-B5 and lncRNAs may reveal the potential targets for antiviral drugs and provide insights for disease affecting the central nervous system.
In this study, RNA sequencing was used to determine the expression of lncRNAs with or without CV-B5 infection in human rhabdomyosarcoma(RD)cell and the human neuroblastoma-derived neuronal cell lineSH-SY5Y cell. And Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were adopted to de ne the potential role of differentially expressed lncRNA in CV-B5-host interactions. We also further conducted the candidate lncRNA-IL12A to elucidate its effect on CV-B5 proliferation. Our results will identify the novel molecular pathways for the treatment of severe HFMD.

Materials And Methods
Cells and CV-B5 virus infection RD and SH-SY5Y cells were cultured in Dulbecco's Modi ed Eagle Medium (DMEM) (SH30243.01, HyClone) supplemented with 10% fetal bovine serum (FBS) (FBSSA500-S, AusGeneX) or 2% FBS (maintenance medium) at 37°C in a humidi ed incubator with 5% CO 2 . The CV-B5 virus was stored at -80℃ in our laboratory. When RD and SH-SY5Y cells were cultured approximately to 80% con uency and then were infected with CV-B5 virus at multiplicity of infection (MOI) of 1. The cytopathic effect (CPE) in RD and SH-SY5Y cells were evaluated at 24-hour post-infection (hpi) by imaging under an electron microscope. At the same time, virus infection was veri ed using the CV-B5 speci c antibody (preserved in our laboratory).

RNA isolation
Total RNA was extracted from virus uninfected and infected RD and SH-SY5Y cells using RNAiso Plus (No.9109, TaKaRa) according to the manufacturer's instructions. RNA quality was evaluated on 1% agarose gels, followed by determination of the concentration and purity by measuring the absorption at 260 nm using a Nano Photometer spectrophotometer. Total RNA obtained from the samples was used for the subsequent steps.

RNA library construction and sequencing
Ribosomal RNA was rst eliminated, and rst strand cDNA and second strand cDNA were synthesized using random hexamer primers [23]. It is worth noting that the dTTP in dNTPs is replaced by dUTP in the synthesis of the second chain, followed by end repair, poly-A 3' ends, ligation of adapters, and screening of DNA fragments in order to construct the RNA-seq library [24]. Library quality was assessed using Qubit2.0, Agilent 2100, and RT-PCR [25], while next-generation high throughput sequencing of samples were performed on the Illumina PE 150 [26,27]. The RNA library was sequenced by Novogene Co., LTD (Beijing, China).

Assessment and classi cation of new lncRNAs
High quality clean data was assembled, and then the novel transcripts were ltered in order to remain with transcripts whose length was greater than 200 bp and had at least two exons. CPC2, CNCI, and PFAM tools were then used to predict the coding ability of the transcripts. Previous studies have reported that an intersection of the tools indicates novel lncRNA [28,29]. Then we calculated expression levels of the transcripts using fragments per kilobase of transcript per million fragments mapped (FPKM) values [30]. Differentially-expressed genes were de ned as p<0.05 and/or fold change>2 times based on their FPKM values between the groups, which were calculated using Cuffdiff soft-ware [31]. Differentiallyexpressed lncRNAs are divided into four classes: antisense, intergenic, sense overlapping, and sense intronic lncRNAs according to their position relative to protein-coding genes [32].

GO annotations and KEGG enrichment
Gene ontology (GO) is a comprehensive database for describing molecular functions, biological processes, and cell composition because it provides the functional annotation and classi cation of molecular functions [33]. On the other hand, the Kyoto Encyclopedia of Genes and Genomes (KEGG) is a resource database that integrates information on genome, chemistry, and system functions [34]. We performed GO and KEGG enrichment analysis on target genes in order to determine the co-expression of differentially-expressed genes.
RT-qPCR con rmation RNA was veri ed using a One Step TB Green PrimeScript RT-PCR kit (RR086A, Takara). Reactions were performed using an Applied Biosystems 7500 Real Time PCR System (Applied Biosystems), with the following cycle conditions: 42℃ for 5 min, 95℃ for 10 s, 40 cycles at 95℃ for 5 s, and 60℃ for 34 s. GAPDH from Homo sapiens was chosen as the endogenous reference gene. All the primers used in the study are listed in Table 1. A signi cantly up-regulated lncRNA-IL12A was only obtained in CV-B5 infected SH-SY5Y cell. Using RNAfold predicted the structure of lncRNA-IL12A, and veri ed its in uence on virus replication.

Results
Validation of the CV-B5 infection models RD and SH-SY5Y cells infected with 1MOI CV-B5 at 24h can cause CPE, including rounding up, aggregation, and even death when compared with control cell (Fig. 1A). Western Bolt analysis also detected the viral speci c protein VP1 in both RD and SH-SY5Y cells (Fig. 1B). Three independent experiments were carried out. These results showed that RD and SH-SH5Y cells were successfully permissive to CV-B5 infection.

Analysis and classi cation of differentially-expression of lncRNAs
To analyze the level of transcripts, RD and SH-SY5Y cells infected with CV-B5 were subjected to Human Genome Array. We obtained 6500 and 5375 novel lncRNAs in RD and SH-SY5Y cells respectively, based on CPC2, CNCI and PFAM analysis ( Fig. 2A and 2B). With regard to the differentially-expressed genes, our transcript contained a total of 1268 lncRNAs in RD cell (Supplementary table 1) and 1603 lncRNAs in SH-SY5Y cell (Supplementary table 2). In total, 1264 lncRNAs with transcripts length ≥200 bp and multi-exon in nature were acquired from the primary lncRNAs in RD cell, of which, 46.2% were lincRNAs, 28.6% were anti-sense lncRNAs, 24.1% were sense overlapping lncRNAs, and 1.0% were sense intronic lncRNAs (Fig.   2C). On the other hand, 1598 lncRNAs were ltered from the primary lncRNAs in SH-SY5Y cell and they contained 48.6% lincRNAs, 34.7% anti-sense lncRNAs, 16.0% sense overlapping lncRNAs, and 0.8% sense intronic lncRNAs (Fig. 2D).
Using Cuffdiff and p<0.05 as the criteria, the differentially-expressed lncRNAs in RD and SH-SY5Y cells in response to CV-B5 infection were shown in clustering analysis (Fig. 3A). Preliminary microarray hybridization in RD cell identi ed 508 up-regulated and 760 down-regulated lncRNAs which greatly changed after CVB5 infection (Fig. 3B). And, 792 signi cantly increased and 811 signi cantly decreased lncRNAs were identi ed in SH-SY5Y cell infected with CV-B5 (Fig. 3C). Among them, some have been associated with neurological disorders, such as TCONS_00385873, TCONS_00087939 and TCONS_00450969, when compared to RD cell. Therefore, we conclude that these lncRNAs may have critical roles in the nervous system.

GO annotation and KEGG pathway analysis
We used the clusterPro ler to explore the gene ontology functional classes and pathways in order to further understand our data. The differentially-expression lncRNAs were used to conduct the GO enrichment and KEGG pathway analysis. The top 20 signi cant GO terms of molecular function, biological process, and cellular component, based on the p value and rich factors, were shown in Fig. 4. The results indicated that the most enriched GO terms in RD cell were single-organism process and membrane part, and cellular process (Fig. 4A). And cell part and protein binding were the most enriched GO terms in SH-SY5Y cell (Fig. 4B).
With regard to the KEGG pathway, the most enriched terms were mainly neuroactive ligand-receptor interaction, spliceosome, and viral myocarditis for CV-B5 infection in RD cell (Fig. 5A). On the other hand, the enriches terms in SH-SY5Y cell were mainly associated with several signaling pathways including ECM-receptor interaction, complement and coagulation cascades and TGF-beta signaling pathway (Fig.  5B).

RT-qPCR for validating the lncRNAs after CV-B5 infection
To further verify the accuracy of the RNA sequencing analysis, the signi cantly up-and down-regulated lncRNAs were selected and validated using RT-qPCR. Seven lncRNAs were selected for further validation.
Fold changes in the CV-B5 infected cells were calculated using the 2 -△△CT method (three independent experiments). Results showed that three lncRNAs were up-regulated and four lncRNAs were downregulated (Fig. 6). These ndings were in accordance with our data analysis, which indicated that fold change results determined by RT-qPCR were lower than those determined using RNA sequencing results.
The analysis of antisense lncRNA-IL12A The function of lncRNAs needs further experimentation because it cannot only be inferred from the RNA sequencing results. LncRNA-IL12A belonged to antisense lncRNA which related the gene Interleukin-12 subunit alpha and located on chromosome 3. In our present study, we selected lncRNA-IL12A in order to illustrate its roles in CV-B5 replication. Firstly, we predicted the secondary structure of lncRNA-IL12A (Fig.  7A). Then we contrasted pcDNA3.1-2Flag-lncRNA-IL12A transfect into 293T and found it does not have the ability to code (Fig. 7B). Next, we overexpressed and silenced its RNA expression in order to determine CV-B5 VP1 expression and our results showed that LncRNA-IL12A inhibits viral replication (Fig. 7CD). Finally, from the KEGG analysis we found the lncRNA-IL12A enriched in the Wnt signaling pathway and the key molecules downstream (C-myc and CyclinD1) siginicantly changes by real-time RT-PCR (Fig. 7E).

Discussion
In our research, we used the next-generation sequencing to investigate the non-coding transcriptome of CV-B5 infected RD and SH-SY5Ycells. The analyses revealed signi cant differences in the patterns of lncRNAs expression in different cells. In additional, we performed GO and KEGG analysis to determine the functions of these lncRNAs. Moreover, we explored the regulatory mechanism of the candidate lncRNA-IL12A on virus replication and found that it inhibited viral replication through Wnt signaling pathway. To the best of our knowledge, this study is the rst comprehensive research that has clearly illustrated the role of lncRNA in the host response to CV-B5 infection.
Since the development of a vaccine for EV-71, CV-B5 has been considered as the emerging neurotropic enteroviral pathogen and it has led to an increasing in the prevalence of severe HFMD [9][10][11]. Therefore it is vital to use a neuron-speci c human cell in order to understand the interactions between CV-B5 and host. In this study, RD and SH-SY5Y cells were selected as experimental models. RD cell (derived from rhabdomyosarcoma) are highly susceptible to CV-B5 and often used to isolate CV-B5 from clinical specimens. Neuroblastoma cell line SH-SY5Y is an excellent model system for neurotropic viruses and is useful for neuropathogenicity studies and it has been previously used Thus, we used the SH-SY5Y cell model to determine whether a barrier to viral entry was partly responsible for neurotropism. In this study, we successfully modeled the virus infection (Fig. 1), infect the RD and SH-SY5Y cells with 1 MOI CV-B5 at 24h, we observed the CPE and also measured the CV-B5 speci c protein.
Unlike mRNAs, lncRNAs have some unique and important characteristics. The RNA sequencing conducted identi ed 1268 differentially expressed lncRNAs in RD cell and 1603 differentially expressed lncRNAs in SH-SY5Y cell after CV-B5 infection. From the Heat map analysis, CV-B5 infected the different cells show that lncRNA expression pro les are signi cantly different (Fig. 3). We then used GO and KEGG pathway analysis to identify the differentially expressed genes. GO analysis that they were mainly enriched single-organism process and membrane part in RD cell and cellular process, cell part and protein binding in SH-SY5Y cell (Fig. 4). And KEGG results indicated that contractile ber, myo bril, and chromatin assembly genes were enriched upon CVB5 infection in RD cell (Fig. 5A). On the other hand, embryonic morphogenesis, ossi cation, and membrane raft genes were enriched upon CV-B5 infection in SH-SY5Y cell, and the coding gene was enriched in the neurotrophin and ErbB signaling pathway (Fig. 5B). Previous studies have reported that the neurotrophin signaling pathway can regulate growth, differentiation, and apoptosis of neurons through p75 neurotrophin receptor and tyrosine kinase receptor [35; 36]. ErbB signaling pathway modulates neural development, cell proliferation, and differentitation and is also associated with neurodevelopment disorders [37]. The above-mentioned results suggest that the CV-B5 can lead to neurogenic diseases.
Similar to CV-B5, EV-71 and CV-A16 also belong to the human enterovirus genus within the Picornaviridae family and they are the main pathogens that cause HFMD and even severe neurological diseases [1].
Previous study had showed the lncRNA expression pro le of CV-A16 and EV-71 infecting RD cell [38; 39].
Our results indicated that CV-B5 infected RD cell had the lowest number of total lncRNAs, 1268, while the EV-71 and CV-A16 infected RD cell had 3231 and 1970 lncRNAs, respectively. Among the total lncRNAs, intergenic lncRNAs accounted for about more than half of EV-71, CV-A16, and CV-B5 infection. Several reports have de ned intergenic lncRNAs as a transcription brake that both promotes and represses the expression of in ammatory genes [40; 41]. Therefore, we speculate that EV-71, CV-A16, and CV-B5 having 50% of intergenic lncRNAs may regulate in ammation. However, antisense lncRNAs accounted for 28.6% of CV-B5 infection compared with 12.7% in EV-71 and 1.54% in CV-A16 infection. Antisense lncRNAs can enhance combination of the promoter and enhancer responsible for regulating gene expression [42]. In our study, lncRNA-IL-12A belongs to antisense lncRNAs, but further veri cation should be conducted to determine whether it in uences the promoter level.  [44]. Cao et al reported that lncRNA-ACOD1 can regulate recombinant GOT2 protein and its metabolites in uence VSV replication [45]. In our study, we have also concluded that lncRNA-IL-12A can regulate CV-B5 infection by Wnt signaling pathway, but further veri cation on protein level should be conducted. Despite the progresses in the interaction between virus and lncRNAs, the speci c functions of these lncRNAs have not been fully elucidated. Therefore, the virus infection models used in the study will provide a new platform for studying the mechanism and regulation of lncRNAs. Compliance with ethical standards Con ict of interest None con ict of interest.    The KEGG pathway analysis of differentially-expressed genes after CV-B5 infection. (A/B) The top 20 enriched pathways associated with differentially-expressed lncRNAs in RD and SH-SY5Y (5Y) cells.

Figure 6
The expression of seven differentially expressed lncRNAs was validated using qRT-PCR and RNA-seq