Transcriptome analysis following EV-A71 and CV-A16 infection in respiratory epithelial cells

doi:10.21203/rs.2.20645/v1

PDF

Research

Transcriptome analysis following EV-A71 and CV-A16 infection in respiratory epithelial cells

https://doi.org/10.21203/rs.2.20645/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 13 Jan, 2020

You are reading this latest preprint version

Background: Enterovirus 71 (EV-A71) and coxsackievirus A16 (CV-A16) are the major pathogens responsible for hand, foot and mouth disease (HFMD), but the detailed mechanism caused by the two viruses remains unclear. Methods: In this study, we aimed to adopt transcriptome sequencing technology to investigate changes in the transcriptome profiles after infection with EV-A71 and CV-A16 in human bronchial epithelial (16HBE) cells.

Results: Through systematic bioinformatics analysis, we then searched for useful clues regarding the pathogenesis of HFMD. As a result, a total of 111 common differentially expressed genes were present in both the EV-A71- and CV-A16-infected groups. A trend analysis of these 111 genes showed that there were 91 genes displaying the same trend in the EV-A71 and CV-A16 infection groups, including 49 upregulated genes and 42 downregulated genes. These 91 genes were further used to conduct GO, pathway, and coexpression network analyses. It was discovered that the enriched GO terms (such as Histone acetylation and positive regulation of phosphorylation) and pathways (such as glycosylphosphatidylinositol (GPI)-anchor biosynthesis and DNA replication) might be closely associated with the pathogenic mechanism of the two viruses, and key genes (such as TBCK and GPC) might be involved in the progression of HFMD. Finally, we randomly selected 10 differentially expressed genes for qRT-PCR to validate the transcriptome sequencing data. The experimental qRT-PCR results were roughly in agreement with the results of transcriptome sequencing.

Conclusions: Collectively, our results provide clues for the pathogenesis of HFMD induced by EV-A71 and CV-A16.

Virology

Enterovirus 71 (EV-A71)

Coxsackievirus A16 (CV-A16)

Hand

foot and mouth disease (HFMD)

Transcriptome sequencing

Enterovirus 71 (EV-A71) and Coxsackievirus A16 (CV-A16) have been identified as the most common etiological agents causing hand, foot and mouth disease (HFMD), which is usually seen in Asia-Pacific regions among infants and young children. These two viruses both belong to the human enterovirus A species of the family Picornaviridae, which have positive sense single-stranded RNA with an approximate length of 7400 nucleotides, consisting of four structural viral proteins (VP1 to VP4) and seven nonstructural proteins (2A to 2C, 3A to 3D)[1]. Clinically, most HFMD cases are usually mild and self-limiting, but some cases progress rapidly and are accompanied by severe neurological complications, such as aseptic meningitis, encephalomyelitis, and acute flaccid paralysis, and can further deteriorate into fatal myocarditis and pneumonia[2]. In recent years, the number of these serious cases and deaths has been increasing and has attracted global concern. Fortunately, an inactivated EV-A71 vaccine with completely independent intellectual property rights has been successfully developed by three vaccine organizations, including Beijing Vigoo Biological Co., Ltd. (Vigoo), Sinovac Biotech Co., Ltd. (Sinovac), and the Institute of Medical Biology, Chinese Academy of Medical Science (CAMS). This vaccine is currently available in the market and has been approved by the China Food and Drug Administration (CFDA) [3]. This vaccine shows high efficacy and a satisfactory safety profile to provide protection against clinical EV-A71-associated HFMD, but it does not provide cross-strain protection against CV-A16-induced HFMD. Furthermore, it has no therapeutic effect on patients who have already been infected with EV-A71[4]. In previous studies, it was demonstrated that HFMD patients with EV-A71 and CV-A16 infections presented significantly different clinical manifestations. EV-A71 infections frequently lead to severe central nervous system (CNS) complications and even death, while CV-A16 infection often results in milder symptoms that resolve within a few weeks[5]. Nevertheless, over the past years, accumulating evidence has revealed that HFMD patients infected with CV-A16 can also develop into a serious stage with neurological complications, and their overall condition and clinical symptoms are basically consistent with severe HFMD patients infected with EV-A71[6]. Data from the China National Center for Disease Control and Prevention shows that there were 2.37 million cases of HFMD in China in 2018, including 36 deaths. Therefore, it is urgently needed to strengthen basic research on the replication and pathogenesis of EV-A71 and CV-A16, which could further provide a theoretical basis for the development of HFMD-specific therapeutic drugs.

