key: cord-032801-b2ncmkjg
authors: Song, Jie; Hu, Yajie; Li, Weiyu; Li, Hui; Zheng, Huiwen; Chen, Yanli; Dong, Shaozhong; Liu, Longding
title: Transcriptome analysis following enterovirus 71 and coxsackievirus A16 infection in respiratory epithelial cells
date: 2020-09-29
journal: Arch Virol
DOI: 10.1007/s00705-020-04821-1
sha: 
doc_id: 32801
cord_uid: b2ncmkjg

Enterovirus 71 (EV-A71) and coxsackievirus A16 (CV-A16) are the major pathogens responsible for hand, foot and mouth disease (HFMD), but the mechanism by which these viruses cause disease remains unclear. In this study, we used transcriptome sequencing technology to investigate changes in the transcriptome profiles after infection with EV-A71 and CV-A16 in human bronchial epithelial (16HBE) cells. Using 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 EV-A71- and CV-A16-infected cells. A trend analysis of these 111 genes showed that 91 of them displayed the same trend in EV-A71 and CV-A16 infection, including 49 upregulated genes and 42 downregulated genes. These 91 genes were further used to conduct GO, pathway, and coexpression network analysis. It was discovered that 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 these 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. Collectively, our results provide clues to the mechanism of pathogenesis of HFMD induced by EV-A71 and CV-A16. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (10.1007/s00705-020-04821-1) contains supplementary material, which is available to authorized users.

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 the Asia-Pacific region in infants and young children. Both of these two viruses belong to the species Enterovirus A, family Picornaviridae. They have a positivesense single-stranded RNA genome with an approximate length of 7400 nucleotides, consisting of four structural viral proteins (VP1 to VP4) and seven nonstructural proteins (2A to 2C and 3A to 3D) [10] . Clinically, most HFMD cases are 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 [22] . In recent years, the number of these serious cases and deaths has been increasing and has caused 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 on the market and has been approved by the China Food and Drug Administration (CFDA) [28] . 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 [17] . In previous studies, it was demonstrated that HFMD patients with EV-A71 and CV-A16 infections presented with significantly different clinical manifestations. EV-A71 infections frequently lead to severe central nervous system (CNS) complications and even death, while CV-A16 infections often result in milder symptoms that resolve within a few weeks [15] . Nevertheless, over the past years, accumulating evidence has revealed that HFMD in patients infected with CV-A16 can also develop into a serious stage with neurological complications, and the overall condition of these patients and their clinical symptoms are generally consistent with those of severe HFMD patients infected with EV-A71 [16] . 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 necessary to strengthen basic research on the replication and pathogenesis of EV-A71 and CV-A16, which could provide a theoretical basis for the development of HFMD-specific therapeutic drugs. RNA sequencing (RNA-Seq) can be used to investigate differences in gene expression at a genome-wide level. RNA-Seq has the advantages of more accurate quantification, higher repeatability, a wider detection range, and more reliable analysis than other methods [18] . In addition to analyzing gene expression levels, RNA-Seq can also be used to discover new transcripts, single nucleotide polymorphisms (SNPs), and splice variants to provide information about allele-specific gene expression [23] . 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 a signaling pathway that might participate in inhibiting viral replication [14] . Yao et al. carried out a transcriptome analysis of EV-A71-infected human rhabdomyosarcoma (RD) cells and found that EV-A71-2A protein could be considered a key inducer that triggered cellular apoptosis and death in RD cells mediated by thioredoxin-interacting protein (TXNIP) [27] . 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 [26] . 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 [9] . Song et al. also studied 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 notably increased after CV-A16 infection [21] . 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 an important route for HFMD transmission; therefore, the interaction between EV-A71/CV-A16 and airway epithelial cells should be investigated [30] . 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 by bioinformatics analysis, which might provide clues about the mechanisms of HFMD and possible molecular targets for HFMD treatment.

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 5 cells per well in 6-well sterile plastic culture plates and grown in a base medium of Roswell Park Memorial Institute (RPMI)-1640 medium containing 10% (vol/vol) fetal bovine serum (FBS), 100 units of penicillin per mL, 100 μg of streptomycin per mL and 2 mM L-glutamine at 37°C in a 95% air and 5% carbon dioxide (CO 2 ) 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 no. EU812515.1) or coxsackievirus A16 (CV-A16; subgenotype B, GenBank no. JN590244.1), which were isolated from an epidemic in Fuyang, China, in 2008 and from an HFMD patient in Guangxi, China in 2010, respectively. Both viruses were incubated with cells for one hour to attain virus attachment, and RPMI-1640 medium containing 1% FBS was then added. The cells were incubated for a further 6 h, 12 h and 24 h for virus propagation. It was observed that the proportion of cells infected with each of these viruses at the 24 hour time point had reached about 90% (Fig. S1 ). At each time point, cells were collected using RPMI-1640 medium and centrifuged twice in phosphate-buffered saline (PBS) at 4°C for 5 min at 3000 rpm. Finally, the harvested cell pellets were snapfrozen 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 in place of the virus.

