Bioinformatics Analysis of Key micro-RNAs and mRNAs under the Hand, Foot, and Mouth Disease Virus Infection

Hand, foot, and mouth disease (HFMD) is a common endemic childhood disease worldwide, particularly in Asia. Generally, it is triggered by two major causative agents: enterovirus 71 (EV71) and coxsackievirus A16 (CA16) (Kaminska et al. 2013; Koh et al. 2016), followed by coxsackievirus A6 (CA6) and coxsackie virus A10 (CA10) (Wang et al. 2018). This disease majorly affects children under five years old (Omana-Capeda et al. 2016). The common symptoms include fever, rashes on the volar regions of the hands and feet, herpangina, and difficulties in eating and drinking. Severe circumstances can exhibit potentially fatal complications involving the nervous system, such as brain stem encephalitis, or cardiopulmonary systems, such as pulmonary edema. Severe cases may show drastic progression and die from complicated symptoms.

There have been limited tools for effective diagnosis of HFMD. Although enterovirus infection triggers various pathological responses, the exact mechanisms underlying HFMD development remain largely unknown. For example, the regulatory roles played by EV71/CA16 infection towards endothelial cells and neural systems remain unclear. Previous studies have explored the disordered signals in different aspects, such as inflammatory profiles in cytokine expression (Teo et al. 2018; Linghua et al. 2019), long non-coding RNA (lncRNA) profiles (Meng et al. 2017), and immune cell changes (Wang et al. 2014). It is conceivable and supported that patients with HFMD may exhibit changes in expression profiles of microRNAs and mRNAs, mostly derived from blood samples (Cui et al. 2011; Hu et al. 2016; Yee et al. 2016; Zhu et al. 2016; Song et al. 2017; Hu et al. 2018; Jia et al. 2018; Li et al. 2018; Mi et al. 2018; Song et al. 2018). In addition, a comprehensive understanding of HFMD-related expression profiles and identification key markers may provide substantial diagnostic and prognostic values. To dig the clinical significance of micro-RNA and mRNA changes under HFMD virus infection, we conducted this bioinformatics analysis to clarify crucial differentially expressed genes (DEGs), functional and pathway enrichment, and the relative regulatory network, using four published datasets. This analysis may provide deeper insight into the mechanism of HFMD pathological development and novel strategies to prevent HFMD outbreaks.

Differentially expressed miRNAs after HFMD virus infection. We first analyzed the GSE85829 and GSE94551 datasets and found some common differentially expressed miRNAs (DE-miRNAs). In 16HBE cells based on the GSE85829 dataset, EV71 infection-induced 25 DE-miRNAs and CA16 infection-induced 13 DE-miRNAs (Fig. 1A and 1B), among which seven common DE-miRNAs were found. In the HUVEC line, according to the GSE94551 dataset, 39 DE-miRNAs were observed after EV71 infection, and 99 DE-miR-NAs were found after CA16 infection (Fig. 2A and 2B), and they shared seven common DE-miRNAs (Fig. 2B). However, no intersection was found between two groups of seven common DE-miRNAs. A Venn diagram in Fig. 2C presented the common DE-miRNAs among four groups of cells. Similarly, miRNAs in serum exosomes were analyzed based on the GSE52780 dataset. A total of 258 DE-miRNAs were found (Fig. 3A), including 85 up-regulated and 173 down-regulated ones; the top ten DE-miRNAs were listed in Fig. 3B. Afterward, we screened those DE-miRNAs acquired three times in the above samples (cell lines or exosomes) in different datasets and regarded them as key miRNAs in HFMD. Ultimately, five key miRNAs were acquired: miR-100-3p (appeared four times), miR-125a-3p (appeared three times), miR-1273g-3p (appeared three times), miR-5585-3p (appeared three times), and miR-671-5p (appeared three times).

Fig. 1.

Fig. 2.

Fig. 3.

GO and KEGG enrichment. Based on the above five key miRNAs, GO functional, and KEGG pathway enrichment were performed using the miRPath v. 3 database. First, GO, and KEGG enrichment was acquired by the intersection targets of the five key miRNAs (only four miRNAs were included in this database) (Fig. 4A and 4B). Further, this was verified by the union of key miRNA related GO terms or KEGG pathways; Heatmaps are shown in Fig. 4C and 4D. Overall, this result was consistent with Fig. 4A and 4B. Parallelly, we applied the mRNA profiles in SH-SY5Y cells infected by EV71 from the GSE45589 dataset and performed the enrichment analysis. The enriched GO functions of up-regulated and down-regulated genes were listed in Fig. 4E and 4F, respectively. A total of 29 up-regulated GO functional terms were identified (meiotic nuclear division, response to radiation and regulation of the microtubule-based process, etc.), accomplished by 43 down-regulated terms (positive regulation of transcription, DNA-templated, negative regulation of the apoptotic process, and anatomical structure morphogenesis, etc.). Ten enriched KEGG pathways were shown in Fig. 4G, including two upregulated ones (cell cycle and spliceosome) and eight down-regulated pathways (cytokine-cytokine receptor interaction, hematopoietic cell lineage, and intestinal immune network for IgA production, etc.). In comparison, there were three common enriched GO terms between miRNA-derived prediction and mRNA-derived (dataset GSE45589) analysis: biosynthetic process, cytosol, and nucleoplasm. Moreover, one common KEGG pathway, namely cell cycle, was shared between miRNA-based and mRNA-based enrichment.

