- Short report
- Open Access
Transcriptional profiling of host cell responses to encephalomyocarditis virus (EMCV)
Virology Journal volume 14, Article number: 45 (2017)
Encephalomyocarditis virus (EMCV) has been discovered on pig farms worldwide and can cause myocarditis in piglets and reproductive failure in sows. However, little is known about the host transcriptional responses to infection and host-pathogen interactions.
In this study, transcription profiling was performed by Illumina RNA-Sequencing (RNA-seq) to identify EMCV induced differentially expressed genes in BHK-21 cells at serial time points (12, 24, and 30 h post infection (hpi)), using mock infected cells as control.
We identified 237, 241, and 207 differentially expressed genes (DEGs) respectively, majority of which were up-regulated. A large number of DEGs clustered into host defense, cellular signaling and metabolism categories. Moreover, short time series expression analysis revealed that 12 hpi was an important time point for expression change, indicating host virus resistance.
This RNA-seq analysis provides the first data for understanding the network of virus host interactions under EMCV infection in vitro, and for identifying host components which involved in the virus infection course.
Encephalomyocarditis virus (EMCV) belongs to the Cardiovirus genus of the Picornaviridae family, and is a non-enveloped, positive single stranded RNA virus . EMCV has been regarded as a worldwide distributed zoonotic pathogen infecting a broad range of host species including rodents, non-human primates, swine, as well as human beings [2–4]. Since the first outbreak in 1958, clinical outbreaks due to EMCV in swine farms have been reported in several countries [3, 5, 6]. The number of EMCV strains isolated from pigs and wild animals in China is increasing in recent years [7–9]. Although pigs are the most severely infected domestic animal species to EMCV, rodents are considered to be the natural reservoir of the virus and play a vital role in transmitting the virus to pigs . EMCV induces acute myocarditis with sudden death in preweaning piglets as well as severe reproductive disorder in pregnant sows [6, 10]. Moreover, EMCV increases the risks of xenotransplantation for human recipients . In rodents, EMCV induces encephalitis, member paralysis, myocarditis and type 1 diabetes [12–14]. EMCV infection was also responsible for the deaths of various animals in zoos distributed in different countries [2, 15, 16].
As the genomic RNA of EMCV is in the cytoplasm, it becomes translated into viral proteins which are indispensable for replication and viral particle formation. EMCV manipulates the host lipid metabolism and induces cell membrane rearrangement to generate replication organelles needed for replication . From virus entry into host cells to release of mature viral particles, viral factors must interact with multiple cellular components of host cells and affect host transcription. However, the mechanism of interaction between host and virus is largely unknown, e.g. cellular signal transduction and host cell pathological responses triggered by EMCV infection. The initial and subsequent host responses to the virus are significant to understand host pathogen interactions and can be investigated at the transcriptional level. Next generation RNA sequencing (RNA-seq) analysis makes the study of impact by virus infection on host cell transcription at whole genome scale efficient. Previous RNA-seq transcriptomic analysis supply massive data on host transcriptional responses to various viruses. Reinhard et al. studied transcriptional host responses to feline immunodeficiency virus and discussed the similarities of that to HIV . Moreover, RNA-seq analysis revealed that HCV infection affected multiple metabolic pathways in host hepatocytes . Transcriptome analysis of host and virus revealed that varicella zoster virus (VSV) altered the differentiation of human keratinocyte, providing insight into the pathogenic mechanisms of VSV . Melina et al. discovered that negative NF-kappa B regulators were significantly up regulated during BVDV infection by RNA-seq, suggesting a possible blocking of this signaling pathway by the virus .
In the present study, we investigated the host responses induced by EMCV at different time points (12, 24 and 30 hpi) using RNA-seq, and analyzed differential gene expression between uninfected and infected cell groups, as well as the dynamic of host gene expression during the virus infection. The results provide the first data for understanding of the complex mechanisms of virus host interactions during EMCV infection.
Materials and methods
Cells and Virus
Baby hamster kidney 21 (BHK-21) cells were maintained in Dulbecco's modified Eagle's medium (DMEM, Gibco) supplemented with 10% (v/v) fetal bovine serum (FBS, HyClone) in a humidified 5% CO2 incubator at 37 °C. EMCV PV21 strain (GenBank No. X74312) from ATCC was propagated in BHK-21 cells for viral challenge.
BHK-21 cells were cultured in 75 cm2 flasks until they reached 80% confluency. Cells were washed with 1% phosphate buffered saline and infected with EMCV at a multiplicity of infection (MOI) of 0.01. After incubating for 1 h at 37 °C, the supernatant was removed and DMEM medium supplemented with 3% FBS was added into each flask and incubated at 37 °C. Control cells were mock infected with FBS free DMEM and treated in the same way in parallel, and then were harvested (designated as 0 hpi in dynamic gene expression analysis during EMCV infection). Subsequently, cells were harvested at 12, 24 and 30 h after EMCV infection for transcriptomic analysis. Three individual replicates were set up for each time point.
Library construction and illumina sequencing
Total RNA was extracted from uninfected and infected BHK-21 cells using TRIzol Reagent (Life technologies, USA) according to the manufacturer's instructions. The quality of RNA was checked by capillary electrophoresis separation on Agilent 2100 Bioanalyzer (Agilent, USA). RNA samples were stored at −80 °C until further use. 1 μg of total RNA (three replicates of each time point were pooled together) were used to generate cDNA libraries. mRNA was enriched by removing rRNAs from the total RNA after DNase I treatment, by which long non-coding RNAs (used in further study) were not depleted. mRNA was chemically fragmented into 140–160 bp fragments and used as templates to synthesize cDNA by priming with random hexamers. The synthesized cDNA fragments were purified and resolved with elution buffer for end reparation as well as adding single nucleotide A (adenine). The cDNA fragments were subsequently connected with sequencing adapters. Suitable fragments were selected by agarose gel electrophoresis as templates of RCR amplification. The quantity and quality of the cDNA libraries were assessed by Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System. Finally, the prepared libraries were sequenced using Illumina HiSeq™ 4000. The RNA-Seq raw data files have been uploaded in NCBI (http://www.ncbi.nlm.nih.gov/bioproject/356035) under the accession number PRJNA356035.
RNA-seq data analysis
RNA-seq raw data was filtered using FASTQC software to exclude reads containing adapters, reads in which unknown bases were more than 10% and low quality reads. Clean reads were aligned and mapped against hamster reference genome (http://www.ncbi.nlm.nih.gov/genome/11998) and EMCV PV21 genome (GenBank: X74312.1) by using BWA. Bowtie was used to map reads to hamster gene reference. To analyze the expression abundance of genes, unique mapped gene reads were normalized by Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) and log transformed using Cufflinks. DEGs were identified using Cuffdiff. Genes with a fold change of larger than 1.5 and with adjusted P values < 0.05 were considered as differentially expressed. The identified DEGs were subjected to further analysis, including Gene Ontology (GO) analysis, pathway enrichment based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and gene expression dynamic analysis conducted by STEM .
Real-time reverse transcription quantitative PCR (RT-qPCR) validation
RT-qPCR was carried out to verify the gene expression by RNA-seq analysis. 1 microgram total RNA from each group (mock infected cells, 12 hpi, 24 hpi, and 30 hpi) was used for reverse transcription using a PrimeScript RT reagent kit with a gDNA Eraser (TaKaRa). Equal amounts of RNA from three replicates were pooled together to subject to reverse transcription. All RT-qPCRs were performed in 10 μL reactions including 1 μL of cDNA product and 5 μL of SYBR Premix Ex TaqII (TaKaRa). The reactions were conducted in triplicate using a CFX96 real-time PCR machine (Bio-Rad, USA) following the temperature protocol: Initial denaturation 95 °C for 2 min; 39 cycles of 95 °C for 10 s, 56 °C for 20 s, 72 °C for 20 s. The expression fold-changes were calculated by the relative 2-△△Ct method and the housekeeping gene ACTB was taken as a reference gene.
Results and discussion
Transcriptome sequencing and differential expression analysis
To depict the global picture of host transcriptomic response to EMCV infection, and to identify host factors involved in the infection course, we performed RNA-seq on Illumina platform using cDNA libraries of EMCV infected and mock infected BHK-21 cells. EMCV infected cells were harvested at 12, 24, and 30 hpi. An average of 115 million reads were generated from each sample, and an average of 114 million reads passed the quality control. 72.6% of clean reads were mapped to the hamster genome in the mock infected group. In 24 hpi group, 65.3% and 7.2% of clean reads were aligned to hamster genome and EMCV genome respectively, indicating viral RNA replication before 24 hpi. In 30 hpi group, 41.5% of clean reads were mapped to EMCV genome, indicating active replication of the viral RNA (Table 1). The amounts of EMCV copies were quantified by real-time PCR, and EMCV infection was also validated. The viral RNA levels increased rapidly during infection (Fig. 1).
The number of DEGs (≥1.5fold, with P value < 0.05) at each time point was shown in Table 2. 237, 241, and 207 genes were found differentially expressed at 12, 24, and 30 hpi respectively (Additional file 1: Figure S1, Table 2). There was little difference between the total number of DEGs at each time point, however, the percentage of down regulated genes reduced drastically as the infection time grew. At 30 hpi, only 7.7% DEGs were down regulated. Even at 12 hpi, majority of the DEGs were up regulated. Moreover, DEGs with a fold change of larger than 2 at 30 hpi were much more than that at other time points (12 hpi and 24 hpi) (Table 2).
DEG data were further subjected to Gene Ontology (GO) and KEGG pathway analysis. GO terms classified the functions of all differentially expressed transcripts, producing 41 functional types at 30 hpi, including immune system process, metabolic process, response to stimulus, and signaling (Fig. 2). We designated three functional categories of interest based on GO functional classification and KEGG pathway enrichment (Figs. 2 and 3): host defense, signaling and metabolism. Representative DEGs involved in host defense at the three time points, as well as the KEGG pathways which the DEGs belong to, were listed in Table 3. We identified several innate immunity pathways responsive to EMCV infection, including RIG-I-like receptor signaling pathway, TGF-beta signaling pathway, and NF-kappa B signaling pathway (Fig. 3, Table 3).
LGP2, which is demonstrated to coordinate with MDA5 , showed a fold decrease at 12 h after EMCV infection. Genomic RNA of EMCV is considered to recognized by cellular RLR receptor MDA5, not RIG-I . Mx dynamin-like GTPases have been shown to have antiviral activity in type I and type III interferon systems. Recent data revealed that human MX2 served as a restriction factor for HIV-1 and other primate lentiviruses . Here we showed that MX2 was down-regulated in EMCV infected BHK-21 cells at 12 hpi. The above results indicated a quick restriction of host antiviral activity by EMCV. On the other hand, antiviral genes, TAK1 and complement C3, were found up-regulated in this study at 24 hpi and 30 hpi, suggesting host defense responses. Tumor necrosis factor receptor superfamily, member 25 (TNFRSF25) transcript was discovered decreased remarkably at 30 hpi. The same TNFRSF25 was also found down-regulated in Japanese encephalitis virus infected mouse spleens . Recent data indicated that knockdown of TNFRSF25 led to an increased anti-apoptotic potential . EMCV is a lytic virus causing necrotic cell death, which counteracts host cell antiviral immunity by inhibiting apoptosis. Down-regulation of TNFRSF25 might contribute to this inhibition. TNFAIP3 (A20) transcript expression was found increased at 24 hpi and 30 hpi, which functions in negative regulation of NF-kappa B signaling (Table 3). A20 was demonstrated to inhibit TNF induced apoptosis . Furthermore, a recent study showed that deficiency of A20 in myeloid cells improved the resistance against influenza infection, by an enhanced innate response . Thus, overexpression of A20 probably counteracted the innate immune responses induced by EMCV infection, however, still needs further investigation. These findings indicated that EMCV also restricted host defense at later infection time points.
Representative DEGs involved in signaling and the KEGG pathways which the DEGs belong to, were listed in Additional file 2: Table S1. Representative DEGs involved in metabolism as well as the KEGG pathways which the DEGs belong to, were listed in Additional file 3: Table S2. Several signaling pathways, including MAPK signaling pathway, PPAR signaling pathway, mTOR signaling pathway, were widely affected subsequent to EMCV infection (Additional file 2: Table S1, Fig. 3). The differential analysis also showed that multiple metabolism pathways, especially lipid metabolism were responsive to EMCV, including Terpenoid backbone biosynthesis, Steroid biosynthesis, Biosynthesis of unsaturated fatty acids, and Glycerophospholipid metabolism (Additional file 3: Table S2, Fig. 3). Replication of EMCV RNA needs to modulate host cell lipid landscape . Following the deep analysis of DEGs, we saw that the transcriptional level of host lipid metabolism were changed by EMCV.
Dynamic host gene expression during EMCV infection
To analyze dynamic gene expression at series time points (0, 12, 24, and 30 hpi) during EMCV infection with RNA-seq DEG data, we performed pattern analysis using STEM (Short Time-series Expression Miner) software (mock infection was designated as 0 hpi). STEM is used to analyze short time series gene expression data. The results revealed four significant expression patterns (P < 0.05): 1. decreasing at 12 hpi and then increasing (red, including IL18, LGP2, A20 and TAK1); 2. increasing at 12 hpi and then decreasing (green, including WNT5A); 3. gradually increasing, part of the genes decreasing at 12 hpi, and then increasing gradually at later time points after infection (blue, including PLA2G2A, TXNIP, C3, TRAF2, PTGS2, CLK1, and ZFP36); 4. decreasing at 12 hpi and then increasing, but back to expression level at 0 hpi (yellow, including MRGPRG) (Additional file 4: Table S3, Fig. 4). The STEM analysis was also supported by real-time PCR later. Several antiviral genes, including LGP2, TAK1, ZFP36 and C3 decreased at 12 hpi, then increased, of which C3 increased gradually at later time after 12 hpi. These results suggested that 12 hpi was an important time point for expression pattern change, probably indicating antagonism between virus and host. EMCV might immediately reduce host antiviral gene expression to restrict antiviral activity, however, the host cells were able to make a response later to counteract with or defeat the restriction activity induced by EMCV, although were not capable of inhibiting the replication of the virus.
Real-time PCR verification of differential expressions
Eighteen genes from the three functional categories were selected for validation of our RNA-seq DEG data by Real-time RT-qPCR (Additional file 5: Table S4). Gene expression profiles of Real-time RT-qPCR were compared with those of RNA-seq. Among the 18 genes, 16 genes exhibited similar expression patterns as compared with RNA-seq data. The other two genes (DDX3X and HES2) didn't have obvious expression changes based on RT-qPCR results. Fourteen of the verified differential expression genes were up regulated, and only two genes, TNFRSF25 and WNT9A down regulated. ZFP36, a RNA binding protein, which targets HIV RNA for degradation , showed the most significant changes in differential expression (Fig. 5). Nevertheless, the specific role of ZFP36 in EMCV infection is unclear. Moreover, short time series gene expression pattern was also verified by real-time PCR (Additional file 6: Figure S2).
Differentially expressed genes
Hours post infection
Kyoto encyclopedia of genes and genomes
Reverse transcription quantitative PCR
Carocci M, Bakkali-Kassimi L. The encephalomyocarditis virus. Virulence. 2012;3(4):351–67.
Canelli E, Luppi A, Lavazza A, Lelli D, Sozzi E, Martin AM, Gelmetti D, Pascotto E, Sandri C, Magnone W, et al. Encephalomyocarditis virus infection in an Italian zoo. Virol J. 2010;7:64.
Murnane TG, Craighead JE, Mondragon H, Shelokov A. Fatal disease of swine due to encephalomyocarditis virus. Science. 1960;131(3399):498–9.
Feng R, Wei J, Zhang H, Fan J, Li X, Wang D, Xie J, Qiao Z, Li M, Bai J, et al. National serosurvey of encephalomyocarditis virus in healthy people and pigs in China. Arch Virol. 2015;160(12):2957–64.
Dea S, Bilodeau R, Sauvageau R, Martineau GP. Outbreaks in Quebec pig farms of respiratory and reproductive problems associated with encephalomyocarditis virus. J Vet Diagn Invest. 1991;3(4):275–82.
Koenen F, Vanderhallen H, Castryck F, Miry C. Epidemiologic, pathogenic and molecular analysis of recent encephalomyocarditis outbreaks in Belgium. Zentralbl Veterinarmed B. 1999;46(4):217–31.
Lin W, Liu Y, Cui S, Liu H. Isolation, molecular characterization, and phylogenetic analysis of porcine encephalomyocarditis virus strain HB10 in China. Infect Genet Evol. 2012;12(6):1324–7.
Feng R, Zhang H, Wei J, Li X, Xie J, Li M, Qiao Z, Feng Y, Ma Z. Isolation, molecular and phylogenetic analysis of encephalomyocarditis virus strain GS01 in China. Infect Genet Evol. 2015;30:19–26.
Liu H, He X, Song X, Xu L, Zhang Y, Zhou G, Zhu W, Chang C, Yin Z, Shi Y, et al. Isolation and molecular and phylogenetic analyses of encephalomyocarditis virus from wild boar in central China. Infect Genet Evol. 2016;40:67–72.
Joo HS, Kim HS, Leman AD. Detection of antibody to encephalomyocarditis virus in mummified or stillborn pigs. Arch Virol. 1988;100(1–2):131–4.
Brewer LA, Lwamba HC, Murtaugh MP, Palmenberg AC, Brown C, Njenga MK. Porcine encephalomyocarditis virus persists in pig myocardium and infects human myocardial cells. J Virol. 2001;75(23):11621–9.
Psalla D, Psychas V, Spyrou V, Billinis C, Papaioannou N, Vlemmas I. Pathogenesis of experimental encephalomyocarditis: a histopathological, immunohistochemical and virological study in mice. J Comp Pathol. 2006;135(2–3):142–5.
LaRue R, Myers S, Brewer L, Shaw DP, Brown C, Seal BS, Njenga MK. A wild-type porcine encephalomyocarditis virus containing a short poly(C) tract is pathogenic to mice, pigs, and cynomolgus macaques. J Virol. 2003;77(17):9136–46.
White LL, Smith RA. D variant of encephalomyocarditis virus (EMCV-D)-induced diabetes following natural killer cell depletion in diabetes-resistant male C57BL/6 J mice. Viral Immunol. 1990;3(1):67–76.
Yeo DS, Lian JE, Fernandez CJ, Lin YN, Liaw JC, Soh ML, Lim EA, Chan KP, Ng ML, Tan HC, et al. A highly divergent Encephalomyocarditis virus isolated from nonhuman primates in Singapore. Virol J. 2013;10:248.
Lamglait B, Joris A, Romey A, Bakkali-Kassimi L, Lemberger K. Fatal Encephalomyocarditis Virus Infection in an African Savanna Elephant (Loxodonta Africana) in a French Zoo. J Zoo Wildl Med. 2015;46(2):393–6.
Dorobantu CM, Albulescu L, Harak C, Feng Q, van Kampen M, Strating JR, Gorbalenya AE, Lohmann V, van der Schaar HM, van Kuppeveld FJ. Modulation of the Host Lipid Landscape to Promote RNA Virus Replication: The Picornavirus Encephalomyocarditis Virus Converges on the Pathway Used by Hepatitis C Virus. PLoS Pathog. 2015;11(9):e1005185.
Ertl R, Klein D. Transcriptional profiling of the host cell response to feline immunodeficiency virus infection. Virol J. 2014;11:52.
Woodhouse SD, Narayan R, Latham S, Lee S, Antrobus R, Gangadharan B, Luo S, Schroth GP, Klenerman P, Zitzmann N. Transcriptome sequencing, microarray, and proteomic analyses reveal cellular and metabolic impact of hepatitis C virus infection in vitro. Hepatology. 2010;52(2):443–53.
Jones M, Dry IR, Frampton D, Singh M, Kanda RK, Yee MB, Kellam P, Hollinshead M, Kinchington PR, O’Toole EA, et al. RNA-seq analysis of host and viral gene expression highlights interaction between varicella zoster virus and keratinocyte differentiation. PLoS Pathog. 2014;10(1):e1003896.
Villalba M, Fredericksen F, Otth C, Olavarria V. Transcriptomic analysis of responses to cytopathic bovine viral diarrhea virus-1 (BVDV-1) infection in MDBK cells. Mol Immunol. 2016;71:192–202.
Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinforma. 2006;7:191.
Bruns AM, Horvath CM. LGP2 synergy with MDA5 in RLR-mediated RNA recognition and antiviral signaling. Cytokine. 2015;74(2):198–206.
Gitlin L, Barchet W, Gilfillan S, Cella M, Beutler B, Flavell RA, Diamond MS, Colonna M. Essential role of mda-5 in type I IFN responses to polyriboinosinic:polyribocytidylic acid and encephalomyocarditis picornavirus. Proc Natl Acad Sci U S A. 2006;103(22):8459–64.
Haller O, Staeheli P, Schwemmle M, Kochs G. Mx GTPases: dynamin-like antiviral machines of innate immunity. Trends Microbiol. 2015;23(3):154–63.
Yang Y, Ye J, Yang X, Jiang R, Chen H, Cao S. Japanese encephalitis virus infection induces changes of mRNA profile of mouse spleen and brain. Virol J. 2011;8:80.
Xu LX, Grimaldo S, Qi JW, Yang GL, Qin TT, Xiao HY, Xiang R, Xiao Z, Li LY, Zhang ZS. Death receptor 3 mediates TNFSF15- and TNFalpha-induced endothelial cell apoptosis. Int J Biochem Cell Biol. 2014;55:109–18.
Won M, Park KA, Byun HS, Sohn KC, Kim YR, Jeon J, Hong JH, Park J, Seok JH, Kim JM, et al. Novel anti-apoptotic mechanism of A20 through targeting ASK1 to suppress TNF-induced JNK activation. Cell Death Differ. 2010;17(12):1830–41.
Maelfait J, Roose K, Bogaert P, Sze M, Saelens X, Pasparakis M, Carpentier I, van Loo G, Beyaert R. A20 (Tnfaip3) deficiency in myeloid cells protects against influenza A virus infection. PLoS Pathog. 2012;8(3):e1002570.
Maeda M, Sawa H, Tobiume M, Tokunaga K, Hasegawa H, Ichinohe T, Sata T, Moriyama M, Hall WW, Kurata T, et al. Tristetraprolin inhibits HIV-1 production by binding to genomic RNA. Microbes Infect. 2006;8(11):2647–56.
This work was supported by the Fundamental Research Funds for the Central Universities of China (31920150026), Talent Recruitment Program Funds of Northwest Minzu University (xbmuyjrc201317), National Natural Science Foundation of China (N0. 31460665), and Program for Changjiang Scholars and Innovative Research Team in University (IRT13091).
Availability of data and materials
The raw RNA-seq data files reported in this paper were deposited to NCBI (http://www.ncbi.nlm.nih.gov/bioproject/356035) under the accession number PRJNA356035.
Ruofei Feng and Jia Wei conceived and designed this research. Jia Wei, Haixia Zhang, Xiangrong Li, Qiongyi Li, Zhongren Ma, Jialin Bai, and Zilin Qiao performed the experiments and analyzed the data. Jia Wei wrote the manuscript. All authors have read and approved the manuscript.
The authors declare that they have no competing interests.
The Venn diagram shows common differential expressing genes at three time points (12 hpi, 24 hpi, and 30 hpi). (TIF 3137 kb)
Representative DEGs involved in signaling at different time points. (DOCX 15 kb)
Representative DEGs involved in metabolism at different time points. (DOCX 15 kb)
Representative genes in each STEM pattern (4 patterns: red, green, blue, and yellow). (DOCX 15 kb)
Selected genes for real-time PCR verification. (DOCX 13 kb)
Real-time PCR verification of temporal host gene expression regulated by EMCV. Mock infection was designated as 0 hpi. The x-axis represents infection time, while y-axis represents normalized fold change of transcripts. (TIF 1941 kb)