Total RNA was extracted using an RNeasy ® Mini Kit (QIA-GEN, USA) as recommended by the manufacturer. The extracted RNA was cleared of contaminating genomic DNA by treatment with RNase-free DNase I (Takara Bio, Japan). The concentration of RNA was measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA), and the RNA quality was determined using an Ultrospec 3000 Pro UV/Visible spectrophotometer (GE Healthcare, UK). Samples with an absorbance ratio (A260/A280) of 1.8 to 2.0 were subjected to RNA integrity analysis using an Agilent ® 2100 Bioanalyzer. Samples with an RNA integrity number (RIN) greater than 7 were used for RNA-seq library construction. Briefly, mRNA was enriched using 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 make first-strand 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 using Phusion High-Fidelity DNA Polymerase (NEB, Beijing, China), and the products were sequenced on an Illumina HiSeq™ 2000 platform (Illumina, USA).

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.bioin forma tics.babra ham.ac.uk/proje cts/ fastq c/). The obtained "clean" reads were mapped against the human reference genome GRCh38 using 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) [8] .

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.

For gene expression analysis, the expression level of each gene was calculated by determining the number of reads and was further normalized by a variation of the FPKM method. To identify genes that were differentially expressed 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 the Benjamini and Hochberg approach to control the false discovery rate (FDR). Genes with a two-fold change and an adjusted P-value < 0.05 were regarded to be differentially expressed.

To find commonalities between EV-A71-and CV-A16-infected 16HBE cells, we constructed a Venn diagram using the 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.

