Characterization of Enterovirus Associated m6A RNA Methylation in Children With Neurological Symptoms: A Prospective Cohort Study

Little is known about the particular changes of N6-methyladenosine (m6A) RNA methylation in enterovirus (EV) infection among children with neurologic symptoms. Here, we determined the characterization of EV associated m6A RNA methylation in this population. A prospective cohort study was conducted from 2018/2 to 2019/12 at the Guangzhou Women and Children’s Medical Center. We included EV infected children with and without neurological symptoms. High-throughput m(6)A-RNA immunoprecipitation sequencing (MeRIP-seq) and RNA-seq analysis were used to evaluate the m6A RNA methylation and transcript expression of cerebrospinal fluid samples. The functional annotation and pathways of differentially methylated m6A genes with synchronously differential expression were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Seven patients were enrolled in the control group, and 13 cases were in the neurological symptoms (NS) group. A total of 3472 differentially expressed genes and 957 m6A modified genes were identified. A conjoint analysis of MeRIP-seq and RNA-seq data found 1064 genes with significant changes in both the m6A modifications and mRNA levels. The different m6A RNA methylation was increased in the transcriptome’s CDS regions but decreased in both the 3′UTRs and stop codon among the NS group. Functional annotation like the “oxidative phosphorylation” gene pathway, “Parkinson’s disease” and GO terms like “respiratory electron transport chain,” “cellular metabolic process,” and “oxidation-reduction process” was enriched in symptomatic patients. Our study elucidated the changes of RNA m6A methylation patterns and related cellular functions and signaling pathways in EV patients with neurologic symptoms.

Introduction

Enteroviruses (EV) are common pathogens that causes an array of diseases, especially in neonates. This includes hand-foot and mouth disease (HFMD), aseptic meningitis, encephalitis, acute flaccid paralysis, and acute flaccid myelitis (Gonzalez et al., 2019). It has several genotypes, among which, EV-A71, CV-A6 and CV-A16 have caused large epidemics in Asia since 1997 (Lee et al., 2014; Wang et al., 2017; Huang et al., 2018). Recently, a multicenter prospective study found that the leading cause of viral encephalitis and meningitis in Chinese children was the human EV (Ai et al., 2017). Enteroviruses infection is known to cause thalamus and medulla oblongata damage, thus induces a persisting neurological sequelae or even death (Gonzalez et al., 2019). Therefore, the treatment for EV infections remained mainly supportive apart from the recent anecdotal use of small molecules like ribavirin and DTriP-22 that inhibit EV-A71 by blocking the polymerase (Li et al., 2008). Although our understanding of the molecular mechanisms involved in EV infections has increased dramatically in the past few years, the exact mechanisms of enterovirus infection still remain unknown, especially with regards to the pathogenesis of neurologic damage.

Similar to proteins and DNA, RNA modifications are involved in many aspects of biological functions, like gene expression, protein translation, cell behaviors, and physiological conditions, of which m6A modification of mRNAs is the major form (Du et al., 2019). It has been reported that m6A modification is associated with many human diseases, for example, obesity (Frayling et al., 2007; Ho et al., 2010) and a significant number of cancers (Jin et al., 2012; Reddy et al., 2013; Lin et al., 2016). Also, it plays an important role in brain development as well as a variety of neurological disorders (Du et al., 2019). Some studies have identified that m6A tagged genes are associated with the development of mental disorders, intellectual disability, schizophrenia, and bipolar disorders (Yoon et al., 2017). Recently, the association between RNA methylation and viral infections have gained popularity among researchers especially with the advent of transcriptome sequencing technology. Lichinchi et al. (2016b) found that the host RNA methyltransferase plays a negative post-transcriptional regulatory role in ZIKV virus. Some reports also described the process of m6A modification in HIV-1 RNA and its mechanisms in which it affects viral gene expression (Kennedy E et al., 2016; Lichinchi et al., 2016a). Despite ongoing progress in the field, little is known about the particular changes of m6A RNA methylation in EV infection among children with clinical neurologic symptoms.

In this study, we used RNA-sequencing technology to profile the transcriptome-wide alterations of m6A RNA methylation in cerebrospinal fluid (CSF) samples of simple enterovirus infected children with and without neurologic symptoms, with the aim of understanding the underlying mechanisms of neurologic damage caused by EV infections. This project used the SingleSiteMod program to perform single-base site analysis. The aim was to find the position of the “GGAC” sequence in the m6A peak region, which will reflect the exact position and number of m6A modifications in a peak region. This result has an important reference value for studying m6A modification. We hope this study will contribute to highlighting the pathogenesis of EV infections and will further provide countermeasures to mitigate disease chains of transmission.

Materials and Methods

Ethics Approval and Consent to Participate

This study was approved by the Ethics Committee of the Guangzhou Women and Children’s Medical Center (No.2017122501 and 2021019A01). All procedures performed involving human participants were in accordance with the ethical standards of the institutional and/or National Research Committee and in line with the 1964 Helsinki declaration. All patients signed an informed consent form upon admission.

