Between November to December 2020, a field experiment was set up at three nearby farms in the Copperbelt Province, Zambia (Latitude − 12.9833⁰S, longitude 28.6333⁰E). The sites were located within a radius of approximately two Km for logistical convenience. Each site had unique in-house mabisi processing conditions; site 1 involved incubation in an outdoor and elevated space; site 2 involved incubation on a wooden platform in an enclosure; and site 3 involved incubation on the floor, but also in an enclosure. The study utilized five substrates: raw cow milk (RCM) serving as a control (Schoustra et al., 2013), ultra-high temperature low-fat milk (LFM), ultra-high temperature full-cream milk (FCM), F100 infant formula (F100), and S26 infant formula (S26), representing novel environments differing in pretreatment and nutritional parameters (Table 1) (Kunda et al., 2015; Nyirenda et al., 2009; Park et al., 2012). The RCM was sourced and homogenized from three local dairy farmers. The F100 was prepared according to a World Health Organization (WHO) guided in-house protocol (World Health Organisation, 1999), while S26, LFM, and FCM were sourced ready for use from a local supermarket. The milk substrates were inoculated with a pooled mabisi sample, representing a source community, in the ratio of 1:10, in a total volume of 1 L. There were three replicates for each milk substrate at each farm site. The preparation was left to ferment for a complete cycle of 48 hrs. Analogous to an evolution experiment, a portion of the fermented products was serially transferred into a fresh set of its respective substrates in the same ratio over time: designated arbitrarily as early propagation phase (1st – 3rd cycle) and late propagation phase (4th – 10th cycle) for a total of 10 cycles (~ 66 bacterial generations). Culture products were sampled in 50 mL clean Falcon tubes at the end of each cycle, placed in a cooler box with ice, and transported to the central laboratory at the National Health Research and Training Institute (NHRTI) (formerly the Tropical Diseases Research Centre) for further analysis and storage at -20°C. A summary of the experimental setup is depicted above (Fig. 1).

[Please insert Table 1 here]

For bacterial community profiling, mabisi samples from the 1st, 3rd, 6^th, and 10th cycles of fermentation were selected. The DNA extraction was performed according to the previously described protocol (Schoustra et al., 2013), with minor modifications. Briefly, 1 mL mabisi sample was centrifuged in a 2 mL screw capped tube at 12,000 rpm for 2 min. The pellets were resuspended in a cell digestion solution containing 64 µl of a 0.5 M EDTA, 160 µl nuclei lysis reagent (Promega), 5 µl RNase (10 mg/mL), 120 µl of lysozyme (10 mg/mL), and 40 µl pronase E (20 mg/mL) reagents. The mixture was incubated automated heating block for 60 min at 37°C while shaking at 350 rpm. The tubes were subjected to bead beating with 2 scoops of sand-sized beads (Sigma, Germany) following an in-house Lactoccocus lysis protocol, while cooling on ice in between. After, 400 µl of ice-cold 0.5 M ammonium acetate (Sigma Aldrich) was added, mixed, and incubated for 15 min at room temperature. The tubes were centrifuged at 13,000 rpm for 4 min. Later, 700 µl of the supernatant was transferred into a sterile 1.5 mL capped cryovial, and an equal volume of molecular grade phenol pH 8.0 (Sigma Aldrich) was added. The tubes were spun at 12,000 rpm for 6 min, and 350 µl of the supernatant was transferred into a new 1.5 mL cryovial tube. An equal volume of chloroform (Sigma Aldrich) was added and centrifuged at 12000 rpm for 2 min to remove the phenol. Thereafter, 300 µl of the supernatant was transferred into a new 1.5 mL cryovial tube to which 400 µl of isopropyl alcohol (Sigma Aldrich) was added. The tubes were placed in a -20°C freezer overnight to precipitate the DNA. After, the tubes were centrifuged for 15 min at 13,000 rpm at 4°C and the DNA pellets were washed with 1ml ice-cold 70% ethanol by carefully inverting the tubes 10 times, followed by centrifugation at 12,000 rpm at 4°C for 10 min. The wash was repeated. The supernatant was carefully decanted, and the tubes were left to dry for 5 min the heating block at 37°C. After, the DNA was eluted by adding 20 µl of low EDTA elution buffer (10 mM Tris, bring to pH 8.0 with HCL; 1 mM EDTA) (Sigma Aldrich). A NanoDrop™ ND-2000 and Qubit ^TM 4 fluorometer (Thermal Fisher Scientific, UK) was used to check for the DNA quantity and quality and stored at -20°C until further needed.

