Sixty-four schoolchildren (8–11 years) belonging from five primary schools in the neighborhoods of Turin (northwest Italy) were enrolled and classified as Sedentary based on Triaxial Accelerometry and Physical Activity Questionnaire for Older Children (PAQ-C-It) (42).

Parents/guardians and teachers furnished written informed consent for participation in the study, adhering to the ethical standards outlined in the 1964 Declaration of Helsinki. The ethics committees on human research of the University of Turin (9 March 2020: Protocol #134691) and Naples (17 January 2020: Protocol #376/19) granted approval for the study.

Anthropometric measurements were reported. The procedures used to take anthropometric measures were standardized as described in Lupo et al. (42): height was measured using a portable stadiometer (Model 214; Seca, Hamburg, Germany), weight was measured using an electronic scale (Model 876; Seca, Hamburg, Germany), and waist circumference was measured while standing, halfway between the lowest rib and the iliac crest, using an Ana elastic meter. Body Mass Index (BMI) was calculated by dividing weight by height squared (kg/m²; Table 1). On average, the BMI of all enrolled schoolchildren fall within a range between the 75th and 85th percentile (43). For the daily dietary intake assessment all participants completed the questionnaire (3-Days food records). Records were processed using Winfood software (Medimatica S.u.r.l., Colonnella, TE, Italy). Statistical analysis was performed through Jamovi (2.3.26 version).

All sedentary enrolled schoolchildren were randomly assigned to 2 different types of exercise: Structured Exercise (Sa), 2 h/w muscle activation, strength and coordination exercises (supervised by a kinesiologist); Daily Mile (Dm), 1 h/w Sa plus 15 min/4 for week of walking/running outdoors (supervised by class teacher) or control group (curricular activities, 2 h/w PA supervised by class teacher, Ct) for 6-months (October, 2021–March, 2022). All interventions were planned within the common weekly volume of physical activity (i.e., 2 h), despite differently distributed according to the type of intervention. All children belonging to a single school class were enrolled in one of the three sub-groups, and no information about results of previous fitness measures (classification in active and no-active children) have been shared with them.

In particular, 22 schoolchildren belonging to the same class were involved in the Sa group and performed the following protocol: two (1 h in a day, and one in another day) school weekly hours of physical activity leaded by an expert (i.e., Sport Sciences Master degreed) according to American College of Sports Medicine guidelines (44). In particular, a single hour of this school physical activity type was distributed in 3 parts, 20 min/each focused on: fundamental movement skills, strength activities (strengthening activities) and aerobic activity. For the first part of exercises, basic body movements such as running, jumping, catching, balance, agility and coordination were proposed (45). The strength activities included all those movements that are generated by muscle contraction against resistance, aiming not only to improve individual strength level, but also anaerobic endurance and muscle size, balance, flexibility, mobility, posture, bone density, strength, self-esteem, prevention of chronic conditions such as diabetes, heart disease and depression, reduction in injury risk and pain management (46). Finally, for the last 20 min of the session, the proposed exercises were characterized by the contraction of large muscle groups, exclusively based on the aerobic energy system, by practicing walking, running, playing games, dancing, all activities recognized as able to increase cardiorespiratory fitness, weight maintenance, strengthens muscles and bones, mental health, social skills, and academic performance (44, 47). The second semi-structured intervention group (Dm) included a class of 22 schoolchildren that performed the following protocol: four weekly days characterized by the practice of The Daily Mile,¹ plus the 1 h of curricular physical activity leaded by the generic teacher. In particular, for The Daily Mile, children used to go outside within the spaces belonging to the school, to jog and/or run independently for 15 min, supervised by the generic teacher. A class of 20 schoolchildren continued to carry out the curricular physical activities under the supervision of their generic teacher, control sub-group (Ct).

Lower-Limb Muscle Power: The standing long-jump test was performed to assess lower-limb power. The test was performed as previously reported (48). Two trials were performed (with a 1 min pause in between), and the best score was considered for statistical analysis. Not all enrolled children carried out the POST-intervention assessment, therefore the final sample size for this test was: Ct = 19, Dm = 21, Sa = 13.

Balance: Balance was tested using the single-leg stance test. The dominant leg of the participants was determined by asking them which was their favorite leg (49). The test was performed as previously reported (48). The time (in seconds) was recorded using a stopwatch with an accuracy of 0.01 s. Not all enrolled children carried out the POST-intervention assessment, therefore the final sample size for this test was: Ct = 19, Dm = 18, Sa = 21.

Cardiorespiratory Fitness: The 20 m shuttle-run test, according to the Leger test soundtrack version (CAEP Quebec Faca), was used to evaluate the cardiorespiratory fitness. The test was performed as previously reported (50). The test score was the last reached stage (converted into the corresponding distance expressed in meters) (51). Not all enrolled children carried out the POST-intervention assessment, therefore the final sample size for this test was: Ct = 19, Dm = 21, Sa = 21.

Sprint Ability: Children’s sprint ability was evaluated by the 20 m linear sprint test (52). Individually, participants started the test with feet behind the startling line (53) and ran the 20 m distance as fast as possible, twice with a 1 min rest between the trials. The time elapsed from the start to the finish line was measured through infrared reflex photoelectric cells (Witty—Wireless Training Timer; Microgate, Udine, Italy), measuring to the nearest 0.01 s. Not all children carried out the POST-intervention assessment, therefore the final sample size for this test was: Ct = 19, Dm = 22, Sa = 21.

The volunteers were instructed not to consume food and refrain from using mouth cleanliness products 1 h prior to saliva retrieval. At minimum, 2 milliliters of unstimulated spit were gathered, chilled, and preserved at −80°C until analysis. DNA extraction from saliva specimens was accomplished utilizing the MagPurix Bacterial DNA Extraction Kit (ZP02006; Zinexts Life Science Corp.) following the manufacturer’s instructions. DNA concentration was measured utilizing the Qubit dsDNA BR and HS assay kit (Life Technologies, CA, United States).

PCR amplification was conducted to target the variable V3-V4 segments of the 16S rRNA bacterial genes. Specific primers embedded with barcodes and potent enzymes were utilized for PCR execution. The PCR primer sequences were as follows: (forward 341F: CCTAYGGGRBGCASCAG; reverse 806R: GGACTACNNGGGTATCTAAT). PCR products ranging from 450 to 500 base pairs were isolated using 2% agarose gel electrophoresis. For library construction, equimolar amounts of PCR products from each sample were combined, end-repaired, A-tailed, and subsequently ligated with Illumina adapters. Library quality control was conducted using Qubit and real-time PCR for quantification, and bioanalyzer for verifying the insert size distribution. Libraries were sequenced on a paired-end Illumina platform, producing 250 bp paired-end raw reads. The raw sequencing data are available in Zenodo (Publication date: July 1st, 2024). ²

We conducted the analyses as reported in Mancini et al. (41). In brief, we used the R package DADA2 and its workflow (54) to infer the Amplicon Sequence Variants (ASVs) and for the taxonomic assignments.

Briefly, we first filtered and trimmed raw sequencing reads to remove low quality bases, adapters and identical reads. Then, the reads were denoised, merged and filtered to remove artifacts (PCR, and PhiX related chimeras). Then, we obtained the ASVs counts and assigned taxonomy level annotations using the SILVA database (version: v138, nr99) (55). The data were structured in objects including the ASVs quantifications, the taxonomy annotations, the sample group information and the phylogenetic tree using the phyloseq and the APE packages (56, 57). Finally, based on the initial DNA concentration, we removed possible contaminant ASVs by using the “prevalence” method of the decontam package (58). For the downstream analyses we used the MicrobiomeAnalystR package (59, 60) and conducted data normalization, measures of diversity and differential abundance estimation. As previously described (61), we filtered and normalized the raw ASV counts based on their low abundance (at least the 20% of a given ASV values should contain at least 4 counts) and low variance (based on Inter-quantile range ± 10%) in order to discard possible sequencing errors. Moreover, ASVs showing constant values in all samples could be excluded from the comparative analyses. Finally, we used the total sum scaling in order to obtain normalized ASV counts. Based on these counts, we conducted the downstream analyses in order to evaluate the alpha-and beta-diversity reporting the species richness and evenness and the diversity of species between the groups. Then we evaluated the species richness from the results of sampling and whether samples were sufficiently sequenced to represent their identities by using the rarefaction analysis. Moreover, we measured the differential abundance with the LEfSe (Linear Discriminant Analysis Effect Size) algorithm (62) for biomarker discovery. Finally, we used PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) to infer the metabolic potential of microbial communities (KEGG pathways). In this analysis, we started from the ASVs belonging to the significant genera and species obtained by the differential abundance analysis (LefSE algorithm).

