Multi-omics integration reveals gut microbiota dysbiosis and metabolic alterations of cerebrospinal fluid in children with epilepsy

The study participants were recruited from the Neurology Department of the Second Affiliated Hospital and Yuying Children's Hospital of Wenzhou Medical University, from January 2022 to December 2022.

A total of 40 epileptic patients were enrolled and then divided into the 17 common epilepsy (CEP) patients and 23 refractory epilepsy (REP) patients including 16 with CSF collected, according to the diagnosis and classification criteria of epilepsy in ILAE from 2017 (). We also enrolled 10 non-epilepsy (NEP) patients matched by age and gender (including 8 with CSF collected), who were ultimately diagnosed with functional headaches (including migraine, cluster headache, emotion-related headache, and sleep-related headache) without positive neurological signs and allocated them to the control group. In addition, lumbar puncture for CSF detection for NEP group was required to rule out both CNS and neurological involvement as determined by the clinicians. The enrollment criteria for the epileptic patients in the study were as follows: (i) age 1–17 years; (ii) stable clinical symptomatology and EEG features from the last 3 months; (iii) recent brain MRI negative for potentially epileptogenic alterations (stroke, tumors, infectious diseases); and (iv) ASMs treatment in both CEP and REP groups given at stable dose from at least 3 months. Exclusion criteria for both the epileptic patients and NEP controls were as follows: (i) clear histories of chronic or allergic diseases; (ii) known inherited metabolic diseases; (iii) treatment with antibiotics, probiotics, or proton pump inhibitors within 3 months before the sample collection; (iv) Definite brain organic lesions caused by sequelae of intracranial infections, mechanical trauma, or spontaneous hemorrhage due to vascular malformations; and (v) treatment with ketogenic diet. According to routine pediatric neurological practice, lumbar puncture is not recommended for children with common epilepsy unless specific clinical concerns arise. Therefore, no CSF samples were collected from CEP patients in this study. [Scheffer et al., 2017]

Fecal samples were collected by the parents using a sampling kit (including a sterile collection tube and a cotton swab) and then stored at −80 °C within half an hour. CSF was collected by professional clinicians via lumbar puncture in a sterile polypropylene tube without any additives, labeled, and then stored at −80 °C. This study was approved by the Independent Ethics Committees and Institutional Review Board of the Second Affiliated Hospital and Yuying Children's Hospital, Wenzhou Medical University, and was conducted according to the ethical principle of the Declaration of Helsinki. Prior to enrollment, written informed consent was procured from the patients or their guardians.

The genomic DNA of Microbial community was extracted from fecal samples using the CTAB/SDS method according to the manufacturer's instructions. DNA concentration and purity were monitored on a NanoDrop 2000 and Qubit 3.0 Spectrophotometer (Thermo Fisher Scientific, Wilmington, United States). The hypervariable V3-V4 region of the 16S rRNA gene was amplified using the specific 341F and 805R primers with the barcode. Sequencing libraries were generated using the NEBNextUltra™ II DNA Library Prep Kit (Cat No. E7645), and their quality was evaluated with the Qubit 3.0 Spectrophotometer and the Agilent Bioanalyzer 2100 system. Finally, the libraries were sequenced on an Illumina NovaSeq platform, generating 250 bp paired-end reads. ®

The raw sequence data were processed using the QIIME 2 pipeline (). With the cutadapt plugin implemented in QIIME2, we removed the adaptor and primer sequences, followed by chimeras, low-quality read ends with a quality score below 35, and identification of amplicon sequence variants (ASVs) using the DADA2 plugin. Taxonomic assignments of ASV representative sequences were conducted based on the SILVA database (version 138) with the RDP Naive Bayesian Classifier algorithm. We removed sequences assigned to mitochondrial and chloroplast for bacteria from the ASV tables for subsequent analysis using the R package microeco (). Subsequent analyses were conducted on taxa with a mean relative abundance of more than 0.01% and present in at least 10% of the samples. [Bolyen et al., 2019] [Liu et al., 2021]