To profile the bacterial community structure and their temporal dynmics following the experimental treatments described earlier, the extracted DNA samples were sent for 16S rDNA sequencing to Novogene, United Kingdom. Polymerase chain reaction (PCR) was used for library preparation using 341F CCTAYGGGRBGCASCAG and 806R GGACTACNNGGGTATCTAAT universal primers targeting the 16S rDNA V3-V4 region, according to the previously described protocol (Schoustra et al., 2013), with modifications. Briefly, the PCR products were purified, end-repaired, A-tailed and ligated with Illumina adapters, and sequenced on the NovaSeq PE250 platform to generate 250 bp paired end raw reads. The barcode and primer sequences were truncated and FLASH (v 1.2.11) applied to merge the reads. Then, the reads underwent filtering, denoising, removal of chimera, and generation of amplicon sequence variants (ASV) with DADA2 R package along with annotation using a publicly available Silver database: silva_nr99_v138.2_toGenus_trainset.fa.gz (Wambua, 2025). Then, the ASVs were normalized by rarefying with random repeated sampling of samples to a minimum sequence reads of 19,534 without replacement, nonbacterial taxa excluded, and a cut-off relative abundance of 0.25% applied in at least each sample to exclude singleton or spurious taxa before downstream analysis (Reitmeier et al., 2021), with R (version 4.5.0) and phyloseq package (version 1.52.0).

Volatile organic compounds (VOC) were also used as a proxy measure for the metabolic activity following propagation of mabisi microbial communities in varied milk substrates. The VOC were analyzed by Headspace-Solid Phase Microextraction Gas Chromatography-Mass Spectrometry (HS-SPME, GC-MS), Trace 1300 Gas Chromatograph (Thermo Fisher), TriPlus RSH autosampler (Thermo Fisher) and an ISQ QD mass spectrometer (Thermo Fisher) using an inhouse protocol, as previously described (Moonga et al., 2021). Briefly, 1mL of mabisi sample was injected into GC-MS labeled vials with caps fastened, and placed on the GC-MS platform to incubate at 60°C for 20 min. The VOCs were allowed to vaporize and adsorb on the SPME fiber (Car/DVB/PDMS/Supelco) at 60°C for 20 min. Next, the extracted VOCs were desorbed for 2 min onto a Stabilwax ®-DA column (30 m length, 0.25 mm ID, 0.5 µm df, Restek), PTV split-less mode at 250°C for 5 minutes with helium gas as carrier at 1.5 ml/min. The GC oven was set at 40°C for 2 min, raised to 240°C with a slope of 10°C/min, and kept at 240°C for 5 min. Mass spectra data were collected over a range of m/z 33–250 in full scan mode with 3.0030 scans/sec. The obtained data were analyzed by Chromeleon ® 7.2 software (Thermo Fisher) using the ICIS algorithm and the NIST main library for signal peak integration and compound hit annotation. The RSI match factor of > 750 was applied. The results were exported to Excel files for statistical analysis. Nine samples, spanning substrate type, propagation phases, and sites, showed no detectable values for all the tested VOCs and were therefore considered technical failures and excluded from downstream analysis.

The fermentation process was monitored by recording changes in pH at the initial time point, after 24 hrs and 48 hrs, respectively, using a digital pH meter (Jenway, UK). The pH probe was disinfected with 70% ethanol and rinsed in sterile water before and in-between measurements of each sample.

The R statistical software (version 4.5.0) was used for all data analysis. The microbiome data were visualized by box plots, histographs, and non-metric multidimensional scaling plots. Furthermore, alpha (Chao1 and Shannon) and beta diversity metrics were applied with Wilcoxon rank-sum and PERMANOVA (ANOSIM) tests, corrected by the Benjamin-Hochberg method, to determine how microbial community diversity varied by study parameters. Species differential analysis was further explored through a linear discriminant effective size (LEfSe) analysis using the microbiomeMarker package (version 1.13.2).