Today, with the rapid development of sequencing technology, RNA sequencing (RNA-Seq) can be used for investigating differences in gene expression at a genome-wide level. RNA-Seq has the characteristics of more accurate quantification, higher repeatability, a wider detection range and more reliable analysis[7]. In addition to analyzing gene expression levels, RNA-Seq can also discover new transcripts, single nucleotide polymorphisms (SNPs), and splice variants to further provide information about allele-specific gene expression[8]. Indeed, a large number of studies have been conducted to analyze the transcriptome expression profiles of EV-A71 and CV-A16 infections. For example, Lui et al. analyzed the transcriptome profiles of EV-A71-infected colorectal cells and found that EV-A71 activated the signaling pathway, which might participate in inhibiting viral replication[9]. Yao et al. carried out a transcriptome analysis of EV-A71-infected Rhabdomyoma (RD) cells and revealed that EV-A71-2A protein could be considered a key inducer that triggered cellular apoptosis and death in RD cells through mediating thioredoxin-interacting protein (TXNIP)[10]. Xu et al. reported that differentially regulated mRNAs were associated with the host cellular pathways that directed cell cycle/proliferation, apoptosis and cytokine/chemokine and immune responses in SH-SY5Y human neuroblastoma cells infected with EV-A71[11]. This finding suggests that the changed mRNAs might be involved in the pathophysiological mechanisms of EV-A71 infection in human neural cells. Jin et al. performed transcriptome sequencing in CV-A16-infected HEK293T cells and found that CV-A16 can upregulate the expression of SCARB2 and ECM receptor[12]. Song et al. also studied the transcriptome changes in CV-A16-infected rhesus monkey peripheral blood mononuclear cells and discovered that inflammatory cytokines, such as IL-6 and IL-18, were obviously increased after CA16 infection[13]. However, no studies have analyzed and compared the pathologic attributes in EV-A71 and CV-A16 infections. Previous studies have demonstrated that the respiratory tract is the most important route for HFMD transmission; therefore, the interaction between EV-A71/CV-A16 and airway epithelial cells should be investigated[14]. In the current study, we aimed to discover significant differentially expressed genes in EV-A71- and CV-A16-infected respiratory epithelial cells through transcriptome sequencing. We then investigated the potential pathogenesis of HFMD induced by EV-A71 and CV-A16 via systematically analyzing these differentially expressed genes with bioinformatics analysis, which might provide clues to further explore the mechanisms of HFMD and seek molecular targets for HFMD treatment.

Virus and cell lines

Monolayers of human bronchial epithelial (16HBE) cells were purchased from Jennino Biological Technology (Guangzhou, China). The cells were seeded at a density of 5 × 10⁵ cells per well in 6-well sterile plastic culture plates and grown in a base media of Roswell Park Memorial Institute (RPMI)-1640 containing 10% (vol/vol) fetal bovine serum (FBS), 100 units/mL penicillin, 100 μg/mL streptomycin and 2 mM L-glutamine at 37°C in a 95% air and 5% carbon dioxide (CO₂) incubator. For transcriptome study, 16HBE cells were infected at a multiplicity of infection (MOI) of 0.01 with enterovirus 71 (EV-A71; subgenotype C4, GenBank: EU812515.1) or coxsackievirus A16 (CV-A16; subgenotype B, GenBank: JN590244.1), which were isolated from an epidemic in Fuyang, China in 2008 and an HFMD patient in Guangxi, China in 2010, respectively. Both viruses were incubated for one hour to attain virus attachment, and RPMI-1640 containing 1% FBS was then added. The cells were incubated for a further 6 h, 12 h and 24 h for virus propagation. At the end of each incubation point, cells were collected using RPMI-1640 and centrifuged twice in phosphate-buffered saline (PBS) at 4°C for 5 min at 3000 rpm/min. Finally, the harvested cell pellets were snap-frozen in liquid nitrogen and stored until RNA extraction. For the control sample, the same process was carried out with the exception that 2 mL of sterile PBS was used to replace the virus.

RNA extraction and RNA sequencing (RNA-seq) library construction

