One hundred and thirty children (8–10 years) belonging to five primary schools in the neighborhoods of Turin (northwest Italy) were enrolled. All information on the aim of the study has been provided to children’s parents/guardians and teachers as previously described (27).

Children meeting any of the following criteria were excluded from the study: (i) recent infections (1 month prior to sample collection), (ii) having disorders affecting diet or physical activity, and (iii) recent usage of either antibiotic, prebiotics and probiotics supplements (1 month prior to data and sample collection). The enrolled children were classified in two groups: active (A) and sedentary (S) on the basis of Physical Activity Questionnaire for Older Children (PAQ-C-It cut-off score of 2.75), the related procedures are detailed in Lupo et al. (27). Parents/guardians and teachers provided written informed consent for participation to the study, according to the ethical standards provided in the 1964 Declaration of Helsinki. Ethics committee on human research of the University of Turin (9 March 2020: Protocol #134691) and Naples (17 January 2020: Protocol #376/19) approved the study. The procedures used to take anthropometric measures were described in Lupo et al. (27); briefly, stature was measured by a portable stadiometer (Model 214; Seca, Hamburg, Germany), body mass was measured by an electronic scale (Model 876; Seca, Hamburg, Germany), participants’ waist circumference was measured in the standing position, midway between the lowest rib and the iliac crest by Ana elastic meter. The Body Mass Index (BMI) was calculated as body mass divided by height squared (kg/m²).

To estimate the daily food intake, all participants filled 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 a one-way ANOVA (Statview software).

The donor was asked not to eat and not to use oral hygiene products 1 h before saliva collection. At least 2 mL of unstimulated saliva was collected, put on ice and stored at − 80°C until the analysis. DNA was extracted from saliva samples using the MagPurix Bacterial DNA Extraction Kit (ZP02006; Zinexts Life Science Corp.) according to the manufacturer’s instructions. DNA was quantified using the Qubit dsDNA BR and HS assay kit (Life Technologies, CA, United States).

PCR amplification was conducted to target the hypervariable V3–V4 regions of the 16S rRNA bacterial genes. Specific primers with barcodes and high-efficiency enzymes were used to perform PCR. The PCR primers were: forward 341F: CCTAYGGGRBGCASCAG; reverse 806R: GGACTACNNGGGTATCTAAT. The PCR products of 450–500 bp were collected with 2% agarose gel electrophoresis. To build library, same amount of PCR products from each sample is pooled, end-repaired, A-tailed and further ligated with Illumina adapters. The library QC was performed with Qubit and real-time PCR for quantification and with bioanalyzer to check the insert size distribution. Libraries were sequenced on a paired-end Illumina platform to generate 250 bp paired-end raw reads. The raw sequencing data are available in Zenodo (https://doi.org/10.5281/zenodo.7920752; Publication date: May 10, 2023).

We used the R platform for statistical analysis and for the data processing. We applied the DADA2 workflow (28) to infer the Amplicon Sequence Variants (ASVs) and for the taxonomic assignments.

In brief, we first filtered and trimmed raw sequencing reads in order to remove low quality bases and adapter contamination. Then, we removed identical reads. Moreover, the reads were denoised, merged filtered to remove artifacts (PCR, and PhiX related chimeras). We obtained the ASVs quantifications and assigned taxonomy annotations (including the Species level) using the SILVA database of non-redundant sequences (version: v138, nr99) (29). The data were structured in objects including the ASVs quantifications, the taxonomy annotations, the sample group data and the phylogenetic tree using the phyloseq and the APE packages (30, 31). Finally, based on the initial DNA concentration, we removed possible contaminant ASVs by using the “prevalence” method of the decontam package (32).

Downstream analyses were performed using the MicrobiomeAnalystR package (33, 34) and included data normalization, measures of diversity and differential abundance estimation.

Briefly, we normalized the ASV counts based on their abundance (low count filter: for any ASV to be retained, at least 20% of its values should contain at least 4 counts) and variance (low variance filter: based on Inter-quantile range ± 10%). This, because ASVs with small counts (in few samples) could represent sequencing errors. Moreover, ASVs that are closely constant in all samples could be excluded from the comparative analyses. Finally, we used the total sum scaling in order to bring all the samples to the same scale.

We evaluated the alpha-diversity by calculating the Abundance-based coverage estimator (ACE) and a nonparametric estimator of species richness (Chao1) indices and by the Fisher metrics (to consider both richness and eveness). The degree to which the species composition changes between the two groups (the beta-diversity) was obtained by PCoA (Principal Coordinates Analysis) of the distances calculated with the un-weighted UniFrac method and the statistical significance assessed by the PERMANOVA test.

We also used the rarefaction curves to evaluate whether the samples were sufficiently sampled and sequenced to represent their species richness. We assessed the statistical significance of comparisons between the two groups of samples under study by using the Mann–Whitney test.

The differential abundance was assessed by the LEfSe (Linear Discriminant Analysis Effect Size) algorithm (35) for biomarker discovery and interpretation of metagenomics data. It involves the Kruskal-Wallis rank sum test to identify features (e.g., Species or Genera) with significant differential abundance in the two groups, followed by linear discriminant analysis (LDA) to evaluate the relevance (the effect size) of the selected features. Different abundant features were considered if the FDR adjusted value of p was less than or equal to 0.05 and if the Log LDA was greater than or less than 0.5. We used Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) to assess the metabolic potential of microbial communities (KEGG pathways). In this analysis, we started from the ASVs belonging to the significant genera obtained by LefSE algorithm.

Anthropometric characteristics and eating habits of the children enrolled in this study are shown in Table 1. No significant differences in individuals’ gender (Chi-square 0.6748; p = 0.411399). anthropometric characteristics such as height, weight, BMI and the waist/height ratio were observed (p > 0.05). Similarly, no significant differences in all diet components analyzed in A and S groups were found (Table 1).

The Illumina sequencing of the hypervariable V3–V4 regions of the 16S rRNA bacterial genes generated 2 × 250 bp paired-end reads. On average, we obtained 169.124 reads per-sample (Supplementary Figure S1A). Overall, the percentage of bases with quality scores above 20 and 30 (Q20 and Q30, respectively) was of 96.45 and 91.22, respectively (Supplementary Figures S1B,C). The percentage of GC nucleotides was of 52.03 (Supplementary Figure S1D). The set of reads was used to run the DADA2 workflow including the filtering and trimming (median = 169.02), denoising of forward (median = 165.00) and reverse (median = 164.80) reads; merging (median = 149.27) and chimeric reads removal (median = 114.43; Supplementary Figure S1E). After merging, the median length of reads was of 424 bp. Overall, starting from filtered reads, we obtained a merging rate of 88.38% and a final rate of read processing (non-chimera over merged reads) of 67.62% (Supplementary Figure S1F). Details on read processing are reported in Supplementary Table S1.

For taxonomic assignments we used the SILVA database of non-redundant sequences (version: v138, nr99). Overall, we could identify a total of 14.197 taxa (ASVs) that were annotated to the seven taxonomic ranks as follows. All the ASVs were taxonomically assigned to the kingdom of bacteria. The 98.06% of the ASVs was annotated at the phylum level (43 phyla), the 96.97% at the class level (94 classes), the 94.93% at the order level (195 orders), the 88.43% at the family level (231 families), the 81.04% at the genus level (404 genera) and the 4.86% was annotated up to the species level (257 species; Supplementary Figure S2A). From the initial set of annotated ASVs, we discarded a total of 68 taxa as possible contaminants. Moreover, as described in Methods we removed low abundant and low variable ASVs to obtain the final set of 472 ASVs that was normalized and used for downstream analyses.

At the genus level, alpha-diversity estimates were significantly different between the two groups, A and S. Indeed, ACE (Figure 1A), Chao1 (Figure 1B) and Fisher (Figure 1C) indices showed p < 0.01 (Mann–Whitney test). The beta-diversity analysis, as measured by unweighted UniFrac distances, showed a significantly different microbial composition between the two groups (r² = 0.026, p = 0.001; Figure 1D), the diversity was also confirmed by using the weighted UniFrac distance metric (r² = 0.018, p = 0.001; Figure 1E). As previously reported (36), unweighted and weighted UniFrac distance measures can be considered as quality-based and quantity-based indexes, respectively. Indeed, we can assess that the observed variation between the two groups was due to the different taxa abundances and to the types of taxa in their microbiome. Moreover, the rarefaction analysis clearly evidenced the capacity to capture the species richness from the results of sampling and sequencing in both groups without any statistically significant difference (p = 0.084; Supplementary Figure S2B).

We evaluated and compared the taxa abundance in the final set of 472 filtered and normalized ASVs. Overall, we identified 8 phyla, 12 classes, 30 orders, 46 families, 84 genera and 96 species.

At the phylum level, on average, the most abundant bacteria were Firmicutes, Bacteroides and Proteobacteria accounting for the 32.08, 26.97% and the 25.58% of the taxa, respectively (Supplementary Figure S3A). The most represented classes were Bacteroides (26.97%), Gammaproteobacteria (25.58%) and Bacilli (21.36%; Supplementary Figure S3B). The most prevalent orders were Bacteroides (26.63%), Lactobacillales (19.58%), Pasteurellales (12.88%) and Burkholderiales (12.12%; Supplementary Figure S3C). Among the most abundant families, we found Prevotellaceae (23.07%), Streptococcaceae (17.66%), Pasteurellaceae (12.88%) and Neisseriaceae (11.77%; Supplementary Figure S3D). At the genus level, we found the Prevotella (19.74%), Streptococcus (17.66%), Haemophilus (11.87%), Neisseria (11.57%) and Veillonella (6.06%; Figure 2A). Finally, the top abundant species that we were able to classify were Prevotella melaninogenica (10.86%), Fusobacterium periodonticum (4.44%), Haemophilus parainfluenzae (2.27%), Rothia mucilaginosa (2.10%) and Veillonella dispar (1.55%; Figure 2B).

As described in Methods, we used the LEfSe algorithm to perform the differential abundance analysis and to identify the taxa that could explain the differences between the two groups A and S. We set the threshold at 5% and the Log LDA cutoff at ± 0.5. Interestingly, using these very stringent criteria, we found Coriobacteriaceae family as more abundant in the A compared to the S group (LDA score = −3.17; FDR = 0.021).

Further, at the genus level, we found that ten genera were responsible for the differences between the two groups. In particular, Agathobacter (LDA score = −3.40; FDR = 0.015), Escherichia–Shigella (LDA score = −3.37; FDR = 7.68×10-04), Collinsella (LDA score = −3.17; FDR = 0.012), Simonsiella (LDA score = −2.95; FDR = 0.044), Eubacterium-yurii group (LDA score = −2.79; FDR = 0.041) and Parabacteroides (LDA score = −2.43; FDR = 0.015) were more abundant in the A group. On the contrary, Mogibacterium (LDA score = 2.71; FDR = 9.21×10–04), Stomatobaculum (LDA score = 3.24; FDR = 0.44), TM7× (also known as Nanosynbacter lyticus, LDA score = 3.90; FDR = 0.045) and Granulicatella (LDA score = 4.14; FDR = 0.045) were more abundant genera in the S group (Figures 3A,B; Supplementary Table S2). Eight species showed significant differences in the LEfSe analysis. Indeed, Prevotella nigrescens (LDA score = −3.76; FDR = 1.5 × 10–03), Collinsella aerofaciens (LDA score = −3.17; FDR = 7.45 × 10–03), Simonsiella muelleri (LDA score = −2.96; FDR = 2.76 × 10–05), Parabacteroides merdae (LDA score = −2.43; FDR = 1.3 × 10–02), were the most represented species in the A group. Conversely, Gemella parahaemolysans, Prevotella aurantiaca (LDA score = −3.9; FDR = 5.27 × 10–04), Prevotella pallens (LDA score = 4.23; FDR = 1.93 × 10–02), Neisseria mucosa (LDA score = 4.43; FDR = 1.31 × 10–02) were more abundant species in the S group (LDA score = 2.94; FDR = 7.45×10-03; Supplementary Table S3; Figures 3C,D).

PICRUSt analysis highlighted the predominance of super pathway of hexitol degradation, L-glutamate degradation VII (to propionate), 2-methylcitrate cycle II, tetrapyrrole biosynthesis I, L-histidine degradation II, superpathway of beta-D-gluconide and D-gluconate degradation, biotin biosynthesis I, and L-arginine biosynthesis pathways activation in saliva from S compared to A children. Conversely, we found a prevalence of superpathway of fatty acid biosynthesis initiation (E. coli) and catechol degradation II (meta-cleavage pathway) in saliva from A respect to S children (Figure 4).