Serial fermentation in milk generates functionally diverse community lineages with different degrees of structure stabilization

ABSTRACT Microbial communities offer considerable potential for tackling environmental challenges by improving the functioning of ecosystems. Top-down community engineering is a promising strategy that could be used to obtain communities of desired function. However, the ecological factors that control the balance between community shaping and propagation are not well understood. Dairy backslopping, which consists of using part of the previous production to inoculate a new one, can be used as a model engineering approach to investigate community dynamics during serial propagations. In this study, 26 raw milk samples were serially propagated 6 times each, giving rise to 26 community lineages. Bacterial community structures were analyzed by metabarcoding, and acidification was recorded by pH monitoring. The results revealed that different types of community lineages could be obtained in terms of taxonomic composition and dynamics. Five lineages reached a repeatable community structure in a few propagation steps, with little variation between the final generations, giving rise to stable acidification kinetics. Moreover, these stabilized communities presented a high variability of structure and diverse acidification properties between community lineages. Besides, the other lineages were characterized by different levels of dynamics leading to parallel or divergent trajectories. The functional properties and dynamics of the communities were mainly related to the relative abundance and the taxonomic composition of lactic acid bacteria within the communities. These findings highlight that short-term schemes of serial fermentation can produce communities with a wide range of dynamics and that the balance between community shaping and propagation is intimately linked to community structure. IMPORTANCE Microbiome applications require approaches for shaping and propagating microbial communities. Shaping allows the selection of communities with desired taxonomic and functional properties, while propagation allows the production of the biomass required to inoculate the engineered communities in the target ecosystem. In top-down community engineering, where communities are obtained from a pool of mixed microorganisms by acting on environmental variables, a major challenge is to master the balance between shaping and propagation. However, the ecological factors that favor high dynamics of community structure and, conversely, those that favor stability during propagation are not well understood. In this work, short-term dairy backslopping was used to investigate the key role of the taxonomic composition and structure of bacterial communities on their dynamics. The results obtained open up interesting prospects for the biotechnological use of microbiomes, particularly in the field of dairy fermentation, to diversify approaches for injecting microbial biodiversity into cheesemaking processes.

IMPORTANCE Microbiome applications require approaches for shaping and propagat ing microbial communities.Shaping allows the selection of communities with desired taxonomic and functional properties, while propagation allows the production of the biomass required to inoculate the engineered communities in the target ecosystem.In top-down community engineering, where communities are obtained from a pool of mixed microorganisms by acting on environmental variables, a major challenge is to master the balance between shaping and propagation.However, the ecological factors that favor high dynamics of community structure and, conversely, those that favor stability during propagation are not well understood.In this work, short-term dairy backslopping was used to investigate the key role of the taxonomic composition and structure of bacterial communities on their dynamics.The results obtained open up interesting prospects for the biotechnological use of microbiomes, particularly in the field of dairy fermentation, to diversify approaches for injecting microbial biodiversity into cheesemaking processes.abundances of species during the propagations (23).However, this study used three samples of raw milk, and the propagations were carried out in culture medium, which does not allow us to take into account the microbial diversity of raw milk and the behavior of the communities when the propagations are carried out in milk.
The aim of this study was to investigate the dynamics of bacterial communities during the early stages of a backslopping process and to analyze the relationship between the dynamics of community structures and their acidification properties.We set up an experiment based on serial fermentation starting from raw milk samples collected from 26 different farms of the same protected designation of origin (PDO) geographical area in the north east of France.Each milk sample was serially fermented six times in mesophilic conditions.Acidification kinetics were recorded by pH monitoring during each fermentation, and bacterial communities were analyzed by metabarcoding.

Sample collection
Raw milk samples, called M1-M26, were collected from the tanks of 26 farms in the north east of France (Grand Est region).Samples M1-M16 were collected on 2 September 2021, and samples M17-M26 were collected on 29 September 2021.Approximately 500 mL was taken directly from the tanks for each sample, put in a sterile container, and kept at a low temperature until arrival at the laboratory.Samples were stored at 4°C overnight.Then, 300 mL was used for centrifugal concentration for DNA extractions, and 25 mL was stored at −80°C until serial fermentation.