We calculated the ACE and Chao1 indices (for accounting species richness) and the Fisher metrics (to consider both richness and evenness).

To evaluate the alpha-diversity, we calculated the ACE and Chao1 indices and the Fisher metrics and assessed the statistical significance of group comparisons between the groups of samples under study by using the Mann–Whitney test. We evaluated diversity of species between the groups (the beta-diversity) by PCoA (Principal Coordinates Analysis) of the UniFrac distances and assessed the statistical significance with the PERMANOVA test. The statistical significance of rarefaction profiles was assessed by using the Mann–Whitney test. We assessed the statistical significance by using the Mann–Whitney test. Finally, we measured the differential abundance with the LEfSe algorithm (62) for biomarker discovery. It involves the Kruskal-Wallis rank sum test to identify features (e.g. Species or Genera) with significant differential abundance in the groups, followed by linear discriminant analysis (LDA) to evaluate the relevance (the effect size) of the selected features. Differentially abundant features were considered if the p-value was less than or equal to 0.05 and if the Log LDA was greater than or less than 0.5. We used PICRUSt to infer the metabolic potential of microbial communities (KEGG pathways). In this analysis, we started from the ASVs belonging to the significant genera and species obtained by the differential abundance analysis (LefSE algorithm). Statistical analysis of variables shown in Tables 1–3 were performed using Jamovi software (version 2.2.5.0). The normality of the data was verified by applying the Shapiro–Wilk test and mean ± standard deviation (SD) for each variable reported in Tables. Comparisons between groups were determined with Two- or One-way ANOVA for all normally distributed variables or with Kruskal-Wallis H test for non-parametric variables. Paired Sample t-test within intervention groups for physical fitness parameters was performed and followed by Effect size (>0.2 and < 0.5 small; >0.5 and < 0.8 moderate; >0.8 large). The level of significance was set at p < 0.05.