Study Subjects and Sampling

A descriptive, prospective cohort study was conducted from January 2018 to December 2019 in the Guangzhou Women and Children’s Medical Center. Detailed demographic, clinical characteristics at admission, biochemistry, CSF, and hematologic indicators of the enrolled patients were extracted from the structured electronic medical records system (EMRS). The earliest value of indicators within 24 h after admission was used. Pediatric patients whose CSF samples were collected within 24 h of admission were enrolled. Patients with either positive throat or anal swabs or those with a positive CSF EV-Rt-PCR result were confirmed as EV infection (Zhou et al., 2015). We included all children with EV infection, including those with and without neurological symptoms. Among the patients with neurology involvement (NS group), convulsions, vomiting, headache or altered state of consciousness occurred showing in Table 1. The patients who had fever of unknown origin (FUO) or infection lesions were not clear with high inflammatory indicators had undergone the CSF test and was divided into control group. We excluded those coinfected with other pathogens. Patients with or without neurologic symptoms were matched for age to reduce bias, and then divided into control and NS group.

TABLE 1

Table 1. Demographic and clinical characteristics of enrolled patients.

Pharyngeal swabs, anal swabs, and CSF samples were collected for detection of EV. Real-time reverse transcription-polymerase chain reaction (RT-PCR) was performed using enterovirus dual fluorescent quantitative RT-PCR kit (Guangdong huayin pharmaceutical technology co. LTD, China). For the clinical routine biochemical testing, we used an automatic analyzer (Automatic Analyzer 7600, Hitachi High-Technologies Corporation). An automatic blood and body fluid analyzer (Sysmex XN550, Japan) was used for the detection of serum and CSF cells.

RNA Library Construction and Sequencing

For transcriptome profiling, the CSF samples were thawed at 37°C and centrifuged for 30 min at 2000 × g at 4°C to clear cells. Then the supernatant was transferred to a new tube and centrifuged at 10,000 × g at 4°C for 45min to remove cell fragment. After that, the supernatant was filtered by 0.45 μm filter membrane, and the filtrate was collected. RNA was then extracted from the filtrate and purified by the Zymo RNA clean and concentrator-5 kit (Zymo Research, cat. no R1013). Purified RNA was cut into 200 nt lengths by hyperthermia and added into the precipitation buffer containing the anti-m6A Antibody (Sigma-Aldrich: ABE572), dynabeads Protein A (Invitrogen™, cat. no 10002D) and dynabeads Protein G (Invitrogen™, cat. no 10004D). After magnetic separation and supernatant removal, RNA enzyme inhibitor was added and incubated at 4°C for 1-3 h, and then washed three times by low and high salt precipitation buffer. RNA in the eluate was purified and used for library generation. The obtained products were sequentially removed by ribosomal RNA and synthesized by SMART technique into the first strand cDNA then amplified by PCR and enriched by library fragments, and then DNA purified magnetic bead library fragments were obtained to detect the library by ultrahigh RNA methylation of m6A. Bioptic Qsep100 Analyzer was then used for quality control of the completed libraries.

Sequencing was performed using NovaSeq PE150 in accordance with the manufacturer’s recommendations. HISAT2 software (version 2.1.0) was used to compare the filtered clean reads with the reference genome of the corresponding species of the sample to obtain unique mapped reads for further analysis. A R/Bioconductor package “exomePeak” was used for peak calling. P < 0.05 was considered to be differentially peaks distribution. We used HOMER software (version 4.10.4) to perform motif analysis on the Peaks¹. SingleSiteMod program was used to determine the position and number of m6A modification in the peak region. The raw sequence reads were deposited into the Genome Sequence Archive (GSA) database (File number: OMIX716).

Conjoint Analysis of MeRIP-seq and RNA-seq Data

To identify the genes with significant changes in both the m6A modification and mRNA levels, conjoint analysis of MeRIP-seq and RNA-seq (MeRIP-seq input library data) data was performed.

The expression level of each gene was calculated by determining the length of the fragments and the read counts mapped to those fragments and were further normalized by a variation of the FPKM method. In our study, the differentially expressed genes were screened and sorted according to log 2FC (log 2FC > 0, up-regulated genes; log 2FC < 0, down-regulated genes), P value, FDR of RNA-sequence difference. The “exomePeak” was used for differential RNA methylation analysis. P < 0.05 was considered to be differentially peaks distribution.

Gene Ontology Annotations and Kyoto Encyclopedia of Genes and Genomes Enrichment

The Gene Ontology (GO) was used to identify the functional annotation and the classification of molecular functions, together with the biological processes and cellular component aspects of the identified m6A methylated RNAs and the differentially m6A methylated RNAs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was applied to identify pathways of differentially expressed genes (DEGs) in this study. Meanwhile, Fisher’s exact test was used to calculate the P-value according to the annotations, and the rich factor was calculated based on the numbers of symbols in the list.

