ORIGINAL RESEARCH article

Front. Neurosci., 07 December 2021

Sec. Neurodegeneration

Volume 15 - 2021 | https://doi.org/10.3389/fnins.2021.791544

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

  • 1. Department of Pediatric Emergency, Guangzhou Women and Children’s Medical Center, Guangzhou Medical University, Guangzhou, China

  • 2. Department of Pediatric Neurology, Guangzhou Women and Children’s Medical Center, Guangzhou Medical University, Guangzhou, China

Abstract

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

GroupsAll (n = 20)Control group (n = 7)NS group (n = 13)P-value
Sex female, n (%)9 (45)2 (28.6)7 (53.8)0.29
Age, months (IQR)18 (2, 26)21 (13, 24)13 (2, 30)0.66
Weight, mean ± SD, kg9.52 ± 3.2510.01 ± 2.289.25 ± 3.730.63
Onset, mean ± SD, days2.80 ± 4.251.8 ± 1.393.3 ± 5.20.45
Fever, n (%)17 (85.00)5 (71.43)12 (92.31)0.22
Convulsions, n (%)11 (55)011 (84.62)0
Vomiting, n (%)2 (10)02 (15.38)0.29
Headache, n (%)1 (5)01 (7.69)0.46
Fatigue/somnolence, n (%)1 (5)01 (7.69)0.46
Acute disorder of consciousness, n (%)1 (5)01 (7.69)0.46
Aspartate aminotransferase, mean ± SD, (U/L)44.68 ± 25.2932.83 ± 3.5451.15 ± 29.18*0.04
CSF microprotein, mean ± SD, (g/L)0.44 ± 0.360.25 ± 0.150.54 ± 0.41*0.04
Lymphocytes, mean ± SD (10^9/L)3.35 ± 1.774.24 ± 2.442.88 ± 1.120.20
Glucose, mean ± SD (mmol/L)5.88 ± 0.666.07 ± 0.775.79 ± 0.620.42
Immunoglobulin E, mean ± SD (IU/ML)72.67 ± 140.43137.33 ± 229.3140.33 ± 56.180.35
Neuroimaging abnormal findings, n (%)5/13(38.46)05/8(62.50)0.07
EEG abnormal findings, n (%)5/9(55.56)1/3(33.33)4/6(66.67)0.63

Demographic and clinical characteristics of enrolled patients.

* The difference between the EV infected patients with and without neurological symptoms were statistically significant, P<0.05.

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

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

Gene nameGene IDRegulationGene biotypeLog2FC#FDRP-value
MT-RNR2ENSG00000210082UpMt-rRNA1.79335500
MT-ATP6ENSG00000198899DownProtein coding−4.576300
RPS23ENSG00000186468DownProtein coding−5.7675900
RPS15ENSG00000115268DownProtein coding−9.855200
MT-CO2ENSG00000198712DownProtein coding−9.6750900
EMP3ENSG00000142227DownProtein coding−6.5623800
HSP90AB1ENSG00000096384DownProtein coding−4.7373900
HSP90AA1ENSG00000080824DownProtein coding−3.4052500
MAP4ENSG00000047849DownProtein coding−10.94300
CASC3ENSG00000108349DownProtein coding−4.7076100
ENO1ENSG00000074800DownProtein coding−10.136200
BCL2L1ENSG00000171552DownProtein coding−10.888500
PTK2ENSG00000169398DownProtein coding−6.6535700
MT-ND2ENSG00000198763UpProtein coding4.86554400
EIF3BENSG00000106263DownProtein coding−5.379441.41E-2991.84E-301
STIP1ENSG00000168439DownProtein coding−6.38861.33E-2991.85E-301
AMOTL1ENSG00000166025DownProtein coding−9.794765.09E-2757.53E-277
FILIP1LENSG00000168386DownProtein coding−13.09642.56E-2574.01E-259
JUNBENSG00000171223DownProtein coding−13.06413.23E-2535.34E-255
RAB10ENSG00000084733DownProtein coding−11.68099.90E-2511.72E-252

Differentially expressed input RNA genes (top 20).

#The log2FC is the log2 conversion value of the difference fold change.

TABLE 3

