Raw data processing phase

After applying the amplicon workflow (filtering, trimming, and taxonomy assignment) to the raw sequencing data from the three datasets using DADA2, we obtained the corresponding ASVs. Table 1 summarizes the technical specifications of each dataset: the 16S rRNA target regions, the characteristics of the raw reads, the trimming/truncation parameters applied, the lengths of the resulting ASVs, and the total number of ASVs generated for each dataset. The trimming parameter applied for PRJDB7295 means that all reads are trimmed of the first 10 bp (trimLeft = 10), forward reads are truncated at 250 bp, and reverse reads are truncated at 220 bp (truncLen = c(250,220)).

Table 1.

Summary of datasets characteristics: NCBI accession numbers, targeted 16S rRNA gene regions, approximate raw read lengths, final ASV lengths resulting from DADA2 processing, trimming/truncation parameters, and the total number of ASVs generated on each dataset.

* Discovery dataset.

Because our methodology focuses on finding a signature (a sequence contained within another) that is present in other databases to validate its effectiveness in distinguishing between groups (cases versus controls) through classification algorithms, we selected PRJNA633365 as the discovery dataset based on the eligibility criteria which states: “the discovery dataset is the one that contains the shortest sequence length after the raw data processing phase” (Rojas-Velazquez et al., Reference Rojas-Velazquez, Kidwai, Kraneveld, Tonda, Oberski, Garssen and Lopez-Rincon2024b). It is important to note that PRJNA633365 16S rRNA V4, while PRJDB7295 and PRJNA562650 are V3-V4. This means that any feature found in PRJNA633365 can be potentially identified in the V4 region of the other two datasets, making the length criterion secondary compared to the region they share. PRJDB7295 and PRJNA562650 were selected as testing datasets to validate the set of features identified in the discovery phase.

Feature selection

After running REFS on the discovery dataset, the resulting reduced set of features contained 16 out of the original 1,227 ASVs. This means that REFS achieved its highest accuracy ( $ >0.90 $ ) with the minimum number of features (16 features), Figure 2A. After running the validation module for the resulting set of features, the Multilayer Perceptron (MLP) classifier had the best performance, with an AUC-ROC of $ 0.93 $ , Figure 2B. In the case of SelectKbest, k = 16, the resulting set of features selected by SelectKbest has 3 in common with those selected by REFS: Clostridium (feature 6), Clostridium sensu stricto 13 (feature 11), and Dysgonomonas (feature 15). The classifier with the best performance was Extra Trees with an AUC-ROC of 0.84. According to Šimundić (Reference Šimundić2009), the AUC-ROC of 0.93 is considered as “excellent” diagnostic accuracy.

Figure 2.

(A) Feature selection phase: the reduced set of features that achieves the highest accuracy is identified by the red line; (B) AUC-ROC for the classifier with the best performance in the validation module applied to the discovery dataset PRJNA633365.

Additionally to the validation module, we executed MaAsLin2 (Mallick et al., Reference Mallick, Rahnavard, McIver, Ma, Zhang, Nguyen, Tickle, Weingart, Ren, Schwager, Chatterjee, Thompson, Wilkinson, Subramanian, Lu, Waldron, Paulson, Franzosa, Bravo and Huttenhower2021) with a linear model on the resulting reduced set and simultaneously performed a Multivariate Analysis of Variance (MANOVA) model (IBM Corporation, 2011). The results from MaAsLin2 indicate that from the 16 bacterial taxa identified with our methodology, 10 are the most significant: g.Dialister, g.Bifidobacterium, g.Clostridium, g.Clostridium sensu stricto 13, g.Dysgonomonas, g.Erysipelatoclostridium, g.Clostridium sensu stricto 1, g.Clostridium innocuum group, g.Streptococcus, f.Peptostreptococcaceae. In the same way, we executed MaAsLin2 on all ASVs in the discovery dataset and used the top 20 as input in the validation module for additional validation. After this process, the resulting AUC-ROC from the classifier with the best performance (Extra Trees) was 0.71. This AUC-ROC is lower compared with the AUC-ROC of REFS (MLP, 0.93) and SelectKbest (Extra Trees, 0.84). The elements at the genus level shared by the set selected by REFS and the top 20 of MaAsLin2 are Bifidobacterium (feature 2), Dialister (feature 1), Erysipelatoclostridium (feature 12), and Streptococcus (feature 9). The MANOVA model applied to the microbiome taxa relative abundance matrix demonstrates strong statistical significance, indicating that the collective profiles of the selected taxa have a substantial effect in distinguishing breast milk and infant formula groups (p < 0.0001). For detailed information, see Supplementary File 3.

