Initial community composition determines the long-term dynamics of a microbial cross-feeding interaction by modulating niche availability

Multi-step substrate consumption pathways can promote microbial biodiversity via cross-feeding. If one cell type preferentially consumes a primary substrate rather than the subsequently formed intermediates, then other cell types can specialize in consuming the intermediates. While this mechanism for promoting biodiversity is established, predicting the long-term persistence of such cross-feeding interactions remains challenging. Under what conditions will the interaction (and thus biodiversity) persist or disappear? To address this question, we propagated co-cultures of two isogenic strains of the bacterium Pseudomonas stutzeri. One completely reduces nitrate to nitrogen gas but preferentially reduces nitrate rather than nitrite (referred to as the generalist), while the other only reduces nitrite to nitrogen gas (referred to as the specialist). We found that the two strains coexist via nitrite cross-feeding when grown together, but the initial ratio of specialist-to-generalist (rS/G) determines the long-term dynamics of the co-culture. Co-cultures with large initial rS/Gs converge to the same rS/G and persist thereafter. Co-cultures with small initial rS/Gs also converge to the same rS/G but then become increasingly dominated by the generalist. The likely cause of these different dynamics is that the initial rS/G determines the initial environment, which in turn determines the initial selection pressures and phenotypes acquired by the generalist. Our results demonstrate that initial community composition controls the long-term dynamics and persistence of a cross-feeding interaction, and is therefore an important factor for community development and for engineering communities to achieve desired outcomes.

and removed non-relevant image areas (e.g., the overgrown central area of the agar plate or the area outside petri dish) and background pixels. E We then identified circular objects within a predefined parameter range (radius range, sensitivity factor) from the reflected light image.
We finally combined F the fluorescent image with G the identified circular objects to quantify frequencies of colonies expressing egfp or echerry. We grew the ancestral generalist and phenotype A of the generalist alone in batch culture in microtiter plates. We used OD 600 data to calibrate the growth yield, maximum growth rate, and intermediate toxicity parameters of our dynamical model. We first calibrated the model with data from the ancestral generalist grown at pH 7.5 (weak nitrite [NO 2 -] toxicity). Our objective was to describe growth under different conditions (nitrite toxicity) or of different phenotypes (phenotype A) with as few parameter changes as possible, rather than to make the best possible model fit to the measured data for each scenario. Solid lines are the measured OD 600 data and dotted lines are model predictions. A, B Ancestral generalist grown at pH 7.5 (weak nitrite toxicity). C, D Ancestral generalist grown at pH 6.5 (strong nitrite toxicity). E, F Phenotype A of the generalist grown at pH 6.5 (strong nitrite toxicity). Log ratio of specialist-to-generalist (r S/G ) Initial nitrate (NO 3 -) concentration (mM) Supplementary Fig. S6: Dynamics of the duration of nitrite (NO 2 -) availability during serial transfer in anoxic ACS medium. We mixed the generalist and specialist togeter at the indicated initial frequencies and serially transferred them with nitrate (NO 3 -) in batch cultures at pH 6.5 (strong nitrite toxicity). Colored symbols are means and colored error bars are one standard error of two biological replicates. Black symbols are means and black error bars are one standard deviation of eight biological replicates (generalist-only controls). ] toxicity). We initiated co-cultures with a specialist frequency of approximately 0.1% and acquired fluorescent photographs of spiral streaks from the co-cultures. The left column images are of plates streaked from the initial specialist/generalist mixture (T0). The center column images are of plates streaked from co-cultures after three transfers (T3). The right column images are of plates streaked from co-cultures after twelve transfers (T12). The first and third rows are reflected light images showing both specialist and generalist colonies while the second and fourth rows are the respective fluorescent light images (image of the same agar plate) where the specialist colonies are fluorescent and the generalist colonies are not. The first and second rows show one replicate co-culture where the specialist is tagged with egfp, while the third and fourth rows show a second replicate co-culture where the specialist is tagged with echerry. We grew individual generalist isolates alone in batch culture, and these generalist isolates displayed five distinct phenotypes. A Phenotype A immediately switches between nitrate (NO 3 -) and nitrite reduction, B phenotype B has a short time delay between nitrate and nitrite reduction, C phenotype C has a long time delay between nitrate and nitrite reduction, D phenotype D has a short time delay between nitrate and nitrite reduction and its cell density declines during that time, and E phenotype E has a long time delay between nitrate and nitrite reduction and its cell density declines during that time. A detailed description of how we delineated the phenotypes is provided in the Supplementary Text. Supplementary Fig. S9: Model simulations for co-cultures of the generalist and specialist initiated at different r S/G s. We varied the initial r S/G while holding the initial population size of A the generalist or B the specialist constant. The simulations are otherwise identical to those presented in Fig. 7.
Supplementary Table S1. Isogenic mutant strains of Pseudomonas stutzeri A1501 used in this study.

