The impact of phage treatment on bacterial community structure is minor compared to antibiotics

Phage treatment is suggested as an alternative to antibiotics; however, there is limited knowledge of how phage treatment impacts resident bacterial community structure. When phages induce bacterial lysis, resources become available to the resident community. Therefore, the density of the target bacterium is essential to consider when investigating the effect of phage treatment. This has never been studied. Thus, we invaded microcosms containing a lake-derived community with Flavobacterium columnare strain Fc7 at no, low or high densities, and treated them with either the bacteriophage FCL-2, the antibiotic Penicillin or kept them untreated (3 × 3 factorial design). The communities were sampled over the course of one week, and bacterial community composition and density were examined by 16S rDNA amplicon sequencing and flow cytometry. We show that phage treatment had minor impacts on the resident community when the host F. columnare Fc7 of the phage was present, as it caused no significant differences in bacterial density α- and β-diversity, successional patterns, and community assembly. However, a significant change was observed in community composition when the phage host was absent, mainly driven by a substantial increase in Aquirufa. In contrast, antibiotics induced significant changes in all community characteristics investigated. The most crucial finding was a bloom of γ-proteobacteria and a shift from selection to ecological drift dominating community assembly. This study investigated whether the amount of a bacterial host impacted the effect of phage treatment on community structure. We conclude that phage treatment did not significantly affect the diversity or composition of the bacterial communities when the phage host was present, but introduced changes when the host was absent. In contrast, antibiotic treatment was highly disturbing to community structure. Moreover, higher amounts of the bacterial host of the phage increased the contribution of stochastic community assembly and resulted in a feast-famine like response in bacterial density in all treatment groups. This finding emphasises that the invader density used in bacterial invasion studies impacts the experimental reproducibility. Overall, this study supports that phage treatment is substantially less disturbing to bacterial communities than antibiotic treatments.

Inspection of the control samples indicated that 127 ASVs were contaminants, and these were removed from the dataset.Additionally, we removed 4 ASVs identified as archaeal and another 127 identified as Chloroplast.Two samples (DNA_66.C.P.2.D3_S17, DNA_29.T.N.5.D1_S5) were of low quality due to low sequencing depth and were removed from the dataset.After the removals, the high-quality dataset contained 7 894 016 sequencing reads and 4630 ASV.The 141 samples had an average sequencing depth of 55 985±17 363 (mean±SD) reads.Rarefaction curves indicated that most samples had flattened off at the minimum sampling depth of 26 448 reads, with an average slope of 0.0084±0.004(Supplementary Figure 2).For α-diversity analysis, the dataset was normalised by scaling each sample to 26 448 reads before each sample was rarefied to 26 447 reads to ensure equal sampling depth.69 ASVs were removed from the dataset due to normalisation.Before analysis, we inspected the normalised 5 of 21 dataset to ensure it reflected the full dataset satisfactorily.First, we inspected the correlation between the normalised and full dataset regarding Hill diversity numbers (N order ) of order 0 and 1 [3] (Supplementary Figure 3a).The correlation coefficient was 1 for N 1 and 0.92 for N 0 .While the linear regression slope between the normalised and full dataset was 0.98 for N 1 , it was only 0.67 for N 0 with increasing deviation from the 1:1 ratio with higher sample ASV richness (N 0 ).These highly rich samples had a higher sequencing depth (data not shown).Thus, while only 69 ASVs were lost, highly diverse and deeply sequenced samples lost more of their rare taxa due to the normalisation method.
Next, we visually inspected the community composition differences through PCoA plots based on the Bray-Curtis and Sørensen similarity (Supplementary Figure 3b and c).For the Bray-Curtis-based ordinations, the sample-clustering was highly similar between the normalised and full-dataset, while it was more condensed for the Sørensen similarity.Thus, the normalised dataset did indeed represent the core-microbiome stronger.However, the datasets were deemed to reflect the overall trends satisfactorily.

Quantification of bacterial density using flow cytometry
The total bacterial density and the proportion of living-dying, and dead cells were estimated throughout the experiment.The total bacterial density was quantified from fixated samples after the experiment was conducted, while the live-dead estimations were conducted daily from fresh samples.For both types of acquisitions, samples were diluted in 0.2µm filtered phosphatebuffered saline (PBS, 1x) to obtain stable sample aquation and dilute background noise present in the sample.