In addition to MANOVA and MaAsLin2, we applied CCREPE (Emma Schwager and Weingart, Reference Emma Schwager and Weingart2021; Schwager and Huttenhower, Reference Schwager and Huttenhower2016) to the selected ASVs to evaluate whether they exhibit meaningful associations rather than being isolated signals. The results confirmed that several pairs have statistically significant relationships after multiple-testing correction, indicating non-random co-occurrence patterns. For example, g.Clostridium sensu stricto 13 (feature 11) and g.Erysipelatoclostridium (feature 13) showed a positive association ( $ p=7.34\times {10}^{-17} $ , $ q=4.03\times {10}^{-14} $ , similarity score = 0.80), g.Bifidobacterium (feature 2) and g.Bifidobacterium (feature 14) also displayed a positive association ( $ p=2.28\times {10}^{-6} $ , $ q=4.17\times {10}^{-4} $ , similarity score = 0.37), and g.Clostridium (feature 6) and g.Klebsiella (feature 16) showed a negative association ( $ p=2.10\times {10}^{-6} $ , $ q=5.77\times {10}^{-4} $ , similarity score = −0.48). These findings confirm that the selected ASVs are embedded in ecological structures, providing biological context and supporting the robustness of the REFS selection. Therefore, the 16 ASVs are not only sufficient for classification but also represent meaningful, non-random relationships that strengthen the interpretability of the results.

The addition of BLAST to the taxonomy assignment allowed us to reach the genus level in 15 of 16 taxa, and in some cases, the species level, so the taxonomy distribution is expressed at the genus level. For taxa that did not have resolution at the genus level, it is expressed as Peptostreptococcaceae;_;_ indicating the family level as the maximum resolution. The visual representation of the taxonomy distribution of the selected ASVs is shown in Figure 3. The taxonomy assignment at all levels, SILVA and BLAST, and the corresponding sequence of the 16 features, as well as which features were found in the testing datasets, are detailed in Supplementary File 1.

Figure 3.

Taxonomy distribution at the genus level. The numbers represent the number of elements contained in each genus-level bacterial taxa.

Testing phase

When searching for the selected 16 ASVs in the testing datasets, we found 13 out of 16 ASVs in PRJDB7295 and 12 out of 16 ASVs in PRJNA562650. After running the validation module on the two testing datasets, the best performing classifiers were the AdaBoost classifier for PRJDB7295 with an AUC-ROC of 0.69 and MLP for PRJNA562650 with an AUC-ROC of 0.92, Figure 4. According to Šimundić (Reference Šimundić2009), the AUC-ROC for PRJDB7295 corresponds to a “sufficient” diagnostic accuracy; for PRJNA562650, the diagnostic accuracy is considered as “excellent.” Although an AUC-ROC value below 0.7 might be considered marginally acceptable (PRJDB7295), it can still indicate a reasonable discriminatory ability for diagnosis (Mandrekar, Reference Mandrekar2010). In the case of SelectKbest, we found 5 out of 16 ASVs in PRJDB7295 and 8 out of 16 ASVs in PRJNA562650, with Extra trees being the classifier with the best performance, with an AUC-ROC of 0.58 and 0.79 respectively. In both testing datasets, two of three common features between REFS and SelectKbest were present.

Figure 4.

AUC-ROC plots corresponding to the classifier with the best performance for the two testing datasets: (A) AdaBoost for PRJDB7295, (B) MLP for PRJNA562650.

We conducted an analysis on the relative abundance of the 16 bacterial taxa by calculating the arithmetic average of the relative abundance of each taxon in the breast-milk and formula-fed groups for each dataset. These averages were visualized in a heatmap to provide an intuitive overview of group-level patterns among the identified taxa, as shown in Rojas-Velazquez et al. (Reference Rojas-Velazquez, Kidwai, Liu, El-Yacoubi, Garssen, Tonda and Lopez-Rincon2025), see Figure 5. The color coding (dark blue = higher in breast milk; light blue = lower) was intended as a descriptive tool to illustrate relative trends, not as a formal statistical test. We acknowledge that microbiome data are compositional and that simple arithmetic means can be misleading for hypothesis testing because of the constant-sum constraint (Gloor et al., Reference Gloor, Macklaim, Pawlowsky-Glahn and Egozcue2017). However, in the context of this manuscript, the heatmap serves as a complementary visualization to highlight the behavior of selected taxa across datasets, rather than to infer significance or effect size.