Public Databases and Analysis

We downloaded 31 human tissues from Genotype-Tissue Expression (GTEX), covering the expression levels of 60498 genes (FPKM). According to the method in a previous report (Gerstberger et al., 2014), the tissue specificity score (S) of each gene was calculated and its value ranged from 0 to 4.95. The higher the value, the higher the tissue specificity. S ≥ 1 was considered tissue specific, while S < 0.3 was considered to be widely expressed. When 0.3 ≤ S < 1, the tissue-specific was at the intermediate level. In this study, the tissue specific genes with S ≥ 1 were screened for display, and the heat map was drawn with FPKM value.

Statistical Analysis

Statistical analysis of the demographic and the clinical data was performed using SPSS version 22.0 (IBM Corp., Armonk, NY, United States). When the data were normally distributed, we expressed it as a mean ± SD; if not, we expressed it as the median and the interquartile range (IQR). P-value of < 0.05 was considered statistically significant.

For analysis of the sequenced data difference, the absolute value of diff. log 2FC ≥ 1, and a P < 0.05 were set as the threshold for DEGs. For the analysis of m6A data difference, the input RNA and m6A RNA reads were aligned to genome reference sequences using HISAT2 software. ExomePeak was used for peak calling and differential RNA methylation analysis. HOMER (version 4.10.4) was used for motif search. FDR ≤ 0.05, P < 0.05, and diff.log 2FC > 0 was considered as hypermethylation, while diff. log2FC < 0 was considered as demethylation. The absolute value of diff. log 2FC ≥ 1, P < 0.05 were set as the threshold for differential m6A methylated RNAs. The global distribution of hypermethylated and hypomethylated m6A peaks was displayed by Circos plot using CIRCOS.

Results

Demographic and Clinical Characteristics of Enrolled Patients

A total of 1015 CSF specimens were collected, in which, 573 cases had no detectable pathogens, 271 cases were non-EV infections, 87 cases had no nucleic acid detection, and 84 cases with EV infection. Among the 84 cases with confirmed EV, 48 cases were coinfected with other pathogens and 36 cases were simple EV infection. In those with simple EV infection, patients without neurologic symptoms (7 cases) were all under 4 years old, so we choose the remaining 13 cases with neurologic symptoms under 4 years old for comparison in order to reduce bias (Figure 1A).

FIGURE 1

Figure 1. (A) Flow chart of patient enrollment. (B) The flowchart of meRIP- Sequencing. (C) FPKM distribution in control group and NS group. Fragments Per Kilobase of transcript sequence per Millions base pairs. (D) m6A RNA Methylation Peaks Distribution of Host Cell RNA Transcripts. 1. m6A peaks distribution on mRNA of EV infected patients without neurological symptoms (control group). 2. m6A peaks distribution on mRNA of EV infected patients with neurological symptoms (NS group). 3. Pie chart of m6A peaks distribution in control group. 4. Pie chart of m6A peaks distribution in NS group. (E) Circos plot showing the distribution of hypermethylated and hypomethylated m6A peaks in the human transcriptome of neurological symptoms group compared with control group. Red mark represents the distribution of m6A peaks in NS group. Green mark represents the distribution of m6A peaks in control group. The line in red or green represents an enrichment fold of m6A methylation peaks in different chromosomes with fold change at Y-axis. Reference genome: GRCh38/hg38.

Demographic and clinical information of the enrolled patients was shown in Table 1 and Supplementary Table 1. There were no statistically significant differences in any of the demographic characteristics between the EV infected patients with neurological symptoms (NS group) and those without (control group). The main neurological symptoms of EV infected patients were convulsions, vomiting, headache, dizziness, fatigue/somnolence, and acute disorder of consciousness (ADOC), of which convulsions was the most common. Serum levels of aspartate aminotransferase (AST) were significantly higher in patients with neurological symptoms (P = 0.044), and those with a higher concentration of microprotein in CSF samples (P = 0.035). The lymphocytes count, immunoglobulin E levels, and glucose concentration were lower in the NS group, but with no statistical significance. In the NS group, eight patients underwent neuroimaging, five of whom had abnormal findings, one had leptomeningeal enhancement and subdural effusion, two patients had abnormal signals located in posterior horn of ventricle or cerebral hemispheres and midbrain, respectively, one had swollen cerebral cortex and basal ganglia, and the last one had enlarged sulci. But generally, most patients had a favorable prognosis and a recovery rate of 69.2%.

Discharge summaries found that, in the control group, 3 cases had hand-foot-and-mouth disease, and 4 cases with herpangina. In the NS group, 4 cases had hand-foot-and-mouth disease or herpangina with febrile convulsions, 5 cases had enterovirus encephalitis, 2 cases had epilepsy with enterovirus infection, one case had enterovirus infection co-infected with bacteria which caused purulent meningitis, and a case with enterovirus infection accompany with Kawasaki disease. All cases in the control group recovered fully, while in the NS group, only 9 cases had a full recovery, 3 cases had clinical improvement and 1 case of death.