The VOC data was plotted with heatmaps and principal component analysis (PCA) for visualization and analyzed by multivariate analysis (ADONIS) with default free permutation of 999 after normalization using a log transformation and median scaling by compound (column). The pH and consistency were summarized by median, lower and upper quartiles, and visualized by violin plots. Kruskal-Wallis and Dunn tests, with the Benjamin-Hochberg adjusted method, were applied whenever appropriate to analyze the statistical differences between experimental parameters. The significance test was set at 0.05 for all analyses.

A total of 4322 bacterial ASVs were obtained. At the phylum level of classification, these taxa were predominantly represented by the Bacillota, and this pattern was consistent with the starting mabisi microbial community (Figure S1). However, at the genus level, and specifically focusing on key lactic acid bacteria (LAB) and acetic acid (AAB) community members known for their important role in dairy ecosystems, the members, including Lactococcus, Acetobacter, Leuconostoc, Lactiplantibacillus, Paucilactobacillus, and Lacticaseibacillus were dominant in all substrate treatments at all farm sites over time (Fig. 2; Fig. S2).

[Please insert Fig. 2 here]

Figure 2 Distribution of propagated mabisi microbial community in different milk substrates over time. The y-axis represents the relative abundance (%) of each taxon in each sample, while the x-axis shows sample identities (Sample ID), faceted by propagation phases: starting mabisi microbial community (O), early phase (early_phase), and late phase (late_phase). Substrates are represented by raw cow milk (RCM), F100 infant formula (F100), S26 infant formula (S26), ultra-high temperature low-fat milk (LFM), and ultra-high temperature full-cream milk (FCM). The legend lists the most abundant genera; Acetobacter (dark teal), Clostridium (brownish orange), Enterobacter (warm terracotta), Enterococcus (soft lavender purple), Klebsiella (magenta), Lacticaseibacillus (mutated brick red), Lactiplantibacillus (bright olive green), Lactococcus (golden yellow), Leuconostoc (dark mustard yellow), Paucilactobacillus (earthly brown), and unassigned - referring to ASVs that could not be annotated in the silva database (medium gray)

To gain further insights into whether propagation of a shared starting mabisi microbial communities in varied substrates exerted diversity changes in the microbial community composition at different farm sites over time, alpha and beta diversity analysis was applied. The alpha diversity was analyzed by Chao1 (richness) and Shannon (richness and evenness) indices, with raw cow milk (RCM) used as a reference (the usual substrate for mabisi microbial community). Relative to the RCM, all substrate treatments did not show significant differences in the resultant community richness (Fig. 3; Table S1) or across sites (Fig. S3; Table S1). Also, the richness did not differ between the early and late propagation phases (Fig. 3; Table S1). However, both the richness and evenness of the propagated microbial community from FCM and S26-based substrate treatments increased significantly (p < 0.05), but not in the F100 and LFM-based treatments (p > 0.05) (Fig. 3; Table S1) relative to RCM. Between the early and late propagation phases or sites, the richness and evenness did not differ (p > 0.05) (Fig. S3; Table S1). Overall, this suggests that substrate variation constituted niche division, shaping differences in the community alpha diversity, while propagation phases or site differences did not.

To assess the dissimilarity in the community composition following the repeated propagation of mabisi microbial communities across substrates and sites over time, a non-metric multidimensional scaling (NMDS) with permutation of multivariate analysis (PERMANOVA) was conducted. The microbial community significantly differed in its composition mediated by the substrate variation (ANOSIM: R = 0.16, p < 0.05), propagation phases (ANOSIM: R = 0.40, p < 0.05), and site differences (ANOSIM: R = 0.03, p < 0.05) (Fig. 4, Table S2).

[Please insert Fig. 3 here]