Raw milk concentration
To concentrate raw milk for the metabarcoding analysis, 300 mL of each milk sample was centrifuged for 30 min at 5,300 × g at 4°C.Liquid and cream were discarded, and the pellets were resuspended in 1.5 mL of phosphatebuffered saline.Then, the samples were put in 2 mL tubes and centrifuged for 5 min at 13,000 × g at 4°C.Supernatants were discarded, and the pellets were stored at −80°C until DNA extraction.

Serial fermentation
Tubes containing 25 mL of frozen milk were thawed and placed in a water bath at 34°C for 24 h resulting in a first fermentation step.Then, 250 µL of each fermented sample was inoculated in 25 mL of fresh commercial ultra-high-temperature (UHT) processed milk (Marque Repère, E. Leclerc, Ivry-sur-Seine, France) and incubated at 34°C.This step was carried out five times to get six fermentations in total for each raw milk sample (Fig. 1).As acidification for steps 2-6 was faster than for step 1, the incubation duration was reduced to 18 h.During each fermentation, the pH was monitored, and at the end of each fermentation, samples were collected, placed in 2 mL tubes, and stored at −80°C until DNA extraction.Three campaigns were carried out: the serial fermentation of the raw milk samples M1-M10, M11-M18, and M19-M26 were started on 13 September 2021, 4 October 2021, and 11 October 2021, respectively.

Acidification kinetics
Acidification was recorded by pH monitoring during each fermentation.Before each pH time course, pH sensors (VWR pHenomenal 111, VWR International, Radnor, PA, USA) were cleaned with deionized water, sterilized by immersion in 70% ethanol, and then rinsed by immersion in sterile deionized water.The pH was measured every 5 min using a multi-channel pH meter (Consort datalogger D291, Consort bvba, Turnhout, Belgium).

Acidification parameter analysis
Acidification parameters (lag time, maximum acidification rate, and final pH) were extracted from each curve using a custom R script (version 4.3.1)(24).
Maximum acidification rate (MAR) was the lowest slope coefficient among the linear regressions obtained on a 20-point sliding window.The resulting linear model was used to calculate the intersection point with a line of null slope and intersecting y axis at the first pH value of each time point.The time point corresponding to the intersection between this line and the linear model for the MAR was used as an estimation of the lag time (corresponding to the time before the pH drop).For the final pH, a linear model was built with the last 20 points of the curves and used to predict the final pH at 18 h (for F2-F6) or 24 h (for F1).

DNA extraction
DNA was extracted from the raw milk pellets and the fermented samples using a FastDNA Spin Kit for Soil (MP Biomedicals, IllkirchGraffenstaden, France) with bead beating, following the manufacturer's recommendations.DNA was stored at −20°C.Contamination controls were performed following the same protocol without DNA, every time extractions were performed.