Supplementary Materials and Methods
Bioreactor experiment. To characterize the growth and nitrogen oxide reduction properties of the generalist, we grew the generalist in 1 L continuously-stirred gastight culture flasks at 34°C (referred to as bioreactors). We first purged 445.5 mL of a completely defined synthetic citrate-asparagine medium (ACS medium) [1] amended with 5 mM or 12 mM sodium nitrate (NaNO3), 10 mg/L gentamycin, and 1 mM sodium bromide (NaBr) with an anoxic nitrogen (N2):hydrogen (H2) (97:3) gas mixture while stirring for 45 minutes. This depleted dissolved oxygen in the medium to below the detection limit (0.3 µM) of our oxygen microsensor (Unisense, Aarhus, Denmark). We next inoculated the medium with 4.5 mL of an overnight culture of the generalist grown in oxic ACS medium. After inoculating the bioreactor, we continuously sampled small aliquots from the bioreactor using a Cavro XCalibur syringe pump (Tecan, Männedorf, Switzerland) and deposited the aliquots into individual wells of a microtiter plate sealed with a silicon seal to prevent evaporation using an XYZ positioning robot (Makeblock, Shenzhen, China). While depositing the aliquots into the microtiter plate, we used the syringe pump to mix the aliquots with glutaraldehyde (final concentration of 0.5% [vol:vol]), thereby fixing the cells and quenching further biological activity. We used these aliquots to quantify cell densities and extracellular nitrate (NO3 -) and nitrite (NO2 -) concentrations as described below.
We used flow cytometry to quantify cell densities during bioreactor operation. We first diluted the glutaraldehyde-fixed cells in 10 mM Tris EDTA (pH 8) buffer supplemented with SYBR Green I dye dissolved in dimethyl sulfoxide (ThermoFisher Scientific, Waltham, MA). We incubated the cells with SYBR Green I for 30 minutes and then counted the number of SYBR Green I-stained cells with a Gallios (Beckman Coulter, Brea, CA) or an Accuri C6 (BD Biosciences, Eysins, Switzerland) flow cytometer. We separated SYBR Green I-stained cells from background signal as described elsewhere [2,3].
We used ion chromatography (IC) to measure extracellular nitrate (NO3 -) and nitrite (NO2 -) concentrations during bioreactor operation. We first filtered the glutaraldehyde-fixed samples through 0.20 µm hydrophilic filters and diluted the filtrates 25-fold in deionized water. We determined extracellular nitrate and nitrite concentrations in the filtered and diluted samples using an IC Flex 930 chromatography system equipped with a Metrosep A Supp 16 -250/4.0 column (Metrohm, Herisau, Switzerland) as we described elsewhere [4]. We used chloride as an internal reference.