Gene nameLocusRegulationGene biotypeFold enrichmentlog2FCP-value
MT-CO2MT: 7644–7972DownProtein coding55−9.6750895770
MT-ND6MT: 14148–14382UpProtein coding944.1042193972.15E-51
SETChr9: 128693653–128693950UpProtein coding82.72.1042193970.014500156
ERICH3Chr1: 74572124-74572365UpProtein coding61.12.9393584953.68E-199
SFXN3Chr10: 101040198–101040378UpProtein coding44.23.9437547258.35E-06
PTMAChr2: 231711309–231711549UpProtein coding3144.0060390031.82E-47
LRRC75A-AS1Chr17: 16441669–16441849UpProcessed transcript39.43.8104881942.68E-09
MT-ATP6MT: 8526–8793DownProtein coding178−4.576302330
MT-RNR1MT: 707–946UpMt rRNA42.51.9588616183.59E-207
PLEKHA4Chr19: 48868452–48868632UpProtein coding19.36.8506453211.44E-37
PRRC2CChr1: 171532482–171532872DownProtein coding1080−3.3661005376.08E-05
DDX21Chr10: 68959968–68960208UpProtein coding1426.1416941039.22E-92
MALAT1Chr11: 65499083–65499293DownLincRNA86−1.9543307523.51E-247
SAMHD1Chr20: 36890438–36890589UpProtein coding6.282.0362009743.75E-11
MT-ND4MT: 10998–11269DownProtein coding92.9−2.9213156958.59E-07
LMNAChr1: 156139788–156139969UpProtein coding1175.7266566041.22E-18
MT-ND1MT: 3783–4172UpProtein coding1171.5821169211.40E-09
HSP90AB1Chr6: 44252017–44252226DownProtein coding57.8−4.7373878340
MT-RNR2MT: 2539–3229UpMt rRNA5541.7933549330
RPL18Chr19: 48617356–48617536UpProtein coding457.9983989412.14E-74

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

Fold Enrichment: Fold Enrichment for m6A methylation; Log2FC: Fold change for gene expression.

TABLE 4

Gene nameLocusRegulationGene biotypeFold enrichmentLog2FC#P-value
SET*Chr9:128693801–128693950UpProtein coding82.72.1042193970.0145002
PTMAChr2:231711309–231711549DownProtein coding3144.0060390031.82E-47
MT-ND6*MT: 14148–14382UpProtein coding944.1042193972.15E-51
MT-CO2MT: 7644–7972UpProtein coding55−9.6750895770
MT-ND4MT: 10998–11269UpProtein coding92.9−2.9213156958.59E-07
MT-ND1*MT: 3783–4172UpProtein coding1171.5821169211.40E-09
MT-ATP6MT: 8526-8793UpProtein coding178−4.576302330
MT-RNR2MT: 2539–3229DownMt rRNA5541.7933549330
MT-RNR1MT: 707–946DownMt rRNA42.51.9588616183.59E-207

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

*Methylated m6A and mRNA were up-regulated at the same time. #The log2FC is the log2 conversion value of the difference foldchange.

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 3

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

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.

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.

Statements

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

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.

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.

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

