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

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.

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 and ECM receptor [12]. Song et al. also studied the transcriptome changes in CV-A16infected 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-A71and 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 5 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 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: 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. 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 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-
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 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  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 RNAseq expression profiles.

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

Conclusions
In

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. Supporting Information Table S1. Sequences of primers used for qRT-PCR assays.      Differentially expressed genes were validated by qRT-PCR. The RNA-seq data and qRT-PCR validation results were compared. This is a list of supplementary files associated with the primary manuscript. Click to download.