Spiral plating and fluorescence photography.
To quantify the ratios of the specialist to generalist (rS/Gs) of co-cultures, we used a spiral plating method in conjunction with fluorescence photography. This approach allowed us to obtain more than 2000 isolated colonies per LB agar plate, which was sufficient for us to accurately quantify rS/G at the lowest observed frequencies of the generalist or specialist and to distinguish the specialist from the generalist.
To perform spiral plating, we used a modified version of a method described elsewhere [5] and presented in Supplementary Fig. S3. We first diluted co-cultures in 0.85% sodium chloride (NaCl) solution by 200-fold. We then attached a fresh LB agar plate to the top of a mini lab-top centrifuge (500 rpm), dipped a sterile toothpick into the diluted co-cultures, and made a slow streak from the center to the rim of the rotating LB agar plate. This resulted in a dense spiral-shaped distribution of colonies on the LB agar plate (Supplementary Fig. S3D). To reduce autofluorescence from the LB agar and improve contrast between colonies and LB agar, we added 5 g of We adjusted the exposure time as needed while keeping the lens aperture set at f/5.6 and sensor sensitivity at ISO 800. We obtained better images once we eliminated exposure to ambient light. We therefore kept the LEDs, LB agar plates and the camera in a black box shielded from ambient light (Supplementary Fig. S3A) and used tethering software for camera control and image acquisition (Canon EOS Utility and digiCamControl). Indirect and even (from all sides) LED illumination were critical for two reasons. First, when illuminated directly, colonies were recorded as two circular objects; one was the colony itself while the other was surface reflection. Surface reflection typically resulted in detecting additional (false positive) objects during automated colony calling. Second, the surface reflection intensity was much higher than the fluorescence signal and could therefore lead to false colony identification during automated analysis.
We analyzed the photographs in Matlab (MathWorks, MA, USA). After background removal ( Supplementary Fig. S3D), round objects (in our case microbial colonies) were detected from the reflected light images with a circular Hough transform implemented with the 'imfindcircles' function. This allowed us to detect colonies without information of their fluorescence, thus avoiding analytical bias originating from the presence/absence of a fluorescence signal (Supplementary Fig. S3E). We took extra precautions to not to move the LB agar plate and the camera during imaging (both the camera and agar plate holder were surface mounted [Supplementary. Fig.   S3A]), which enabled us to retain spatial information -the colony center in the reflected light image therefore corresponds to the colony center in the fluorescence images. In the next step, we queried fluorescence images ( Supplementary Fig. S3F) for fluorescence intensity in the colony center area (Supplementary Fig. S3G; red colonies expressed the echerry fluorescent protein-encoding gene). For example, if the captured fluorescent light intensity in an image exceeded the threshold for the area corresponding to the center of a colony, then this colony would be considered fluorescent. As we used only a small area around the colony center for this analysis (81 pixels in 3456 x 3456 pixel images), colony size had no effect on colony binning. If light intensity was below a threshold for both red and green fluorescent proteins, then it was considered as tag-free. In some cases, two merged colonies appeared as one Delineating different individual-level phenotypes of the generalist. One explanation for the dynamics of the specialist is that the generalist diversifies its phenotype with regards to its nitrite (NO2 -) utilization properties. Such diversification of the generalist could include phenotypes that occupy or deplete temporally available niche opportunities (i.e., nitrite availability) for the specialist. To test for diversification in nitrite utilization properties by the generalist, we isolated 420 individual colonies of the generalist from one evolved co-culture series (6 x 70 colonies) serially transferred at pH 6.5 (strong nitrite toxicity), and from four (out of eight) control generalist-only control cultures (4 x 70 colonies). We measured growth properties of those colonies in anoxic ACS medium containing nitrate (NO3 -) as the growth limiting substrate and the pH set to 6.5 (strong nitrite toxicity). Briefly, we transferred each of the 70 colonies into a separate well of a single 96-well microtiter plate. As controls, we also transferred seven colonies of the ancestral generalist into wells of the same 96-well plate. We first archived the plate and then incubated it in anoxic conditions and collected OD600 measurements over time as described in the Materials and Methods section of the main text.
After completing the plate incubation and data acquisition, we first pooled all the recorded growth curves into one collection and identified distinct growth patterns using K-means clustering. We first performed K-means clustering along two properties; the duration of nitrite (NO2 -). availability and the maximum observed OD600 at the end of growth with nitrate (NO3 -). A subset of the generalist isolates did not show biphasic growth, and we could therefore not distinguish growth with nitrate from growth with nitrite. For these isolates, we set the maximum OD600 to that observed after the end of growth with both nitrate and nitrite. In total, this approach identified three main clusters based on the time lag between nitrate and nitrite reduction. One cluster had no time lag (Supplementary Fig. S8A; phenotype A), a second cluster had a short time lag (Supplementary Fig. S8B and S8D; phenotypes B and D), and a third cluster had a long time lag (Supplementary Fig. S8C and S8E; phenotypes C and E).
After this first round of K-means clustering, we performed a second round of K-means clustering. In this second round, we excluded the cluster that had no observable time lag between nitrate (NO3 -) and nitrite (NO2 -) reduction (Supplementary Fig. S8A; phenotype A) from the dataset and separated the remaining isolates based on the difference in the OD600 at the beginning and end of the time lag period divided by the duration of the time lag period, and based on the ratio of the OD600 at the end of growth with nitrate and at the end of growth with both nitrate and nitrite. We performed this clustering because the OD600 was clearly lower at the end of the time lag period than at the beginning for some isolates (suggesting cell death, aggregation, or deposition to the edges of the microtiter plate wells), and the difference in the OD600 was proportional to the length of the time lag period. This clustering round effectively separated those phenotypes that did not show a reduction in OD600 during the time lag period (Supplementary Fig. S8B and S8C; phenotypes B and C) and those that did (Supplementary Fig. S8D and S8E; phenotypes D and E).
When we combined the results from both K-means clustering rounds, we distinguished four phenotypes. Phenotype B has a short time lag between nitrate (NO3 -) and nitrite (NO2 -) reduction and a consistent OD600 over that time ( Supplementary Fig. S8B). Phenotype C has a long time lag and a consistent OD600 over that time (Supplementary Fig. S8C). Phenotype D has a short time lag and a decreasing OD600 over that time (Supplementary Fig. S8D). Finally, phenotype E has a long lag phase and a decreasing OD600 over that time (Supplementary Fig. S8E).