Anthropometric measurements of the schoolchildren participating in this study, were acquired before and after the training intervention (PRE-POST-) and reported in Table 1. No significant difference in the gender frequency of individuals was found [Χ² (2) = 1.14, p = 0.57]. Basal differences (PRE-) found in height [F (2, 63) = 19.4; p < 0.001 one-way ANOVA followed by Tukey post-hoc test revealed a statistically significant difference (p < 0.01) in Ct vs. Sa and Dm vs. Sa] can be attributed to the school class to which the children belong. Weight, BMI, Waist/Height parameters did not vary among groups in baseline conditions (Krustal-Wallis One-way ANOVA; p > 0.05). Furthermore, two-way ANOVA revealed that there was not a statistically significant interaction between time (PRE-and POST-) and anthropometric measurements (relatives p-value and Tukey post-hoc were reported in Table 1).

On the same time, One-way ANOVA revealed no significant differences in all diet components analyzed among Sa, Dm and Ct groups (Table 2).

The schoolchildren were assessed for physical fitness both before and after the 6-mo training intervention (PRE-and POST-). Significant PRE-POST improvements were observed for lower limb power (standing long-jump) and sprint (linear 20 m) only in Sa group. Conversely, balance (single-leg stance) and endurance (20 m shuttle-run test) resulted significantly improved in both Sa and Dm groups; slight significant improvement only in balance test was found in Ct (Table 3). Two-way ANOVA (Type of intervention × Time PRE/POST with Physical test as dependent variable) resulted not significant for 20 m shuttle-run, Single-leg stance and Sprint (p > 0.05). Differently, significant interaction between time (PRE-and POST-) and type of intervention for Standing long-jump test [F (2, 50) = 4.91; p = 0.011] was found.

The Illumina sequencing of the hypervariable V3-V4 regions of the 16S rRNA bacterial genes generated 2×250 bp paired-end reads. Overall, the percentage of bases with quality scores above 20 and 30 (Q20 and Q30, respectively) was of 97.47 and of 92.93, respectively (Supplementary Figure S1A). The percentage of GC nucleotides was of 51.85 (Supplementary Figure S1B).

On average, we obtained 172,066 reads per-sample (min = 150,914; max = 189,035). The raw reads were filtered and trimmed (median = 171,879; min = 150,769; max = 188,851), denoised on the forward (median = 170,557; min = 149,387; max = 187,291) and reverse (median = 170,671; min = 149,387; max = 187,291) directions, merged (median = 149,274; min = 112,706; max = 170,971) and chimeric reads removal (median = 114,432; min = 87,964; max = 136,831). After merging, the median length of reads was of 425 bp (min = 421; max = 427; Supplementary Figures S1C,D).