PCR amplification
A nested PCR was performed for all raw milk samples as bacterial concentrations were low.First, the 16S rRNA gene (~1,450 bp) was amplified using Taq polymerase (Appli gene, IllkirchGraffenstaden, France), and the universal bacterial primers W02 (5′-GNTAC CTTGTTACGACTT-3′) and W18 (5'-GAGTTTGATCMTGGCTCAG-3′) as described previously (25,26).The PCR program was initiated with incubation at 94°C for 3 min, followed by 17 cycles of 94°C for 30 s, 50°C for 30 s, and 72°C for 1 min 30 s, followed by a final extension at 72°C for 10 min.Then, the product obtained with this first PCR was used as a matrix for the next PCR, targeting the variable region V3-V4 of the 16S rRNA gene (~510 bp) using MTP taq polymerase (Sigma, Saint-Louis, MO, USA) and the primers MSQ-16SV3F (5′-TACGGRAGGCWGCAG-3′) and PCR1R_460 (5′-TTACCAGGGTATCTAATCCT -3′), as described previously (27).The PCR program was initiated with incubation at 94°C for 1 min, followed by 30 cycles of 94°C for 1 min, 65°C for 1 min, and 72°C for 1 min, followed by a final extension at 72°C for 10 min.For all fermented samples, only the second PCR for V3-V4 amplification was performed using the extracted DNA as matrix.

Sequencing and data analysis
The 16S amplicon sequencing was performed by the GeT platform (GenoToul, Toulouse, France) using the Illumina Miseq technology, for the 26 raw milk samples and every fermented sample (182 total), as well as the contamination controls for the DNA extraction.
Further analysis was performed using custom R scripts (R version 4.3.1)(24).Alpha and beta diversity analyses were conducted using tools from the "vegan" package (33).
Results were expressed as means ± standard error of the mean (SEM).Pairwise comparisons were performed with the Tukey honestly significant difference (HSD) test following analysis of variance (ANOVA) or the Wilcoxon rank sum exact test with the Bonferroni P value adjustment method following a Kruskal-Wallis rank sum test when the data did not meet the application conditions.Variance comparisons were carried out using an F test.Correlation tests were carried out using Pearson's product-moment correlation test.Simple linear regressions were performed to further investigate the relationship between quantitative variables.In the figures, different superscript letters indicate a significant difference (P < 0.05).

Acidification dynamics during serial fermentation
Twenty-six milk samples were serially fermented 6 times each, giving rise to 26 commun ity lineages.A lineage corresponds to one single raw milk sample and the resulting six communities obtained by serial fermentation.During each fermentation step, the pH was monitored over time (Fig. S1).Three parameters of interest were extracted from each pH time course: the time at which pH starts to drop (lag time), the maximum acidification rate (MAR), and the pH at the end of the fermentation (final pH).
The analysis of the parameters showed that the lag time was always longer for the first fermentation compared to the following steps (Fig. S2A, P = 6.0 × 10 −14 , Wilcoxon rank sum exact test [WRSET]).It took 14.3 h on average before the pH drop occurred in the F1 and from 2.1 to 2.7 h for the following steps.The speed of pH drop was lower for fermentation F1 (−0.27 pH units per hour) compared to the other fermentations (from −0.36 to −0.39, Fig. S2B, P < 9.0 × 10 −4 , WRSET).These results show that the first fermentation was markedly different from the others.
In GP1, the lag time, the MAR, and the final pH did not significantly change over the course of fermentation (Fig. 2B through D, P > 0.05 Pearson's product-moment correla tion [PPMC]).Interestingly, even if the final pH in GP1 was stable from F2 to F6, it was significantly different among the lineages: the final pH was 4.13 for M15 on average, 4.33 for M3 and M26, 4.40 for M18, and 4.86 for M19 (all lineages were significantly different-P < 0.022-except lineages M3 and M26 according to the Tukey HSD test following ANOVA).
Overall, these results show that three groups can be delineated according to their acidification patterns from the second fermentation step onwards.The group GP1 was characterized by stable acidification parameters but contained community lineages characterized by different final pH values.In contrast, the other two groups GP2 and GP3 were characterized by different levels of kinetics variability during the serial propaga tions.

Bacterial community dynamics during serial fermentation
Overall, bacterial community structures of the raw milk samples and of serial fermenta tion revealed a high variability in terms of both taxonomic composition and dynamics (Fig. 3).The majority of initial raw milk samples (17/26) contained more than 50% relative abundance of the Pseudomonas genus (with species P. yamanorum for 16 samples and P. gessardii for 1 sample, T0 Fig. 3).One sample (M6) was dominated by Streptococcus porcinus, whereas other samples were more diversified (M18, M21, M20, M22, M4, M17, M19, M16, with a richness index of more than 75 and a Shannon index above 2.5, Fig. S3A  and C).No significant differences in richness and Shannon index were observed between the starting raw milk samples of the three groups (Fig. S3B and D).
After one fermentation (F1), the richness markedly decreased for all lineages from 60.9 ± 3.9 in T0 to less than 23.0 (Fig. S4A, P < 1.6 × 10 −8 , WRSET).The richness in FIG 3 Bacterial composition of the community lineages at the genus or species level.T0: raw milk samples, F1-F6: fermented milk samples.Genera representing less than 5% of the total relative abundance for each fermentation series have been aggregated into the class "Other" (gray).The barplot series are presented according to groups GP1, GP2, and GP3. the last fermentation step was 13.2 ± 0.7 for all 26 communities.The first fermentation also exhibited a decreased Shannon index (from 1.9 ± 0.2 to 1.2 ± 0.1), although less noticeably (Fig. S4B, only the comparison of T0 with F2 was significant with P = 1.3 × 10 −2 , WRSET, otherwise P > 0.19).The poor significance is likely due to a higher variability of the Shannon index between raw milk samples compared to propagated communities (Fig. S4B, P < 2.2 × 10 −16 , F test to compare the variances between T0 and fermented samples).These results showed that serial fermentation induces a decrease in bacterial alpha diversity.
Among the F1 communities, 17/26 were dominated by ASV of the phylum Bacillota (Fig. S5), and more specifically by LAB ASV from the genera Lactococcus or Streptococcus (with species L. lactis and S. thermophilus, respectively, Fig. 3).The other remaining F1 communities were dominated by Gram-negative bacteria of the phylum Pseudomona dota (Escherichia genus in 7 F1 communities, Klebsiella genus for one and Enterobacter genus for another one).In all cases, the communities ended up being dominated by Bacillota from step F3 (Fig. S5).However, the Pseudomonadota remained present in a large proportion until fermentation F6 in several fermented milk samples (e.g., M6, M7, M9).Interestingly, the abundance of Lactobacillus genus sensu lato was very low in F1 (<1%) and started to increase to at least 25% (up to 63%) in 9 F3 communities.They ended up being dominant (>50%) at F6 for 11 lineages where the main species was L. delbrueckii.
The bacterial communities from GP1 were characterized by minor changes in community structure during serial fermentation (Fig. 3).Lactococcus (M3, M26) or Streptococcus ASV (M15, M18) were dominant in GP1 (>75% relative abundance), except for M19 in which proportions of Lactococcus and Klebsiella were similar during the fermentation steps, with Klebsiella even becoming dominant in the last steps (Fig. 3).The GP2 communities showed more changes than GP1 communities and were dominated by multiple LAB (Lactococcus, Streptococcus, Lactobacillus, and Enterococcus, Fig. 3).GP3 showed broader changes in the communities: while Lactococcus ASV were dominant during the first propagation steps, they were subsequently replaced by Lactobacillus ASV sensu lato during the course of serial fermentation (Fig. 3).The relative abundance of Lactobacillus sensu lato was significantly higher in GP3 compared to GP2 (P = 6.2 × 10 −10 , WRSET) and in GP2 compared to GP1 (P = 2.2 × 10 −4 , WRSET).These results revealed that community structures from GP1, GP2, and GP3 exhibit different patterns of change.
The analysis of richness and Shannon index showed that GP1 was less diversified than GP2 and GP3 (Fig. S6B andE).In addition, the groups appeared to differ in the range of variation of these metrics over the course of serial fermentation (Fig. S6A and D).Therefore, the sums of the differences between adjacent fermentation steps (n + 1 and n) were calculated to report specifically on these variations in alpha diversity.The results show that Shannon diversity fluctuated significantly more in GP3 compared to GP1 and GP2 (Fig. S6F, P < 0.0098, WRSET).On the other hand, no significant differences were observed between the groups, regarding the variation in richness (Fig. S6C).
The non-metric multidimensional scaling (NMDS) ordination based on Bray-Curtis dissimilarity revealed that GP1 community lineages were stable but segregated, showing a high inter-lineage variability (Fig. 4D).GP2 lineages showed longer ecological trajectories than GP1, and their directions were divergent.In GP3, communities were highly dynamic and were characterized by longer trajectories.Interestingly, the direction was parallel for the majority of the GP3 lineages, showing that the change in taxonomic composition was similar in those communities compared to the other groups.These results show that the groups are characterized by different typologies of lineage trajectories and segregation.
A possible correlation was investigated between the acidification properties of the communities and Bray-Curtis dissimilarity.For this purpose, the difference in pH metrics between adjacent steps from F2 to F6 was calculated for each lineage.A significant correlation was highlighted between the delta of the final pH and Bray-Curtis dissimilar ity: the higher the pH variation, the greater the Bray-Curtis dissimilarity (Fig. 4C, r = 0.62, P = 2.5 × 10 −12 , PPMC).No significant correlation was obtained with the delta of lag time or MAR.This result reveals that the stability of community structures is related to functional stability.

Relationship between acidification properties and community structure
LAB are known to be responsible for acidification by producing lactic acid from lactose in milk.As they were highly abundant in the communities, we focused on the impact of LAB on acidification properties.The proportion of LAB only slightly impacted the acidification lag time when they represented less than 88% of relative abundance within the communities (P = 0.10 for the linear model, Fig. 5A).However, when the communities contained more than 88% LAB, the lag time significantly increased (slope = 0.090, P = 4.4 × 10 −3 for the linear model).This increase in lag time mostly concerned communities dominated by Lactobacillus sensu lato ASV (3.19 ± 0.16 h on average, Fig. 5D) compared to the other communities (from 1.95 ± 0.04 to 2.14 ± 0.05 h, P < 2.3 × 10 −5 , WRSET).
In addition, when the relative abundance of ASV corresponding to LAB increased, the speed of pH drop increased and the final pH decreased (P = 5.8 × 10 −6 for MAR, P < 2.2 × 10 −16 for the final pH, Fig. 5B and C).Accordingly, when the communities were dominated by Pseudomonadota ASV, the decrease of pH was slower compared to communities dominated by LAB (Fig. 5E).The MAR also depended on the LAB genus: the decrease in pH was significantly slower when the communities were dominated by Lactobacillus sensu lato ASV (−0.35 ± 0.01 pH units per hour) compared to Lactococcus (−0.39 ± 0.01 pH units per hour, P = 3.5 × 10 −2 , WRSET, Fig. 5E).On the other hand, a high variability in the speed of pH drop was observed when the communities were dominated by Streptococcus ASV (Fig. 5E).The final pH of fermentation was strongly related to the dominant ASV and was the highest for communities dominated by Pseudomonadota ASV (4.86 ± 0.04, P < 7.9 × 10 −5 , WRSET, Fig. 5F).Conversely, the final pH was the lowest when Lactobacillus sensu lato ASV (3.95 ± 0.02) were dominant (P < 3.8 × 10 −8 , WRSET, Fig. 5F).The communities dominated by the Lactococcus or Streptococcus ASV were characterized by intermediate final pH (4.41 ± 0.03 and 4.34 ± 0.04, respectively, Fig. 5F).
These results show that the acidification properties are linked to the proportion of LAB in the community and more specifically to the dominant LAB genus.

DISCUSSION
We have investigated the changes in bacterial community structures and their acidifica tion properties during experimental backslopping in milk.To do so, 26 raw milk samples were incubated, and the resulting communities were serially propagated in sterile milk.Communities were analyzed by metabarcoding, and acidification was recorded by pH monitoring.
The raw milk samples investigated in this study showed high bacterial diversity.The microorganisms present in raw milk can have different origins including the animal's teat, the dairy equipment, and the animal's surroundings.Although hygiene practices induce a strong decrease in microbiological load, leading to a stable low microbial load of approximately 10 3 -10 4 CFU/mL (37), raw milk microbiota is still characterized by high diversity (38)(39)(40)(41).The majority of raw milk samples were enriched in Pseudomonas, which are psychrotrophic bacteria.After dispersal, i.e., after transfer of microorganisms from the environment into the milk (42), milk samples were stored at 4°C for up to 72 h until fermentation.This step could have led to the enrichment of milk in psychro trophic bacteria as previously described (39,41).After the first fermentation step, a drastic reduction in bacterial diversity was observed.This reduction likely resulted from the strong selective pressure that applies during fermentation and only allows certain microorganisms such as LAB to persist (42).
Although the lag time for the first fermentation step was long (approximately 14 h), which is in agreement with previous reports (43), the successful spontaneous fermentation of all raw milk samples implies that they were all initially colonized by microorganisms able to acidify milk.In line with this, metabarcoding analyses revealed that all were colonized by LAB, which can be considered the main contributors to milk fermentation.Thus, even though LAB were present in very low relative abundance in raw milk, they were able to assure milk acidification, suggesting that adaptation to the milk matrix is more determinant than the relative abundance of bacteria initially present in raw milk.Besides LAB, several, but not all, fermented milk samples were colonized by other bacteria including species of the phylum Pseudomonadota.The impact of these bacteria on cheese is poorly described in the literature.Gram-negative bacteria such as Escherichia coli were exclusively considered spoilage bacteria for years, as members of this phylum are unwanted because of their animal intestine origin.However, some Gram-negative bacteria such as Psychrobacter celer or Hafnia alvei could have a beneficial impact on cheese by contributing to aroma development and protection against pathogens (44,45).Top-down community engineering through raw milk spontaneous fermentation could also be an interesting source of beneficial bacteria other than LAB, including Gram-negative bacteria.
Although richness decreased in our experimental engineering setup, it was still high compared to the communities previously described and obtained by advanced backslopping (22).Indeed, Swiss hard cheese starter cultures were colonized by two LAB species: Streptococcus thermophilus and Lactobacillus delbrueckii subsp.lactis.By contrast, most of the fermented milk samples obtained in our study were colonized by approximately 13 different ASV on average, including different species.This suggests that early backslopping, starting from raw milk, could be better suited to capturing microbial diversity and could be interesting for the design of starter cultures.
In this study, three groups based on the change in lag time, MAR, and final pH of the acidification curves were distinguished.GP1 showed a reproducible profile from F2 to F6, GP2 was more variable, and GP3 showed higher variations with intersecting curves.Metabarcoding analyses revealed that changes in community structures mirrored the acidification patterns: GP1 showed stable communities from F2 to F6, mainly domina ted by one LAB genus (Lactococcus or Streptococcus), GP2 showed moderate variations (with more diversified communities), and GP3 showed the highest variations (with a dominance of Lactococcus in the first steps and a dominance of Lactobacillus sensu lato in the last steps).Milk is a medium with high concentrations of nutrients.It was previously shown that serial propagation of communities in nutrient-rich media leads to communi ties with low diversity and high shaping (i.e., low stability) (46).In nutrient-rich media, bacteria greatly alter the chemical environment, causing more negative interactions between species.Strong species interactions could explain the reduction in bacterial diversity during milk fermentation performed in our study.However, besides the lineages characterized by low stability (in GP2 and GP3), some lineages were stabilized in a few propagation steps.The study of Ratzke et al. was conducted following a bottomup design approach where 15 different soil bacteria were mixed to build synthetic communities (46).The serial propagations conducted in our study involved 26 milk samples, each colonized by an average of 60 ASV.The higher diversity sampled in our study could have allowed us to highlight cases where community stabilization is possible and gives rise to the GP1 communities.The GP1 communities were less diversified than GP2 and GP3 communities.This suggests that conservative propagation during backslopping is facilitated when communities have a relatively low diversity, which would be consistent with data obtained in previous studies (22,43).
Community lineages with a stable final pH also presented reproducible metabarcod ing profiles.Furthermore, although GP1 communities are stabilized, they displayed not only interesting structural diversity (as shown by NMDS) but also functional inter-lineage diversity.Indeed, the final pH of GP1 communities differed significantly between the lineages in this group.Thus, in contrast to previous work (21,22), which showed a high degree of functional redundancy between communities obtained by advanced backslopping, here, we show that it is possible to obtain community lineages exhibit ing dissimilar community structures at the interspecific level, and presenting different functions.In cheesemaking, the pH varies strongly depending on the technology (47).The possibility of obtaining communities with different acidification properties opens up very interesting prospects for innovation in the field of dairy processing.
Our results showed that the dominance of LAB over Pseudomonadota led to both a lower MAR and a lower final pH.More precisely, these parameters depended on the dominant LAB genus Lactobacillus sensu lato, Lactococcus, or Streptococcus.The highest lag times were associated with Lactobacillus sensu lato, in comparison with the other LAB or the Pseudomonadota.In GP3, the ecological trajectories of the lineages were the longest and tended to be parallel, reflecting the microbial community succession in this group.Lactococcus was dominant in the first fermentation steps and then decreased while Lactobacillus increased and became dominant in the final steps.Interestingly, the level of dynamics was associated with the presence of Lactobacillus sensu lato: the higher the relative abundance of Lactobacillus sensu lato in the lineages, the higher the dynamics.This suggests that Lactobacillus sensu lato has a strong influence on the balance between community shaping and propagation.This role of lactobacilli on community dynamics could be due to their organic acid production and to their resistance to acid stress.Lactobacillus sensu lato can actually produce more acid and can be more resistant to acid stress than Lactococcus (7,48,49).In addition, Herreros et al. explained that lactobacilli metabolize lactose more slowly than lactococci (50).
The metabarcoding approach used in this study allowed us to investigate the microbial diversity at the genus and/or species level and is, therefore, not suitable for unraveling strain-level diversity.However, it has been shown that the microbiota of fermented dairy products can be characterized by significant strain-level diversity (22,51,52).It is, thus, reasonable to hypothesize that intraspecific diversity may have played a role in community assembly and functioning in our study.Similar to undefined starter cultures, where strain-level diversity leads to stable communities (22), the serially propagated communities characterized by structural stability at the genus/species level could be subject to intraspecific dynamics, where the relative abundance of different strains of the same species evolves during the course of propagations.Alternatively, strain-level diversity could have had a strong impact on community assembly as previously shown for model ripening cheese communities (52).Furthermore, although several lineages of communities exhibited similar structures, they could be character ized by slightly different functions.Indeed, even if M15 and M18 were stable and showed similar community structures (dominated by S. thermophilus), their final pH were significantly different (4.13 and 4.40), which could be explained by strain-level variabil ity.Further investigations are, therefore, needed to determine the role that strain-level variability may play in the assembly and function of cultivated communities.

Conclusion
In conclusion, this work has shown the wide diversity of community behavior in terms of dynamics, diversity, and trajectories, resulting in significant functional impacts in the framework of a top-down community engineering approach.The use of dairy backslopping at early stages is a powerful model to investigate the impact of microbial diversity and interactions on the balance between community shaping and propagation.Furthermore, this work on the top-down engineering of complex microbial communities offers very interesting prospects in all areas where expectations are high in terms of restoring microbial diversity in ecosystems.

FIG 2
FIG 2 Analysis of acidification kinetics.(A) Examples of acidification kinetics showing the different patterns in the three groups: M3, M2, and M1 for GP1, GP2, and GP3, respectively.The acidification parameters (means and standard errors) lag time, MAR, and the final pH are represented in panels B, C, and D, respectively.

FIG 4
FIG 4 Beta diversity analysis.(A) Mean Bray-Curtis dissimilarities between adjacent fermentation steps from F2 to F6. Error bars represent the standard errors of the means.(B) Sums of Bray-Curtis dissimilarities between adjacent steps for each group of communities.Different superscript letters indicate a significant difference (P < 0.05; pairwise comparisons using the Wilcoxon rank sum exact test with the Bonferroni adjustment method for the P value).(C) Bray-Curtis dissimilarities and delta of final pH between adjacent steps from F2 to F6. (D) Non-metric multi-dimensional scaling ordination representing ecological trajectories for each group.

(FIG 5
FIG 5 Relationship between the acidification parameters and the LAB proportion and composition.The lag time, the MAR, and the final pH as a function of the proportion of LAB in communities (A, B, C, respectively), or the dominant taxa in communities (D, E, and F, respectively).Different superscript letters indicate a significant difference (P < 0.05; pairwise comparisons using the Wilcoxon rank sum exact test with the Bonferroni adjustment method for the P value).