Dynamics of the Bacterial Community Associated with Phaeodactylum tricornutum Cultures

Abstract

The pennate diatom Phaeodactylum tricornutum is a model organism able to synthesize industrially-relevant molecules. Commercial-scale cultivation currently requires large monocultures, prone to bio-contamination. However, little is known about the identity of the invading organisms. To reduce the complexity of natural systems, we systematically investigated the microbiome of non-axenic P. tricornutum cultures from a culture collection in reproducible experiments. The results revealed a dynamic bacterial community that developed differently in “complete” and “minimal” media conditions. In complete media, we observed an accelerated “culture crash”, indicating a more stable culture in minimal media. The identification of only four bacterial families as major players within the microbiome suggests specific roles depending on environmental conditions. From our results we propose a network of putative interactions between P. tricornutum and these main bacterial factions. We demonstrate that, even with rather sparse data, a mathematical model can be reconstructed that qualitatively reproduces the observed population dynamics, thus indicating that our hypotheses regarding the molecular interactions are in agreement with experimental data. Whereas the model in its current state is only qualitative, we argue that it serves as a starting point to develop quantitative and predictive mathematical models, which may guide experimental efforts to synthetically construct and monitor stable communities required for robust upscaling strategies.

1. Introduction

Microalgae are photosynthesis-driven cells able to store light energy by converting carbon dioxide into carbohydrates, lipids, proteins, and other cellular components with potential biofuel, food, feed, and pharmaceutical and nutraceutical applications [1]. Novel applications also include the use of microalgae as an alternative sustainable development tool [2]. One such microalgae is the pennate diatom Phaeodactylum tricornutum that is able to synthesize a number of industrially relevant molecules applicable in: aquaculture as feed in e.g., bivalve, echinoderm, crustacean and fish hatcheries [3,4]; as biomass for biofuels [5,6]; pharmaceuticals and nutraceuticals [5,7,8,9]; and nanotechnology [10], and bioremediation industries [11]. To fully exploit the industrial potential of microalgal-derived products, substantial quantities of microalgal biomass is required, preferably obtained while maintaining low production costs. This is achieved by implementation of large-scale cultivation methods such as open raceway ponds and photobioreactors. The majority of conventional cultivation methods rely on keeping monocultures of the desired species, particularly if the final product is a bioactive molecule for human consumption [12]. Photobioreactors are closed systems that allow for the production of monoseptic cultures, fully isolated from potential contamination if cultivation protocols are followed correctly [13]. However, high operational costs of photobioreactors might not be sustainable. Another option is open raceway ponds, which are simple open-air cultivation systems that have been in use since the 1950s [1]. They are highly susceptible to contamination, and unless the desired species is a halophile or thermophile [14], it is hard to maintain monocultures. Irrespective of the cultivation method, the establishment of unwanted organisms such as amoeba, ciliates, rotifers, bacteria, viruses, and other photosynthetic organisms in microalgal cultures, is a serious obstacle for large-scale microalgae cultivation [15,16]. Although much research is carried out in the field of microalgal culture upscaling, very little is known about the true identity and characteristics of these invading organisms, responsible for microalgal culture “crashes” which lead to loss of biomass, and therefore, loss of revenue.

Microalgae are not found in monoculture in nature and it is not surprising that imposing such an artificial environment results in unstable large-scale cultures. By understanding rather than attempting to push out these micro-invaders, potential alternatives such as “synthetic ecology” as novel scaling up techniques should be explored [17]. This concept is based on the Competitive Exclusion Principle, or Gause’s Law, which states that two species competing for the same natural resource scarcely ever occupy a similar niche [18,19]. By “synthesizing” a community of organisms that fills every niche in the ecosystem of the microalgal culture and support, rather than harm, the growth of the phototroph, we would automatically optimize the utilization of nutrients and prevent the establishment of other potentially harmful organisms [17]. In order for synthetic ecology to be a legitimate contender as a novel scaling up technique, a greater understanding of species-specific interactions is required, starting with the bacterial faction, which are present in all of the Earths’ biomes [20], and arguably the key players in maintaining balance within a system. Theoretical ecology employs mathematical models to study the emergent patterns in ecosystems dynamics [21]. Because of the many industrial applications of microbial communities, current research has shown great interest in improving our understanding of such systems [22]. In particular, mathematical models and interdisciplinary approaches are fundamental to understanding the crucial underlying mechanisms that regulate community dynamics [23,24]. Since the same system can be inspected at different spatio-temporal scales and at a different degree of complexity, it is important to select the most suitable method to describe the biological phenomena under study in mathematical terms [25]. The first ecosystem models at the population-scale date back to the 1920s with the well-known Lotka-Volterra (LV) predator-prey model [26,27]. Since then, LV models have been extensively used to represent cooperation/competition population dynamics with a system of ordinary differential equations (ODEs). In generalized LV models (gLV) the system includes an arbitrary number of co-existing organisms and they directly represent positive/negative pairwise interactions as fixed parameters [28]. Today, a gLV model can be developed by inferring a co-occurrence network from a time series of metagenomics data [29]. This however requires a reasonable number of time resolved metagenomics data and will provide information on direct, one-on-one interactions only.