Figure 5.

Heatmap for a graphic representation of the relative abundance of the selected taxa (increase or decrease) in breast-milk samples compared to formula samples. Genera are displayed from left to right following the top-to-bottom order as presented in Supplementary File 1.

The results indicate that taxa relative abundance varies across datasets, likely due to sequence quality and variations in sequencing equipment, a phenomenon known as the batch effect (Rincon et al., Reference Rincon, Kraneveld and Tonda2020), the geographical origins of the datasets, differences in age at which faecal samples were collected, formula composition, and different diets. For example, three of the identified taxa are consistently increased or decreased across all three datasets where they are observed: Clostridium, Peptostreptococcaceae;_;_, and Erysipelatoclostridium (taxa 6, 10, and 12). In contrast with the four Bifidobacterium (taxa 2, 4, 5, and 14), which actually show a variation in relative abundance with respect to breastfed versus formula-fed, depending on the dataset. This phenomenon is common in this type of data, so REFS incorporates a nested k-fold cross-validation to mitigate these variances in data and avoid biased results and overfitting (Vabalas et al., Reference Vabalas, Gowen, Poliakoff and Casson2019). By design, our methodology prioritizes the selection of potential biomarkers that are consistently predictive across different datasets, ensuring that the selected ASVs are robust to the identified technical and biological noise. For example, while regional dietary changes in China, such as the introduction of solid foods after 6 months, influence the composition of the microbiota (Brink et al., Reference Brink, Mercer, Piccolo, Chintapalli, Elolimy, Bowlin, Matazel, Pack, Adams, Shankar, Badger, Andres and Yeruva2020), the ensemble-based architecture in a nested k-fold cross-validation of REFS allows the model to keep discriminatory efficiency by focusing on a stable biological signature. Thus, the differences in validation accuracy do not merely reflect data inconsistency but rather demonstrate the model’s ability to generalize findings despite the variability in microbiome datasets.

MicrobiomeAnalyst

We used MicrobiomeAnalyst (https://microbiomeanalyst.ca) to examine the relationship between the bacterial taxa, diet, and lifestyle using the host-intrinsic taxon functionality. Host-intrinsic taxon analysis revealed that feeding breast milk is associated with a lower abundance of Peptostreptococcaceae;_;_ and Erysipelatoclostridium (taxa 10 and 12), which is consistent with the heatmap results, Figure 5. Additionally, a higher abundance of Clostridium and Enterobacteriaceae was observed in the host-intrinsic taxon analysis, whereas one of the Clostridium (taxa 6) from the 16 selected bacterial taxa showed higher abundance in the breast-fed group in Figure 5. However, one of the Clostridium (taxa 7) from the 16 selected bacterial taxa and Enterobacteriaceae showed a mixed response in Figure 5. From the analysis, fructose and glucose are associated with an increased abundance of Clostridium (taxa 6, 7, and 8), Klebsiella (taxa 16), and Enterococcus (taxa 3), which is consistent with the heatmap results, Figure 5. Goran et al. (Reference Goran, Martin, Alderete, Fujiwara and Fields2017) show that fructose and lactose are naturally present in breast milk and detected in infants’ gut at 6 months of age.

We also analyzed with the host-intrinsic taxon functionality the selected bacterial taxa to explore their connection with different diseases. The results of the analysis indicated associations with chronic diseases, including cardiovascular diseases, type 2 diabetes, and asthma. In cardiovascular diseases, there was an increased abundance of Erysipelatoclostridium, Bifidobacterium, Peptostreptococcaceae;_;_, and Streptococcus. In allergic diseases, such as eczema (allergic dermatitis), a decreased abundance of Bifidobacterium, Dialister, and Streptococcus was observed. For asthma (children under 1 year of age at risk), a decreased abundance of Bifidobacterium and Streptococcus was also observed.