General Features of m6A RNA Methylation

The flowchart of meRIP-sequencing was shown in Figure 1B. The expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) showed no difference in total gene expression between the two groups (Figure 1C). A total of 957 m6A methylated genes were obtained from the two groups. The distribution of m6A RNA methylation in EV infected children was illustrated in Figure 1D. The m6A peaks of all enrolled patients mainly occurred in the four non-overlapping transcript segments, including coding sequence (CDS), 3′ untranslated regions (3′UTR), stop codon and start codon. Metagene analysis showed that NS patients had an increased level of m6A in the CDS regions of the transcriptome but correspondingly decreased in the 3′UTRs and the stop codon. At the same time, we compared the enrichment folds of m6A methylation peaks in different chromosomes between the NS group and the control group (Figure 1E). No difference in the methylation enrichment on Y chromosome between control group and NS group was observed, while the m6A methylation distribution on other chromosomes was more common in control group than that of the NS group (Figure 1E). The sequence logo of m6A motif was shown in Supplementary Figure 1.

Conjoint Analysis of MeRIP-seq and RNA-seq Data

We also further investigated changes in the expression pattern of the m6A methylated genes in EV infected patients with neurological symptoms and the control group. A total of 3472 DEGs were identified. The top 20 DEGs according to the absolute log2FC were summarized in Table 2 and the all DEGs were presented in Supplementary Table 2. The DEGs with m6A methylation was indicated in Table 3. Differentially methylated m6A peaks with synchronously differential expression were summarized in Table 4 and Supplementary Figure 3. The m6A methylated genes with synchronously differential expression were identified between the two groups, including PLEKHA4 (Supplementary Figure 2).

TABLE 2

Table 2. Differentially expressed input RNA genes (top 20).

TABLE 3

Table 3. Top 20 of input differentially expressed Genes with m6A methylation in NS group.

TABLE 4

Table 4. Differentially methylated m6A Peaks with synchronously differential expression in NS group.

Additionally, 31 human tissues, covering the expression levels of 60498 genes were downloaded from the GTEX database among which, the tissue specificity scores of m6A abundant genes in the NS group were shown in Supplementary Figure 3. The heatmap exhibits that PLEKHA4 was a tissue specific gene in nerve and MAP2 was a tissue specific gene in brain (Supplementary Figure 3).

Gene Ontology Annotation and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis

We performed the GO and KEGG pathway analysis to identify the functions of these genes with significant changes in m6A modifications using 15 significant enrichment GO terms in biological process (BP), molecular function (MF) and cellular components (CC) as shown in Figure 2. The GO molecular functional enrichment results showed that the genes with m6A modifications were mainly related to poly(A) RNA binding in all EV infected children especially in those with neurological symptoms (Figures 2A,B). We also annotated the m6A into biological process, and found that, the NS group were mainly involved in the respiratory electron transposition, cellular metabolic process and translational initiation (Figure 2B). Lastly, the cellular component analysis found that, the genes were enriched in the respiratory chain, nuclear speck, and nucleoplasm in NS group, while in the control group, the GO functional classification analysis exhibited distinct results, as shown in Figures 2A,B. Additionally, we performed GO annotation analysis for the genes with significant changes in both the m6A modification and the mRNA levels (Figure 3). Compare with control group (Figure 3A), the GO terms result showed that those genes of NS group were most related to respiratory electron transposition, cellular metabolic and oxidation-reduction in biological process. In molecular function, those genes were mainly involved in oxidoreductase activity, NADH dehydrogenation and site-specificity. We also annotated those genes into cellular component, the results showed that they were enriched in respiratory chain, mitochondrial membrane and mitochondrion (Figure 3A). Furthermore, according to the expression level on mRNA, the differently expressed m6A methylated genes were divided into up-regulation and down-regulation genes, GO annotation of those genes were also performed, respectively, as presented in Figures 3B,C. The result of up-regulation genes was similar with that of genes with significant changes in both m6A modification and mRNA levels (Figures 3A,B), while the down-regulation genes showed distinct results (Figure 3C).

FIGURE 2

Figure 2. GO Enrichment Analysis of m6A RNA Methylation Genes. (A) The top 15 significant enrichment GO terms of genes with m6A modification in control group. (B) The top 15 significant enrichment GO terms of genes with m6A modification in neurological symptoms group. NS: neurological symptoms. BP: biological process. MF: molecular function. CC: cellular component.

FIGURE 3

Figure 3. GO Enrichment Analysis of genes with significant changes in both m6A modification and mRNA levels (A) The top 15 significant enrichment GO terms of genes with significant changes in both the m6A modifications and mRNA levels in NS group (NS vs. control group). (B) The top 15 significant enrichment GO terms of differently m6A methylated genes with significant upregulation of mRNA levels in NS group (NS vs. control group). (C) The top 15 significant enrichment GO terms of differently m6A methylated genes with significant downregulation of mRNA levels in NS group (NS vs. control group). NS: neurological symptoms, Red: significant difference, blue: no significant difference. BP: biological process. MF: molecular function. CC: cellular component.