Diatoms and bacteria have co-evolved for more than 200 million years [30], and their co-existence is most likely based on a “biological barter trade system”, where substances such as trace metals, vitamins, and nutrients (nitrate, phosphate, silicate, carbon) are exchanged. In this work, we built on previous research that investigated algae-bacterial interactions including Provasoli’s work from 1958 where he suggested that bacteria can enhance the growth of algae [31], and subsequent species-specific studies that further corroborated his initial idea [32,33,34,35]. We first characterized the relative composition of the bacterial community in non-axenic P. tricornutum cultivated in the presence and absence of trace metals, vitamins and sodium metasilicate at different time points. Secondly, using critical peer-reviewed literature we defined the most likely functional roles of the bacterial factions and constructed a putative interaction network. Lastly, from the derived putative network of interactions, we built an ODE model with modified Verhulst equations [36] for microbial growth that included the direct effect of nutrient availability. Mortality rates were also introduced as dependent on specific bactericidal substances. The qualitative mathematical model, with parameters fitted to the available experimental data, served as a proof-of-concept that data as obtained here is sufficient to reconstruct a theoretical model that (a) reproduces the experimental observations, thus demonstrating consistency of our assumptions, and (b) allows for testing different hypotheses regarding the nature of the metabolic interactions underlying the ecosystem dynamics. It therefore represents a starting point to gain a deeper understanding of the principles of microbial community dynamics by an iterative experimental and theoretical approach.

2. Materials and Methods

2.1. Strains and Culture Conditions

All P. tricornutum strains were obtained from the Culture Collection of Algae and Protozoa (CCAP) based in Oban, Scotland [37]. All cultures were obtained non-axenic. Based on previous experimental evidence [38], the P. tricornutum strain CCAP1052/1B displayed optimal growth in 5L cultures. P. tricornutum was cultured in Guillard’s medium for diatoms (F/2 + Si) in filtered natural seawater chemically sterilised using sodium hypochlorite and sodium thiosulphate pentahydrate. P. tricornutum was grown in two media conditions; (1) complete F/2 medium containing sources of nitrogen (NaNO${}_3$) and phosphorus (NaH${}_2$PO${}_4$·2H${}_2$O), as well as trace metals and vitamins with the addition of sodium metasilicate, as per Guillard and Ryther 1962 [39] and Guillard 1975 [40], and (2) minimal media which contained just sources of nitrogen (NaNO${}_3$) and phosphorus (NaH${}_2$PO${}_4$·2H${}_2$O) at the same concentration as in the F/2 medium recipe. All cultures were grown in hanging 5L polyethylene bags with a “V” shaped bottom prepared using a heat sealer (Supplementary Figure S1). All cultures had a modified aeration system provided by a 10 mL pipette attached to the main pressurised air supply via 0.2 $\mathsf{\mu }$m sterile air filters. Cultures were kept at 18–20 °C and 24 h light at an average of 132.3 $\mathsf{\mu }$mol m^$-2$ s^$-1$. All cultures, irrespective of media condition, were inoculated with 250 mL from the same 5L stock culture of actively growing non-axenic P. tricornutum.

2.2. Growth Measurements

Growth was monitored every 24 to 48 h using a light microscope and carrying out cell counts of each culture in quadruplicate. During the cell counts the ratios of the four different morphotypes (oval, fusiform, triradiate and cruciform) were recorded, and descriptions of each culture noted. Samples of each culture were subsequently taken using a sterile 10 mL syringe and placed in 50 mL Falcon centrifuge tubes and placed in a −20 °C freezer.

2.3. Genomic DNA Extraction

All samples from days 1, 8, 15, and 22 were thawed in a water bath set at 25 °C. As per de Gouvion Saint Cyr et al. [41], samples were centrifuged for 5 min at 2000 x g to gather the P. tricornutum in the pellet while particles such as debris, other organisms, bacteria, and soluble substances remain in the supernatant. Because the bacteria might be attached to the P. tricornutum cells in the pellet, the pellet was washed with deionised water and then centrifuged for 5 min at 2000× g. This was repeated twice to ensure that majority of the bacteria attached to the pellet were released and were included in the community analysis. Genomic DNA extraction was carried out in the Aquaculture and Fisheries Development Centre, University College Cork. Mo Bio’s PowerWater^® DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA, catalogue No. 14900-100-NF) was utilized to carry out the genomic DNA extraction. Presence of gDNA was detected by running a 1% agarose-ethidium bromide gel with 72 wells. The samples were sent on dry ice to Heinrich Heine University, Düsseldorf, for the V6 16S sequencing.

