Influence of metabolic guilds on a temporal scale in an experimental fermented food derived microbial community

Abstract The influence of community diversity, which can be measured at the level of metabolic guilds, on community function is a central question in ecology. Particularly, the long-term temporal dynamic between a community's function and its diversity remains unclear. We investigated the influence of metabolic guild diversity on associated community function by propagating natural microbial communities from a traditionally fermented milk beverage diluted to various levels. Specifically, we assessed the influence of less abundant microbial types, such as yeast, on community functionality and bacterial community compositions over repeated propagation cycles amounting to ∼100 generations. The starting richness of metabolic guilds had a repeatable effect on bacterial community compositions, metabolic profiles, and acidity. The influence of a single metabolic guild, yeast in our study, played a dramatic role on function, but interestingly not on long-term species sorting trajectories of the remaining bacterial community. Our results together suggest an unexpected niche division between yeast and bacterial communities and evidence ecological selection on the microbial communities in our system.


Introduction
What factors affect dynamics of diversity in natural communities and how this links to community function has long been a central question in ecology (Cardinale et al. 2012, Gonzalez et al. 2020 ).Defining community function and community diversity at the le v el of species can take many forms.Community function has been estimated by quantifying parameters such as community pr oductivity, ov er all metabolic output, and stability or resistance a gainst inv asion (Aubr ee et al. 2020 ).Comm unity div ersity ca ptures the living organisms present and their relative abundances.This diversity can be expressed at several levels of taxonomic richness; at the br oader le v el of functional types or metabolic guilds (i.e.capable of performing the same metabolic processes (Allen et al. 2023, Reynolds et al. 2023 )), then genus or species, or more narr owl y at intr aspecific div ersity of genotypes within a species.A gener al positiv e association between higher div ersity and incr eased measur es of comm unity function has been shown (Balv aner a et al. 2006 ), and it is typically interpreted to result from functional complementarity (Tilman 1999, Cardinale et al. 2007 ) of collections of species with differ ent nic hes (Cardinale et al. 2012 ).
Ecological pr ocesses suc h as species sorting-the sorting of variation at the level of species along an ecological or evolutionary timescale-may alter comm unity div ersity.For instance, when a particular species or type is better adapted to the new selecti ve condition, the y may incr ease in r elativ e abundance (Vellend 2010 ).The response to selection may depend on initial species diversity-both in the number of species present and their relati ve abundance-or di versity at the level of metabolic guilds.In this way, species sorting may affect ov er all comm unity functioning, when species or metabolic guilds that increase or decrease in abundance are k e y to overall community function.While bodies of theory exist on the theme of species sorting in communities (Loeuille and Leibold 2008 ), experimental tests are relatively few (Langenheder and Székely 2011, Cairns et al. 2018, 2020, Groenenboom et al. 2022 ).Since most experimental studies on the influence of community or metabolic guild diversity on community function focus on short-term ecological timescales of one or a few gener ations (Wa gg et al. 2014, Sier ocinski et al. 2018, Aubr ee et al. 2020, Gonzalez et al. 2020 ), r elativ el y little is known on how this influence changes over multiple generations (Fiegna et al. 2015 ).
As an experimental model system, the microbial communities of fermented foods provide a po w erful method to study the effects of selection on species or metabolic guild diversity and on community function (Wolfe and Dutton 2015, Wolfe 2018, Aleksee v a et al. 2021, Conacher et al. 2021 ).In these communities, function is often quantified through changes in pH and through metabolic output measured as volatile compound production, contributing to aroma and taste, which can be further connected to known biochemical pathways (Krömer et al. 2004, Smid et al. 2005, De Filippis et al. 2016 ).Metabolic guild diversity can be traced to Figure 1.Experimental design.Dilution a ppr oac h to cr eate initial comm unities (t0), pr ogr essiv el y r emoving r ar e types, arriving at four le v els of initial species diversity (full, medium, and low diversity, and one synthetic community).The synthetic community was created of five single isolates, thus creating the lo w est diversity community.From each diversity level eight replicates were then propagated for 16 transfers of 1% to sterile milk.Image created in BioRender.Abbreviations: med = medium, t0 = transfer 0 (i.e. starting inoculum).
the le v el of micr obial gr oups known to be r esponsible for fermentation, such as lactic acid bacteria, acetic acid bacteria, and alcohol-producing yeast.Metabolic guilds or functional types can be considered equivalent to the classic ecological classification of species into guilds , lifeforms , or strategies (Louca et al. 2018 ).This is the basis of ours and others' terminology of metabolic guilds to r efer to gr oups of micr obes that can perform the same ecological functions (Allen et al. 2023, Reynolds et al. 2023 ).The well-studied microbial metabolic guilds, or fermentative types of fermented foods are therefore relevant communities for studying fundamental questions in ecology.
Specificall y, a tr aditionall y fermented milk be v er a ge fr om Zambia, Mabisi (Moonga et al. 2020 ;Sc houstr a et al. 2013 ), pr ovides a model system to study the relationship between metabolic guild diversity and functioning of the community, and how this may c hange ov er r epeated cycles of sequential pr opa gation.The moder ate div ersity of Mabisi, composed of six to ten dominant lactic and acetic acid bacterial species, plus numerous other low abundance types (other bacterial species , yeasts , and viruses), further facilitates experimental design and analysis (Moonga et al. 2020, Sc houstr a et al. 2013 ).The defined and measurable functional properties of Mabisi, including metabolite profiles and acidity, fr ame inv estigations of the influence between metabolic guild diversity and community function.Four fermentative types ar e natur all y pr esent at v arious le v els of r elativ e abundance.These include: 1 alcohol producers; 2 alcohol consumers producing acetic acid; 3 homofermentative lactic acid producers (only lactic acid produced); and four heterofermentative lactic acid producers (lactic acid, ethanol, acetic acid, and carbon dioxide produced) (Gänzle 2015 ).The presence and ratios of these four fermentative types, or metabolic guilds, is expected to affect comm unity metabolic pr ofiles due to the distinctiv e metabolic ca pabilities each community member possesses.
Her e, we pr esent an experimental test of pr edictions on community functioning and species sorting using microbial com-munities differing in diversity of metabolic guilds , species , and genotypes, whic h wer e r epeatedl y pr opa gated for 16 cycles ( ∼100 generations) in milk in a laboratory environment.Our design was inspired by a prediction of altered community function at v arious le v els of species div ersity, and consequent metabolic guild diversity.Upon diluting a natural, "full" community (10 0 ) to medium (10 −4 ) and low (10 −9 ) diversity levels, then by using a synthetic community of five isolates from the community (Fig. 1 ), we pr ogr essiv el y eliminated r ar e types.Specificall y, y east w as a priori determined to be eliminated in low and synthetic communities as visualized b y ey e and calculated fr om pr e vious abundance estimates (Moonga et al. 2019, Sc houstr a et al. 2013 ).At transfers 1, 5, and 17 we then measured metabolic output as a proxy for community function, and bacterial diversity by full 16S rRNA gene amplicon sequencing.At e v ery second tr ansfer, pH was also measured as indication of community function since acidity str ongl y impacts both consumer pr efer ences (Ott et al . 2000) and stability against pathogen invasion (Mpofu et al. 2016).With this we aimed to answer three main r esearc h questions: 1. does altering initial metabolic guild diversity, with the associated loss of low abundance microbial types, influence community function?2. and if so, are such shifts in community function stable over r epeated pr opa gation?3. does r emoving metabolic guilds influence bacterial species sorting trajectories over repeated cycles of pr opa gation?