The KEGG pathway analysis also revealed various top 20 pathway enrichments as shown in Figure 4. In the control group, the protein processing in the endoplasmic reticulum, seienoamino acid metabolism, and the olfactory transduction pathway were significantly enriched (Figure 4A), while in the NS group, the enriched pathways were tight junction, ascorbate and aldarate metabolism, pancreatic secretion, and cGMP-PKG signaling pathways (Figure 4B). Furthermore, GO terms like respiratory electron transport chain, cellular metabolic process and oxidation-reduction process were all enriched as revealed by the analysis of the differentially methylated genes in NS group (Figure 4C). In addition, the KEGG pathway analysis based on the differentially methylated genes between the two groups unveiled that oxidative phosphorylation and Parkinson’s disease pathway were also enriched (Figure 4D). Additionally, the enriched pathways of all genes with significant changes in both the m6A modifications and mRNA levels in NS group were shown in Figures 4E–G.

FIGURE 4

Figure 4. KEGG Pathway Analysis of m6A RNA Methylation Genes. (A) The top 20 enriched pathway of genes with m6A modification in control group. (B) The top 20 enriched pathway of genes with m6A modification in NS group. (C) The top 20 enriched GO terms of the differentially methylated genes in NS group compared with those in control group. (D) The top 9 enriched pathways of the differentially methylated genes NS group compared with those in control group. (E) The enriched pathways of all genes with significant changes in both the m6A modifications and mRNA levels in NS group (NS vs. Control group). (F) The enriched pathways of differently m6A methylated genes with significant upregulation of mRNA levels in NS group (NS vs. Control group). (G) The enriched pathways of differently m6A methylated genes with significant downregulation of mRNA levels in NS group (NS vs. Control group).

Discussion

Globally, EV infection has caused several notable outbreaks in pediatric population due its virulence and mostly in Asian countries. Despite the vast majority of EV infections are mild and self-limiting, some children will develop central nervous complications with serious sequelae and even death. Vaccines developed in China confers protection (Mao et al., 2016), but it is not effective in symptomatic patients. The aim of our study was to explore the underlying mechanism of EV related neurological damage by using RNA-sequencing technology to profile the transcriptome alterations of m6A RNA methylation in CSF samples of simple EV infected children with and without neurologic symptoms. We believe this study will contribute to highlighting the pathogenesis of EV infections and will further provide countermeasures to mitigate disease chains of transmission.

In our cohort, the most prevalent neurological symptoms were convulsions, vomiting, headaches, somnolence, and acute disorder of consciousness (ADOC). We previously found that somnolence in EV patients was correlated with high mortality (Yang et al., 2018). Similarly, previous studies have also shown that symptoms like convulsion, dyspnea, cyanosis, and vomiting were also associated with increased high risk of death from severe HFMD (Long et al., 2016). CVA6-associated severe cases were characterized by high fever with intermittent twitching, while EV71-associated severe cases were characterized by poor mental condition, loss of pupillary reflex, and vomiting (Li et al., 2016; Yang et al., 2017, 2018). However, all our enrolled patients in our study had a favorable prognosis, possibly due to our small sample size. It has been reported that leukocytosis and increased CSF protein is more common in encephalitis than in febrile convulsion cases of HFMD (Xu et al., 2019). In our study, we found a similar pattern of elevated microprotein levels in EV infected patients with neurological symptoms. Contrary to a previously published study on EV-A71 infection associated pulmonary edema (Crabol et al., 2017), there was a low glucose levels and high AST level in the NS group. There is still no consensus on the risk factors associated with severe Infection (Ooi et al., 2009; Crabol et al., 2017; Qin et al., 2019; Xu et al., 2019), hence further well conducted, large sample size studies are needed.

Previous studies have also described the changing patterns of m6A deposition in response to heat shock and Flaviviridae infections, like dengue virus (DENV), Zika virus (ZIKV), West Nile virus (WNV), and hepatitis C virus (HCV) (Zhou et al., 2015; Lichinchi et al., 2016b; Gokhale et al., 2020). Our data indicated that there was a differentially expressed pattern of m6A RNA methylation in EV infected children with neurologic symptoms. We observed changes in the distribution of m6A, with increased peaks in the CDS regions but with a decreased trend in the 3′UTRs and stop codon regions of the children with neurological symptoms. This indicates that EV infection might affect gene translation, gene splicing and mRNA stability by changing the deposition of m6A, leading to the development of severe neurological symptoms and sequalae. Nine m6A methylated genes with synchronously differential expression were identified by conjoint analysis of MeRIP-seq and RNA-seq data. Among which, PLEKHA4, a tissue-specific gene in the nerve, was upregulated and hypermethylated in patients with neurological symptoms. It’s known that PLEKHA4 is involved in autism spectrum disorder (Hashimoto et al., 2016). Previous study demonstrated that PLEKHA4 was a signaling strength modulator in Wnt signaling pathway, which controls key cell fate decisions in the development of multicellular eukaryotes (Shami Shah et al., 2019). One study about Melanoma found that, PLEKHA4 was required for melanoma proliferation and survival, PLEKHA4 knockdown could attenuate tumor growth and enhance the effect of clinically used inhibitor encorafenib (Shami Shah et al., 2021). We speculated that PLEKHA4 and related Wnt signaling may also involve with EV-induced neuron damage. Whether the targeting therapy of PLEKHA4 may help reduce the neuron damage in EV infection needs further study. We infer that those other m6A methylated genes with synchronously differential expression might also participate in EV-induced neuron damage.