Total RNA was extracted using an RNeasy^® Mini Kit (Qiagen, USA) as per the methods recommended by the manufacturer. The extracted RNA was cleared of contaminating genomic DNA by RNase-free DNase I (Takara Bio, Japan) treatment. The concentration of RNA was measured with a NanoDrop 2000 (Thermo Scientific, USA), and the RNA quality was determined by an Ultrospec 3000 Pro UV/Visible spectrophotometer (GE Healthcare, UK), where samples with an absorbance ratio (A260/A280) of 1.8 to 2.0 were considered for RNA integrity analysis with an Agilent^® 2100 Bioanalyzer. Samples with an RNA integrity number (RIN) greater than 7 were sent for RNA-seq library construction. Brieﬂy, mRNA was enriched by Oligo (dT) beads, and poly(A)-containing mRNA was then purified using Dynabeads (Life Technologies, USA) and further fragmented into smaller pieces with fragmentation buffer using RNase III and an Ion adaptor. Next, the RNA fragments were reverse-transcribed and amplified to form first-stranded complementary DNA (cDNA) with random primers, followed by second-strand cDNA synthesis. The second-strand cDNA was further purified, adenylated at the 3’ ends, and ligated with sequencing adaptor. These fragments (250~300 bp in size) were subjected to PCR amplification with Phusion High-Fidelity DNA polymerase (NEB, Beijing, China), and the products were sequenced on an Illumina HiSeq™ 2000 platform (Illumina, USA).

Computational analysis of RNA-seq data

Transcriptome assembly