To address variations in sequence depth, we used the rarefy_even_depth method in the "phyloseq" package to rarefy the ASV table to the minimum sequence depth (). Rarefaction curves were produced for individual samples to evaluate the depth of sequencing by simulating the resampling process based on the microeco package. Three indices of alpha diversity (i.e., Shannon, Pielou_evenness and Richness) were evaluated by diversity functions from the R package vegan. The alpha diversity, community composition and ternary plot were visualized using the "ggplot2" packages in R. Beta diversity matrices were calculated using Bray–Curtis distances, and PCoA plots were generated from Bray–Curtis similarity matrices and visualized using "phyloseq" package. Group significance was determined by analysis of similarities (ANOSIM). [Mcmurdie and Holmes, 2013]

MaAsLin 3 was used to identify microbial features associated with diagnostic groups through multivariable linear modeling (). Age, sex, and sequencing depth were included as covariates, and associations with both abundance and prevalence were evaluated using False Discovery Rate (FDR) correction (q < 0.05). The gene functions related to the microbial community based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) were predicted by PICRUSt2 (). Ecologically relevant functional annotation of the gut microbiota was predicted by the trans_func class in the microeco package using the Functional Annotation of Prokaryotic Taxa (FAPROTAX) database (). The prediction of Microbial phenotypic characteristics was analyzed by EasyAmplicon tool based on the BugBase database (). [Nickols et al., 2024] [Douglas et al., 2020] [Louca et al., 2016] [Liu et al., 2023]

Random forest models were constructed using differential microbial genera (p < 0.05) identified by Wilcoxon rank-sum tests via the trans_diff function from the microeco package in R, to classify samples according to group labels. Model performance was evaluated using Leave-One-Out Cross-Validation (LOOCV). Probabilistic predictions from the LOOCV procedure were used to compute the Receiver Operating Characteristic (ROC) curve and the corresponding Area Under the Curve (AUC) using the pROC package, in order to evaluate the diagnostic effectiveness of the model. In addition, the Precision-Recall (PR) curve and its corresponding AUC were calculated using the PRROC package to assess performance under class-imbalanced conditions. ROC and PR curves were visualized based on predicted probabilities from LOOCV.

To assess the statistical significance of the model's predictive performance, permutation tests were performed for both AUC and classification accuracy. Specifically, group labels were randomly shuffled 1,000 times, and for each permutation, the random forest model was re-trained using the same feature set and evaluated under the identical LOOCV strategy. The resulting null distributions of AUC and accuracy were then compared to the corresponding values from the original, non-permuted model. Empirical-values were calculated as the proportion of permuted AUC or accuracy values that were greater than or equal to those observed in the original model. For visualization, histograms of AUC and classification accuracy from the 1,000 label-shuffled models were generated using the ggplot2 package, depicting the null distributions. Vertical dashed lines representing the original (non-permuted) AUC and accuracy values were overlaid to highlight their deviation from random expectations. p

For CSF metabolomics analysis, 100 μL of CSF samples were placed in EP tubes, resuspended in prechilled 80% methanol with vigorous vortexing, incubated on ice for 5 min, and centrifuged at 15,000 g, 4 °C for 20 min. A portion of the supernatant was diluted to a final concentration of 53% methanol using LC-MS grade water. The samples were then transferred to fresh Eppendorf tubes, centrifuged again at 15,000 g, 4 °C for 20 min, and the supernatant was finally injected into the LC-MS/MS system for analysis. UHPLC-MS/MS analyses were performed using a Vanquish UHPLC system (ThermoFisher, Germany) coupled with an Orbitrap Q Exactive™ HF-X mass spectrometer (Thermo Fisher, Germany).

The raw data files generated by UHPLC-MS/MS were processed using Compound Discoverer 3.3 (CD3.3, ThermoFisher) for peak alignment, peak picking, and quantitation of each metabolite. Main parameters included peak area correction with the first quality control (QC), actual mass tolerance of 5 ppm, signal intensity tolerance of 30%, and minimum intensity. Peak intensities were normalized to total spectral intensity and used to predict molecular formulas based on additive ions, molecular ion peaks, and fragment ions. Peaks were matched with the mzCloud (), mzVault (), and MassList databases for accurate qualitative and relative quantitative results. For non-normally distributed data, relative peak areas were obtained by standardizing according to the formula: sample raw quantitation value/(sum of sample metabolite quantitation value/sum of QC1 sample metabolite quantitation value). Compounds with CVs of relative peak areas in QC samples >30% were removed, and the identification and relative quantification of metabolites were finalized. https://www.mzcloud.org/↗ https://www.mzcloud.org/↗