We carried out the GO and KEEG pathways analysis of m6A methylated RNAs which showed that, the GO terms including respiratory electron transport chain, cellular metabolic process, and oxidation-reduction process pathways were all enriched in the EV infected patients with neuropsychiatric symptoms, citing a possible correlation with the underlying mechanism of central nervous system damage in EV infections. Recently, an experiment of SH-SY5Y cells with EV71 infection also found a dysregulated transcriptomes profile, which revealed that the enriched GO term “Nervous system development” and enriched pathway “CCKR signaling map” might be important contributors to EV71-induced neuropathological mechanism (Hu et al., 2020). When we compared the difference in signaling pathways of m6A peaks between the children with and without neurological symptoms, we found that, the oxidative phosphorylation and Parkinson’s disease pathway were both enriched, a finding that is consistent with previous study that concluded the association between acute parkinsonism and Coxsackie virus (Jang et al., 2009). Also, another study found the presence of virus-like particles and enterovirus antigen in the brainstem neurons in patients with Parkinson’s disease (Dourmashkin et al., 2018). But the association of EV with aged-onset idiopathic Parkinson’s disease has never been established. One recent study shows that the acute course of EV71-infected patients with neurological symptoms, such as parkinsonian-like “limb tremor” did not persist as a neurological sequela (Yang et al., 2017). Underlying mechanisms need to be elucidated in further studies.

Limitations of our study include the lack of distinct identification of EV serotype among the 100 serotypes of human EV (Khetsuriani et al., 2006). Secondly, this study involves a small number of patients from the same hospital and the results may lack general applicability and extrapolation.

In summary, we described the m6A RNA methylation transcriptomic patterns of EV by high-throughput RNA-sequencing analysis in EV infected children under four years old. We found a variable m6A distribution, pathways, and expression in children with neurological symptoms, which were mostly associated with respiratory chains and cell metabolism. These specific genes with aberrant expression might play an important role in the process and the pathogenesis of CNS lesions caused by EV. Further studies are required to explore the more succinct underlying mechanisms of EV infection.

Data Availability Statement