We assigned taxa using the SILVA database of non-redundant sequences (version: v138, nr99). Overall, we could identify a total of 4,759 taxa (ASVs) that were taxonomically assigned to the kingdom of bacteria. The 97.92% of the ASVs was annotated at the phylum level (n = 13 phyla), the 97.75% at the class level (n = 21 classes), the 97.18% at the order level (n = 56 orders), the 92.54% at the family level (n = 96 families), the 82.92% at the genus level (n = 215 genera) and the 6.51% was annotated up to the species level (n = 172 species; Supplementary Figure S1E). From the initial set of annotated ASVs, we discarded a total of 31 taxa as possible contaminants. Moreover, as described in Methods we removed low abundant and low variable ASVs to obtain the final set of 3,503 ASVs that was normalized and used for down-stream analyses.

We estimated the alpha-diversity at the genus level. Overall, we observed slight, no significant differences among the groups (Ct, n = 19; Dm, n = 20; Sa, n = 22) indeed, ACE (Figure 1A), Chao1 (Figure 1B) indexes are not significant, p > 0.05 (Mann–Whitney test). In particular, the SA group showed higher richness and evenness (Fisher metrics; Figure 1C) compared to the Control Group (Ct, p = 0.048).

The beta-diversity analysis, evaluated by weighted UniFrac distance metric did not evidentiate significant differences in microbial composition among the groups (R² = 0.0414, p = 0.147, PERMANOVA test; Figure 1D). Differently, the microbial composition resulted significant different when measured by unweighted UniFrac distances (R² = 0.04229, p = 0.024, PERMANOVA test; Figure 1E). Moreover, the rarefaction analysis clearly showed that we were able to capture in depth the species richness from the results of sampling and sequencing in the groups without any statistically significant difference between the groups under study (p > 0.8; Mann–Whitney test) Indeed, we can assess that the variation observed between the two groups was due to the different taxa abundance and to the types of taxa present in their saliva microbiome.

We evaluated and compared the taxa abundance in the final set of 3,503 filtered and normalized ASVs (see above). Overall, we identified n = 13 phyla, n = 21 classes, n = 56 orders, n = 96 families, n = 215 genera and n = 172 species.

At the phylum level, on average, the most abundant bacteria were the Proteobacteria, the Firmicutes, the Bacteroidota accounting for the 34.27, 32.27%, and the 19.93%, respectively (Figure 2A). The most represented classes were the Gammaproteobacteria (34.27%), the Bacilli (21.99%), and the Bacteroidia (19.93%; Figure 2B). The most prevalent orders were the Lactobacillales (19.52%), the Bacteroidales (19.45%), the Pasteurellales (12.3%) and the Pseudomonadales (11.72%; Figure 2C). Among the most abundant families, Streptococcaceae (17.37%), the Prevotellaceae (16.35%), the Pasteurellaceae (12.3%) and the Pseudomonadaceae (11.55%; Figure 2D) were found. At the genera level, the Streptococcus (17.36%), Prevotella (13.72%), Pseudomonas (11.55%), Haemophilus (11.13%), and Neisseria (9.63%; Figure 2E) are the most aboundant. Finally, the top abundant species that we were able to classify were the Prevotella melaninogenica (7.28%), Fusobacterium periodonticum (4.69%), Haemophilus parainfluenzae (1.73%), Veillonella dispar (1.53%) and the Rothia mucilaginosa (1.46%; Figure 2F).

Taken together, the diversity and the overall abundance estimates are in line with our recent results (41). Indeed, we can assess that the data resulting from the sampling, sequencing and analyses protocols are reliable. Moreover, the variation observed among the groups, although slight and probably due to the small sample size, is due to the different taxa abundances and to the types of taxa in the analyzed microbiomes.

As described in Methods we used the LEfSe algorithm to perform the differential abundance analysis and to find the taxa that could explain the differences between the groups under study. We set the significance threshold at 5% and the Log LDA cutoffs at ±0.5. With these stringent filters, we observed the first significant differences at the Phylum level. Indeed, the Bacteroidota (LDA score = 5.61; p = 0.013) and the Desulfobacterota (LDA score = 2.8; p = 0.002) were more abundant in the Sa group as compared to the others. Conversely, the Phylum of the Proteobacteria was enriched in the Ct group (LDA score = 5.54; p = 0.042; Figure 3). At the Class level, the Bacteroidia (LDA score = 5.61; p = 0.013) and the Desulfovibrionia (LDA score = 2.8; p = 0.003) were more abundant in the Sa group whereas the Gammaproteobacteria was enriched in the Ct group (LDA score = 5.54; p = 0.042; Figure 3). In details, five Orders (four in the Sa and one in the Ct group, respectively) and eight Families (six in the Sa and two in the Ct group, respectively) were differentially abundant after the LEfSe analysis (Figure 3).