Fig. 4.

Common differential mRNAs and PPI network. Using TarBase V8 in DIANA tools, we acquired 1,520 potential targets (mRNA) from the five key DE-miR-NAs, among which the159 DE-mRNAs also included 11 DEGs in the GSE45589 dataset: MACF1, MARS, SF3B3, SMARCC1, BRMS1L, SMC1A, SPHK2, LIG1, CSF3, CYR61, and FGFR1OP (Fig. 5A). Theoretically, these genes were the most likely to be influenced by HFMD virus infection. GO functional analysis showed three terms might be enriched according to these DEGs: positive regulation of cell proliferation, anatomical structure morphogenesis, and ATP binding (Fig. 5B). These common DEGs showed a PPI network mainly connected by SMC1A, SMARCC1, SF3B3, LIG1, and BRMS1L (Fig. 5C), and this network locates at a core place in the PPI network constructed by the 159 DE-mRNAs in the GSE45589 dataset (Fig. 5D, the isolated nodes were removed). Together, changes in five key miRNAs and 11 key mRNAs may play crucial roles in HFMD virus-induced pathological changes and can be used as diagnostic markers for the HFMD.

Fig. 5.

In this study, we used five datasets to identify key RNA members in HFMD development. We were interested in five key miRNAs, 11 mRNAs, and several important GO and KEGG enrichment after filtering progressively. Our results might provide some theoretical perspectives about HFMD development and a potential strategy in its early warning.

At the miRNA level, several potentially useful markers have been proposed in clinical diagnosis. A survey in Singapore reported a 6-miRNA scoring model which predicts HFMD with an overall accuracy of 85.11% in the training set and 92.86% in the blinded test set, and circulating salivary miRNA hsa-miR-221 (downregulated in that work) was regarded as a highly validated marker (Mi et al. 2018). Song et al. have applied rhesus monkey peripheral blood mononuclear cells to search DE-miRNAs, and they identified 13 novel DE-miRNAs with 2501 targets (Song et al. 2018). Zhu et al. (2016) performed the microarray examination and noticed 27 DE-miRNAs (15 up-regulated and 12 downregulated) associated with CA16 and EV71 infection. There were some other important findings of specific miRNAs. MiR-1303 has been known to promote CNS lesions following CA16 infections by targeting MMP9 (Song et al. 2018). EV71 can evade the immune surveillance system to proliferate by activating miR-21 (Feng et al. 2017), antagonize the antiviral activity of host STAT3 and IL-6R through miR-124 (Chang et al. 2017), and induce autophagy by regulating miR-30a to promote viral replication (Fu et al. 2015). So far, there have been very few direct reports about the relationship between the five key miRNAs and HFMD. Only one study had surveyed the miRNA expression profile in the exosome of HFMD patients (Jia et al. 2014), and it revealed that the expression level of three miRNAs (miR-671-5p, miR-16-5p, and miR-150-3p) was significantly different between mild HFMD, extremely severe SHFMD, and the healthy controls. We also noticed that miR-671-5p was among the key miRNAs in HFMD.

The PPI network implied that five targets, SMC1A, SMARCC1, SF3B3, LIG1, and BRMS1L, might play the most crucial roles during HFMD progression. However, none of them has been paid enough attention to date, and they are worth more concerns in further researches. Taken different datasets together, we found three common enriched GO terms between miRNA-derived prediction and mRNA-derived analysis: biosynthetic process, cytosol, and nucleoplasm; and a common KEGG pathway, cell cycle, was noticed (Fig. 4). These functions and pathways suggest that HFMD viruses strongly drive the host proliferation. This fact could be also be supported by the GO functional enrichment constructed based on the 11 key mRNAs (Fig. 5B), which exhibited that positive regulation of cell proliferation was the most enriched functional term. This process may contribute to virus amplification and also be a homeostasis response to fight against virus invasion, particularly for epithelial cells. However, the definite mechanism needs more evidence to unravel.

However, this study has some limitations. First, when we probed the key roles, common DEGs were screened between EV71 and CA16 infection. However, despite belonging to the same genus, Enterovirus, these two viruses do not necessarily have similar effects. For example, Chinese scholars identified that miR-4516 presented down-regulation in EV71 infection and upregulation in CA16 infection, and it was an important regulator of intercellular junctions by targeting PVRL1 (Hu et al. 2016). Liu et al. had analyzed miRNA profiles and acquired diverse outcomes induced by EV71 and CA16 infection (Song et al. 2017). The inconsistency was also shown in Fig. 2B. There were two common up-regulated miRNAs, two common down-regulated miRNAs, and three inconsistent ones (miR-502-5p, miR-503-5p, and miR-542-3p). Besides, the direct regulatory relationship between five key miRNAs and 11 key mRNAs had not been validated in the present study. Our further efforts would focus on the construction of diagnostic and prognostic models based on the large real-world sample using these miRNA and mRNA factors.

Our results shed light on the potentially crucial roles of the five miRNAs, 11 coding genes, and several functions and pathways in HFMD. A combination of these roles may benefit the early diagnosis and treatment of HFMD.