Preparing communities (t0)
A fermented Mabisi sample containing its full microbial community was serially diluted up to a factor of 10 −10 using 7.5 ml culture with 67.5 ml UHT full fat milk (Jumbo brand, Houdbare Volle Melk).Only one replicate dilution series was performed, creating one culture bottle for each dilution level.The entire volumes of these dilutions were then fermented for 72 h at 28 • C, then 0.75 ml Yeast and whey production was ther efor e seen only in medium and full communities.
A synthetic community of five individual isolates was also created as a fourth div ersity le v el ( Table S1 ).A Mabisi community that had been diluted to 10 −8 and undergone two 72 h fermentation cycles at 28 • C was diluted and gr own aer obicall y at 28 • C on MRS (de Man, Rogosa, Shar pe) a gar.Initiall y eight unique a ppearing morphotypes (labelled A-H) were selected from the agar plate for in vestigation.T he community was also plated on PC A (plate count agar) and M17 agars, but the MRS provided the clearest and most diverse collection of colony types.Colonies were grown for 5 days in MRS broth and arc hiv ed at −80 • C (1 ml culture + 0.5 ml 85% gl ycer ol).The eight colonies were assessed for their growth in MRS broth, morphotype on MRS agar, API metabolism test (BioMérieux), and ability to acidify milk in isolation.Using the mentioned assessments, a final five colonies (A, B, C, F, and G) were chosen based upon being the most dissimilar (confident at least thr ee wer e unique).Eac h colon y was str eaked twice on MRS a gar and inoculated 5 days growth at 28 • C in MRS broth (1.5 ml broth in 24 well plate).The cultures' optical density was measured (600 nm wav elength r eading) and then combined in volumes containing equal cell densities of each.This mixed culture formed the synthetic T 0 diversity community and was arc hiv ed at −80 • C (1 ml culture + 0.5 ml 85% gl ycer ol).Later Sanger sequencing of the full 16S rRNA gene (27F, 1492R primers) with Blast search in the NCBI database was unable to decipher the colonies to the species le v el (Table 1 ), as top identifications were all a > 99% match.Ho w ever, the 16S rRNA gene sequences indicated that at least three unique types were present-two Acetobacter and one Lactobacillus .