Furthermore, we observed that 13 genera were responsible for the differences between the groups (Figures 4A,B). In particular, the Prevotella (LDA score = 5.56; p = 0.018), the Dubosiella (LDA score = 3.68; p = 0.042) and the Family XIII AD3011 group (LDA score = 3.28; p = 0.010) were the most abundant genera in the Sa group. Conversely, the Neisseria (LDA score = 5.34; p = 0.014) and the Abiotrophia (LDA score = 3.98; p = 0.047) were more abundant in the Ct group.

Four species showed significant results in the LEfSe analysis (Figures 4C,D). Indeed, the Prevotella melaninogenica (LDA score = 5.3; p = 0.041) and the Prevotella nanceiensis (LDA score = 4.4; p = 0.012) were more abundant in the Sa group. The Abiotrophia defectiva (LDA score = 3.98; p = 0.047) was enriched in the Ct group. Of note, the Gemella sanguinis (LDA score = 3.92; p = 0.029) was the only enriched species in the Dm group.

We used PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) to predict the activity of metabolic pathways (KEGG) starting from the significant genera and species obtained by LefSE analysis between the Ct and Sa groups. At the genus level, PICRUSt analysis highlighted the prevalence of the superpathway of fucose and rhamnose degradation (LDA score = 3.59; p = 0.037), the aerobic respiration I (cytochrome c; LDA score = 3.59; p = 0.014), the superpathway of phylloquinol biosynthesis (LDA score = 3.57; p = 0.046), the superpathway of menaquinol-7 biosynthesis (LDA score = 3.55; p = 0.046), the ppGpp biosynthesis (LDA score = 3.55; p = 0.012), the biotin biosynthesis II (LDA score = 3.54; p = 0.013), the 1,4-dihydroxy-2-naphthoate biosynthesis I (LDA score = 3.53; p = 0.046), the 2-methylcitrate cycle II (LDA score = 3.53; p = 0.014), the arginine, ornithine and proline interconversion (LDA score = 3.52; p = 0.029), the anaerobic heme biosynthesis II (LDA score = 3.52; p = 0.012) for the Sa compared to the Ct group. By contrast, the Ct group showed an enrichment of the ubiquinol-7 biosynthesis (LDA score = 2.06; p = 0.016), the superpathway of L-isoleucine biosynthesis I (LDA score = 2.4; p = 0.014), the NAD salvage pathway II (LDA score = 2.4; p = 0.018), the superpathway of menaquinol-8 biosynthesis I (LDA score = 2.41; p = 0.046), the pentose phosphate pathway (LDA score = 2.43; p = 0.011), the myo-inositol degradation I (LDA score = 2.73; p = 0.02), the superpathway of menaquinol-12 biosynthesis (LDA score = 2.96; p = 0.046; Figure 5A).

At the species level, we found that only the Sa showed enriched pathways, compared to Ct group. Indeed, we found enrichment of the superpathway of N-acetylglucosamine, N-acetylmannosamine and N-acetylneuraminate degradation (LDA score = 2.54; p = 0.011), the L-methionine biosynthesis I (LDA score = 2.5; p = 0.017), the L-methionine biosynthesis III (LDA score = 2.37; p = 0.014), the superpathway of S-adenosyl-L-methionine biosynthesis (LDA score = 2.31; p = 0.023), the L-lysine biosynthesis I (LDA score = 2.3; p = 0.014), the homolactic fermentation (LDA score = 2.25; p = 0.048), the glycogen biosynthesis I (from ADP-D-Glucose; LDA score = 2.18; p = 0.011; Figure 5B).

Our study presents several strengths: the design and implementation of a structured and controlled training model in primary schoolchildren allowed us to obtain information, for the first time, on the effect of this type of exercise on the salivary microbiota in children. The homogeneity of age-matched schoolchildren belonging from the same geographical area, with similar eating habits, do not assuming drugs in days prior the analysis and with a good oral hygiene, enabled us to highlight the effect of physical exercise on the composition of the salivary microbiota, minimizing the influence of confounding variables (63, 81, 100).

However, this study presents some limitations: first, the limited number of children tested. Due to the small sample size, we were unable to perform a gender-specific analysis of saliva microbiota composition, which may introduce bias, as the extent to which gender influences the diversity and abundance of the salivary microbiota remains unclear (84, 101). We aim in future studies, to increase the sample size and provide a gender-based analysis of saliva microbiota composition.