2.4. Barcoded 16S-V6-Next Generation Sequencing

Ion Torrent™ barcoded Next Generation Sequencing protocol (Thermo Fischer Scientific Inc., Waltham, MA, USA) was used to sequence the bacterial gDNA [42,43]. Amplification of the V6 hyper variable region of 16S rRNA with forward and reverse primers (Supplementary Table S1) was carried out. Ion Reporter™ software (Thermo Fischer Scientific Inc., Waltham, MA, USA) assembled all the raw sequencing data and sorted all the reads using the unique sample-specific barcode sequences and removed them from the reads. The outcome was raw FASTQ files which were ready for analysis using bioinformatics tools.

2.5. Bioinformatics Analysis

A total of 87,077,374 reads were identified. The smallest sample was just over 1 million reads; the largest sample was just under 10 million reads. The sequencing data was subjected to a pipeline adapted and modified from Pylro et al. [44]. Primers were trimmed with fastq-mcf (version 1.04.807) [45], the resulting sequences were quality filterted and clustered into OTUs with usearch (version 8.0.1517; 32Bit-opensource) [46,47]. Taxonomy assignment was done by QIIME (version 1.9.0) [48] with the implemented uclust classifier based on 97% sequence identity to the reference 16S sequences from SILVA 111 database [49]. Statistical analyses were performed in R [50]. The complete protocol containing all processing steps is available on GitHub (see Supplementary Materials).

2.6. Mathematical Model

Starting from our understanding of the organism-to-organism interactions, we developed a dynamic model consisting of 13 ordinary differential equations (ODEs) and including 56 (55 free) parameters (see Appendix B.1). The model was built from the following working hypotheses:

(1)

the growth rate $\gamma$ of each population followed a standard Verhulst equation [36] parametrized with a carrying capacity and scaled by Monod-type terms [51] that describe the dependency on (micro)nutrients. These terms are in practice positive scaling factors <1.

(2)

the mortality rate of each population was inversely proportional to (1 + $\gamma$), to account for the fact that cells during replication (high growth rate) were healthier;

(3)

additional contributions to population mortality was given by the presence in the environment of noxious elements like bactericidal substances;

(4)

changes in metabolite concentrations are in general directly proportional to the growth $\gamma$ of the consumers and producers;

(5)

in the event of micronutrient scarcity (Iron and Vitamins in our model), P. tricornutum will secrete more organic carbons favored by those bacteria able to provide the needed micronutrients.

The initial conditions for simulations are different between minimal and complete media:

(1)

the initial quantity of Iron and Vitamins is 10 times higher in complete media;

(2)

the initial quantity of P. tricornutum biomass is matched to the first data point.