Total bacterial density
Each sample was first diluted to reduce occurrences of coincidence events.Next, SYBR green II stock (10 000x) was diluted in 0.2µm filtered PBS (1x) to 200x.This solution was used to stain samples to a final concentration of 2x.Stained samples were incubated for 15 minutes in the dark at 37°C before analysis.Five samples were stained at a time and data were acquired before staining more samples to ensure a stable fluorescent signal.Data were collected using the blue laser (488 nm) with detection in BL1 (530±30 nm) and BL3 (>630 nm) using a BL1 threshold of 400-800 (depending on the sample).Instrument voltages were as follows; FSC 320V, SSC 340V, BL1 420 V and BL3 400V and YL2 620V.Data was acquired by running 110µL of the sample at a 100 µL/min flow rate.Each sample was vortexed before data acquisition.0.2µm filtered PBS and 0.2µm filtered lake water were used as negative controls, and the Fc7 as a positive control.
All samples were gated similarly with minor modifications to each sample gate to ensure high-quality data (Supplementary Figure 4).The gating strategy was as follows: First, a stable signal collection was identified in the time-histogram plot, and the stable events were gated.
Then all co-coincidence events were excluded by removing events with a higher BL1-A (area) than BL1-H (height) signal.Lastly, background noise was identified by comparing 0.2µm filtered lake water with non-filtered lake water and was gated out.The remaining events were identified as the bacterial community.Supplementary Figure 4: Gating strategy for total bacterial density quantification.A stable signal was first selected in the "stb" gate.Then all coincidence events were removed by only keeping events in the "no-coin" gate.We kept events in the "not-cigar" gate based on the negative controls and then gated out the bacterial populations in the SSC-H vs BL1-H scatterplot.a) 0.2µm filtered PBS solution was used as a negative control and for diluting the bacterial communities.b) 0.2µm filtered lake water and c) The bacterial community in the microcosms (here showing reactor replicate one from the control group added low amounts of F. columnare Fc7 (T-C-1) at 0 DPI).

Living-, dying-and dead density quantification
For quantifying the living-, dying-and dead-bacterial population, we stained 1mL of the sample within one hour of sampling daily from replicate 3 in each experimental group.SYBR green I (10 000x, Invitrogen) and propidium iodide (PI, 20 mM, Invitrogen) were diluted in 0.2µm filtered PBS (1x) to a final concentration of 100x and 0.4mM.Each sample was stained with this solution giving a final dye concentration of SYBR green I (1x) and PI (0.4µM).SYBR green I and PI bind to DNA and emit a green and red fluorescent signal if exited at 488 nm, respectively.SYBR green I can enter all cells, while PI only enters membrane-compromised cells (i.e.dying-or dead cells).
Samples were incubated for 15 minutes in the dark at 37°C with the stain before acquiring data.Data were collected using the blue laser (488 nm) with detection in BL1 (530±30 nm), BL3 (>630 nm) and YL-2 (620±15 nm) using a BL1 threshold of 400-800 (depending on the sample).