Metabolite identification confidence was classified in accordance with the Metabolomics Standards Initiative (MSI). Due to the untargeted nature of the study, most annotated metabolites correspond to MSI Level 2, indicating identification based on high-resolution mass spectra, retention time alignment, and MS/MS spectral matching against reference libraries. MSI Level 1 confirmation, which requires validation using authentic chemical standards, was not feasible for all compounds.

These metabolites were annotated through the KEGG database () and the HMDB database (). Principal components analysis (PCA) and Partial least squares discriminant analysis (PLS-DA) were conducted using the MetaboAnalyst 6.0 (). In the PLS-DA analysis, metabolites with Variable Importance in Projection (VIP) >1 were considered as important variables for classification. Metabolic pathway analyses associated with these metabolites with VIP > 1 were annotated using MetaboAnalyst 6.0. We applied univariate analysis (-test) to calculate the statistical significance (-value). Metabolites passing the threshold of-values < 0.05 and fold change ≥1.5 or ≤ 2/3 were regarded as statistically different substances. Volcano plots were used to filter metabolites of interest which based on log(fold change) and -log(-value) of metabolites by ggplot2 in R. Multivariate ROC analyses on a combination of selected features using random forest analysis were performed to calculate the AUC and analyze the prediction model using MetaboAnalyst 6.0. https://www.genome.jp/kegg/pathway.html↗ https://hmdb.ca/metabolites↗ [Pang et al., 2024] t p p p 2 10

Spearman correlation analysis was employed to examine the intricate relationship between the microbiome and metabolome. Differential metabolites and key microbiota identified via random forest analysis were selected to compute correlation coefficients and statistical significance using the cal_cor function from the microeco package with visualization through a heatmap. The-value < 0.05 was considered statistically significant for correlation. p

A total of 40 patients with epilepsy including 17 children with CEP and 23 children with REP, along with 10 NEP controls were enrolled in the study. The NEP cohort included children presenting with functional headache syndromes—such as migraine, emotion-related headache, and sleep-related headache—who exhibited no signs of organic or neurological disease. These participants underwent lumbar puncture as part of standard clinical evaluation to exclude central nervous system involvement. No significant differences were observed among the three groups in terms of age, gender, body mass index, and dietary habits. All clinical data in the cohort are summarized inand. Table 1 Supplementary Table S1

To characterize gut microbiota structure among participants, fifty stool samples were collected from the three groups and analyzed via Illumina MiSeq sequencing of the V3-V4 region of the bacterial 16S rRNA gene. After quality filtering, a total of 3,625,932 high-quality readings were obtained, with an average of 72,519 ± 8,891 sequences per sample (). After denoising and filtering, these sequences were classified into a total of 8,890 ASVs (). Supplementary Table S2 Supplementary Table S3

Rarefaction curve analysis of richness, an alpha diversity index, demonstrated sufficient phylogenetic coverage, as curves reached saturation with increasing sequencing depth (). Community composition analysis revealed that bacterial ASVs were affiliated with 54 known phyla, with the top nine phyla represented in a bar graph (;). Of these, five dominant phyla (>1% abundance across all groups) were identified, including, and(). Further examination showed that while, andwere primarily observed in REP,was predominantly present in CEP (,), although differences in the average relative abundance of these taxa across groups were not statistically significant (). Hierarchical clustering analysis revealed that the CEP cohort clustered together with REP, while distinctly separating from NEP (). Figure 1A Figure 1B Supplementary Table S4 Figure 1B Figures 1B C Supplementary Figure S1 Figure 1D Firmicutes, Bacteroidota, Proteobacteria, Actinobacteriota Campylobacterota Desulfobacterota, Actinobacteriota Campylobacterota Fusobacteriota

Gut microbial diversity analysis showed no significant differences across groups in three alpha diversity indices: Pielou_evenness, Richness, and Shannon index (). However, Principal Coordinates Analysis (PCoA) and ANOSIM analysis based on Bray–Curtis distance indicated significant differences in the gut microbial community composition (= 0.091,= 0.001;). To characterize microbial features associated with diagnostic groups, we applied multivariable association modeling using MaAsLin 3, enabling robust identification of taxa differentially abundant across the REP, CEP, and NEP cohorts while accounting for covariates and data compositionality. The analysis identified 13 taxa significantly associated with diagnostic groups (< 0.05), including eight enriched in REP and five negatively associated with CEP, as detailed in. Notably, members of the Clostridia lineage—including, and—exhibited strong positive associations with the REP group. Taxa belonging to theclade—namely, and—were negatively associated with CEP. The top 13 group-discriminatory taxa identified by MaAsLin 3 are depicted in, highlighting distinct microbial signatures characteristic of REP and CEP cohorts. Figure 1E Figure 1F Supplementary Table S5 Figure 1G R p q Clostridiaceae, Clostridiales, Clostridium_sensu_stricto_1 Erysipelotrichales Synergistia Synergistia, Synergistaceae, Synergistales Synergistota

PICRUSt2-based functional predictions identified 212 KEGG pathways from filtered ASVs. While PCA analysis indicated no significant differences in the overall composition of predicted functions among the three groups (), ALDEx2's Kruskal–Wallace test revealed statistically significant differences in the relative abundance of eight KEGG pathways spanning five biological processes: translation, immune disease, biosynthesis of other secondary metabolites, carbohydrate metabolism, and glycan biosynthesis and metabolism (). Figure 2A Figure 2B

To further assess potential biotic nutrient cycling mechanisms, we employed the FAPROTAX database, linking microbial taxonomy to functional traits. FAPROTAX analysis identified 38 functional groups associated with carbon (C), nitrogen (N), and sulfur (S) cycling. Among these, fermentation emerged as the most dominant functional category, followed by anaerobic chemoheterotrophy and animal parasites or symbionts (left panel in). ANOVA analysis revealed significant differences in the relative abundance of 21 functional groups (right panel in). Figure 2C Figure 2C

Additionally, bacterial phenotypic characteristics were inferred using the BugBase database, predicting microbial functions across nine phenotypic categories, including pathogenicity, aerobic and anaerobic metabolism, facultative anaerobic capability, mobile elements, biofilm formation, Gram-positive/Gram-negative classification, and oxidative-stress tolerance. BugBase predictions indicated that the REP and NEP groups were enriched in aerobic taxa compared to CEP (). Conversely, the CEP group exhibited significantly higher abundance of anaerobic taxa relative to REP (). Figure 2D Figure 2E

To explore whether gut microbial composition could be utilized for the diagnosis of epilepsy in children, differential microbial genera were identified using the Wilcoxon rank-sum test via the trans_diff function in the microeco R package. Genera with FDR-adjusted-values < 0.05 were retained as microbial signatures for subsequent classification modeling. p

There were eight genera with significant differences between CEP and NEP groups (FDR adjusted< 0.05, Wilcoxon rank sum test;). These genera belonged to, and, with CEP enriched in genera fromand, whereas NEP showed a higher abundance of genera fromand. Specifically,, andwere more prevalent in CEP, while, andwere dominant in NEP (). These differential genera were used as input features for a random forest classification model, which was evaluated through LOOCV. The model demonstrated high discriminatory performance, achieving an AUC of 0.906 for the ROC curve (), an AUC of 0.942 for the PR curve (), and an overall classification accuracy of 0.852. Statistical robustness was confirmed via permutation testing (1,000 iterations), with empirical-values of 0 for both AUC () and accuracy (), indicating that the observed classification performance was highly unlikely to result from random chance. p Proteobacteria, Firmicutes Actinobacteriota Proteobacteria Firmicutes Proteobacteria Actinobacteriota Klebsiella, Megasphaera, Romboutsia Erysipelatoclostridium Neisseria, Corynebacterium, Haemophilus Rothia p Figure 3A Figure 3A Figure 3B Figure 3E Figure 3C Figure 3D

Further analysis identified 14 microbial genera displaying significant differences between REP and NEP groups (FDR adjusted< 0.05, Wilcoxon rank sum test;). These genera were distributed across, and, with REP characterized by an increased presence ofand, whereas NEP exhibited greater representation of, and. Notably, genera such as, andwere enriched in the REP group, while the other the 10 genera such as, andwere more abundant in the NEP group. A random forest classification model constructed using these genera also demonstrated strong predictive capability, yielding an AUC of 0.913 for the ROC curve (), an AUC of 0.965 for the PR curve (), and an accuracy of 0.818. Permutation tests reinforced model reliability, with empirical-values of 0 for AUC () and 0.002 for accuracy (), confirming the robustness of the classification approach. p Firmicutes, Bacteroidota, Proteobacteria, Actinobacteriota Fusobacteriota Firmicutes Bacteroidota Firmicutes, Proteobacteria, Actinobacteriota Fusobacteriota Allobaculum, Alistipes, Dubosiella Ileibacterium Gemella, Neisseria Ruminococcus p Figure 3F Figure 3G Figure 3J Figure 3H Figure 3I

To assess the generalizability of our microbial classification models, we applied the CEP vs. NEP and REP vs. NEP random forest classifiers to an independent pediatric epilepsy 16S rRNA dataset derived from a European cohort (), representing a geographically and ethnically distinct population from our original study. This external cohort includes medication-sensitive (MS) and medication-resistant (MR) patients, as well as healthy controls (HC). The CEP-trained model demonstrated strong discriminatory ability when tested on MS vs HC samples, achieving an AUC of 0.838 and an accuracy of 0.812 (). Permutation testing with 1,000 iterations yielded empirical-values of 0 for both metrics, confirming the statistical robustness of the classification (,). Similarly, the REP-trained model achieved an AUC of 0.796 and an accuracy of 0.723 for distinguishing MR from HC samples, with permutation-values of 0 for AUC and 0.004 for accuracy (–), further supporting the reliability of the model's performance in an external cohort. Feature overlap analysis revealed that 7 of the 8 genera used in the CEP classifier () and 9 of the 14 genera used in the REP classifier () were present in the external European dataset. These retained genera were among the top-ranked features in the original models, suggesting that conserved microbial signals contributed to the preserved classification performance. [Riva et al., 2025] Supplementary Figure S2A Supplementary Figures S2B C Supplementary Figures S2D F Supplementary Figure S2G Supplementary Figure S2H p p

These findings highlight the potential of genus-level microbial signatures as effective diagnostic markers for epilepsy in children. The high classification accuracy and statistical validation suggest that microbial composition could be leveraged for reliable disease stratification, providing valuable insights into the microbial characteristics distinguishing different clinical conditions.

To further explore metabolic dysregulation associated with epilepsy, we conducted untargeted metabolomic profiling of CSF samples from 16 children with REP and 8 children with NEP using LC-MS. Following peak detection, alignment, and deconvolution, a total of 1,082 unique putative metabolites were detected in both positive and negative ionization modes. After filtering exogenous compounds and duplicates, 389 endogenous metabolites were retained based on HMDB annotations (). These metabolites were categorized into 31 secondary categories and further grouped into six primary classifications based on KEGG pathway annotation (). The majority of metabolites were involved in metabolic pathways, particularly global and overview maps, amino acid metabolism, and carbohydrate metabolism. HMDB-based classification indicated that organic acids and derivatives (29.56%), lipids and lipid-like molecules (20.82%), and organoheterocyclic compounds (15.94%), followed by benzenoids (11.05%), collectively accounted for 77.37% of the annotated metabolites, reflecting a substantial representation of central metabolic intermediates and bioactive compounds (). Supplementary Table S6 Supplementary Figure S3 Figure 4A

Unsupervised principal component analysis (PCA) revealed clear metabolic separation between REP and NEP groups, with the first two principal components explaining 68.6% (PC1) and 10.8% (PC2) of the total variance, respectively. Statistical significance of the observed separation was confirmed via permutation testing (= 999,= 0.015), supporting the existence of distinct metabolic profiles between the two groups (). To further assess class discrimination and identify key metabolites, a PLS-DA model was constructed, demonstrating clear intergroup separation (). Ten-fold cross-validation yielded a classification accuracy of 84.7% (= 0.584,= 0.403;), while ROC analysis returned an AUC of 0.805 (), indicating robust explanatory and predictive power. Model reliability was further supported by permutation testing (= 1,000), which showed that the observed classification accuracy significantly (0.847) exceeded that of randomly permuted models (= 0.044;). n p R Q n p Figure 4B Figure 4C Figure 4C Supplementary Figure S4 Figure 4D 2 2

Based on VIP scores derived from the PLS-DA model, 23 key discriminatory metabolites were identified (VIP > 1.0), including citric acid, taurine, L-glutamic acid, and acetyl-L-carnitine, all of which are involved in energy metabolism and neurotransmitter pathways (). Metabolic pathway analysis using KEGG revealed 19 significantly impacted pathways, with the six most affected (impact score > 0.05) including taurine and hypotaurine metabolism, alanine, aspartate and glutamate metabolism, histidine metabolism, arginine biosynthesis, citrate cycle (TCA cycle), and vitamin B6 metabolism (and). These pathways are functionally linked to neuronal excitability, neurotransmission, and mitochondrial activity. Figure 4E Figure 4F Supplementary Table S7

To visualize metabolite expression differences between children with REP and NEP controls, an independent-test was conducted. This univariate analysis identified 40 significantly altered metabolites (< 0.05), including six upregulated metabolites (fold change > 1.5) and 34 downregulated metabolites (fold change < 2/3) in the REP group (). To further refine feature selection for diagnostic modeling, significantly altered metabolites from the-test were intersected with those exhibiting VIP scores > 1 in the PLS-DA mode. This approach yielded six shared metabolites for model construction (), comprising three elevated and three reduced metabolites in REP compared to NEP (). t p t Figure 5A Figure 5B Figure 5C

These six discriminatory metabolites were used to construct a classification model employing a random forest algorithm with LOOCV to evaluate their predictive capacity (). The model exhibited promising diagnostic potential, achieving an AUC of 0.875 for ROC curve (), an AUC of 0.943 for PR curve (), and an overall classification accuracy of 0.833. Model reliability was further assessed through permutation testing (= 1,000), where sample labels were randomly shuffled while preserving model architecture. The original model consistently outperformed the permuted models, yielding empirical-values of 0.006 for both AUC and accuracy (,), supporting reliability and non-random predictive ability of the model. These findings highlight the diagnostic relevance of metabolite-based classification. Figure 5D Figure 5E Figure 5F Figures 5G H n p

To investigate potential associations between metabolite alterations and gut microbiota composition, Spearman's rank correlation analysis was performed. This analysis examined correlations between six shared metabolites and the relative abundance of 14 differential bacterial genera (identified by the Wilcoxon rank-sum test) in paired fecal and CSF samples from 16 children with REP and 8 with NEP. Applying a selection criterion requiring at least one significant correlation (p < 0.05) for both metabolites and bacterial genera, four metabolites exhibited significant associations with 10 bacterial genera (). Correlation analysis showed that 5-Methyluridine, alpha-Ketoisocaproic acid (KIC), and alpha-Ketoisovaleric acid (KIV) had strong positive correlations with three bacterial genera while exhibiting negative associations with seven others, though some of these negative correlations were not statistically significant (). Notably, these three metabolites were more abundant in the REP group compared to NEP (). Conversely, Acetyl-L-carnitine exhibited negative associations with three bacterial genera and positive correlations with seven, although two of these positive correlations were not statistically significant, and its abundance was significantly higher in NEP (,). Figure 6A Figure 6A Figure 5C Figures 5C 6A

To assess the combined discriminatory power of these microbial and metabolic features, a random forest model was constructed using the 10ten bacterial genera and four significantly correlated metabolites, employing LOOCV for evaluation. Feature importance rankings () identified, and Acetyl-L-carnitine as the top three contributors to classification performance. The model demonstrated strong discriminatory capability, as evidenced by an ROC curve with an AUC of 0.953 (). Model robustness was further evaluated via permutation testing (1,000 iterations). The AUC permutation test yielded a-value of 0 (0/1,000), indicating that the observed AUC was highly unlikely to occur by chance (). Similarly, accuracy distribution analysis showed the original model achieving an accuracy of 0.875 (), significantly outperforming permuted models (= 0.001, 1/1,000). The PR analysis further reinforced classification performance, with an AUC of 0.976 (), reflecting high precision and recall in distinguishing REP from NEP. These findings highlight the strong interplay between gut microbiota and CSF metabolites, suggesting that their combined signatures may serve as effective biomarkers for distinguishing REP from NEP. The robustness of the classification model further supports the potential utility of microbial-metabolite interactions in disease stratification. Figure 6B Figure 6C Figure 6D Figure 6E Figure 6F Corynebacterium, Ruminococcus p p