Genome analyses.
We sequenced the genomes of representative generalist isolates for each phenotype. Briefly, we streaked each generalist isolate onto LB agar plates, picked single colonies, and grew the generalist isolates overnight in oxic LB medium until reaching stationary phase. We then extracted and purified genomic DNA using the Wizard Genomic DNA purification kit (Promega, Madison, WI) and confirmed the quality of the DNA with gel electrophoresis. We then sent the purified genomic DNA to GATC Biotech AG (Konstanz, Germany) for sequencing with an Illumina HiSeq 4000 sequencer (Illumina, San Diego, CA). The sequencing was conducted with 150 cycles of paired-end sequencing and a target depth-of-coverage of 300-500x. Primary data and quality control analyses were performed using FastQC (Illumina, San Diego, CA).
We reported the complete set of parameters used for quality control elsewhere [6,7,8].
We identified putative genetic differences between each generalist isolate serially transferred at pH 6.5 (strong nitrite [NO2 -] toxicity) and its respective ancestor, including synonymous and non-synonymous mutations, insertions, deletions, and multiplications, in collaboration with the ETH Genetic Diversity Centre (Zürich, Switzerland). PRINSEQ-lite v0.20.4 was used to quality filter the raw reads, remove duplicate reads, and trim ambiguous base pairs [9]. breseq v.0.24rc5 was used to identify putative genetic differences relative to reference genomes [10]. The procedures and parameters used for sequence analyses are identical to those reported in our previous studies [6,7,8].
Mathematical models and simulations. We used adapted versions of models formulated and described in detail elsewhere that simulate primary substrate use and intermediate substrate production and use by multiple cell-types growing within the same completely-mixed physical compartment [4,11]. Growth of individual We applied the model to simulate the invasion of the rare phenotypes (including evolved generalists) during serial transfer in completely mixed batch reactors, as we performed in our experiments. Thus, the final community composition at the end of one round of batch growth was the initial community composition in the subsequent round of batch growth. We adjusted the initial cell density of the subsequent round of batch growth to account for the appropriate dilution; in our case, we diluted the community 100-fold prior to inoculating into the subsequent round of batch growth. We solved the system of differential equations with the 'odeint' function in the SciPy package in Anaconda Python distribution (version 2020.02). The complete model is as follows: Eq. 1. Yieldie is the yield coefficient for the evolved generalist growing with nitrite, and Yieldis is the yield coefficient for the specialist growing with nitrite. ) is the Andrews inhibition term that has been used previously to account for nitrite inhibition [4]. The term (1 + We could describe the growth of the evolved generalist by substantially increasing the term Km_se (0.5) compared to the term Km_sg (2 x 10 -6 ) in Equations 9 and 10. This approach simulates a scenario where the ability of the generalist to transport nitrate (NO3 -) from the bulk medium into the cell is impaired, thus reducing the accumulation of nitrite (NO2 -) and preventing its toxic effects.

Eq. 4.
Eq. 9. period. Phenotype E, on the other hand, has the lowest maximum OD600 measurements but a higher OD600 when growing with nitrate and a higher growth rate when growing with nitrite (NO2 -) ( Supplementary Fig. S8). Of the four phenotypes with observed time lags between nitrate and nitrite reduction, phenotypes B and D spend less time in the lag phase. (Supplementary Fig. S8). They are second only to phenotype A with regards to the maximum OD600 and have the highest growth rates at low cell densities (not significant for phenotype D). On the other hand, phenotype B has lower growth rates with nitrite than phenotypes D and E isolates. Overall, our data suggest many possible tradeoffs may exist between different growth properties with nitrate and nitrite. Genome sequencing of individual-level phenotypes of the generalist. We sequenced whole genomes of several randomly chosen representatives of each phenotype ( Supplementary Fig. S8). Of the 43 sequenced isolates (including several ancestral isolates), 19 have genetic differences when compared to their respective ancestor; 15 have one genetic difference, three have two genetic differences, and one has three genetic differences (Supplementary Table S2). All isolates belonging to phenotypes A, D, and E have at least one genetic difference, while only one isolate displaying phenotype B has a genetic difference. All isolates displaying phenotype D and phenotype E have genetic differences in genes associated with flagella/motility, either encoding structural proteins, or proteins involved in flagella regulation and folding (Supplementary Table S2).
Approximately 85% of the analyzed ancestral isolates display phenotypes classified as phenotype C, 15 % display phenotypes classified as phenotype B, and none display phenotypes classified as phenotypes A, D, or E (Supplementary Table S2). We sequenced six isolates that display phenotype B from serially transferred populations and three ancestral isolates that also display phenotype B. Aside from one isolate that displays phenotype B with a synonymous CAG>CAA mutation (Supplementary Table   S2), we could not detect any genetic differences in that phenotype when compared to