The data reported in this paper have been deposited in the OMIX, China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (https://ngdc.cncb.ac.cn/omix: accession no.OMIX716).

Ethics Statement

The studies involving human participants were reviewed and approved by Ethics Committee of the Guangzhou Women and Children’s Medical Center. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author Contributions

DZ, YS, and DH performed and analyzed the experiments and wrote the manuscript. SL and GL carried out the data collection and data analysis. PL and SY conceived and coordinated the study and designed and revised the manuscript. All authors reviewed the results and approved the final version of the manuscript.

Funding

This work was supported by Guangdong Natural Science Foundation (No. 2020A1515010014); National Natural Science Foundation of China (No. 81801206); and Innovative Project of Children’s Research Institute, Guangzhou Women and Children’s Medical Center, China (Nos. Pre-NSFC-2019-002 and Pre-NSFC-2019-007).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Thanks all the enrolled patients and the Clinical Biological Resource Bank of Guangzhou Women and Children’s Medical Center for providing the clinical samples used in this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2021.791544/full#supplementary-material

Footnotes

1. ^ http://homer.ucsd.edu/homer/ngs/peakMotifs.html

References

Ai, J., Xie, Z., Liu, G., Chen, Z., Yang, Y., Li, Y., et al. (2017). Etiology and prognosis of acute viral encephalitis and meningitis in Chinese children: a multicentre prospective study. BMC Infect. Dis. 17:494. doi: 10.1186/s12879-017-2572-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Crabol, Y., Pean, P., Mey, C., Duong, V., Richner, B., Laurent, D., et al. (2017). A prospective, comparative study of severe neurological and uncomplicated hand, foot and mouth forms of paediatric enterovirus 71 infections. Int. J. Infect. Dis. 59, 69–76. doi: 10.1016/j.ijid.2017.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Dourmashkin, R. R., McCall, S. A., Dourmashkin, N., and Hannah, M. J. (2018). Virus-like particles and enterovirus antigen found in the brainstem neurons of Parkinson’s disease. F1000Res. 7:302. doi: 10.12688/f1000research.13626.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, K., Zhang, L., Lee, T., and Sun, T. (2019). m(6)A RNA methylation controls neural development and is involved in human diseases. Mol. Neurobiol. 56, 1596–1606.

Google Scholar

Frayling, T. M., Timpson, N. J., Weedon, M. N., Zeggini, E., Freathy, R. M., Lindgren, C. M., et al. (2007). A common variant in the FTO gene is associated with body mass index and predisposes to childhood and adult obesity. Science 316, 889–894. doi: 10.1126/science.1141634

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstberger, S., Hafner, M., and Tuschl, T. (2014). A census of human RNA-binding proteins. Nat. Rev. Genet. 15, 829–845. doi: 10.1038/nrg3813

PubMed Abstract | CrossRef Full Text | Google Scholar

Gokhale, N. S., McIntyre, A. B. R., Mattocks, M. D., Holley, C. L., Lazear, H. M., Mason, C. E., et al. (2020). Altered m6A modification of specific cellular transcripts affects flaviviridae infection. Mol. Cell 77, 542.e–555.e. doi: 10.1016/j.molcel.2019.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez, G., Carr, M. J., Kobayashi, M., Hanaoka, N., and Fujimoto, T. (2019). Enterovirus-associated hand-foot and mouth disease and neurological complications in japan and the rest of the world. Int. J. Mol. Sci. 20:5201. doi: 10.3390/ijms20205201

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashimoto, R., Nakazawa, T., Tsurusaki, Y., Yasuda, Y., Nagayasu, K., Matsumura, K., et al. (2016). Whole-exome sequencing and neurite outgrowth analysis in autism spectrum disorder. J. Hum. Genet. 61, 199–206. doi: 10.1038/jhg.2015.141

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho, A. J., Stein, J. L., Hua, X., Lee, S., Hibar, D. P., Leow, A. D., et al. (2010). A commonly carried allele of the obesity-related FTO gene is associated with reduced brain volume in the healthy elderly. Proc. Natl. Acad. Sci. U.S.A. 107, 8404–8409. doi: 10.1073/pnas.0910878107

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Y., Xu, Y., Huang, Z., Deng, Z., Fan, J., Yang, R., et al. (2020). Transcriptome sequencing analysis of SH-SY5Y cells infected with EV71 reveals the potential neuropathic mechanisms. Virus Res. 282:197945. doi: 10.1016/j.virusres.2020.197945

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Liao, Q., Ooi, M. H., Cowling, B. J., Chang, Z., Wu, P., et al. (2018). Epidemiology of recurrent hand, foot and mouth disease, China, 2008-2015. Emerg. Infect. Dis. 24, 432–442.

Google Scholar

Jang, H., Boltz, D. A., Webster, R. G., and Smeyne, R. J. (2009). Viral parkinsonism. Biochim. Biophys. Acta 1792, 714–721.

Google Scholar

Jin, D. I., Lee, S. W., Han, M. E., Kim, H.-J., Seo, S.-A., Hur, G.-Y., et al. (2012). Expression and roles of Wilms’ tumor 1-associating protein in glioblastoma. Cancer Sci. 103, 2102–2109. doi: 10.1111/cas.12022

PubMed Abstract | CrossRef Full Text | Google Scholar

Kennedy E, M., Bogerd H, P., Kornepati A, V. R., Kang, D., Ghoshal, D., Marshall, J. B., et al. (2016). Posttranscriptional m 6 A editing of HIV-1 mRNAs enhances viral gene expression. Cell Host Microbe 19, 675–685. doi: 10.1016/j.chom.2016.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Khetsuriani, N., Lamonte, A., Oberste, M. S., and Pallansch, M. (2006). Neonatal enterovirus infections reported to the national enterovirus surveillance system in the United States, 1983-2003. Pediatr. Infect. Dis. J. 25, 889–893. doi: 10.1097/01.inf.0000237798.07462.32

CrossRef Full Text | Google Scholar

Lee, K. Y., Lee, Y. J., Kim, T. H., Cheon, D. S., and Nam, S. O. (2014). Clinico-radiological spectrum in enterovirus 71 infection involving the central nervous system in children. J. Clin. Neurosci. 21, 416–420. doi: 10.1016/j.jocn.2013.04.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Sun, Y., Du, Y., Yan, Y., Huo, D., Liu, Y., et al. (2016). Characterization of coxsackievirus A6- and enterovirus 71-associated hand foot and mouth disease in Beijing, China, from 2013 to 2015. Front. Microbiol. 7:391. doi: 10.3389/fmicb.2016.00391

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z. H., Li, C. M., Ling, P., Shen, F.-H., Chen, S.-H., Liu, C.-C., et al. (2008). Ribavirin reduces mortality in enterovirus 71-infected mice by decreasing viral replication. J. Infect. Dis. 197, 854–857. doi: 10.1086/527326

PubMed Abstract | CrossRef Full Text | Google Scholar

Lichinchi, G., Gao, S., Saletore, Y., Gonzalez, G. M., Bansal, V., Wang, Y., et al. (2016a). Dynamics of the human and viral m(6)A RNA methylomes during HIV-1 infection of T cells. Nat. Microbiol. 1:16011. doi: 10.1038/nmicrobiol.2016.11

PubMed Abstract | CrossRef Full Text | Google Scholar

Lichinchi, G., Zhao, B. S., Wu, Y., Lu, Z., Qin, Y., He, C., et al. (2016b). Dynamics of human and viral RNA methylation during zika virus infection. Cell Host Microbe 20, 666–673. doi: 10.1016/j.chom.2016.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, S., Choe, J., Du, P., Triboulet, R., and Richard, G. I. (2016). The m 6 A Methyltransferase METTL3 promotes translation in human cancer cells. Mol. Cell 62, 335–345. doi: 10.1016/j.molcel.2016.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, L., Gao, L. D., Hu, S. X., Luo, K.-W., Chen, Z.-H., Ronsmans, C., et al. (2016). Risk factors for death in children with severe hand, foot, and mouth disease in Hunan. China. Infect Dis 48, 744–748. doi: 10.1080/23744235.2016.1185801

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, Q.-Y., Wang, Y., Bian, L., Xu, M., and Liang, Z. (2016). EV71 vaccine, a new tool to control outbreaks of hand, foot and mouth disease (HFMD). Expert Rev. Vaccines 5, 599–606. doi: 10.1586/14760584.2016.1138862

PubMed Abstract | CrossRef Full Text | Google Scholar

Ooi, M. H., Wong, S. C., Mohan, A., Podin, Y., Perera, D., and Clear, D. (2009). Identification and validation of clinical predictors for the risk of neurological involvement in children with hand, foot, and mouth disease in Sarawak. BMC Infect. Dis. 9:3. doi: 10.1186/1471-2334-9-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, L., Dang, D., Wang, X., Zhang, R., Feng, H., Ren, J., et al. (2019). Identification of immune and metabolic predictors of severe hand-foot-mouth disease. PLoS One 14:e0216993. doi: 10.1371/journal.pone.0216993

PubMed Abstract | CrossRef Full Text | Google Scholar

Reddy, S. M., Sadim, M., Helenowski, I. B., Li, J., Yi, N., Kaklamani, V. G., et al. (2013). Clinical and genetic predictors of weight gain in patients diagnosed with breast cancer. Br. J. Cancer 109, 872–881. doi: 10.1038/bjc.2013.441

PubMed Abstract | CrossRef Full Text | Google Scholar

Shami Shah, A., Batrouni, A. G., Kim, D., Punyala, A., Cao, W., Han, C., et al. (2019). PLEKHA4/kramer attenuates dishevelled ubiquitination to modulate WNT and planar cell polarity signaling. Cell Rep. 27, 2157–2170 e8. doi: 10.1016/j.celrep.2019.04.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Shami Shah, A., Cao, X., White, A. C., and Baskin, J. M. (2021). PLEKHA4 promotes Wnt/beta-catenin signaling-mediated G1-S transition and proliferation in melanoma. Cancer Res. 81, 2029–2043. doi: 10.1158/0008-5472.CAN-20-2584

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Hu, T., Sun, D., Carr, M. J., Xing, W., Li, S., et al. (2017). Epidemiological characteristics of hand, foot, and mouth disease in Shandong, China, 2009-2016. Sci Rep. 7:8900. doi: 10.1038/s41598-017-09196-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, Y., Li, S., Cai, C., Liu, J., Wang, Y., Jiang, Y., et al. (2019). Characterization of inflammatory cytokine profiles in cerebrospinal fluid of hand, foot, and mouth disease children with enterovirus 71-related encephalitis in Hangzhou, Zhejiang, China. Medicine 98:e18464. doi: 10.1097/MD.0000000000018464

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S. D., Li, P. Q., Huang, Y. G., Li, W., Ma, L. Z., Wu, L., et al. (2018). Factors associated with fatal outcome of children with enterovirus A71 infection: a case series. Epidemiol. Infect. 146, 788–798. doi: 10.1017/S0950268818000468

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S. D., Li, P. Q., Li, Y. M., Li, W., Lai, W.-Y., Zhu, C.-P., et al. (2017). Clinical manifestations of severe enterovirus 71 infection and early assessment in a Southern China population. BMC Infect. Dis. 17:153. doi: 10.1186/s12879-017-2228-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoon, K. J., Ringeling, F. R., Vissers, C., Jacob, F., Pokrass, M., Jimenez-Cyrus, D., et al. (2017). Temporal control of mammalian cortical neurogenesis by m(6)A methylation. Cell 171:e17. doi: 10.1016/j.cell.2017.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Wan, J., Gao, X., Zhang, X., Jaffrey, S. R., and Qian, S. B. (2015). Dynamic m(6)A mRNA methylation directs translational control of heat shock response. Nature 526, 591–594. doi: 10.1038/nature15377

PubMed Abstract | CrossRef Full Text | Google Scholar