References

  • 1

    AiJ.XieZ.LiuG.ChenZ.YangY.LiY.et al (2017). Etiology and prognosis of acute viral encephalitis and meningitis in Chinese children: a multicentre prospective study.BMC Infect. Dis.17:494. 10.1186/s12879-017-2572-9

  • 2

    CrabolY.PeanP.MeyC.DuongV.RichnerB.LaurentD.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.596976. 10.1016/j.ijid.2017.04.005

  • 3

    DourmashkinR. R.McCallS. A.DourmashkinN.HannahM. J. (2018). Virus-like particles and enterovirus antigen found in the brainstem neurons of Parkinson’s disease.F1000Res.7:302. 10.12688/f1000research.13626.2

  • 4

    DuK.ZhangL.LeeT.SunT. (2019). m(6)A RNA methylation controls neural development and is involved in human diseases.Mol. Neurobiol.5615961606.

  • 5

    FraylingT. M.TimpsonN. J.WeedonM. N.ZegginiE.FreathyR. M.LindgrenC. M.et al (2007). A common variant in the FTO gene is associated with body mass index and predisposes to childhood and adult obesity.Science316889894. 10.1126/science.1141634

  • 6

    GerstbergerS.HafnerM.TuschlT. (2014). A census of human RNA-binding proteins.Nat. Rev. Genet.15829845. 10.1038/nrg3813

  • 7

    GokhaleN. S.McIntyreA. B. R.MattocksM. D.HolleyC. L.LazearH. M.MasonC. E.et al (2020). Altered m6A modification of specific cellular transcripts affects flaviviridae infection.Mol. Cell77542.e555.e. 10.1016/j.molcel.2019.11.007

  • 8

    GonzalezG.CarrM. J.KobayashiM.HanaokaN.FujimotoT. (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. 10.3390/ijms20205201

  • 9

    HashimotoR.NakazawaT.TsurusakiY.YasudaY.NagayasuK.MatsumuraK.et al (2016). Whole-exome sequencing and neurite outgrowth analysis in autism spectrum disorder.J. Hum. Genet.61199206. 10.1038/jhg.2015.141

  • 10

    HoA. J.SteinJ. L.HuaX.LeeS.HibarD. P.LeowA. 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.10784048409. 10.1073/pnas.0910878107

  • 11

    HuY.XuY.HuangZ.DengZ.FanJ.YangR.et al (2020). Transcriptome sequencing analysis of SH-SY5Y cells infected with EV71 reveals the potential neuropathic mechanisms.Virus Res.282:197945. 10.1016/j.virusres.2020.197945

  • 12

    HuangJ.LiaoQ.OoiM. H.CowlingB. J.ChangZ.WuP.et al (2018). Epidemiology of recurrent hand, foot and mouth disease, China, 2008-2015.Emerg. Infect. Dis.24432442.

  • 13

    JangH.BoltzD. A.WebsterR. G.SmeyneR. J. (2009). Viral parkinsonism.Biochim. Biophys. Acta1792714721.

  • 14

    JinD. I.LeeS. W.HanM. E.KimH.-J.SeoS.-A.HurG.-Y.et al (2012). Expression and roles of Wilms’ tumor 1-associating protein in glioblastoma.Cancer Sci.10321022109. 10.1111/cas.12022

  • 15

    Kennedy EM.Bogerd HP.Kornepati AV. R.KangD.GhoshalD.MarshallJ. B.et al (2016). Posttranscriptional m 6 A editing of HIV-1 mRNAs enhances viral gene expression.Cell Host Microbe19675685. 10.1016/j.chom.2016.04.002

  • 16

    KhetsurianiN.LamonteA.ObersteM. S.PallanschM. (2006). Neonatal enterovirus infections reported to the national enterovirus surveillance system in the United States, 1983-2003.Pediatr. Infect. Dis. J.25889893. 10.1097/01.inf.0000237798.07462.32

  • 17

    LeeK. Y.LeeY. J.KimT. H.CheonD. S.NamS. O. (2014). Clinico-radiological spectrum in enterovirus 71 infection involving the central nervous system in children.J. Clin. Neurosci.21416420. 10.1016/j.jocn.2013.04.032

  • 18

    LiJ.SunY.DuY.YanY.HuoD.LiuY.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. 10.3389/fmicb.2016.00391

  • 19

    LiZ. H.LiC. M.LingP.ShenF.-H.ChenS.-H.LiuC.-C.et al (2008). Ribavirin reduces mortality in enterovirus 71-infected mice by decreasing viral replication.J. Infect. Dis.197854857. 10.1086/527326

  • 20

    LichinchiG.GaoS.SaletoreY.GonzalezG. M.BansalV.WangY.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. 10.1038/nmicrobiol.2016.11

  • 21

    LichinchiG.ZhaoB. S.WuY.LuZ.QinY.HeC.et al (2016b). Dynamics of human and viral RNA methylation during zika virus infection.Cell Host Microbe20666673. 10.1016/j.chom.2016.10.002

  • 22

    LinS.ChoeJ.DuP.TribouletR.RichardG. I. (2016). The m 6 A Methyltransferase METTL3 promotes translation in human cancer cells.Mol. Cell62335345. 10.1016/j.molcel.2016.03.021

  • 23

    LongL.GaoL. D.HuS. X.LuoK.-W.ChenZ.-H.RonsmansC.et al (2016). Risk factors for death in children with severe hand, foot, and mouth disease in Hunan.China. Infect Dis48744748. 10.1080/23744235.2016.1185801

  • 24

    MaoQ.-Y.WangY.BianL.XuM.LiangZ. (2016). EV71 vaccine, a new tool to control outbreaks of hand, foot and mouth disease (HFMD).Expert Rev. Vaccines5599606. 10.1586/14760584.2016.1138862

  • 25

    OoiM. H.WongS. C.MohanA.PodinY.PereraD.ClearD. (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. 10.1186/1471-2334-9-3

  • 26

    QinL.DangD.WangX.ZhangR.FengH.RenJ.et al (2019). Identification of immune and metabolic predictors of severe hand-foot-mouth disease.PLoS One14:e0216993. 10.1371/journal.pone.0216993

  • 27

    ReddyS. M.SadimM.HelenowskiI. B.LiJ.YiN.KaklamaniV. G.et al (2013). Clinical and genetic predictors of weight gain in patients diagnosed with breast cancer.Br. J. Cancer109872881. 10.1038/bjc.2013.441

  • 28

    Shami ShahA.BatrouniA. G.KimD.PunyalaA.CaoW.HanC.et al (2019). PLEKHA4/kramer attenuates dishevelled ubiquitination to modulate WNT and planar cell polarity signaling.Cell Rep.2721572170 e8. 10.1016/j.celrep.2019.04.060

  • 29

    Shami ShahA.CaoX.WhiteA. C.BaskinJ. M. (2021). PLEKHA4 promotes Wnt/beta-catenin signaling-mediated G1-S transition and proliferation in melanoma.Cancer Res.8120292043. 10.1158/0008-5472.CAN-20-2584

  • 30

    WangJ.HuT.SunD.CarrM. J.XingW.LiS.et al (2017). Epidemiological characteristics of hand, foot, and mouth disease in Shandong, China, 2009-2016.Sci Rep.7:8900. 10.1038/s41598-017-09196-z

  • 31

    XuY.LiS.CaiC.LiuJ.WangY.JiangY.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.Medicine98:e18464. 10.1097/MD.0000000000018464

  • 32

    YangS. D.LiP. Q.HuangY. G.LiW.MaL. Z.WuL.et al (2018). Factors associated with fatal outcome of children with enterovirus A71 infection: a case series.Epidemiol. Infect.146788798. 10.1017/S0950268818000468

  • 33

    YangS. D.LiP. Q.LiY. M.LiW.LaiW.-Y.ZhuC.-P.et al (2017). Clinical manifestations of severe enterovirus 71 infection and early assessment in a Southern China population.BMC Infect. Dis.17:153. 10.1186/s12879-017-2228-9

  • 34

    YoonK. J.RingelingF. R.VissersC.JacobF.PokrassM.Jimenez-CyrusD.et al (2017). Temporal control of mammalian cortical neurogenesis by m(6)A methylation.Cell171:e17. 10.1016/j.cell.2017.09.003

  • 35

    ZhouJ.WanJ.GaoX.ZhangX.JaffreyS. R.QianS. B. (2015). Dynamic m(6)A mRNA methylation directs translational control of heat shock response.Nature526591594. 10.1038/nature15377

Summary

Keywords

enterovirus, neurological symptoms, children, N6-methyladenosine (m6A), m(6)A-RNA immunoprecipitation sequencing

Citation

Zhu D, Song Y, Hu D, Li S, Liu G, Li P and Yang S (2021) Characterization of Enterovirus Associated m6A RNA Methylation in Children With Neurological Symptoms: A Prospective Cohort Study. Front. Neurosci. 15:791544. doi: 10.3389/fnins.2021.791544

Received

08 October 2021

Accepted

08 November 2021

Published

07 December 2021

Volume

15 - 2021

Edited by

Fangfang Qi, Sun Yat-sen University, China

Reviewed by

Yan Xia He, Shenzhen Childen’s Hospital, China; Xiaoyuan Lin, Guangdong Academy of Medical Sciences, China

Updates

Copyright

*Correspondence: Peiqing Li, Sida Yang,

†These authors have contributed equally to this work and share first authorship

This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neuroscience

Disclaimer

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.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics