We enrolled children (aged 3-16 years) admitted to ASST Santi Paolo e Carlo Hospital in Milan (Italy) after their first seizures. Inclusion criteria were non-lesional focal or generalized epilepsy and neurotypical development. Exclusion criteria were chronic or acute intestinal diseases, special diets, and treatments with antibiotics or probiotics within the three months before enrollment.

We collected fecal samples at the time of enrollment and at 4 (T4) and 12 (T12) months after the introduction of ASMs. As a control group, we included 7 mentally and physically healthy children, age- and sex-matched, who were not on any medications. A 3-day dietary survey was completed by the caregivers at the time of enrollment. The diary included three consecutive days, one of which was during the weekend. Children were asked to maintain their usual eating pattern, and caregivers were instructed to record everything including drinks, the cooking method, and the amount of food. Dietary food records were processed using commercially available software (MètaDieta, METEDA srl, Italy). To monitor the children’s gastrointestinal function, caregivers were asked to describe and categorize feces according to the Bristol Stool Form Scale (BSFS, Lewis and Heaton, 1997).

We created a dedicated database that included information about gender, family history for epilepsy and/or febrile seizures, epileptic features, and neurologic examination results. Regarding the epileptic phenotype, we evaluated the type of seizures at onset and at last follow-up according to the ILAE 2017 classification, drug therapy, and electroencephalographic (EEG) pattern at onset and at last follow-up.

Psychomotor and cognitive development was evaluated by formal neuropsychological testing (if available) or best clinical assessment (based on developmental milestones and academic achievement). Brain MRI and metabolic findings were also included in the database.

The study was approved by the Local Ethics Committee (protocol number 2016/ST/199, 28 July 2016). Written informed consent was obtained from the parents and/or legal guardians of the enrolled patients/healthy subjects.

Fecal DNA extraction was performed using the Spin stool DNA kit (Stratec Molecular, Berlin, Germany), according to the manufacturer’s instructions. The V3–V4 hypervariable regions of the bacterial 16S rRNA gene were amplified with a two-step barcoding approach according to the Illumina 16S Metagenomic Sequencing Library Preparation (Illumina, San Diego, CA, USA). Briefly, DNA samples were amplified with dual-index primers using a Nextera XT DNA Library Preparation Kit (Illumina), while library concentration and quantification were determined using a KAPA Library Quantification Kit (Kapa Biosystems, Woburn, MA, USA) and Agilent 2100 Bioanalyzer System (Agilent, Santa Clara, CA, USA), respectively. The libraries were pooled and sequenced with a MiSeq platform (Illumina) for 2 × 250 base paired-end reads and a total of 2.5 Gbases raw reads were obtained.

The obtained 16S rRNA gene paired sequences were merged using Pandaseq (release 2.5; Masella et al., 2012). Reads were filtered by trimming stretches of 3 or more low-quality bases (quality < 3) and discarding the trimmed sequences whenever they were shorter than 75% of the original one. Bioinformatic analyses were conducted using the QIIME pipeline (release 1.9.0; Caporaso et al., 2010), clustering filtered reads into Operational Taxonomic Unit (OTUs) at 97% identity level and discarding singletons as well as possible chimeras. Taxonomic assignment was performed via the RDP classifier (Wang et al., 2007) against the SILVA database (release 132; Pruesse et al., 2007) from phylum to genus level. Alpha-diversity was computed through the QIIME pipeline using the Chao1, the number of OTUs (“observed species” index), Shannon diversity, and Faith’s Phylogenetic Diversity whole tree (PD whole tree) metrics. To compare the microbial community structure of the subjects for the beta-diversity analysis, weighted and unweighted UniFrac distances were used. A functional prediction analysis of the bacterial metabolic pathways has been performed through the PICRUSt software (v 1.0.1) (Langille et al., 2013) and KEGG pathways database (Kanehisa, 2004).

Statistical evaluation among alpha-diversity indices was performed by a non-parametric Monte Carlo-based test in the QIIME pipeline. The Permanova test (adonis function) in the R package vegan (version 2.0-10; Oksanen et al., 2013) was used to determine differences between the dataset cohorts through the beta-diversity analysis. Taxonomic and functional differences were computed and established through the non-parametric Mann-Whitney U-test; the Bonferroni correction has been applied to statistical results. Correlation analysis between nutritional and taxonomic data was conducted using Pearson’s rank correlation coefficients. All analyses have been conducted using R (version 3.6.3 via RStudio, version 1.2.1335).

P-values below 0.05 were considered significant among all comparisons and analyses.

We enrolled 8 drug-naive children (6 females, 75%; 2 males, 25%) with epilepsy onset between 3 and 13 years of age. Mean age at epilepsy onset was 8.9 ± 4.3 years. As control group, we collected fecal samples from 7 healthy subjects (healthy control, HC), who were age- and sex-matched (4 females, 66%; 3 males, 33%; mean age 8.0 ± 4.2 years; 5 delivered vaginally and 2 via cesarean-section).

Epilepsy diagnosis was Focal epilepsy in 4 subjects, Childhood Absence Epilepsy in 3 subjects, and Adolescence Absence Epilepsy in one. Brain MRIs were unremarkable in all patients and EEGs were consistent with the epilepsy diagnosis. In children with Focal epilepsy, the etiology was undetermined.

Twenty-two fecal samples were collected over the two years of study, with a three time-point scheme: the samples collected at seizure onset and before starting ASMs were named “drug-naive” (DN); the samples obtained 4 months after initiating treatment were labeled as “drug-therapy 4” (DT4); 6 out of 8 patients provided an additional sample after 12 months (labeled “drug-therapy 12”, DT12). Drug therapy was initiated after at least two seizures and titrated according to the clinical practice. All the patients were drug-responders and did not experience epileptic seizures during the experimental year of sampling and observation. Metabolic testing was negative for all patients. All clinical data are summarized in Table S1.

According to the Bristol Stool Form Scale (BSFC), none of the children was neither severely constipated nor suffering from diarrhea. Bristol scale was established at DN and remained unchanged at DT4 and DT12 for 7/8 patients; one patient reported a shift from type 2 to type 3 after therapy assumption.

None of the caregivers reported children gastrointestinal discomfort during the 12-month follow-up.

The dietary survey showed no differences in intake of macronutrients. Nutritional parameters are detailed in the relative supplementary table (Supplementary Material, Table S2).

Initially, possible differences between healthy controls and children experiencing seizures were investigated (HC vs DN). Data showed a reduced trend of bacterial abundances and biodiversity across the Chao1 and Shannon alpha-diversity metrics (Figure 1A). Beta-diversity reported a separation as well, in both Unifrac weighted (Figure 1B) and unweighted matrices (Figure 1C), even though it was not significant.

We then considered possible changes in the microbiota composition of children with epilepsy before and after the ASM introduction. A decreasing trend was observed in the alpha-diversity analysis, for all metrics, showing how the children’s gut microbial biodiversity was gradually and constantly reduced over the months of pharmacological treatment (Figure 2A).

The trend along samplings was observed also in the beta-diversity analysis among groups (Figure 2B). After 12 months of ASM treatment (DT12), we observed that the gut microbial communities of patients showed similarities to those of the HC subjects.

Considering each patient signature across the months of observation, we observed an overall shift along the first principal coordinate of the beta-diversity metric for most samples: 5/8 patients (62.5%) have their DN time-point at lower PCoA1 values, while 4/6 patients (66.7%) have the DT12 samples at the higher values (Figure 2C). Interestingly, all the patients who were on valproic acid (VPA) at DT12 were found to be closer to one another at the highest portion of the PCoA1 axis.

At phylum level (Figure 3A), Proteobacteria was found significantly increased at DT4 compared to HC (0.5% versus 2.9%, p=0.0137). On the contrary, the relative abundance of Verrucomicrobia increased substantially in the DN group with statistical significance when compared to HC (p=0.028) and to DT4 (p=0.012).

Among the main abundant families (Figure 3B), Lachnospiraceae and Akkermansiaceae were found significantly altered between the cohorts: Lachnospiraceae decreased steadily in the epileptic patients, from DN to DT12 (respectively, 14.9% and 11.2% compared to HC, 19.1%; p=0.025 HC vs DT12); Akkermansiaceae, instead, spiked at DN (8.7%, p=0.028 vs HC, 2.1%; p=0.012 vs DT4, 0.4%).

At genus level (Figure 3C), the peak of Akkermansia in patients at seizure onset was confirmed from the higher phylogenetic level (8.7% vs 2.1% in HC, p=0.0462). After 4 months of drug intake, this relative abundance dropped rapidly (0.5% at DT4) and, after 12 months, returned to an amount comparable to the abundance found in HC (3.2% at DT12). On the other hand, Faecalibacterium relative abundance was decreased in DN (4.7% vs 11% in HC), gradually returning to a normal level after the drug therapy (8.6% at DT4, 11.1% at DT12). ‘Prevotella 9’ showed the opposite trend: it was nearly absent before the therapy (1.0% in DN patients) and increased after 4 months of drug therapy (4.7% at DT4, 5.9% at DT12) to reach the level seen in HC (4%).

Changes in taxonomic composition across time and during ASM therapy in epileptic patients are detailed in Table S3 and significant bacterial groups are summarized in Figure 3D.

The Bristol score was 3 for all the patients (6/8) but two, who reported a mild constipation (Bristol score = 2) at DN, a condition that recovered at DT12 for one child.

Grouping together samples with BSFC 2 and with BSFC 3, we highlighted two bacterial genus that significantly differ. In particular, ‘Prevotella 9’, given the high abundance in one patient with BSFC 2 (14.77% compared to 0.01 in BSFC 3) was found significantly more abundant (p=0.004). ‘Coprococcus 2’, as well, was found more abundant in patients with Bristol score 2 (3.30 vs 0.82 in BSFC 3, p=0.021). Taxonomic abundances depending on the stool type are shown in Figure S1.

To investigate the ratio between Gram-positive and Gram-negative bacteria, the first 100 genera identified in the cohorts, corresponding to 99.99% of the cumulative abundance of the total bacterial genera, were analyzed. Children with epilepsy showed a higher abundance of Gram-negative bacteria, even after ASM therapy, compared with HC (Supplementary Material, Figure S2).

We used PICRUSt to predict possible pathways enriched or depleted in bacterial communities of children with epilepsy; 6,909 KO genes were analyzed. At KEGG functional level 3, the analysis predicted a significant enrichment in genes encoding enzymes for the “Bacterial secretion system” pathways in the microbiota of DN patients compared to HC (0.65% vs 0.56%, respectively; adjusted p-value=0.034), whereas genes involved in “Energy metabolism” pathways were increased at the beginning of ASM therapy (0.92% in DT4 compared to 0.85% in DN, adjusted p-value=0.029). Although not significantly, the predictive metabolic pathways showed a slight increase in the abundance of genes encoding enzymes involved in LPS-related functions in epileptic patients at all time points compared to that of the control group. In particular, the genes encoding for “Lipopolysaccharide biosynthesis proteins” were 0.32% in HC and 0.38%, 0.35%, 0.36% in DN, DT4, and DT12, respectively, and genes involved in “Lipopolysaccharide biosynthesis pathways” were 0.22% in HC and 0.27%, 0.24%, and 0.26% in children with epilepsy (DN, DT4, and DT12, respectively). Predicted gene abundances are reported in the Supplementary Figure S3.

Although dietary macronutrient intakes were found not statistically different among HC and DN subjects, the correlation analysis between the diet and the most abundant bacterial taxa revealed some divergent relationships in the two groups (Figure 4).

The Christensenellaceae R-7 group was negatively correlated in the HC group with carbohydrates (p=0.030), energy intake (p=0.040), and proteins (p=0.016), whereas was positively correlated in the DN subjects with carbohydrates (p=0.006) and energy intake (p=0.008). Akkermansia was only found positively correlated to the total fiber (p=0.012) intake in HC, as no dietary correlations were observed in DN subjects. On the other hand, Dialister was only found negatively correlated with total lipids (p=0.035) and energy intake (p=0.039) only in the DN group. The Lachnospiraceae NK4A136 group, which was found profoundly depleted in the DN subjects, was positively correlated with lipids (p=0.021), fibers (p=0.013), carbohydrates (p=0.001), proteins (p=0.006), and energy intake (p=0.004) in HC, whereas no correlations were shown for children with epilepsy.