The raw RNA-seq paired-end reads were filtered to remove the “dirty” reads, i.e., those containing sequencing adapters, reads with 10% > Q < 20% bases, and reads of low quality (reads with ambiguous bases “N”) using the Fast QC package (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The obtained “clean” reads were mapped against the human reference genome GRCh38 with HISAT2 software. The fragments per kilobase of transcript sequence per million base pair (FPKM) values of each gene were determined by the length of the gene and read count mapped to the gene with Cufflinks (version 1.3.0)[15].

Principal component analysis (PCA)

PCA was employed for quality control, to identify problems with the experimental design, to find mislabeled samples and to visualize variations between the expression analysis samples by using the data from clean reads.

Identification of differentially expressed genes

For gene expression analysis, the expression level of each gene was calculated by quantifying the number of reads and was further normalized by a variation of the FPKM method. To identify differentially expressed genes between the two samples, Cufflinks was applied to calculate the T-statistic and the p-values for each gene. We calculated the expression ratios of 6 h/0 h, 12 h/0 h and 24 h/0 h as fold changes. All P values were adjusted by Benjamini and Hochberg’s approach to control the false discovery rate (FDR). Genes with a 2-fold change and an adjusted P value < 0.05 were regarded to be differentially expressed.

Unsupervised hierarchical clustering analysis of the common differentially expressed genes

To find commonalities between the EV-A71- and CV-A16-infected 16HBE cells, we constructed a Venn diagram using the obtained common differentially expressed genes in each group. The common differentially expressed genes were subject to unsupervised hierarchical clustering using Euclidean distance and average linkage to measure the cluster similarity/dissimilarity. A trend analysis was further utilized to identify commonalities between EV-A71 and CV-A16 infection. The common differentially expressed genes in the EV-A71 and CV-A16 infections that had the same trend were identified as genes with similar changes.

Functional group analysis

The Database for Annotation, Visualization and Integrated Discovery (http://david.abcc.ncifcrf.gov/), which utilizes Gene Ontology (GO) to check the biological process, molecular function and cellular components of common differentially expressed genes, was applied in the current study. In addition, the Biocarta and Reactome database (http://www.genome.jp/kegg/), which uses the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, was applied to identify pathways of common differentially expressed genes in this study. The FDR-corrected P-value threshold was set at 0.05, which denotes the significance of GO term enrichment and pathway correlation.

Coexpression network construction

The coexpression network construction was based on the GenMANIA algorithm by using the common differentially expressed genes. The construction of coexpression networks is conducive to finding potential mechanisms associated with differentially expressed genes.

Validation of differentially expressed genes by quantitative RT-PCR (qRT-PCR) analysis

To validate the transcriptome results, we designed 10 pairs of primers using Primer Premier 5.0 to perform qRT-PCR analysis targeting 7 upregulated (BLM, GBP3, HDAC9, NNMT, PNISR, RNF6 and TAF1) and 3 downregulated genes (ANO1, IRX2 and PCDH7). Briefly, reverse transcription was first performed using total RNA (1 μg) isolated as described above with the Prime Script^® RT reagent kit (Takara, Japan) according to the manufacturer’s protocol. Next, qRT-PCR was carried out using 1 μl of cDNA template, 0.2 μl of gene-specific primers (Supplementary Table S1), 5 μl of SYBR^®Premix Ex Taq™ (2×), 0.2 μl of ROX Reference Dye II (50×), and water up to 10 μl. The qRT-PCR reaction program was as follows: predenaturation at 95 ℃ for 30 s, followed by 45 cycles of denaturation at 95 ℃ for 10 s, annealing at 55 ℃ for 15 s, and extension at 72 ℃ for 30 s on a 7500 Fast Realtime PCR system (Applied Biosystems, USA). Relative expression levels of the chosen genes were determined by normalizing the data for the target transcripts against the glyceraldehyde-3-phosphate dehydrogenase (GAPDH) transcript data according to the 2⁻^△△^Ctmethod. Three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed.

Overview of the RNA-Seq data from 16HBE cells in response to EV-A71 and CV-A16 infection

Using the Illumina HiSeq™ 2000 platform, sequencing of the seven cellular samples generated approximately 233.68 million raw reads, from which 209.58 million clean reads were obtained, with an average of 29.94 million clean reads per sample. Through rRNA trimming and alignment to the human reference genome GRCh38, 198.68 million reads (94.51%) were mapped to the genome. The detailed data of raw reads, clean reads, rRNA-trimmed reads and mapping reads in each group are shown in Table 1. The mapping reads were used for PCA to assess the discrete degree analysis of each group. The results showed a clear separation of groups between EV-A71 and CV-A16 infection (Fig. 1A). The virus-infected groups were significantly removed from the control group. Thus, these data showed differences between the infection groups.

Table 1

Overview of sequencing data information.
Group	Raw reads	Clean reads	rRNA trimed	Mapping reads
0 h	32,355,870	29,757,717	29,706,403	28,803,947
EV71-6 h	31,914,208	28,895,872	28,833,270	27,439,676
EV71-12 h	41,780,160	37,000,177	36,918,228	34,682,529
EV71-24 h	30,695,202	27,546,778	27,484,897	26,206,126
CA16-6 h	34,101,512	30,206,999	30,142,248	26,206,126
CA16-12 h	31,542,638	28,262,661	28,191,333	28,713,177
CA16-24 h	31,296,916	27,918,538	27,828,716	26,636,815

Expression profile of differentially expressed transcripts in EV-A71- and CV-A16-infected 16HBE cells

To investigate the transcriptome responses to EV-A71 and CV-A16 infection in 16HBE cells, the differentially expressed transcripts based on the criterion of a 2-fold change and an adjusted P-value <0.05” were screened. The distribution of gene expression levels is shown in Fig. S1 using the data of log2 (FPKM). A total of 19,425 differentially expressed genes were found to be significantly differentially expressed, with 8,903 upregulated and 10,522 downregulated differentially expressed genes. The detailed numbers of upregulated and downregulated genes in each group are listed in Fig. 1B. A Venn diagram was constructed to find the common differentially expressed genes in all groups, and it showed that 111 differentially expressed genes cooccurred in each group (Fig. 1C). These 111 common differentially expressed genes were further applied to unsupervised hierarchical clustering analysis. As illustrated in Fig. 1D, by clustering the 111 differentially expressed genes detected in all infected samples, the EV-A71 and CV-A16 infections could be perfectly separated.

Genes analysis of the same trends using GO, Pathways and a Coexpression network

To further explore commonalities between the EV-A71- and CV-A16-infected 16HBE cells, trend analysis of the 111 common differentially expressed genes was carried out to reveal the differentially expressed genes that had the same change trends. The results showed that 49 upregulated and 52 downregulated genes presented a similar trend in both the EV-A71- and CV-A16-infected groups (Fig. 2). To investigate the biological processes that contribute to changes in the transcripts during the development of EV-A71 and CV-A16 infection, these upregulated and downregulated genes were separately utilized to analyze the GO, related pathways and coexpression network. Regarding the GO BP terms, the upregulated genes were enriched in 9 terms (Fig. 3A) and the downregulated genes were enriched in 5 terms (Fig. 3D). The upregulated genes were remarkably enriched in 5 MF terms (Fig. 3B), and the downregulated genes were markedly enriched in 3 MF terms (Fig. 3E). Regarding the CC terms, 5 terms were enriched in the upregulated genes (Fig. 3C) and 5 CC terms were enriched in the downregulated genes (Fig. 3F). Additionally, KEGG pathway enrichment analysis for the significantly differentially expressed genes was used to understand the pathways related to the genes. Our data showed that only 2 pathways were associated with the upregulated genes (Fig. 4A) and 1 was related to the downregulated genes (Fig. 4B). Ultimately, the construction of a coexpression network was implemented to further explain the molecular interactions (Fig. 5).

Validation of transcriptome sequencing results by qRT-PCR

Ten differentially expressed genes from the RNA-seq data were randomly chosen for validation by qRT-PCR (Fig. 6). The results demonstrated that the expression of BLM, GBP3, HDAC9, NNMT, PNISR, RNF6 and TAF1 was upregulated, whereas the expression of ANO1, IRX2 and PCDH7 was downregulated. These findings were consistent with the RNA-seq expression profiles.

Over the last 20 years, EV-A71 and CV-A16 have become important emerging viruses that pose a threat to global public health, especially in children under five years old. They both cause HFMD outbreaks with many CNS-complicated cases and deaths in different parts of the world, particularly in the Asia-Pacific region[16]. Vaccination is considered to be one of the most effective ways to protect against EV-A71 and CV-A16 infection. Therefore, researchers have been working on vaccine development in recent decades, and there are currently three vaccine organizations in China that have developed an inactivated whole EV-A71 vaccine that is safe and has good efficacy for protecting against EV-A71-associated HFMD in children[17]. However, there are no vaccines against HFMD caused by other enteroviruses, including CV-A16. Hence, it is urgent to explore the pathogenesis of HFMD triggered by EV-A71 and CV-A16 infection Transcriptome sequencing technology is able to identify all transcripts, even when lacking detailed genetic information or a reference genome[18]. Therefore, in this study, we adopted this technology to uncover potential information about the pathogenic mechanisms of HFMD induced by EV-A71 and CV-A16 infection.

It was discovered that there was still a significant difference between the infections caused by EV-A71 and CV-A16, because the PCA data clearly showed that each group of EV-A71 infection and each group of CV-A16 infection were individually gathered together. Applying a Venn diagram enabled us to find 111 common differentially expressed genes that appear after EV-A71 and CV-A16 infection. These 111 common differentially expressed genes were used to perform hierarchical cluster analysis, and the heatmap result showed that these genes were mainly clustered into two categories—either significantly upregulated or significantly downregulated— suggesting that EV-A71 and CV-A16 have largely similar transcriptome-level changes after infection. Thus, these common genes might be closely related to the pathogenesis of HFMD. Next, we analyzed the expression trends of the 111 differentially expressed genes and found that they were all classified into 10 expression trends, including 5 upregulated trends after EV-A71 and CV-A16 infection, 4 downregulated trends after EV-A71 and CV-A16 infection, and 1 opposite trend after EV-A71 and CV-A16 infection. The differentially expressed genes that showed the same trends were the genes we chiefly focused on and were applied in the subsequent GO, pathway, and co-expression network analysis. The upregulated genes were enriched in 9 GO-BPs and mainly included Positive regulation of transcription, DNA-templated, DNA-templated transcription, initiation, Histone acetylation, Regulation of endocytosis, GPI anchor biosynthetic process, Protein polyubiquitination, Preassembly of GPI anchor in ER membrane, Cellular response to DNA Damage stimulus and DNA replication. The downregulated genes were enriched in 5 GO-BPs, mainly including Cellular response to heat, Positive regulation of phosphorylation, Branching involved in labyrinthine layer morphogenesis, Positive regulation of epithelial cell proliferation involved in lung morphogenesis, and Specification of loop of Henle identity. Previous studies have reported that the above GO-BPs are intimately associated with viral infections, e.g., Histone acetylation. Growing evidence has proven that histone deacetylase 6 (HDAC6) plays a very important role in natural immunity[19]. For example, in RNA virus-infected hosts, HDAC6 is able to bind to RIG-I and catalyze RIG-I deacetylation, which is essential for RIG-I to recognize double-stranded RNA. In addition, HDAC6 can also promote IFN production by catalyzing the deacetylation of β-catenin, which further hinders viral replication[20]. Hence, these studies imply that the enriched “Histone acetylation” might participate in the infection and replication process of EV-A71 and CV-A16. Moreover, the enriched GO-BP “Positive regulation of phosphorylation” was also observed to be involved in viral infections. For instance, Epstein–Barr virus (EBV)-encoded BGLF4 kinase could directly downregulate the NF-κB signaling pathway by phosphorylating coactivator UXT[21]. Furthermore, EBV nuclear antigen 1 (EBNA1) can stimulate its own nuclear entry by phosphorylating S385 in the nuclear localization signal[22]. Thus, these studies indirectly indicate that the enriched “Positive regulation of phosphorylation” might be a potential mechanism through which HFMD is induced by EV-A71 and CV-A16 infection. KEGG pathway enrichment analysis for these dysregulated genes is useful for revealing related pathways and molecular interactions among genes. The upregulated genes were markedly enriched in two pathways, namely, Glycosylphosphatidylinositol (GPI)-anchor biosynthesis and DNA replication, while the downregulated genes were only markedly enriched in one pathway, namely, Biosynthesis of antibiotics. Among these pathways, the GPI-anchor is considered to play an important role in the infection and pathogenesis of many viruses. For example, Enk et al. found that the herpes simplex virus (HSV)-1-encoded miR H8 could target the GPI anchoring pathway, which further reduced the expression levels of several immune-modulating proteins and finally enhanced viral spread and enabled evasion of natural killer cell elimination[23]. Amet et al. demonstrated that the deficiency of GPI-anchor could attenuate the production of infectious HIV-1 and render virions sensitive to complement attack[24]. Therefore, changes in the “GPI-anchor” pathway might promote the spread of EV-A71 and CV-A16 virus. In addition, it has been previously reported that EV-A71 can affect the DNA replication of host cells through its nonstructural protein 3D and then block it in S phase[25]. Moreover, another enterovirus (CV-A6) can also block host cells at G0/G1 with its nonstructural proteins 3C and 3D by affecting DNA replication[26]. Thus, it is clear that the enriched “DNA replication” pathway might contribute to the development of HFMD caused by EV-A71 and CV-A16. Regarding the “Biosynthesis of antibiotics” pathway, no reports have shown it to be associated with viral infection. However, since this pathway appeared in EV-A71 and CV-A16 infection, it might indicate that this pathway may be a new research direction. Ultimately, the coexpression network analysis for these dysregulated genes was carried out to seek key genes that regulate the pathogenesis of HFMD during EV-A71 and CV-A16 infection. In the coexpression network of upregulated genes, TBCK1 is located at a key node position. A previous study confirmed that mutations in the TBCK1 gene could lead to neurological diseases such as neuronal cerebello-lipidosis and neurodegeneration[27]. In addition, TBCK1 might play an important role in cell proliferation, cell growth and actin organization by modulating the mTOR pathway[28]. As mTOR is a pivotal pathway of autophagy activation, we speculated that TBCK1 may be involved in the induction of autophagy induced by EV-A71 and CV-A16 infection. In the coexpression network of downregulated genes, GPC4/6 is located at a key node position. However, GPC family proteins are mainly involved in tumorigenesis via regulating the Hedgehog signaling and Wnt signaling pathways[29, 30]. There was no indication that GPC had any function in immunity or virus infection.

In conclusion, transcriptome sequencing technology and bioinformatics approaches were employed to identify the differentially expressed genes in 16HBE cells in response to EV71 and CA16 infection. GO and pathway enrichment analysis, combined with the construction of coexpression networks can greatly contribute to a better understanding of the genes that are involved in EV71 and CA16 infection. However, the present study had several limitations. This study involved in vitro experiments, which might not completely reflect the in vivo setting; therefore, further in vivo investigations are necessary to provide more profound insights into the related host–pathogen interactions and pathogenesis.

EV-A71: Enterovirus 71; CV-A16: coxsackievirus A16; HFMD: hand, foot and mouth disease; 16HBE: human bronchial epithelial; Vigoo: Beijing Vigoo Biological Co., Ltd; Sinovac: Sinovac Biotech Co., Ltd; CAMS: Chinese Academy of Medical Science; CFDA: China Food and Drug Administration; CNS: central nervous system; RNA-Seq: RNA sequencing; SNPs: single nucleotide polymorphisms; RD cells: Rhabdomyoma cells; TXNIP: thioredoxin-interacting protein; RPMI: Roswell Park Memorial Institute; FBS: fetal bovine serum; MOI: multiplicity of infection; FDR: false discovery rate; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; HDAC6: histone deacetylase 6; EBV: Epstein–Barr virus; EBNA1: EBV nuclear antigen 1; HSV-1: herpes simplex virus-1.

Ethics Approval and consent to participate

Not applicable.

Consent for Publication

Not applicable.

Availability of data and material

The sequencing data were deposited in the National Center for Biotechnology Information’s Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/) under accession number GSE66757. All other data generated or analyzed in this study are included in this article.

Competing interests

The authors declare that they have no competing interests.

Funding

This study was funded by the CAMS Innovation Fund for Medical Sciences (2016-I2M-1-014), the National Natural Sciences Foundation of China (31700153 and 81373142), the Fundamental Research Funds for the Central Universities and PUMC Youth Fund (3332019004), the Yunnan Applied Basic Research Project (2018FB026), the Medical Reserve Talents of Yunnan Province Health and Family Planning (H-2017034) and Top young talents of Yunnan province ten thousand talents plan. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Authors' contributions

Designed experiment: LDL, SZD; Conducted the experiment: JS, WY L, YJH, HL; Analyzed data: HWZ, YLC; Wrote the manuscript: JS and YJH; All authors read and approved the final manuscript.

Acknowledgements

Not applicable.

Jin Y, Zhang R, Wu W, Duan G: Innate Immunity Evasion by Enteroviruses Linked to Epidemic Hand-Foot-Mouth Disease. Front Microbiol 2018, 9:2422.
Song J, Hu Y, Li H, Huang X, Zheng H, Hu Y, Wang J, Jiang X, Li J, Yang Z, et al: miR-1303 regulates BBB permeability and promotes CNS lesions following CA16 infections by directly targeting MMP9. Emerg Microbes Infect 2018, 7:155.
Yi EJ, Shin YJ, Kim JH, Kim TG, Chang SY: Enterovirus 71 infection and vaccines. Clin Exp Vaccine Res 2017, 6:4-14.
Mao QY, Wang Y, Bian L, Xu M, Liang Z: EV71 vaccine, a new tool to control outbreaks of hand, foot and mouth disease (HFMD). Expert Rev Vaccines 2016, 15:599-606.
Luo Z, Su R, Wang W, Liang Y, Zeng X, Shereen MA, Bashir N, Zhang Q, Zhao L, Wu K, et al: EV71 infection induces neurodegeneration via activating TLR7 signaling and IL-6 production. PLoS Pathog 2019, 15:e1008142.
Mao Q, Wang Y, Yao X, Bian L, Wu X, Xu M, Liang Z: Coxsackievirus A16: epidemiology, diagnosis, and vaccine. Hum Vaccin Immunother 2014, 10:360-367.
Marco-Puche G, Lois S, Benitez J, Trivino JC: RNA-Seq Perspectives to Improve Clinical Diagnosis. Front Genet 2019, 10:1152.
Stark R, Grzelak M, Hadfield J: RNA sequencing: the teenage years. Nat Rev Genet 2019, 20:631-656.
Lui YL, Tan TL, Timms P, Hafner LM, Tan KH, Tan EL: Elucidating the host-pathogen interaction between human colorectal cells and invading Enterovirus 71 using transcriptomics profiling. FEBS Open Bio 2014, 4:426-431.
Yao C, Hu K, Xi C, Li N, Wei Y: Transcriptomic analysis of cells in response to EV71 infection and 2A(pro) as a trigger for apoptosis via TXNIP gene. Genes Genomics 2019, 41:343-357.
Xu LJ, Jiang T, Zhang FJ, Han JF, Liu J, Zhao H, Li XF, Liu RJ, Deng YQ, Wu XY, et al: Global transcriptomic analysis of human neuroblastoma cells in response to enterovirus type 71 infection. PLoS One 2013, 8:e65948.
Jin J, Li R, Jiang C, Zhang R, Ge X, Liang F, Sheng X, Dai W, Chen M, Wu J, et al: Transcriptome analysis reveals dynamic changes in coxsackievirus A16 infected HEK 293T cells. BMC Genomics 2017, 18:933.
Song J, Hu Y, Hu Y, Wang J, Zhang X, Wang L, Guo L, Wang Y, Ning R, Liao Y, et al: Global gene expression analysis of peripheral blood mononuclear cells in rhesus monkey infants with CA16 infection-induced HFMD. Virus Res 2016, 214:1-10.
Zhao T, Zhang Z, Zhang Y, Feng M, Fan S, Wang L, Liu L, Wang X, Wang Q, Zhang X, et al: Dynamic Interaction of Enterovirus 71 and Dendritic Cells in Infected Neonatal Rhesus Macaques. Front Cell Infect Microbiol 2017, 7:171.
Ghosh S, Chan CK: Analysis of RNA-Seq Data Using TopHat and Cufflinks. Methods Mol Biol 2016, 1374:339-361.
Lee JT, Yen TY, Shih WL, Lu CY, Liu DP, Huang YC, Chang LY, Huang LM, Lin TY: Enterovirus 71 seroepidemiology in Taiwan in 2017 and comparison of those rates in 1997, 1999 and 2007. PLoS One 2019, 14:e0224110.
Liang Z, Wang J: EV71 vaccine, an invaluable gift for children. Clin Transl Immunology 2014, 3:e28.
Mellis R, Chandler N, Chitty LS: Next-generation sequencing and the impact on prenatal diagnosis. Expert Rev Mol Diagn 2018, 18:689-699.
Choi SJ, Lee HC, Kim JH, Park SY, Kim TH, Lee WK, Jang DJ, Yoon JE, Choi YI, Kim S, et al: HDAC6 regulates cellular viral RNA sensing by deacetylation of RIG-I. EMBO J 2016, 35:429-442.
Zhu J, Coyne CB, Sarkar SN: PKC alpha regulates Sendai virus-mediated interferon induction through HDAC6 and beta-catenin. EMBO J 2011, 30:4838-4849.
Chang LS, Wang JT, Doong SL, Lee CP, Chang CW, Tsai CH, Yeh SW, Hsieh CY, Chen MR: Epstein-Barr virus BGLF4 kinase downregulates NF-kappaB transactivation through phosphorylation of coactivator UXT. J Virol 2012, 86:12176-12186.
Sun D, Luthra P, Xu P, Yoon H, He B: Identification of a phosphorylation site within the P protein important for mRNA transcription and growth of parainfluenza virus 5. J Virol 2011, 85:8376-8385.
Enk J, Levi A, Weisblum Y, Yamin R, Charpak-Amikam Y, Wolf DG, Mandelboim O: HSV1 MicroRNA Modulation of GPI Anchoring and Downstream Immune Evasion. Cell Rep 2016, 17:949-956.
Amet T, Lan J, Shepherd N, Yang K, Byrd D, Xing Y, Yu Q: Glycosylphosphatidylinositol Anchor Deficiency Attenuates the Production of Infectious HIV-1 and Renders Virions Sensitive to Complement Attack. AIDS Res Hum Retroviruses 2016, 32:1100-1112.
Yu J, Zhang L, Ren P, Zhong T, Li Z, Wang Z, Li J, Liu X, Zhao K, Zhang W, Yu XF: Enterovirus 71 mediates cell cycle arrest in S phase through non-structural protein 3D. Cell Cycle 2015, 14:425-436.
Wang Z, Wang Y, Wang S, Meng X, Song F, Huo W, Zhang S, Chang J, Li J, Zheng B, et al: Coxsackievirus A6 Induces Cell Cycle Arrest in G0/G1 Phase for Viral Production. Front Cell Infect Microbiol 2018, 8:279.
Beck-Wodl S, Harzer K, Sturm M, Buchert R, Riess O, Mennel HD, Latta E, Pagenstecher A, Keber U: Homozygous TBC1 domain-containing kinase (TBCK) mutation causes a novel lysosomal storage disease - a new type of neuronal ceroid lipofuscinosis (CLN15)? Acta Neuropathol Commun 2018, 6:145.
Ortiz-Gonzalez XR, Tintos-Hernandez JA, Keller K, Li X, Foley AR, Bharucha-Goebel DX, Kessler SK, Yum SW, Crino PB, He M, et al: Homozygous boricua TBCK mutation causes neurodegeneration and aberrant autophagy. Ann Neurol 2018, 83:153-165.
Filmus J, Capurro M: The role of glypicans in Hedgehog signaling. Matrix Biol 2014, 35:248-252.
Capurro M, Martin T, Shi W, Filmus J: Glypican-3 binds to Frizzled and plays a direct role in the stimulation of canonical Wnt signaling. J Cell Sci 2014, 127:1565-1575.

Table S1. Sequences of primers used for qRT-PCR assays.

Fig. S1 Distribution of gene expression levels in each group.

PDF

Version 1

posted 13 Jan, 2020

You are reading this latest preprint version

Transcriptome analysis following EV-A71 and CV-A16 infection in respiratory epithelial cells

Status:

Version 1

Abstract

Figures

Background

Materials And Methods

Virus and cell lines

RNA extraction and RNA sequencing (RNA-seq) library construction

Computational analysis of RNA-seq data

Transcriptome assembly

Principal component analysis (PCA)

Identification of differentially expressed genes

Unsupervised hierarchical clustering analysis of the common differentially expressed genes

Functional group analysis

Coexpression network construction

Validation of differentially expressed genes by quantitative RT-PCR (qRT-PCR) analysis

Results

Overview of the RNA-Seq data from 16HBE cells in response to EV-A71 and CV-A16 infection

Expression profile of differentially expressed transcripts in EV-A71- and CV-A16-infected 16HBE cells

Genes analysis of the same trends using GO, Pathways and a Coexpression network

Validation of transcriptome sequencing results by qRT-PCR

Discussion

Conclusions

List Of Abbreviations

Declarations

References

Supporting Information

Supplementary Files

Status:

Version 1