Figure 3 Alpha diversity of mabisi microbial community by substrate treatment over time. Chao1 (richness) and Shannon (incorporating both richness and evenness) diversity metrics were applied to assess microbial community alpha diversity following propagation of a shared starting mabisi microbial community in varied milk substrates. The y-axis represents the diversity indices (Chao1 and Shannon), while the x-axis represents propagation phases: starting mabisi microbial community (O), early phase (early_phase), and late phase (late_phase). Box plots display the median, interquartile range, whiskers, and outliers. The legend indicates the substrate treatment: F100: F100 infant formula (greenish-blue), FCM: ultra-high temperature full-cream milk (redish-orange), LFM: ultra-high temperature low-fat milk (bluish-purple), O: starting mabisi microbial community (pinkish-purple), RCM: raw cow milk (yellow-green), and S26: S26 infant formula (mustard yellow)

[Please insert Fig. 4 here]

Figure 4 Dissimilarity in microbial communities by non-metric multidimensional scaling (NMDS). The NMDS was applied to show the shifts in microbial community diversity following the propagation of mabisi microbiota through varied milk substrates over time. Substrates are distinguished by colors - F100: F100 infant formula (teal), FCM: ultra-high temperature full-cream milk (orange), LFM: ultra-high temperature low-fat milk (purple), O: starting mabisi microbial community (magenta), RCM: raw cow milk (green), S26: S26 infant formula (yellow). Propagation phases are indicated by shape – early_phase: early phase (filled circle), late_phase: late phase (filled triangle point up), and O: starting mabisi microbial community (filled square)

Further analysis by linear discriminant effect size (LEfSe), focusing on LAB and AAB that are known for their key roles in dairy ecosystems, revealed that there were significant differences in the relative enrichment of the mabisi microbial community between substrate treatment and propagation phases. Notably, LDA analysis showed that members of the Lactococcus (LDA = 4.9, p < 0.05 ) were significantly enriched in FCM, Acetobacter (LDA = 4.7, p < 0.05) were significantly enriched in LFM, while Lactiplantibacillus and Leuconostoc (LDA = 4.3, p < 0.05) were significantly enriched in S26 treatment (Fig. 5, Table S3). Similarly, members of the Lactococcus (LDA = 4.9, p < 0.05 ) were significantly enriched during the early phase, while Acetobacter (LDA = 4.5, p < 0.05 ), Lactiplantibacillus (LDA = 4.3, p < 0.05 ), Paucilactobacillus (LDA = 4.3, p < 0.05 ), Leuconostoc (LDA = 4.2, p < 0.05 ), Lacticaseibacillus (LDA = 4.2, p < 0.05 ) and Lactobacillus (LDA = 3.7, p < 0.05 ) were members which were significantly enriched during the late phases of propagation (Fig. 5, Table S3). However, only taxa other than those belonging to LAB and AAB were differentially enriched in farm site 1 (Fig. S4; Table S3). The differential enrichment outcome supports the observed community diversity changes.

[Please insert Fig. 5 here]

Figure 5 Linear discriminant analysis effect size (LEfSe) to identify differentially abundant microbial taxa following the propagation of mabisi in varied milk substrates over time. The x-axis shows the linear discriminant analysis (LDA) score (log 10 transformed effect size) of enriched taxonomic features according to substrate variation (Fig. 5a) and propagation phase (Fig. 5b). The y-axis lists the enriched taxa at genus level. The legends show the enriched groups: substrates - F100: F100 infant formula (dark teal), FCM: ultra-heat temperature full-cream milk (brownish orange), LFM: ultra-heat temperature low-fat milk (warm terracotta), O: starting mabisi microbial community (soft lavender purple), RCM: raw cow milk (magenta), S26: S26 infant formula (mutated brick red); propagation phases – early_phase: early phase (dark teal), late_phase: late phase (brownish orange), and O: starting mabisi microbial community (warm terracotta)

To understand how the changes in the microbial community diversity influenced the community-level functionality following the repeated propagation of mabisi in varied substrates at different sites over time, volatile organic compounds (VOCs), pH, and consistency were analyzed as proxy indicators.

The obtained VOCs belong to the chemical classes of aldehydes, esters, carboxylic acids, alcohols, and ketones (Table S4). The VOC profiles exhibited differences between samples analyzed before and after the propagation of mabisi across substrates and farm sites over time. Particularly, samples analyzed before propagation displayed VOCs that are associated with photo oxidation and undesirable flavors, including octanoic acid, nonanoic acid, ethyl-9 decenoate, and 1-octen-3-one, among others (Fig. S5). Therefore, these specific samples were removed from downstream analysis.