Statistical analysis workflow in R Bacterial density polynomial model
To test for differences in bacterial density between the control and the treatments, we fitted a third-degree polynomial mixed effect model with the log 10 transformed bacterial density as the response variable and treatment, sampling day, propagule pressure and the interactions between these as the explanatory variables.We used lmer() from the lme4 package (version 1.
To test for differences in ASV richness between the control and the treatments, we fitted a second-degree mixed effect polynomial model with richness as the response variable and treatment, sampling day, level of F. columnare Fc7 added and the interactions between these as the explanatory variables using lm().Sampling day was included as a random effect.As for bacterial density, emmeans() and emmip() were used to obtain the estimated marginal means.

Bacterial composition
Bacterial composition was evaluated based on both relative and absolute ASV abundances.To obtain absolute abundances, the ASV table was first transformed into compositional counts.Then each sample was scaled by the quantified bacterial density, and the ASV table was rounded to the closest integer.Finally, bacterial composition bar plots were generated using functions from the package microshades (version 1.10, github.com/KarstensLab/microshades).
ASV1 was identified to be the added F. columnare Fc7.To test for statistical differences in the relative abundance of F. columnare Fc7 between treatment groups within each level of F. columnare Fc7 added, we used kruskal.test()and pairwise.wilcox.test()from the stats package.

Community composition comparisons
We estimated β-diversity using the Bray-Curtis and Sørensen similarity by averaging 100 similarity matrixes generated from randomly subsampled datasets using the function avgdist() from vegan [5,7].The datasets were subsampled to 26 448 reads.Principal coordinate analysis (PCoA) was used to visualise the variation in β-diversity and conducted using ordinate() and plot_ordination() from phyloseq.Within each sampling day and level of F. columnare Fc7 added, a permutational analysis of variance (PERMANOVA) with 999 permutations tested for differences in centroid location and dispersion between samples using the function adonis2() from vegan.
Differences in group dispersions were tested using betadisper() from vegan.We performed 100 PERMANOVA and dispersion tests for each subset and averaged the test statistics.
To identify genera that changed in abundance due to the treatments, we performed differential abundance tests on samples taken at day 7.We used the packages corncob [8] (version 0.3.1),DESeq2 [9] (version 1.38.1) and ANCOMBC [10] (version 2.0.1) for differential abundance testing as different tools vary in which taxa are identified [11].As input to all tools, we filtered the absolute abundance scaled ASV table to only contain ASVs with a total absolute abundance of over 2500 ASVs/mL and a prevalence of over 5%.We used the function differentialTest() from corncob with the Wald significance test without bootstrapping separately for each propagule pressure (n=3 groups x 3 tests).For DeSeq2 and ANCOMB, the data were grouped according to treatment type and propagule pressure (n=9 groups).From DeSeq2, we used the function DeSeq() with parametric dispersion negative binomial GLM fitting and Wald significance testing on median ratio normalised count data.From ANCOMBC, we used the function ancombc2() with pairwise differential group testing and Holm p-value correction.We investigated the ASVs identified to be significant by all three methods.Because minor differences in abundance can be statistically, but maybe not biologically significant, we focused on ASVs with over a 5-fold difference in absolute abundance between the control and the two treatments.

Community assembly
We quantified the change in replicate microcosm similarity per day to investigate the community assembly [12].First, the Bray-Curtis and Sørensen similarity was quantified for each pair of replicate microcosms on each sampling day.Next, a mixed linear regression was performed with similarity as the response variable and sampling day, treatment type, amounts of added F. columnare Fc7 and the interaction between these as the explanatory variables.Sampling day was added as a random variable.Finally, the model's coefficient estimates and 95% confidence intervals were collected using coefficient() and confint() from stat.
Of particular interest was the temporal slope rate for each treatment.In the framework of Gundersen et al., 2021, the temporal slope is interpreted as a similarity rate.Positive similarity rates indicate that the community composition between two replicates became more similar over time, indicative of selection dominating community assembly.Negative rates, on the other hand, reflect that replicates became less similar over time, indicating that drift structured the community assembly.We quantified this similarity rate of change for each experimental group.

Supplementary Table 2:
The ratio between the estimated marginal mean in bacterial density for phage or antibiotic to the control treatment and the ratios standard error (SE), the degrees of freedom for each test (df)

Supplementary Figure 2 :
Rarefaction curves for all samples coloured according to the amount of F. columnare Fc7 added.The samples are organised based on the treatment type and sampling day.The solid vertical line indicates the sequencing depth chosen for normalisation at 26 448 reads.

Supplementary Figure 3 :
Comparisons of the normalised and full dataset in terms of a) Hill diversity of order 0 and 1 (N 0 =richness, N 1 =exponential Shannon), PCoA ordination based on b) Bray-Curtis similarity and c) Sørensen similarity.The normalised dataset contained 4561 ASVs and had a total loss of 69 ASVs compared to the full dataset.Colours indicate the treatment type and shape the sampling day.Abbreviations; Treatment: None = no treatment, Phage = phage treatment (FCL-2), AB = antibiotics (Penicillin).

Supplementary Figure 5 :
Samples were acquired by running 110µL sample at a 100 µL/min flow rate.Each sample was vortexed before data acquisition.0.2µm filtered PBS was used as a negative control.The gating strategy was similar to the total bacterial density.The live, dying, and dead populations were identified in the BL1-H vs BL3-H scatterplot (Supplementary Figure 5).Example of live-, dying-and dead-population estimates in the microcosms added high amounts of F. columnare Fc7 at 0, 1, 5 and 7 DPI.The live-, dying-and dead-population was quantified in the BL1-H vs BL3-H scatterplot.Living = yellow gate, dying = blue, red = dead.Percentages indicate the proportion of events in the gate.The sample dilution factor is indicated above the plots.Note that both phage-and antibiotic treatment increased the proportion of dying and dead cells compared to control after 1 DPI.
1-  34)  to create the model and emmeans() from the emmeans package (version 1.8.2-090007, github.com/rvlenth/emmeans) to obtain estimated marginal means for the density ratio difference between control and either phage-or antibiotic treatment.To estimate the conditional R 2 we used the r.squaredGLMM() function from the MuMIn package (version 1.47.5).P-values obtained from emmeans() were Dunnett adjusted within each sampling day, level of F. columnare Fc7 added and treatment comparison (n=42 comparisons).95% confidence intervals were obtained from emmip() and used for plotting the temporal change in bacterial density.
, 95% confidence interval (lower-upper CI), the t-ratio value statistics, and the corresponding p-value (Student t-test based).P-valueswere Dunnett adjusted within each day and propagule pressure level (n=2 tests).