The Database for Annotation, Visualization and Integrated Discovery (http://david .abcc.ncifc rf.gov/), which utilizes Gene Ontology (GO) to identify 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.genom e.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 was based on the Gen-MANIA algorithm by using the common differentially expressed genes. The construction of coexpression networks is conducive to finding potential mechanisms associated with differentially expressed genes.

To validate the transcriptome results, we designed 10 pairs of primers using Primer Premier 5.0 to perform qRT-PCR analysis targeting seven upregulated (BLM, GBP3, HDAC9, NNMT, PNISR, RNF6 and TAF1) and three downregulated genes (ANO1, IRX2 and PCDH7). Briefly, reverse transcription was first performed using total RNA (1 μg) isolated as described above with a 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 °C for 30 s, followed by 45 cycles of denaturation at 95 °C for 10 s, annealing at 55 °C for 15 s, and extension at 72 °C for 30 s on a 7500 Fast Real-Time 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 by the 2 -ΔΔCt method. Three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed.

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 sequence GRCh38, 198.68 million reads (94.51%) were mapped to the genome. The detailed data on 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 shifted from the control group. Thus, these data showed differences between the infection groups.

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 twofold change and an adjusted P-value < 0.05 were screened. The distribution of gene expression levels is shown in Fig. S2 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 number of upregulated and downregulated genes in each group is 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.

To further explore commonalities between the EV-A71and CV-A16-infected 16HBE cells, trend analysis of the 111 common differentially expressed genes was carried out to reveal differentially expressed genes that had the same change trends. The results showed that 49 upregulated and 52 downregulated genes showed 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 used separately to analyze the GO, related pathways, and coexpression network. Regarding the GO BP terms, the upregulated genes were enriched in nine terms (Fig. 3) , and the downregulated genes were enriched in five terms (Fig. 4) . The upregulated genes were markedly enriched in five MF terms (Fig. 3) , and the downregulated genes were markedly enriched in three MF terms (Fig. 4) .

Regarding the CC terms, five terms were enriched in the upregulated genes (Fig. 3) , and five CC terms were enriched in the downregulated genes (Fig. 4) . Additionally, KEGG pathway enrichment analysis for the significantly differentially expressed genes was used to identify the pathways related to the genes. Our data showed that only two pathways were associated with the upregulated genes (Fig. 5A) and one was related to the downregulated genes (Fig. 5B) . Ultimately, the construction of a coexpression network was implemented to identify the molecular interactions (Fig. 6 ).

Ten differentially expressed genes from the RNA-seq data were randomly chosen for validation by qRT-PCR (Fig. 7) . 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 [11] . Vaccination is considered to be one of the most effective ways to protect against EV-A71 and CV-A16 infections. 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 [13] . 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 detailed genetic information or a reference genome is lacking [19] . There have been many reports on transcriptome analysis of the effect of EV-A71

and CV-A16 infection on different sensitive cells. For example, transcriptome analysis revealed differentially expressed genes, including SCARB2, miR-3605-5p, and enriched ECM-receptor interaction and circadian rhythm pathways involved in the pathogenesis of HFMD in CV-A16-infected HEK 293T cells [9] . Characterization of critical functions of lncRNAs and mRNAs in RD cells and mouse skeletal muscle infected by EV-A71 using RNA-Seq demonstrated that lncRNAs may participate in EV-A71 infection-induced pathogenesis by regulating immune responses, protein binding, cellular component biogenesis, and metabolism [12] . Data from transcriptomics profiling in EV-A71-infected human colorectal cells showed altered expression of human genes involved in critical pathways, including the immune response and the stress response, which provided valuable insights into the host-pathogen interaction between human colorectal cells and EV-A71 [14] . These studies indicated that the changes in the transcriptome of cells infected with EV-A71 and CV-A16 may be different when different cells are examined, which is of great significance for exploring potential pathogenic mechanisms. However, at present, there are no relevant studies on transcriptome comparison between EV-A71 and CV-A16 infections in 16HBE cells. Hence, in this study, we adopted this technology to obtain information about the pathogenic mechanisms of HFMD induced by EV-A71 and CV-A16 infection. It was discovered that there was a significant difference between infections caused by EV-A71 and CV-A16, because the PCA data clearly showed that the groups for EV-A71 infection and CV-A16 infection were individually gathered together. Applying a Venn diagram enabled us to find 111 common genes whose expression changed 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 infections result in largely similar transcriptome-level changes. Thus, these common genes might be involved in 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 five upregulated trends after EV-A71 and CV-A16 infection, four downregulated trends after EV-A71 and CV-A16 infection, and one opposite trend after EV-A71 and CV-A16 infection. The differentially expressed genes that showed the same trends were the genes that we chiefly focused on and were applied in the subsequent GO, pathway, and co-expression network analysis. The upregulated genes were enriched in nine GO-BPs and mainly included positive regulation of transcription, 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 five 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 shown that the above GO-BPs are intimately associated with viral infections. There is growing evidence that histone deacetylase 6 (HDAC6) plays a very important role in natural immunity [5] . 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 [31] . 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 [4] . Furthermore, EBV nuclear antigen 1 (EBNA1) can stimulate its own nuclear entry by phosphorylating S385 in the nuclear localization signal [24] . 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 identifying 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 elimination by natural killer cells [6] . Amet et al. demonstrated that a deficiency of GPI-anchor could inhibit the production of infectious HIV-1 and render virions sensitive to complement attack [1] . Therefore, changes in the "GPI-anchor" pathway might promote the spread of EV-A71 and CV-A16. In addition, it has been reported that EV-A71 can affect DNA replication in host cells through its nonstructural protein 3D and then block the cells in S phase [29] . Moreover, another enterovirus (CV-A6) can also block host cells at G0/G1 through its nonstructural proteins 3C and 3D by affecting DNA replication [25] . Thus, it is clear that the enriched "DNA replication" pathway might contribute to the development (C) Venn diagram displaying the number of EV-A71-specific, CV-A16-specific, and EV-A71/ CV-A16-common genes at different time points after EV-A71 and CV-A16 infection. (D) Hierarchical heatmap showing the 111 overlapping differentially expressed genes in each group. Red indicates genes that were expressed at higher levels, and green indicates genes that were expressed at lower levels. Each column represents one group, and each horizontal line refers to one gene. The cutoff values of log fold change >2 or < -2 and a false discovery rate of <0.05 were applied ◂ Fig. 2 Trend analysis of the differentially expressed genes in 16HBE in response to EV-A71 and CV-A16 infection at different time points postinfection. The genes that were upregulated in both the EV-A71 and CV-A16 infections are indicated in red. The genes that were downregulated in both the EV-A71 and CV-A16 infections are indicated in blue. The genes that showed opposite changes in the EV-A71 and CV-A16 infections are indicated in black 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, more research on this pathway might be warranted. Coexpression network analysis for these dysregulated genes was carried out to identify 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 [2] . In addition, TBCK1 might play an important role in cell proliferation, cell growth, and actin organization by modulating the mTOR pathway [20] . As mTOR is a pivotal pathway of autophagy activation, we speculate 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 [3, 7] . 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 EV-A71 and CV-A16 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 EV-A71 and CV-A16 infection. However, the present study had several limitations. This study involved in vitro experiments, which might not 

Conflict of interest The authors declare that they have no conflicts of interest.

Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.

Glycosylphosphatidylinositol anchor deficiency attenuates the production of infectious HIV-1 and renders virions sensitive to complement attack

Homozygous TBC1 domain-containing kinase (TBCK) mutation causes a novel lysosomal storage disease-a new type of neuronal ceroid lipofuscinosis (CLN15)?

Glypican-3 binds to Frizzled and plays a direct role in the stimulation of canonical Wnt signaling

Epstein-Barr virus BGLF4 kinase downregulates NF-kappaB transactivation through phosphorylation of coactivator UXT

HDAC6 regulates cellular viral RNA sensing by deacetylation of RIG-I

HSV1 microrna modulation of GPI anchoring and downstream immune evasion

The role of glypicans in Hedgehog signaling

Analysis of RNA-Seq data using TopHat and CUFFLINKS

Transcriptome analysis reveals dynamic changes in coxsackievirus A16 infected HEK 293T cells

Innate immunity evasion by enteroviruses linked to epidemic hand-foot-mouth disease

Enterovirus 71 seroepidemiology in Taiwan in 2017 and comparison of those rates in 1997

Characterization of critical functions of long non-coding RNAs and mRNAs in rhabdomyosarcoma cells and mouse skeletal muscle infected by enterovirus 71 using RNA-Seq

EV71 vaccine, an invaluable gift for children

Elucidating the host-pathogen interaction between human colorectal cells and invading Enterovirus 71 using transcriptomics profiling

EV71 infection induces neurodegeneration via activating TLR7 signaling and IL-6 production

Coxsackievirus A16: epidemiology, diagnosis, and vaccine

Differential gene expression confirmed by qRT-PCR. The RNA-seq data and qRT-PCR validation results were compared

EV71 vaccine, a new tool to control outbreaks of hand, foot and mouth disease (HFMD)

RNA-Seq perspectives to improve clinical diagnosis

Next-generation sequencing and the impact on prenatal diagnosis

Homozygous boricua TBCK mutation causes neurodegeneration and aberrant autophagy

Global gene expression analysis of peripheral blood mononuclear cells in rhesus monkey infants with CA16 infectioninduced HFMD

regulates BBB permeability and promotes CNS lesions following CA16 infections by directly targeting MMP9

RNA sequencing: the teenage years

Identification of a phosphorylation site within the P protein important for mRNA transcription and growth of parainfluenza virus 5

Coxsackievirus A6 induces cell cycle arrest in G0/G1 phase for viral production

Global transcriptomic analysis of human neuroblastoma cells in response to enterovirus type 71 infection

Transcriptomic analysis of cells in response to EV71 infection and 2A(pro) as a trigger for apoptosis via TXNIP gene

Enterovirus 71 infection and vaccines

Enterovirus 71 mediates cell cycle arrest in S phase through non-structural protein 3D

Dynamic interaction of enterovirus 71 and dendritic cells in infected neonatal rhesus macaques

PKC alpha regulates Sendai virus-mediated interferon induction through HDAC6 and betacatenin