Transferring and arc hi ving
For each diversity level (full, medium, low, and synthetic) eight r eplicate lines wer e inoculated using their r espectiv e T 0 fr ozen cultures (1.5 ml total, with glycerol) in 75 ml of milk (total 32 evolution lines, plus one uninoculated milk as a negative control).Samples underwent repeated 72-h fermentation cycles at 28ºC for 16 transfers ( ∼100 generations).For logistical purposes, after every second transfer final cultures were inoculated into fresh milk and stored at 4ºC for 24 h before being moved to 28ºC.Ther efor e, a "true" c ycle w as completed after e v ery second tr ansfer (i.e .7 da ys).Fresh samples of final products were archived at transfers 1, 2, 3, 5, 11, and 17, with the follo wing ar c hiv ed: with gl ycer ol at −80 ºC (1.27 ml culture + 0.63 ml 85% gl ycer ol), and without gl ycer ol at −20 ºC for DNA analysis and GC-MS analysis ( ∼9 ml).Please note that arc hiv e labelling was done as follows: t5 = final pr oduct fr om t4 (i.e. the Mabisi culture used to inoculate at transfer 5), t17 = final pr oduct fr om t16.Due to an error in transferring, the M8 line (medium div ersity, r eplicate 8) was lost early in the experiment, thus it is excluded from all analyses.
Data were first normalized by compound using the calculation log 2 ( x median ) , where x is the quantification peak count for a given compound in a given sample (i.e.area under compound peak), and the median ion count across all samples for that compound.Due to the viscous nature of Mabisi, especially the synthetic community samples at later time points, accuracy and reliability of volumes were questionable.For example, pipetting the highly viscous and thick synthetic diversity samples was very difficult, creating pr obable v ariation in sampling volumes.To be conserv ativ e, data was ther efor e further standar dized b y sample using the equation x − μ / σ , where μ is the mean quantification peak count and σ the SD (following normalization by compound) across compounds for a given sample.We realized that this removed all quantitative information about total aroma intensity between samples or div ersity le v els and left information onl y about r elativ e compound peak heights.Ho w e v er, it was the justifiable a ppr oac h considering the likely inaccuracies of sample volumes.

DN A extr action
Ada pted fr om Gr oenenboom et al. 2020 and Sc houstr a et al. 2013(Gr oenenboom et al. 2020 ;Sc houstr a et al. 2013 ).DNA extr action was completed for all eight replicates, of all four communities at transfers t01, t05, and t17.A sample of 1.8 ml of fermented milk was spun down (2 min, 12 000 RPM), then the supernatant and curd r emov ed with a sterile scoopula.Cells wer e r e-suspended in a mix of 64 μL EDTA (0.5 M, pH 8), 160 μL Nucleic Lysis Solution, 5 μL RNAse, 120 μL lysozyme (10 mg/mL), and 40 μL pronase E (10 mg/mL), and incubated for 60 min at 37 • C with agitation of 350 RPM.Cells were dislodged by manually flicking the tube occasionally during incubation.This was to impr ov e mixing with suspension mixture.Bead beading was then performed for 3 min (1 min, 5 min r est, r epeated thr ee times) with sand sized beads, then 400 μL ice-cold ammonium acetate (5 M) added, and the mixtur e immediatel y cooled on ice for 15 min.The mixture was spun down (13 000 r/m x g, 4 min) and 700 μL of supernatant tr ansferr ed to a new 1.5 ml tube.Equal volume of phenol (700 uL) was added, the tube vortexed, and then spun down (6 min, 12 000 RPM, 4 • C).Total volume of 300 μL of supernatant was transferred to a new tube, 300 μL chloroform added, then was vortexed and its content spun down (2 min, 12 000 RPM). Fr om her e, 300 μL of supernatant was tr ansferr ed to another new tube, 400uL of 2isopropanol added and vortexed.This mixture was left at −20 C overnight to precipitate.
Following ov ernight pr ecipitation, the tube w as spun do wn (13 000 r/m, 4 • C, 15 min), then the supernatant car efull y pour ed out so that the DNA pellet remained.A total volume of 1 ml of 70% cold ethanol was added to the tube and then spun down (10 min, 12 000 r/m, 4 • C).The supernatant was car efull y pour ed out and the DN A pellet w ashed again with 1 ml 70% cold ethanol, spun down, and supernatant poured out.The DNA pellet was left to dry at 37 • C for 5 min, then dissolved in 20 μL of TE buffer (pH 8.0).A brief incubation ( < 30 min) at 37 • C impr ov ed dissolving the DNA pellet.Extracted DNA samples were stored at −20 • C. Extracted DN A w as measured on Qubit using dsDNA High Sensitivity kit, and a new diluted sample made with DNA concentration 0.5 ng/ μL.
DNA concentrations in these steps were always measured on Qubit 2.0 fluorometer using dsDNA High Sensitivity Assay kit (Thermo Fisher).

Determination of number of PCR cycles
An a ppr oximate a ppr opriate number of PCR c ycles w as first determined on just three samples (1-F1, 1-L1, and 1-S1), using 13, 16, 18, and 25 cycles .Primers , amounts of r ea gents, and PCR settings were as described below in the section "Tailed PCR reaction".End pr oducts wer e visualized on 1% a gar ose gel and the minimal number of cycles decided based upon the fewest number of cycles where a sufficient band was seen (determined to be 23-30 cycles).The first step in Nanopore sequencing was a PCR reaction using Nanopore specific tailed primers .T he specific number of cycles used for each sample is seen in Table S2 .Two positive contr ols wer e included-the ZymoBIOMICS Micr obial Comm unity DN A Standar d D6305, and a pr e viousl y sequenced Mabisi sample (Groenenboom et al. 2020 ).The tailed primer PCR reaction was as follows:

Tailed PCR reaction
Tailed reaction reagents: 1 uL-DNA [0.5 ng/ μL] 12.5 uL-Phusion High Fidelity PCR 2X master mix (Ther-moFisher) 1.25 uL-forw ar d tailed primer [10uM] 1.25 uL-r e v erse tailed primer [10uM] 9 uL-MilliQ water Tailed cycle conditions: The tailed PCR reaction was performed another two times, resulting in thr ee separ ate tailed primer PCR products per sample.Placement of PCR tubes in machine was adjusted for each PCR to avoid edge effects.Each amplified DNA sample was all visualized on 1% a gar ose gel to confirm successful amplification, then 8 μL of each PCR reaction were combined.A total volume of 24 μL of amplified DNA per samples was used for the clean-up.

PCR clean-up
For each sample, the 24 μL of amplified DN A w as cleaned with 24 μL of homemade SPRI beads (i.e.1:1 ratio) (1 ml Ser a-Ma g SpeedBeads (Cytiv a, Marlbor ough, MA, USA) cleaned and dissolved in 50 ml end volume containing 2.5 M NaCL, 20 mM PEG, 10 mM Tris-HCL and 1 mM EDTA) and eluted into 20 μL of MilliQ w ater.DN A concentration of cleaned amplicons was measured using Qubit 2.0 Fluorometer.A new dilution of 15 μL of the cleaned, amplified 16S rRNA gene PCR product was made into a new diluted sample with DNA concentration 0.5 nM.
In 47 μL of milliQ water, 1 μg of the barcoded, pooled, cleaned libr ary was pr epar ed.Fr om her e, the libr ary was r epair ed, endpr epped, and ada ptor ligated according to the Oxford Nanopore Technologies PCR barcoding ( 96 The pr epar ed libr ary w as loaded on a SpotON Flo w Cell (FLO-MIN106D) on a MinION MK111775 sequencing device (Oxford Nanopor e Tec hnologies).Basecalling was performed shortl y after using Guppy software version 6.2.4 + a11ce76.

Bioinformatics
To assign 16S rRNA gene sequences to taxonomic identity, we first downloaded the SILVA r efer ence database (Quast et al. 2013 ).As w e w er e inter ested in v ariation abov e the species le v el, a custom database was produced using vsearch to cluster 16S rRNA gene sequences at a 95% similarity threshold (Rognes et al. 2016 ).Reads were aligned to this database using minimap2 (Li 2018 ), and the best matching cluster was assigned as the taxonomic identity for eac h r ead.Details of sequence compositions for main clusters ar e found in Table S3 Excel file .Clusters to which less than five reads matc hed wer e assumed to be tr ace contamination and discarded.
T he o v er all fr equency of clusters with < 1% abundance is similar across diversity levels, as evidenced by the amount of "missing reads" to complete 1.00 abundance on bar plot (i.e.blank space at top of bar plot).Since full communities and synthetic communities have similar frequency of the rarest clusters, they are likely present due to general contamination and not true diversity being r emov ed thr ough our anal yses.Unexpected "external" r eads wer e identified (purple and green colouring) but are at reasonable levels as expected from low level laboratory contamination amplified by repeated PCRs during Nanopore sequencing.
We used the SILVA database r efer ence sequences to identify types of bacteria.Our bioinformatic analyses were able to identify Lactobacillus and Limosilactobacillus types to the species le v el, with all clusters containing nearly 100% of sequences with same species identification ( Table S3 Excel file ).This was not true for Acetobacter types, where clusters included sequences matching to se v er al Acetobacter species and these species matched to multiple clusters .For example , sequences labelled as Acetobacter lovaniensis are included in Acetobacter A-C clusters .Hence , for the purpose of our study we only classified into general genus le v el gr oups for Lactobacillus , Limosilactobacillus , and Acetobacters .Our bioinformatic analysis could classify Lactobacillus and Limosilactobacillus clusters as homo or heter ofermentativ e, r espectiv el y.
The very rare Limosilactobacillus type detected in only some full and medium communities, is the sole identified heterofermentative type-Limosilactobacillus fermentum .The yeast present in the Mabisi communities used in this study was isolated and identified as Geotrichum candidum using amplification and sequencing of the ITS (internal transcribed spacer) region (top match GenBank: MK967716.1).Yeast was visibly seen growing at the air interface only in the full diversity and medium diversity communities.

Statistical analyses
All analyses were performed in R version 4.3.1 (R Core Team 2023 ).Significance was defined for all analyses as P -value < 0.05.Results outputs and details of statistical analysis are found in Supplementary Material .

Measures of bacterial di v ersity
The "v egan" R pac ka ge functions "specn umber" and "di versity" were used to calculate bacterial species richness and Shannon div ersity, r espectiv el y.A tw o-w ay ANOVA of comm unity * tr ansfer was performed separ atel y for ric hness and Shannon div ersity v alues using the "lm" function fr om "lme4" R pac ka ge, type 3 partial sum of squares (Bates et al. 2015 ).Post-hoc Tuk e y's pairwise comparisons with adjusted P -values were then performed using the "emmeans" function fr om "v egan" R pac ka ge at 95% confidence intervals (Oksanen et al. 2022 ).

Acidity
A tw o-w ay ANOVA of comm unity * tr ansfer was performed on pH values using the "lm" function from "lme4" R package was used, type 3 partial sum of squares (Bates et al. 2015 ).Post-hoc Tuk e y's pairwise comparisons with adjusted P -v alues wer e then performed using the "emmeans" function from "vegan" R package at 95% confidence intervals (Oksanen et al. 2022 ).

Community compositions
A PERMANOVA was performed on the full data set of community compositions to e v aluate the effect of community (i.e.div ersity le v el), time point, and their inter action on comm unity composition.The "adonis" function from "vegan" R package was used, with "bray" method (Oksanen et al. 2022 ).Post-hoc comparisons between communities were made separately for each transfer subset.The "pairwise.perm.manova"function with "Pillai test" was used from the "RVAideMemoire" R package (Herve 2023 ), follo w ed b y Bonferr oni corr ection of P -v alues using "p.adjust" function.
Differences in dispersion of communities was tested to confirm that the significant effects from PERMANOVA were because group centr oids differ ed and not because gr oup v ariances differ ed.This was performed using the "betadisper" function of "v egan" R pac kage (Oksanen et al. 2022 ).

Bacterial species sorting trajectories are similar following dilution of di v ersity
Figure 2 shows relative abundance of bacterial species in each of the communities at transfer 1, 5, and 17.Species profiles of the bacterial communities at transfer 1 show a clear signature of the pr ogr essiv e serial dilution of diversity into full diversity, medium di versity, low di versity, and synthetic communities: Lactobacillus A type (light turquoise) has a higher r elativ e abundance in the Lactobacillus clusters are shown in blue colour shades, Acetobacter in yellow-orange, Limosilactobacillus in dark navy blue, and low abundance contaminant types in green or purple.Clusters with < 1% abundance were not identified nor plotted.Some replicate populations are missing due to insufficient DNA concentr ations, r esulting either fr om poor DNA extr action or libr ary pr epar ation.Medium comm unity, r eplicate eight was lost early in pr opa gation and r emov ed fr om all anal yses.Abbr e viations: med = medium, synt = synthetic; t01 = tr ansfer 1, t05 = tr ansfer 5, and t17 = tr ansfer 17. medium div ersity comm unities (mean = 0.21) compared to the full div ersity comm unities (mean = 0.08), and is the dominant type in low diversity communities (mean = 0.52).The compositions of bacterial types at transfer 1 significantly differ ( Table S13 : P < 0.05 for all pairwise comparisons of community compositions with Bonferr oni corr ection), yet the bacterial species richness is similar across full, medium, and low diversity communities ( Fig. S5 : mean bacterial species richness transfer 1, full = 4.5, medium = 4.43, low = 4, synthetic = 3).T hus , while our dilution a ppr oac h significantl y c hanges initial bacterial community composition, it is not evident with the standard diversity measure of species richness.Ho w ever, our dilution approach did eliminate metabolic guilds (Table 1 , based from Fig. 2 ), and communities did differ in their initial metabolic guild diversity due to serial dilution of the full community.Most notably is the loss of yeast in the low and synthetic communities.
The increase in relative abundance seen at transfer 1 of Lactobacillus A type diminishes o ver time , with Lactobacillus B type (medium blue) dominating in all replicates and all initial dilution treatments by transfer 17 (mean relative abundance Lactobacillus B : full = 0.87, medium = 0.86, low = 0.79, and synthetic = 0.70).Ho w e v er, the abundance of Lactobacillus A cluster remains compar ativ el y slightl y higher in low div ersity comm unities at the end of pr opa gation (tr ansfer 17 mean r elativ e abundance Lactobacillus A : full = 0.04, medium = 0.05, low = 0.11, and synthetic = 0.00).This Lactobacillus A cluster was not included in the synthetic comm unity but sur prisingl y a ppears in one synthetic comm unity replicate at transfer 17, likely due to cross contamination during the serial pr opa gation.
Bacterial community composition significantly differed per div ersity le v el and time point ( Table S12 : PERMANOVA on community compositions , Bra y method used: effect of initial diversity (p = 0.001, F = 70.9,DF = 3), transfer (p = 0.001, F = 55.1,DF = 2), diversity x transfer interaction (p = 0.001, F = 27.8,DF = 6)).We further conclude that significant results are due to true differences in group centroids and not dispersion within groups ( Table S11 : P > 0.05 for ANOVA community effect of distances to group centroid at all time points).Deriv ed fr om the same data presented in Fig. 4 , NMDS plots of bacterial community compositions for transfer 1, 5, and 17 per initial div ersity le v el (Fig. 3 ) demonstrate convergence of bacterial community composition across full, medium, and low div ersity comm unities ov er time.Furthermor e, final r elativ e r atios of Acetobacter A, Acetobacter B, Lactobacillus A, and Lactobacillus B ar e similar acr oss all comm unities at tr ansfer 17, a part fr om synthetic diversity ones ( Fig. S3 , Table S13 : P < 0.05 for pairwise comparisons at transfer 17 of community compositions for synthetic versus full, medium, and low).
Acetobacters are at highest relative abundance in the low diversity comm unities earl y in pr opa gation (mean r elativ e abundance Acetobacters transfer 1, full = 0.12 (SD = 0.07), medium = 0.08 (SD = 0.05), low = 0.23 (SD = 0.08), synthetic = 0.11 (SD = 0.08)), but synthetic comm unities demonstr ate the highest pr oportions of Aceto- Although r emov al of metabolic guilds acr oss our starting communities affected initial and final metabolic function (see the "Results section" below), it inter estingl y did not hav e a lar ge effect on the outcome of species sorting trajectories of the remaining bacterial comm unity ov er 16 cycles of pr opa gation.This is especially true when focusing on the relative ratios of the major bacterial metabolic guilds-acetic acid versus homofermentative lactic acid bacteria, which were surprisingly uniform across full diversity , medium diversity , and low diversity communities after 16 cycles of pr opa gation.If r eplicate comm unities of higher div ersity div er ged mor e so in their compositions o vertime , then our results would have aligned with the theory that "diversity begets diversity" via gr eater nic he construction (Madi et al. 2020 ).Re v ersel y, our results could have exhibited replicate communities of lo w er di versity di verging from one another, which would support the concept that greater div ersity r estricts possible tr ajectories due to less available niche space (van Moorsel et al. 2021 ).We observe neither outcome .T he minimal effect of div ersity on sorting tr ajectories in our experiment could be explained by high functional redundancies (Allison and Martiny 2008 ) in the major bacterial types (i.e.acetic acid and lactic acid bacteria).

Yeast dri v es di vision of function shown in metabolic profiles
Metabolic profiles of replicate communities at four levels of diversity in metabolic guilds show that two groupings persist across the three time points analysed, with metabolic profiles of full and medium div ersity comm unities clustering together, v ersus low and synthetic (Fig. 4 ).The divide exists along principal component 1 (PC1) but not the second component (PC2).PC1 explains between 53.1% (transfer 5) to 65.5% (transfer 1) of the variation, while PC2 explains between 11.9% (transfer 1) and 19.6% (transfer 17).A division in metabolic profiles between full and medium communities is observed along PC2 at transfer 1 (Fig. 4 A), disappearing at transfer 5 (Fig. 4 B).The metabolic pr ofiles of comm unities at the four div ersity le v els inconsistentl y exhibit r elativ el y lar ger or smaller spread in data points across the time points but replicates of the synthetic comm unity (pur ple ellipses) ar guabl y show closer clustering in metabolic profiles throughout the experiment.
Although at minority abundances in typical Mabisi communities (Sc houstr a et al. 2013 ), yeasts contribute a unique metabolism since they are the sole alcohol fermentative type and are known to produce other distinct compounds such as methyl-esters.We thus pr edicted that pr esence of yeast in the full and medium diversity comm unities or, conv ersel y, its absence in the low diversity and synthetic communities, would significantly influence community function.Initial metabolic pr ofiles acr oss all time points provide support for this expectation, where alcohol and ester compounds typically linked to yeast metabolism are found.The presence of yeast in only the full and medium diversity communities is reflected by presence of ethanol, acetaldehyde, and hexanoic acidethyl ester among the main metabolites detected.Acetaldehydes and esters are a known products of yeast metabolism (Liu andPilone 2000 , Dzialo et al. 2017 ).An additional compound with vector contributions to w ar ds the full and medium community gr ouping is 1-butanol-3-methyl, whic h is a br eakdown pr oduct of leucine via the Ehrlich pathway found in various yeast species (Szudera-Ko ńczal et al. 2020 ).The heter ofermentativ e type Limosilactobacillus is also absent in low diversity and synthetic communities, whose presence is expected to alter metabolite profile; howe v er, since this ov erla ps with the presence/absence of yeast in the communities, we are unable to disentangle which metabolite c hanges ar e specific to the pr esence or absence of heter ofermentative types.
All div ersity le v els contained the two dominant fermentative types-homofermentative lactic acid bacteria ( Lactobacillus A and B ) and acetic acid bacteria (Acetobacter A and B )-which our data suggest are the core bacterial members for community function at the le v el of metabolite production, in addition to yeast.While yeast appear as a core microbial community member in our study's Mabisi sample, inter estingl y, not all natur al Mabisi pr oducts contain yeast (Moonga et al. 2020 ;Sc houstr a et al. 2013 ).Mi-cr obial comm unity pr ofiles of Mabisi samples v ary acr oss r egions and processors in Zambia.There is past and ongoing r esearc h to link micr obial comm unity compositions to pr ocessing methods (Moonga et al. 2020 ) and consumer pr efer ences.
Furthermor e, metabolic pr ofiles ov erla p betw een lo w diversity and synthetic communities (Fig. 4 ), yet their bacterial compositions differ with loss of intra-species genotype diversity and Lactobacillus A in synthetic communities (Fig. 2 , Table 1 ).Hence, in Mabisi comm unities, intr aspecies bacterial genotype compositions do not appear to strongly influence metabolic profiles, suggesting functional redundancies in bacterial metabolic capacities.Others have found low functional redundancy in microbial communities for particular functions, for example methane production (Sierocinski et al. 2018 ), where there is considerable impact by loss of any species.Our results align regarding yeast but not for species or genotypes of the abundant lactic acid and acetic acid bacteria.Ho w e v er, we measur ed br oader, mor e gener alist functions of ov er all ar oma pr ofiles and acidity.High functional r edundancy of Mabisi micr obial comm unities compar ed to envir onmental systems such as feedstock fermenters (Sierocinski et al. 2018, Wagg et al. 2019 ) or soil (Wagg et al. 2014(Wagg et al. , 2019 ) ) is likely due to immense microbial diversity in these other systems and their specific measures of function.

Functional redundancy exhibited for acidity
Figure 5 shows pH of cultures after a full cycle of gr owth ov er the course of the experiment.The ov er all pH differ ences acr oss div ersity tr eatments r anged onl y betw een ∼pH 3.4-3.8,y et with a significant effect of comm unity, tr ansfer, and their inter action on acidity ( Table S8 : P < 2.2e-16 for all effects).Full diversity communities and medium div ersity comm unities maintained similar pH thr oughout pr opa gation ( Table S9 : Tuk e y's pairwise comparison full vs. medium P > 0.05 for all time points).Low diversity communities maintained a higher pH until a drop between transfer 15 and 17 where they converge to a comparable value as the full and medium communities ( Table S9 : Tuk e y's pairwise comparison low vs. full and low vs. medium both P < 0.05 for time points 1 through 15, then both P > 0.05 for transfer 17).Whereas the pH in synthetic communities was significantly higher than all others by transfer 17 ( Table S9 : P < 0.0001 for Tuk e y's pairwise comparisons of synthetic vs. full, medium, and low at tr ansfer 17).The pH dr opped significantl y between tr ansfer 1 and tr ansfer 17 for full, medium, and low div ersity tr eatments ( Table S10 : Tukey's pairwise comparison transfer 1 vs. transfer t17 P < 0.05 for full, medium, and low communities).
High-functional redundancy is observed in the bacterial part of the microbial community (Wagg et al. 2019, Gralka et al. 2020 ), thus we predicted that general acidification properties of lactic acid and acetic acid bacteria would be maintained in diluted comm unities, r egardless of the loss of yeast or r ar e genotypes (Wagg et al. 2019, White et al. 2020 ).We found support for this prediction; ho w e v er, a v ery slightl y higher pH was observ ed for the synthetic, hence lo w est div ersity, comm unity.This modest increase suggests that the presence of yeast and r ar e bacterial genotypes does not impact community acidification.A dynamic jump in pH in synthetic diversity communities between transfers 1 and 3 is not surprising since this community was more naïve in its member structure and substrate (these isolates had been pre-cultured in MRS broth).The range in pH between diversity levels is minimal ( ∼pH 3.4-3.85)but may still r epr esent a selective force on the microbial communities.Mabisi acidity increases quickly during fermentation, with the most notable pH drop between 24 and Figure 5. Loss of diversity initially increases pH, which decreases over continual propagation.Measures of pH of full diversity, medium diversity, low div ersity, and synthetic comm unities ov er 16 r ounds of serial pr opa gation.Mean pH acr oss r eplicates of the same comm unity div ersity le v el measur ed at e v ery second tr ansfer.Points show mean v alue for the eight r eplicate comm unities (exce pt medium di v ersity tr eatment with se v en r eplicates), vertical bars show SEM.Abbreviations: med = medium and synt = synthetic.48 h (Groenenboom et al. 2020 ).Our observed trend of continually decreasing pH may be explained by the preferential propagation of, and thus selection for, persisting abundant types in the acidic conditions created after 72 h of fermentation.Increased acidity is a common outcome of bacterial community domestication and r epeated bac k-slop pr opa gation of fermentation starter cultur es (Bachmann et al. 2011, Spu ś 2016, van Kerrebroeck et al. 2016 ).It would be interesting to test how m uc h lower the pH would reduce with further pr opa gation and when or if pH stabilisation would be r eac hed.Evolving comm unities could be impr oving r esource use, consequentl y excr eting mor e metabolites due to larger supported population sizes, further lo w ering the pH over ad ditional re peated transfers.

Conclusions
In this study, we assessed how pr ogr essiv el y diluting a natur al micr obial comm unity alter ed comm unity function and how this function, as well as bacterial community composition, would subsequentl y c hange on an ecological time scale due to selection upon repeated cycles of pr opa gation.By exploring the relationship between community function and metabolic guild diversity over time, we observed repeatable changes of replicate lineages for metabolic profiles and acidity related to starting diversity le v els of metabolic guilds .T hese results show a clear division in metabolic profiles with full diversity and medium diversity communities on one hand, and low diversity and synthetic communities together on the other hand.This division was sustained throughout repeated rounds of propagation.Further, we found that changes in bacterial community composition (i.e.bacterial species sorting tr ajectories) ov er r epeated cycles of pr opa ga-tion gener all y r esulted in a conv er gence of bacterial communities to the same composition, irr espectiv e of initial metabolic guild diversity.
Most surprising was the seeming lack of influence of yeast presence or absence on bacterial compositions after 16 cycles of propagation.T he con vergence of bacterial community compositions, regardless of the presence of yeast communities, suggests yeast, and bacteria exist in sufficiently unique niches of resource use.Evidence for strong division in resources and hence lack of influence of yeast on bacterial community composition in a fermented food was surprising and contr asts a gainst pr e vious findings showing metabolic associations between the two (Mendes et al. 2013, Suharja et al. 2014, P onomaro va et al. 2017, Blasche et al. 2021, Xu et al. 2021 ).Ho w e v er, these inv estigations mostl y focus on Sacc harom yces cerevisiae and specific Lactobacilli species or strains; the yeast in our system is identified as G. candidum .If Lactobacilli versus Acetobacters had differing metabolic associations with y east, w e w ould expect shifts in their r elativ e abundances following the r emov al of G. candidum, but this was not observed.We can thus hypothesize that in the Mabisi system, yeast do not affect bacterial growth by consuming end products of bacterial metabolism such as lactate, thus are simply secondary in the metabolic route.In addition to more cross-feeding like interactions, yeast and bacteria could be in resource competition, in which case the removal of y east w ould open niche space for certain metabolic guilds or species .T hose whose niche space previousl y ov erla pped mor e with y east, or w ere poor competitors, would expectedly establish at higher abundances in the absence of yeast; this is not what we observ e. Ov er all, our r esults support that amongst bacterial metabolic types in our community, they are all in similar, at most very w eak, resour ce competition with yeast since bacterial community compositions converge across all comm unity div ersities r egar dless of y east's presence.A next step to further elucidate the influence of yeast on community function in Mabisi is to add isolated G. candidum to the low and synthetic comm unities, then compar e metabolic pr ofiles and acidity.The r e v erse dir ection could also be taken, by eliminating yeasts from the full and medium communities with fungicide, which would maintain r ar e bacterial types.
We inter pr et the observ ed r epeatability during the pr opagation c ycles betw een r eplicate linea ges and conv er gence in bacterial communities to be evidence of ecological selection (Vellend 2010(Vellend , 2016 ) in our study system.In evolutionary biology r esearc h, r epeated or conv er gent r atios of genotypes under a common environment is interpreted as evidence for adaptive selection (Hughes 1999 ); we take the same inter pr etation for our results but at the ecological le v el of r atios of bacterial metabolic guilds.Curr entl y, we can onl y hypothesize about the sources of ecological selection in our study, with possibilities including our chosen laboratory conditions (temperature, fermentation time and associated acidity, supermarket milk content) and unspecified biotic species interactions.While prokaryotes and bacteriophage are known members of the microbial community of our natural experimental system (Schoustra et al. 2013 ), we have focussed our analyses on species sorting in the bacterial part of the microbial community since bacteria are present in all communities regardless of starting diversity (i.e.dilution) tr eatment.Ther e is gr owing e vidence of highl y r epeatable species sorting trajectories in microbial community assembly, both in varied environmental samples (Goldford et al. 2018, Diaz-Colunga et al. 2022 ) and communities of isolated strains (Cairns et al. 2020 ).Our r esearc h similarl y finds conv er gence ov ertime in community compositions but instead after disturbance of initial diversity.
Studying functionality responses on longer time scales after an erosion of metabolic guilds , species , and genotypic diversity r emains an under explor ed topic that we made initial steps here to explore .T he ecological processes of species sorting in a community can be paralleled to evolutionary dynamics of genotypic changes within a species (Vellend 2010(Vellend , 2016 ) ). Understanding how c hanges acr oss hier arc hical le v els-fr om genes to comm unities, to ecosystems-influence one another is undoubtedly complicated and it merits further steps to elucidate.In this study, we focused on ecological processes without investigating genotypic changes.Ho w ever, w e foresee exciting future avenues of research to unr av el the r ole of ecology v ersus e volution in altering communities and their functionality.Creating a synthetic community comprised of four fermentative types (i.e.metabolic guilds) that grow in a defined media, combined with whole genome sequencing and/or metagenomics could allow identification of novel genotypes, and is an exciting avenue to better explore ecologicalevolutionary dynamics.

Figure
FigureStarting diversity transiently alters bacterial community composition .Bacterial community compositions at transfers 1, 5, and 17.Lactobacillus clusters are shown in blue colour shades, Acetobacter in yellow-orange, Limosilactobacillus in dark navy blue, and low abundance contaminant types in green or purple.Clusters with < 1% abundance were not identified nor plotted.Some replicate populations are missing due to insufficient DNA concentr ations, r esulting either fr om poor DNA extr action or libr ary pr epar ation.Medium comm unity, r eplicate eight was lost early in pr opa gation and r emov ed fr om all anal yses.Abbr e viations: med = medium, synt = synthetic; t01 = tr ansfer 1, t05 = tr ansfer 5, and t17 = tr ansfer 17.

Figure 4 .
Figure 4. Diversity of metabolic profiles persist through propagation .Compounds associated with yeast metabolism primarily impact divisions in metabolic profile .PCA analysis of metabolic profiles compiled from 19 compounds using GC-MS analyses at transfer 1 (A), 5 (B), and 17 (C).Dots r epr esent individual replicates.Ellipses are 95% confidence interv als.Dir ection of v ectors indicate the compound's contribution to either principal component, whereas the vector length indicates the amount of variation explained by the two plotted principal components.Yeast associated compounds: ethanol, acetaldehyde, hexanoic acid-ethyl ester, 1-butanol 3-methyl.Abbr e viations: med = medium, synt = synthetic.

Table 1 .
Presence/absence of functional types (metabolic guilds) between diversity level treatments.