Transcriptome analysis following enterovirus 71 and coxsackievirus A16 infection in respiratory epithelial cells.

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.

Introduction

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 positive-sense 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.

Materials and methods

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 × 105 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 (CO2) 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 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 in place of the virus.

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

Total RNA was extracted using an RNeasy® Mini Kit (QIAGEN, 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).

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 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].

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 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.

Unsupervised hierarchical clustering analysis of the common differentially expressed genes

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.

Functional group analysis

The Database for Annotation, Visualization and Integrated Discovery (http://david.abcc.ncifcrf.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.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

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 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.

Results

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 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.

Table 1

Overview of sequencing data information

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

Fig. 1

(A) Principal component analysis (PCA) plot of the differentially expressed genes for determining the degree of aggregation and dispersion between groups. (B) The detailed number of upregulated and downregulated differentially expressed genes in each group. (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

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 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.

Identification of genes with 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 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).

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

Fig. 3

GO functional analysis of differentially expressed genes. The abnormal expression levels of upregulated genes were analyzed. The results are summarized in three main categories: biological process, molecular function, and cellular component

Fig. 4

GO functional analysis of differentially expressed genes. The abnormal expression levels of downregulated genes were analyzed. The results are summarized in three main categories: biological process, molecular function, and cellular component

Fig. 5

KEGG pathway enrichment of the differentially expressed genes. (A) KEGG analysis for upregulated genes. (B) KEGG analysis for downregulated genes

Fig. 6

Coexpression networks constructed using the upregulated (A) and downregulated (B) differentially expressed genes. Red circles represent the upregulated differentially expressed genes, dark blue circles represent the downregulated differentially expressed genes, and yellow circles represent the predicted coexpressed genes. The size of the circle reflects the fold change in gene expression. The connections between the genes are shown as different colored lines representing different connection relationships

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. 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.

Fig. 7

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

Discussion

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 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 in tumorigenesis via regulation of the Hedgehog signaling and Wnt signaling pathways [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 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.

Below is the link to the electronic supplementary material.

Fig. S1 The proportion of EV-A71- and CV-A16-infected cells measured using an immunofluorescence assay at 24 h

Fig. S2 Distribution of gene expression levels in each group

Supplementary material 3 (DOCX 14 kb)