The parameters were fitted separately in minimal and complete media using a genetic algorithm [52] which was run in different steps to optimise the fit of P. tricornutum growth and/or the bacteria relative abundances to the experimental data in evolving system conditions (see Appendix B.2 and Supplementary Material 2). The model was written in Python (Python Software Foundation, https://www.python.org/) and is available on GitHub with instructions and scripts for running (see Supplementary Materials).

3. Results

3.1. Characteristics of Phaeodactylum tricornutum Growth

The media composition was shown to have a significant effect on the growth characteristics of P. tricornutum. P. tricornutum cultivated in minimal media exhibited a statistically significantly ($p=0.042$, unpaired Wilcoxon signed rank) higher cell density ($11.2\times {10}^6$ cells/mL) when compared to cultivation in complete media ($9.3\times {10}^6$ cells/mL). The growth rates during the exponential phase in both cultures were ${\mu }_{\mathrm{complete}}=0.43\pm 0.07$ d${}^{-1}$ and ${\mu }_{\mathrm{minimal}}=0.51\pm 0.04$ d${}^{-1}$ respectively. In contrast, the death rates when the cultures “crash” are ${\delta }_{\mathrm{complete}}=0.09\pm 0.02$ d${}^{-1}$ and ${\delta }_{\mathrm{minimal}}=0.08\pm 0.04$ d${}^{-1}$ respectively.

3.2. Bacterial Community Profile of Phaeodactylum tricornutum Cultures

Bacterial gDNA analysis showed that most of the operational taxonomic units (OTUs) could be assigned to the genera level (Supplementary Figure S2). Of the 9727 OTUs identified, 8109 corresponded to known sequences in the SILVA database (v.118) [49]. The OTU abundance at the phylum level showed that 99.97% of all OTUs belonged to Proteobacteria, Bacteroidetes, Actinobacteria and Firmicutes (Figure 1a). A comparison of the number of individual reads to the number of unique OTUs showed that the high number of reads per phyla was not the result of a single OTU (Supplementary Figure S3). OTUs with hits to known 16S P. tricornutum sequences were discarded.

Rarefaction curves were used to evaluate the alpha diversity in the different media conditions as well as at the different time points (Supplementary Figure S4). Species richness in both minimal and complete media was ∼3000. Species richness over time remained between ∼2400 and ∼2600, with reduced species richness (∼1300) on day 8 (both minimal and complete media) possibly due to elevated levels of 16S P. tricornutum chloroplast reads which had to be omitted. Greatest species richness (∼3000) was shown on day 22. All datasets showed a diminished increase in the number of unique species as the sample size increased, confirming adequate species richness in all culture conditions. To compare the species composition between the different samples (days/media) we used a non-metric multidimensional scaling (NMDS) function based on generalized UniFrac distances [53]. This allowed us to characterize the relationship between the particular samples on a visual level by displaying the information contained in the distance matrix. Therefore, similar samples would be placed together in an N-dimensional space. Here we observed a clear gradient of similarity between the bacterial samples from the different time points. The ordination based on the sampling day indicated that the bacterial community was dynamic with a clear divergence visible between day 1 and the other three sampling days. Interestingly the similarity between the different time samples showed the evolving processes of the community over time (overlaps between day 8 with day 15 and day 15 with day 22) and the recovery to the original one (overlap day 22 and day 1) (Figure 1b).

The existence of one dominant family at each investigated time point was a particularly interesting observation. In minimal media (Figure 2a), the lag phase of P. tricornutum growth was dominated by Pseudoalteromonadaceae (85%). However, during the log phase, a wide diversity of bacterial families was observed, with members of the Alteromonadaceae family (21%) beginning to dominate. During the stationary phase, a clear dominance of Alteromonadaceae species (55%) in the community was observed. The decline phase, however, showed the Pseudomonadaceae (39%) as a dominant family, with Pseudoalteromonadaceae species (37%) increasing in abundance again. In complete media (Figure 2b), the lag phase was also dominated by Pseudoalteromonadaceae (63%). During the log phase, 50% of the community was composed of members of the Flavobacteriaceae family, with the other 50% distributed among a number of different families. Flavobacteriaceae (46%) remained high in abundance during the stationary phase, with Pseudoalteromonadaceae species (44%) beginning to increase in abundance again. As for minimal media, Pseudoalteromonadaceae (57%) showed a clear dominance of the community during the decline phase.

An adapted version of PermanovaG was used to carry out permutational multivariate analysis of variance using multiple distance matrices which were previously calculated based on the generalized UniFrac distance [53]. The significance for the test was assessed by 5000 permutations. The results of the PermanovaG tests support the NMDS ordination, confirming a statistically significant effect in the bacterial community profile at the different sampling points and in the two media conditions whereas no significant effect was found in the experimental replicates (Figure A1).

3.3. Effect of Temporal Evolution and Media Composition on the Bacterial Community Profile

We compared the bacterial community profiles over time and in the different media conditions at the family level to avoid diluting the signal of the less abundant genera. Supplementary Figures S5 and S6 show no dynamical difference within the genera that cannot be observed at the family level. By investigating the bacterial community dynamics at the family level, we also included taxonomical information that is unavailable at the genus level. The families over-represented in all samples were Pseudoalteromonadaceae, Alteromonadaceae, Flavobacteriaceae and Pseudomonadaceae. Figure 2 illustrates the temporal evolution of the bacterial community in both minimal and complete media with a unique composition at each time point. A remarkable feature is that at all investigated time points there exist one or two dominant families.

In complete media, members of the Pseudoalteromonadaceae family were highly abundant when P. tricornutum cell densities were low (63% and 57% on day 1 and day 22, respectively). Flavobacteriaceae species dominated (50%) when the P. tricornutum culture was growing exponentially (day 8). Day 15, when P. tricornutum cell densities were at their highest, showed co-dominance of both Flavobacteriaceae (46%) and Pseudoalteromonadaceae (44%).

In minimal media, members of the Pseudoalteromonadaceae family were highly abundant when P. tricornutum cell densities were low. On day 22 Pseudomonadaceae (39%) and Pseudoalteromonadaceae (37%) were both overrepresented. When the P. tricornutum culture was in its exponential growth phase (day 8), a cluster of Families dominated; namely Alteromonadaceae (21%), Pseudoalteromonadaceae (20%), Pseudomonadaceae (18%), Halomonadaceae (15%) and Flavobacteriaceae (9%). When the cell density of P. tricornutum peaked (day 15), the Alteromonadaceae species took over (55%).

The bacterial communities within the two media conditions on day 1 were more closely related than the communities on days 8 and 15 (see Supplementary Table S2 for generalized UniFrac distances). As the cultures begin to “crash” (day 22), the bacterial communities in the two media conditions increased in similarity again. In general, the main families identified showed a distinct pattern of disappearance and regeneration within the bacterial community. In the complete media, Pseudoalteromonadaceae species started at 63% (day 1), dropped in abundance to 7% (day 8) then recovered to 57% (day 22). Flavobacteriaceae species, in complete media, started at 4.5% (day 1), increased in abundance to 50% (day 8), and then fell back to 8% (day 22). In the minimal media, Alteromonadaceae species had an abundance of only 1% (day 1), peaked at 55% (day 15), and decreased down to 18% (day 22).

3.4. Network of Putative Interactions between Phaeodactylum tricornutum and Identified Bacterial Families

The putative roles of each of the dominant families are illustrated in Figure 3. Based on an extensive literature review, five metabolites were identified as playing a crucial role in the interactions between P. tricornutum and the identified bacterial families. These are: bactericidal metabolites; iron; vitamins; dissolved organic carbons; dissolved organic phosphates.

• Bactericidal metabolites. Several species of the Pseudoalteromonadaceae family have been reported to possess bactericidal effects [54]. This ability to suppress the growth of competing bacteria could explain the dominance of Pseudoalteromonadaceae in almost all cultures irrespective of media composition. P. tricornutum also demonstrates bactericidal properties by excreting fatty acids (such as eicosapentaenoic acid or EPA), nucleotides, peptides, and pigment derivatives [55].
• Iron. Iron acquisition is essential for biological processes such as photosynthesis, respiration and nitrogen fixation. Bacteria produce and excrete siderophores, which scavenge iron. Diatoms are not known to produce siderophores, but genome sequence analyses identified the presence of a gene orthologue of a bacterial ferrichrome binding protein that suggests the possibility of iron (III)-siderophore utilization by P. tricornutum [56,57]. Furthermore, it was shown that P. tricornutum was able to uptake siderophores ferrioxamines B and E [58].
• Vitamins. Prokaryotes are thought to be the main producers of B vitamins [59,60]. Although P. tricornutum does not require cobalamin, thiamine and biotin [61], production of organic compounds such as EPA can be considerably enhanced by the bioavailability of co-factors such as cobalamin [62]. This provides the basis for potential mutualistic interactions. For example, Alteromonadales, dominant in our cultures, are thought to be capable of producing B vitamins [63].
• Dissolved Organic Carbon (DOC). It is estimated that up to 50% of carbon fixed via phytoplankton-mediated photosynthesis is utilized by marine bacteria [64], mainly as DOC compounds, defined as the organic material <$0.7\mathsf{\mu }$m in size [65]. DOC from diatoms originates either from live cells or recently lysed or grazed cells, which determines the type of DOCs available, and therefore are likely to influence the bacterial consortia associated with the diatom [30].
• Dissolved Organic Phosphate (DOP). Both diatoms and bacteria primarily utilize orthophosphate as a source of phosphorus. However, to access phosphate from DOP compounds, both diatoms and bacteria developed mechanisms to release orthophosphate (PO${}_4^{3-}$) from DOP. The mechanism is not species-specific, which consequently means the “free” orthophosphates can be acquired by any organism [66].

3.5. Mathematical Model Simulations

Based on the network of putative interactions between diatoms, bacteria, and the environment, we constructed a dynamic mathematical model, based on generalized Verhulst growth-laws [36] extended with Monod-type terms [51] to reflect the dependencies on metabolites (see Materials and Methods and Appendix B.1). Figure 4 presents results of the model simulations after the model parameters were fitted to the data in minimal and complete media conditions, respectively (see Appendix B.1). Experimental data are superimposed. The top panel shows biomasses of the five organisms (data available only for the diatom), the bottom panel shows relative bacteria abundance versus time (individual biomass divided by total bacterial biomass). Because of the qualitative nature of the model, units are arbitrary. The figures show that the model is able to reproduce the main features of the bacterial community dynamics, such as the disappearance and return of Pseudoalteromonadaceae in complete media and the peak of Alteromonadaceae at the end of the exponential growth phase of P. tricornutum in minimal media. Supplementary Figure S7b and S7d show the dynamics of metabolite concentrations, for which no data are available.

Due to the large number of free parameters, the fit was certainly not unique. Supplementary Material 2 presents additional checks we performed on the parameter fitting procedures. The parameter space could be in principle reduced to 43 free parameters, but this did not change the results. It was not possible to find a unique parameter set valid in both minimal and complete media conditions. This is however consistent with the fact that the model was not constructed to capture effects like metabolic re-adjustments, something that would be observed e.g., as a different parameter value for growth rate or metabolite consumption. With the data available, it was not possible to make any quantitative statements about the actual interaction parameters, neither could it be assumed that simulation results are of general validity. Despite these limitations, the model did represent a possible configuration of diatom-bacteria-environment interactions, which was in agreement with the experimentally observed bacterial dynamics.

In order to test the stability of the bacterial community and how it supports the growth of P. tricornutum we ran simulations using the same set of parameters (either minimal or complete media conditions) and varied the initial community composition removing one bacteria per simulation. In complete media the simulated growth of P. tricornutum still fit the experimental data, rather independently from the bacterial community (Figure 5). Also under axenic conditions, diatom growth was predicted to be largely unperturbed. This situation was different in minimal media. While under these conditions diatom growth was also unaffected upon removal of the three bacterial families Pseudoalteromonadaceae, Flavobacteriaceae and Pseudomonadaceae, removing Alteromonadaceae from the community resulted in the total absence of P. tricornutum growth (Figure 6). This behavior was expected from the hypothesized central role of Alteromonadaceae in supplying the diatom with micronutrients. Surprisingly, the removal of a single bacteria from the community in both media conditions still gave, in general, a good fit of the (recomputed) relative abundances, except when removing Alteromonadaceae in minimal media (as direct consequence of what stated previously) and when removing Pseudoalteromonadaceae in both media conditions. This hinted to a relevant role of Pseudoalteromonadaceae in regulating the community composition through its predatory strategy of releasing bactericidal substances. Finally, the community composition at the last time point was overall better captured. This suggests that the mathematical model consistently captured the general interactions leading one bacterial family to dominate over the others on the long term.

4. Discussion

In nature, P. tricornutum does not exist as an isolated entity. In fact, it is part of a complex ecosystem whose complete network of interactions with both its environment and other organisms remains poorly understood. Microbial ecosystems are of high interest for a wide range of applications in fields, such as medicine, renewable energy, and agriculture. Within the scope of this project, the complexity of a natural, variable system was reduced by investigating the batch growth of non-axenic laboratory strains of P. tricornutum from a culture collection. The cultivation method we developed was designed to compromise between highly controlled small-scale laboratory conditions and a large-scale industrial set-up. The bacterial community was characterized and the community dynamics investigated in two conditions: minimal and complete media. The data was then complimented with an extensive study of existing, peer-reviewed literature to identify the putative role of the dominant bacterial families associated with P. tricornutum. We then validated the derived network of interactions by developing a mathematical model which could reproduce the observed dynamics. The presented approach, integrating experiments, bioinformatics and mathematical modeling, illustrates a possible way towards the development of a monitoring pipeline for non-axenic microalgae cultures.

4.1. Experimental Observation of the Dynamics of the Bacterial Community Associated to Phaeodactylum tricornutum

The growth dynamics of P. tricornutum in the two media conditions showed an accelerated “culture crash” in the complete media compared to the minimal media, which indicated a more stable culture in the minimal media (Figure 2). This also suggested that non-axenic cultures of P. tricornutum might not require expensive trace metals and vitamins for optimal growth under the conditions provide, an observation crucial to the large-scale industrial cultivation of P. tricornutum as this would drastically decrease the production costs. Simultaneously, the dynamics of the bacterial community revealed that the community in the minimal media increased in complexity over time. The link between ecosystem complexity and stability based on theoretical and experimental data has been debated by ecologists for over half a century [68,69,70,71]. Our observations are in agreement with more recent hypotheses indicating that diversity generally increases the stability of an ecosystem [72].

4.2. Literature-Based Assessment of the Putative Role of Each Bacterial Family

The bioinformatics analysis of bacterial gDNA abundance showed clear dominance of four bacterial families: Pseudoalteromonadaceae, Alteromonadaceae, Flavobacteriaceae and Pseudomonadaceae. These bacteria were over-represented in all samples and their relative abundances showed different temporal dynamics among the two P. tricornutum growth conditions. In order to understand the putative functional role of these bacteria an extensive study of peer-reviewed literature was carried out.

Pseudoalteromonadaceae. Members of Pseudoalteromonadaceae family have been isolated from coastal, open and deep-sea waters, sediments, marine invertebrates, as well as marine fish and algae [73]. The Pseudoalteromonadaceae family has three genera, namely Pseudoalteromonas, Algicola and Psychrosphaera [74]. Several species of Pseudoalteromonadaceae are reported to possess antibiotic properties with bactericidal effects [54]. For example, concentrated supernatant of a marine bacterium Pseudoalteromonas sp. strain A28 contained various enzymes including proteases, DNAses, cellulases, and amylases, capable of causing the lysis of the diatom Skeletonema costatum [75]. Species of Pseudoalteromonadaceae are also capable of producing cold-adapted enzymes [76,77,78,79,80,81]. Pseudoalteromonadaceae species can produce extracellular polymeric substances allowing them to colonise surfaces, enhancing nutrient uptake whilst limiting diffusion of particular substances across the cell membrane [82]. The ability of Pseudoalteromonadaceae species to suppress the growth of competing bacteria could explain the dominance of Pseudoalteromonadaceae in almost all cultures irrespective of media composition, particularly when P. tricornutum abundance is limited (Figure 2, days 1 and 22). P. tricornutum on the other hand, may protect other bacterial community members from the bacteriolytic ability of Pseudoalteromonadaceae by producing specific antibacterial compounds themselves. Desbois et al. showed that P. tricornutum excreted bacteriolytic fatty acids such as eicosapentaenoic acid (EPA; 20:5n-3), nucleotides, peptides, and pigment derivatives that can eliminate unwanted competition for nutrients such as organic phosphates from certain bacteria [55].

Alteromonadaceae. The Alteromonadaceae family consists of 16 (yet annotated) named genera found predominantly in marine environments [74]. Members of this family were isolated from nutrient-rich environments such as coastal, open, and deep-sea waters, sediments, marine invertebrates and vertebrates, algae, and temperate and Antarctic marine environments [83]. They are able to utilize a vast array of compounds as carbon sources; from glucose to glycerol [74]. Members of this family are known siderophore producers [57,84,85]. Greek for “iron carrier”, siderophores are a group of iron scavengers that act by chelating iron (III) that are produced and excreted by bacteria, and some cyanobacteria, which then reuptake the siderophores with bound iron (III) via outer-membrane transporters that are siderophore-specific [86]. Most bioactive trace metals, including iron, exist at nanomolar to picomolar concentrations in our oceans, approximately one-millionth of the intracellular concentration in diatoms [87,88]. No trace metals, including iron (III), were provided to minimal media cultures. However, natural seawater may contain minute traces of bioactive trace metals. The high abundance of Alteromonadaceae in the minimal media suggests a potential supportive role in sequestering traces of iron (III) that may be present in the sterile natural seawater to the P. tricornutum (Figure 2). This is further supported by the very low level of Alteromonadaceae in the complete media (11% in complete media compared to 55% in minimal media, both on day 15) where the culture has been supplied with 11.7 $\mathsf{\mu }$M of iron (III) chloride hexahydrate.

Flavobacteriaceae. Flavobacteriaceae are members of the Bacteroidetes phylum and include over 120 genera found in soil, sediments and seawater (see [89] for further references). Flavobacteriaceae belong within the Cytophaga-Flavobacterium cluster which has been shown to account for more than 10% of the total bacterial community in coastal and offshore waters [90,91,92]. Members of Flavobacteriaceae can proficiently degrade various biopolymers such as cellulose, chitin and pectin [93,94]. They were shown to be omnipresent during phytoplankton blooms, and their preference for consuming more complex polymers rather than monomers suggests an active role in the processing of organic matter during these blooms [95,96]. Although the exact mechanisms behind them are not perfectly understood, algal blooms are a consequence of exponential growth of phytoplankton [97]. In this respect, the phase of exponential growth of P. tricornutum in complete media, when our results showed highest abundance of Flavobacteriaceae, is the artificial equivalent of an algal bloom of P. tricornutum (Figure 2). In the minimal media, the abundance of Flavobacteriaceae remains very low; at its maximum on day 8 it only accounts for 9% of the total bacterial community. Members of the Flavobacteriaceae family could be more demanding than other bacteria that require lower nutrient levels to thrive.

Pseudomonadaceae. Pseudomonadaceae are an extraordinarily diverse family of bacteria found in almost all habitats on Earth; in soils, freshwater as well as marine environments, as well as plant and animal-associated pathogens [98]. Species from the Pseudomonas genus are the best studied of the Pseudomonadaceae family, whose sheer genetic diversity explains the ability to thrive in such a wide range of environments [99]. Marine isolates from the Pseudomonas genus have been shown to produce a wide range of bioactive compounds, many of which exhibit antibacterial as well as antiviral properties (see [100] for further references). Our results, indeed show an elevated level of Pseudomonadaceae OTUs evident on day 22 of the complete media cultures, and on days 8 and 22 of the minimal media cultures. The increased presence of Pseudomonadaceae when the P. tricornutum culture has “crashed” could be attributed to its ability to produce antibacterial compounds allowing members of this family to begin to thrive in the community through inhibition of its competitors. Given its exceptional genetic diversity, and thus, its metabolic versatility, allows for members of Pseudomonadaceae to be truly saprophytic; providing a hypothetical explanation of its abundance we could measure when the P. tricornutum cultures crash (Figure 2, day 22 in both media conditions).

4.3. Putative Network of Interactions and Validation with a Qualitative Mathematical Model

The literature review work revealed interesting insights into the possible metabolic exchanges going on and allowed to infer interaction links among P. tricornutum and its associated bacterial community. We critically considered which metabolites were most relevant for survival (organic carbons for the bacteria, iron, vitamins and phosphates for the diatom) and which ones could play a role in competition and predation among the microbes (bactericidal metabolites). From these considerations we designed a putative network of interactions that was then translated into a mathematical model. In particular we chose, besides the five microbes’ biomasses, a total of eight metabolites as the variables that directly and specifically influence the interactions among the different organisms. These were: four different possible sources of organic carbons, each preferred by a different bacterial family [30]; two bactericidal substances, one produced by Pseudoalteromonadaceae and affecting all other bacteria [54], the other produced by P. tricornutum and targeting specifically Pseudoalteromonadaceae [55]; vitamins produced by Alteromonadaceae and needed by P. tricornutum [63]; bio-available iron that is chelated by Alteromonadaceae and efficiently absorbed by P. tricornutum [57]. For the scope of the model, we ignored other free iron forms that can be uptaken by all bacteria as well as phosphates that are not species-specific and are present in both minimal and complete media.

Direct metabolic exchanges are known to be central in microbial community interactions [101], but usually population dynamics models like gLV [28] do not include this information. This work, therefore, modified the standard formulation of the Verhulst equation [36] for bacterial growth to include organism-to-organism interactions depending on the production/consumption of metabolites, modeled as Monod-type terms [51]. Nutrients availability can indeed drastically change the “metabolic state” of an organism, inducing a reprogramming of resource allocation to face nutrient scarcity [102]. This was shown at the gene expression level for example for Escherichia coli [103] and Shewanella oneidensis [104] grown in minimal and rich media condition. Micronutrients can as well affect microbial gene expression, as is for example the case with vitamin B${}_{12}$, whose presence induces the expression of the cofactor-dependent methionine synthase enzyme METH, while in its absence the cofactor-independent methionine synthase enzyme METE is expressed [105]. An ODE model at the population level cannot, of course, capture mechanisms such as metabolic shifts caused by changes in the environment such as the supplementation of minimal or complete media [25]. Therefore, we did not expect to find a unique set of parameters for the model in the two conditions. However, the parameters fitted to the data of P. tricornutum with four bacterial families still provide good fits in simulations with altered community composition. Even though the parameter values could not directly be interpreted biologically, we could use the simulation results on metabolites dynamics (an information absent in data) to speculate about the reason for the lower cell count of P. tricornutum in complete media with respect to minimal media (Supplementary Material 3). In minimal media, Alteromonadaceae maintained constant iron levels and the fitted values for parameters characterizing P. tricornutum’s sensitivity to micronutrients levels were significantly lower. This would suggest a key role of Alteromonadaceae in supporting P. tricornutum growth combined with the diatom’s adaptation to scarce micronutrients availability.

Considering the limited information that can be extracted from the current experimental data available, the model we proposed is purely qualitative and provides a proof-of-concept that a quantitative model can, in principle, be constructed if dedicated experiments are designed for calibration. The current qualitative model provides therefore a preliminary validation of our putative network of interactions, and serves as motivation for further research bringing the model to a quantitative, predictive level. Indeed, starting with systematic measurements of model parameters in co-cultivation experiments, the simulations can gain predictive power and become a powerful tool towards the goal of synthetic community design and control.

5. Conclusions and Outlook

This study demonstrated that the bacterial community associated with non-axenic laboratory strains of P. tricornutum is not randomly assembled but follows dynamics that can be reproduced. We postulate that a role within the community can be filled by a number of bacterial species capable of carrying out a certain function (guilds) rather by one specific species of bacteria. Which bacteria fill the role is dependent upon the environmental characteristics and the prevailing needs of the community as a whole at any given time. Unfilled niches will be seized by bacteria with the ideal metabolic functionality. The absence of certain micronutrients creates a new niche that can be filled by a certain unique bacterial faction. Further work is necessary to explore the hypotheses postulated and to further develop the qualitative mathematical model to understand the specific community roles and the ecological niches. In the context of fundamental research, one approach would be to carry out systematic time-resolved omics studies, which provide a holistic view of the genes (genomics) and metabolites (metabolomics) in a specific biological sample in a non-targeted and non-biased manner [106], and use them to develop an “expanded gLV” mathematical model where the species specific interaction terms depend on the metabolite concentrations. This would allow to derive a network of interactions independent of a priori hypotheses, and thus represent a significant step forward in understanding community dynamics based on metabolic exchanges. In the context of industrial scale-up, systematic co-culture experiments with culturable members of the bacterial families of interest, chosen based on desired functional roles, could be used to parametrize a mathematical model like the one we presented and develop it into a powerful predictive tool for culture monitoring. For example, samplings assessing the community composition can be used to predict the harvesting point and avoid “culture crash”. The development of novel co-cultivation strategies for scale-up is extremely relevant for pharma- and nutraceutical, as well as animal feed industries. Therefore there will be increasing interest in further research into co-cultivation approaches and in general in the field of synthetic ecology.