Further analysis of propagated samples through principal component analysis (PCA) and PERMANOVA revealed that VOCs were significantly driven by both the substrate treatment, which separated along the PC1 axis explaining 37.4% of the variation, and the propagation phase, which separated along the PC2 axis explaining 22.4% of the variation (p < 0.05) (Fig. 6; Table S5). Specifically, the VOCs in the F100 and S26 treatments, exhibited higher levels of Ethyl acetate, Hexanal, and 1 Butanol 3-Methyl, which separated from those belonging to LFM and FCM, which showed higher levels of 2-Heptanone, 2-Nonanone, Butanoic acid, and 2-Butanone, Octanoic acid, and Hexanoic acid, whereas, those from RCM appeared to contribute shared VOC from either group (Fig. S6). However, the farm site did not impact any significant separation of VOCs (p > 0.05) (Fig. S7; Table S5).

[Please insert Fig. 6 here]

Figure 6 Volatile organic compounds from mabisi samples after propagation of mabisi microbiota in varied milk substrates over time. Each point represents an individual sample. Legends show substrates distinguished by color – starter: starting mabisi microbial community (orange), RCM: raw cow milk (sky blue), F100: F100 infant formula (green), S26: S26 infant formula (light pink), LFM: ultra-high temperature low-fat UHT (dark blue), and FCM: ultra-high temperature full-cream milk (dark orange). Propagation phases are represented by ellipses at 95% confidence interval – starter: starting mabisi microbial community (solid ellipse not formed due to inadequate data points), early_phase: early phase (dotted ellipse) and late_phase: late phase (dashed ellipse)

Following propagation, the pH dropped in all substrates to a range between 3.98 and 4.51 (Fig. 7a) from a range between 5.69 to 6.77 before propagation (Fig. S8a). When compared to the RCM reference substrate, the median pH for LFM, FCM, and S26-based substrate treatments was significantly lower (p < 0.05), while that for F100 was not significantly different following the propagation process (p > 0.05) (Fig. 7a; Table S6). Further analysis revealed that the pH was significantly lower for propagated mabisi in the late propagation phase compared to the early phase (p < 0.05) (Fig. S8a; Table S7), and was also significantly lower for samples from farm 2 compared to farm 1 and 3 (p < 0.05) (Fig. S8b; Table S8). Thus, the pH was driven by differences in the substrates, propagation, and farm sites.

In addition, the substrate treatments also influenced the consistency outcomes. The consistency ranged between 1.5 to 4.11cm for all substrates following propagation of mabisi in varied substrates, with higher values observed in S26 and lower values observed in F100 treatments (Fig. 7b). Relative to the raw cow milk treatments, the median consistency was significantly lower for F100, FCM, and LFM (p < 0.05) while that for S26-based treatments was not statistically different (p > 0.05) (Fig. 7b; Table S9). However, consistency was not significantly different between early and late propagation phases (p > 0.05) (Fig. 7b; Table S9) or between sites (p > 0.05) (Fig. S8c; Table S9). These outcomes suggest that substrate variation was the key determinant for driving consistency patterns.

[Please insert Fig. 7 here]

Figure 7 The outcome of pH (Fig. 7a) and consistency (Fig. 7b) following repeated propagation of mabisi in varied milk substrates over time. The pH values are shown as medians, with interquartile ranges for each substrate treatment per propagation phase. The x-axis shows the propagation phases: early phase (early_phase) and late phase (late_phase), while the y-axis shows the pH and consistency (cm) measurements. Legends show the substrate treatment - F100: F100 infant formula (light teal), FCM: ultra-heat temperature full-cream milk (light yellow), LFM: ultra-high temperature low-fat milk (lavender), RCM: raw cow milk (coral), and S26: S26 infant formula (light blue). The dotted horizontal line represents the median value of pH (Fig. 7a) and consistency (Fig. 7b) for the reference substrate (RCM)

The datasets generated during the current study are available online in the 4TU.ResearchData repository: https://doi.org/10.4121/82fa6380-c40c-4c40-b5ec-933b0253e6ce