Prokaryotes dominate the biosphere and form diverse communities disrupted by invasion. Invaders and remaining community members experience resource surfeit, competition, and selective pressures. Little is known about invasion in natural microbial communities. We examined invasion by chemotaxis in a meso-tube system at taxonomic, functional, and genomic levels as communities sank, rose, and formed a chemotactic band that migrated for metres. The band velocity increased as the community migrated despite non-motile bacterial hitchhikers and up to 10⁶ viruses/ml. Migrating communities left complex residual communities in their wake, showing dynamic taxonomic composition and adaptation through increased migration-associated genes. Approximately 500 species migrated together, competing for dominance. This system offers a superior method for studying band and residual community dynamics, bacterial hitchhiking, viral transport, gene evolution, and survival strategies, revealing cohesive communities that persist over extended distances. Our methods and results provide an experimental foundation for investigating microbial invasion in multiple ecological settings.
Article
Co-migration of hundreds of species over metres drives selection and promotes non-motile hitchhikers
https://doi.org/10.21203/rs.3.rs-5040873/v1
This work is licensed under a CC BY 4.0 License
Version 1
posted
You are reading this latest preprint version
Bacteria live in complex, heterogeneous communities that develop inter- and intraspecific interactions1,2. The heterogeneous physical, chemical and taxonomic environmental structure requires bacteria to compete for resources, space and nutrients3–6. However, intense competition and continually changing environments mean some bacteria have few resources7. As a result, bacteria have evolved many mechanisms to increase their chances of survival through sporulation, dormancy, competition and migration. Migration provides the opportunity to find new nutrient sources, avoid toxins, elude predation by viruses, escape nutrient-depleted environments, and reduce competition8,9.
The microcapillary assay is valuable for studying chemotaxis10–12. Nutrient-filled microcapillary tubes are inoculated with a single bacterial species at one end, which collectively migrates away from the inoculation site as a sharp, visible band at a constant speed, generating a self-sustained signal gradient as they consume nutrients. This method and its adaptations have been fundamental to our modern understanding of chemotaxis13–15. However, this research has focused on single, isolated bacterial species, which does not accurately represent the diversity of natural environments. Other experimental systems used to study chemotaxis include agar plates6,16,17 and microfluidic devices18,19. Regardless, these systems are limited to short distances and one or a few microbial species.
Here, we examine bacterial community migration dynamics at the genus, species, and gene levels by quantifying the extent of collective bacterial migration. Key objectives included determining whether a species would exclude all others during continuous migration if continual migration confirmed the ‘eat to escape phages’ hypothesis17, assessing how quickly non-motile community members were left behind, and detecting gene selection over one week. To achieve these goals, we performed experiments with mixed communities of bacteria originating from sewage and monitored their speed, abundance, taxa and changes in dominant genes. We observed nomadic communities, a subset of a microbial community that emerges and migrates as a chemotactic band in the direction of nutrients. These nomads included viruses and diverse microbial taxa, selected by their functional capabilities. Our research reveals that the collaborative dynamics among complex communities are key to accelerating migration speed and sustaining functionally diverse microbial communities throughout migration despite changes in community functional potential.
Chemotactic bands form and migrate within a modified tube system
We devised a modified chemotaxis tube assay using transparent vinyl tubing filled with agar, nutrients and methylene blue as a redox indicator. Three variations of this assay were performed to investigate microbial communities' formation and migration, including the residual community left behind.
Band formation experiments were completed using 35 cm-long lengths of tubing inoculated at one end with sewage (Fig. 1a). As methylene blue was reduced, the media changed from blue to colourless. Fifteen hours after inoculation, colourless patches formed along the bottom of the tube. These extended 6 ± 1 cm (standard error of the mean (SEM)) down the length of the tube (Supplementary Table 1). After 18 hours, large agar sections became colourless, with the bottom third of the tube being completely clear (Fig. 1a). At 18 hours, clear areas began to expand upwards into the upper two-thirds of the previously patchy areas towards the top of the tubing (Extended Data Fig. 1). The colourless sections continued to migrate 12.5 ± 2.7 cm down each tube. By 21 hours, the bacteria had travelled 13.4 ± 2.6 cm down the tube length (Supplementary Table 1) and formed a sharp chemotactic band 1–3 mm wide. As the band migrated, the agar behind it became colourless (Fig. 1b).
To investigate the long-term migration of the established communities, the bands were extracted and transferred to 1m-long tubes (Fig. 1a). The bands reformed after approximately 20 hours. They were left to migrate along the tubes for 1 week, with additional transfers made to fresh tubes if the band reached the end of the tube.
Experiments were repeated 3 times with 5 replicate tubes for each experiment and showed initial speeds of 0.53 (± 0.09), 0.32 ± (0.05) and 0.47 ± (0.00) (SEM) cm/hr. These increased to 0.73 (± 0.04), 0.84 (± 0.03) and 1.0 (± 0.1) cm/hr (Fig. 1c-e, Supplementary Table 2). Speed increases were polytonic, and fluctuations were present throughout migration.
The abundance of bacteria and virus-like particle (VLP) populations during community migration and after one week of migration were analysed via flow cytometry (Supplementary Fig. 1). In the inoculum, there were 2.25 ± 1.49 x106 bacterial particles/mL and 4.35 ± 2.04 x106 VLPs/mL, in the chemotactic band formed after 21 hours, there were 9.67 ± 6.01 x106 bacterial particles/mL and 4.62 ± 7.02 x106 VLPs/mL and after one week of migration the chemotactic band contained 6.51 ± 0.65 x106 bacterial particles/mL and 1.24 ± 1.47 x106 VLPs/mL (Supplementary Table 3, Extended Data Fig. 2).
A separate tube assay investigated the residual community left behind as the band migrated forward. After the band reached the end of the tube, community samples were taken along the tube at 3 cm intervals.
Migrating communities shift and support hundreds of prokaryotic species
We investigated community richness during band formation, migration, and the residual community. After formation, the 15 chemotactic bands contained between 355 and 643 bacterial genera and between 1296 and 2741 bacterial species. This decreased to between 187 and 282 genera and 514 to 806 species after 1 week of migration (Supplementary Fig. 2–3, Supplementary Table 4).
The residual community showed this decrease in species richness was not gradual. Within the residual community, there were 2,938 species and 968 genera in the inoculum, which rapidly decreased over the first ten samples before stabilising and remaining at approximately 500 species for all preceding residual communities (Fig. 2a, Supplementary Fig. 4–5, Supplementary Table 5). Despite this decrease, 4,688 unique species were identified across the residual communities indicating species not detectable in the initial community appeared in later samples. In the first three residual communities, a smooth relationship is observed in the rank abundance curves (Fig. 2b). As the community begins to stabilise, a break emerges, indicating a separation between a few highly abundant species and hundreds of rare species (Supplementary Table 5). After stabilisation, the dominant species in the residual communities fluctuated in abundance, sometimes disappearing below detection limits and then reappearing later (Fig. 2c). However, the rare species, including non-motile species, were more consistent in their abundance than the dominant species (Fig. 2d).
Taxonomic sorting during community establishment
During community formation, dynamic changes were observed in the abundance of bacterial genera (Fig. 3a, Supplementary Table 6). Over 70% of the average relative abundance in the inoculum comprised 363 different genera (via 16S sequencing). Three hours into community formation, this proportion decreased to 0.023 of the average relative abundance. Acinetobacter sp. contributed to the highest average relative abundance 3 and 6 hours into community formation. Aeromonas sp. had the highest relative abundance for the remainder of community formation, with the highest being 0.32 relative abundance after 12 hours of community formation (Fig. 3a). The presence of the genera Lactococcus spp., Acinetobacter spp. and Klebsiella spp. at 18 hours, indicated non-motile taxa in the established band formed more than 5 cm from the inoculation point. A Bray Curtis dissimilarity principal coordinates analysis (PCoA) showed that 63% of the total sample variation could be explained by 2 principal coordinates (Extended Data Fig. 3). The initial samples all grouped together with later time points diverging. The two, 6-hour samples and one 12-hour sample grouped together. All remaining samples grouped together, with 12 and 9 hours having the most variation
Taxonomic changes during community migration
After community formation, incremental changes are observed in the taxonomic composition of the residual community as the chemotactic band migrates. This is shown in the PCoA of the relative abundance of bacterial genera (Fig. 3b). As residual communities progress down the tube, samples gradually move along the first principal component that explains 58.26% of the total community variation. This principal component strongly correlates with the distance of the community from the inoculation site, with a Pearson correlation coefficient of 0.905. Fluctuations among high-abundance species cause these incremental changes, while species with low abundance are consistently present (Fig. 2c-d). After these incremental changes, a single bacterial species, Lelliottia amingena dominates the chemotactic band with a relative abundance of 0.861. Similar incremental changes are also observed in the PCoA of the viral component of the residual communities. The first principal component captured 48.42% of the total variation and strongly correlated with the distance of the community from the inoculation site with a Pearson correlation coefficient of 0.956 (Fig. 3b).
Communities were compared at band formation and after one week of migration. After migration, a bacterial species with low relative abundance in the established community became numerically dominant (Fig. 3c). This numerically dominant species was not the same in each replicate community. Leclercia adecarboxylata dominated five replicate communities with relative abundances of 0.916, 0.877, 0.347, 0.912 and 0.908. The remaining five tubes were numerically dominated by Enterobacter sp. RHBSTW-00994, Enterobacter kobei, Enterobacter roggenkampii, Rahnella aceris and Lelliottia amnigena with relative abundances of 0.851, 0.588, 0.912, 0.883 and 0.928 respectively. In each tube, the numerically dominant species in the final community was not numerically dominant in the initial community before (Extended Data Fig. 4).
PCoA was used to assess differences in the overall composition of bacterial genera during migration. Significant separation was observed between the bacterial taxonomic profiles at the species level after community formation and after long-term migration (PERMANOVA, P = 0.0001) (Extended Data Fig. 5a). After band formation, communities form a tight group, while samples collected after long-term migration experience high dispersion (PERMDISP, P = 0.0033). Significant separation was observed between the viral taxonomic profiles at the species level after migration (Extended Data Fig. 5b, PERMANOVA P = 0.0001), coinciding with increased dispersion after long-term migration (PERMDISP, P = 0.0001). No significant differences were observed in the microbial (PERMANOVA, P = 0.8037) or viral (PERMANOVA, P = 0.8265) components of communities arising from the same sewage inoculum.
Non-motile taxa were observed during community formation and migration and in the residual community. Non-motile species included Acinetobacter, Klebsiella, Lactococcus, Raultella, Streptococcus, Tolumonas and Veillonella (Fig. 3a, 3c, Supplementary Fig. 6).
Functional shifts occur during long-term migration.
We investigated the changes in the functional potential of chemotactic bands after migration and in the residual community. PCoA on the level 3 subsystem profiles was used to investigate functional changes before and after long-term migration. Significant separation was observed in the functional profiles after migration (PERMANOVA, P = 0.0001); and in the dispersion of samples (PERMDISP, P = 0.0038)(Extended Data Fig. 5c). Differential expression analysis determined which functions were altered during migration. The volcano plot shows that 92 level 3 subsystems experience a significant difference (Wald test with Benjamini-Hochberg correction, P < 10− 15) of at least 0.5-fold from the community post-formation, indicating many functions decrease in abundance (Fig. 4a). These decreased functions are largely involved in metabolic processes (Supplementary Table 7). Alternatively, eight functions experience a significant increase (P < 10− 15) of at least 0.5-fold during migration. These enriched functions include the two-partner secretion pathway (P = 6.96 x 10− 17) and bacterial chemotaxis (P = 1.27 x 10− 20) as well as the metabolic functions malonate decarboxylase (P = 1.65 x 10− 117), pyrimidine utilisation (P = 7.08 x 10− 22), alkyl phosphate utilisation (P = 8.49 x 10− 17) (Fig. 4b).
As chemotaxis functions were enriched during migration and are essential for the movement of the migrating band, we further examined how specific chemotaxis functions within the Bacterial Chemotaxis subsystem changed during band migration. We observed that methyl-accepting chemotaxis protein genes significantly increased in abundance during migration (Welch’s t-test), including the serine chemoreceptor protein (MCP I)(P = 1.22 x 10− 3), aspartate chemoreceptor protein (MCP II)(P = 2.13 x 10− 2), ribose and galactose chemoreceptor protein (MCP III)(P = 8.05 x 10− 4) and dipeptide chemoreceptor protein (MCP IV)(P = 5.60 x 10− 3). Alternatively, chemotaxis protein methyltransferase cheR decreased in abundance during band migration (P = 1.40 x 10− 2) (Fig. 4c).
The PCoA of the residual community demonstrates that incremental changes occur in the functional composition of the communities left behind with the band separate from the residual communities (Fig. 4d). These functional patterns mirror the changes observed for the bacterial and viral profiles with the first principal component correlated with distance from the inoculation site with a Pearson correlation coefficient of 0.875.
Metagenome-assembled genomes reveal ecological shifts.
Across residual community samples, the cross-assembly of 29 metagenomes recovered eleven metagenome-assembled genomes (MAGs) with completeness > 70% and contamination < 15% (Supplementary Table 8). Two separate bins, S1C52 (completeness = 87.59%, contamination = 9.99%) and S1C53 (completeness = 77.75%, contamination = 5.29%) were annotated as L. adecarboxylata. The NCBI taxonomy classifies them both as the same species, whereas the Genome Taxonomy Database toolkit classifies them as different species belonging to the same genera (Leclercia adecarboxylata and Leclercia adecarboxylata_C)20. Regardless, these genomes are similar, on average, 90.6% identical with an alignment coverage of 74.3% (Extended Data Fig. 6). The abundance of these two bacteria varied throughout the length of the experiment. We used read coverage as a proxy for genome abundance, and initially, S1C52 was more abundant, but S1C53 became dominant after 8cm of migration. In the last two samples, before the band formed, S1C52 resumed community dominance (Fig. 5a). S1C52 and S1C53 have 224 SEED subsystems in common; S1C53 encodes 59 SEED subsystems not present in S1C52, while S1C53 encodes 26 SEED subsystems not present in S1C53. Specifically, S1C53 contains 15 genes in the motility and chemotaxis SEED class absent from S1C52 including genes encoding the flagellar basal body, flagellar biosynthesis, and flagellar hook proteins (Fig. 5b, Supplementary Fig. 7).
Two high-quality MAGs were also obtained for two known non-motile bacterial species, Lactoccocus raffinolactis (completeness = 94.66%, contamination = 9.77%) and Lactococcus chungangensis (completeness = 89.95%, contamination = 4.26%). These MAGs showed no motility or chemotaxis genes (Fig. 5b, Supplementary Fig. 8). Two MAGs contained chemotaxis genes but no flagellar genes, Tissierellales GPF-1 sp. (completeness = 99.05%, contamination = 4.83%) and Peptostreptococcus sp. (completeness = 99.02%, contamination = 1.79%).
Overview. The demonstration of group bacterial migration along microcapillaries catalysed the investigation of bacterial chemotaxis and cell signalling21,22 as critical infection and biosphere processes22. Chemotaxis is crucial for pathogen invasion, resulting in commonly associated maladies, such as gastric, bladder, and colon cancer23. In the biosphere, it is associated with enhanced global oxygen production and nutrient cycling24. This makes laboratory chemotaxis systems important for understanding bacterial behaviour. The meso-tube system we present provides opportunities to study mixed microbial community resilience, adaptations and selection over comparatively long times and distances. Here, co-migrating communities were maintained in the form of visible, accelerating bands composed of complex communities that we term nomadic communities.
Band formation. The classic microcapillary chemotaxis experiments begin with a pre-formed band where the cells must swim. In our meso-tube system, the band is an emergent property after sinking through an agar matrix. Sinking organic microparticles draw bacteria in by acting as micro-patch resources and down by their negative buoyancy25 and the negatively buoyant cells add to the sinking26,27. Bacteria use oxygen as an electron acceptor and the formed micro-patches reduce methylene blue resulting in a patchy appearance. Eventually, 1–3 mm wide tube-filling bands formed and migrated down the tubes. In these bands, oxygen is evenly depleted, and the bands migrate.
Acceleration. For our multi-species migrating bands, we expect the interspecific differences in chemotaxis and metabolic abilities to cause poorly performing species to be lost out of the back of the band. At the same time, we expect the good performers to advance28,29. This self-organisation of bacteria by their chemotactic abilities may explain the accelerations, but details of how this happens remain unclear. Within the meso-tube system, the matrix intermittently confines and organises the forward movement, reducing sideways deviation30. Our results suggest that the intermittent self-organising events are increasingly auto-correlated, implying an increasing hysteresis as the band progresses. The observed self-organisation may also drive species dominance change in the nomadic community.
Taxa loss. The residual and band communities diverged in their taxa abundances and dynamics. Residual community taxa decreased from 2,938 species to 422 species, and the band community decreased from 2,000 species in the inoculant to a range of 514 to 806 species at the ends of the tubes. The residual communities’ rare taxa were at a near-constant abundance while the dominant species fluctuated (Fig. 3C, D). Rare biosphere communities persist31, but the maintenance mechanisms have been difficult to study. We demonstrate rare biosphere persistence and the ability to study it over extended periods.
Hitchhikers. Acinetobacter, Streptomyces, Klebsiella, and Lactococcus, are considered non-motile in agar environments because they lack flagella and cannot move significantly through the matrix. Specifically, Acinetobacter's movement is constrained by its embedded nature in the agar and insufficient motility speed, relying on twitching motility, which is too slow and ineffective in agar, as it requires solid surfaces to function properly 32. The presence of these non-motile bacteria from the beginning to the end of the tubes can only be explained by hitchhiking33, whereby non-motile bacteria directly or indirectly utilise the motility and chemotaxis of other microbes to invade environments. These behaviours are often beneficial for both partners. For example, when non-motile bacteria encode an antibiotic-resistance mechanism, they become ‘cargo’, degrading antibiotics that allow both partners to occupy otherwise toxic environments34,35.
Escape from phage. The “eat to escape” (ETE) hypothesis proposes that chemotaxis is a response to escape bacteriophage (phage) infection, and that nutrient uptake is secondary. However, metagenomic analysis shows a range of phages persist throughout long-term community migration (Fig. 3C). Further, the detection of VLPs in chemotactic bands via flow cytometric analysis suggests that lytic phages remain present as communities migrate. Without a direct, and thus far undiscovered, mechanism that bacteria use to detect phage presence and drive chemotactic machinery, our data support the alternate hypothesis that bacteria are migrating to obtain greater nutrients while avoiding predation is a secondary benefit. Therefore, we refer to this as the “weak eat to escape” hypothesis. The continued presence of phages may play a crucial role in shaping migratory populations by impairing the chemotactic machinery of infected bacteria, preying on susceptible cells, and driving rapid bacterial evolution through gene transfer.
Functional depletion. Species loss does not necessarily mean functional loss as it might in plants and animals, as microbial communities contain extensive functional redundancy36. Continual chemotaxis places strong selection pressure on competing species and functional groups. This competition, combined with essential nutrients provided by the media, has caused a depletion or loss of some functions. For example, the depletion of the B12 synthesis gene is consistent with B12 presence in the growth medium, and the depletion of the attachment fimbriae gene is consistent with the need for continual movement (Fig. 4c). Whether these genes are eventually lost from this contained artificial ecosystem and how long that takes will be an informative future experiment37.
Functions enriched. Given the imposed requirement of continuous migration, we expected the community to be enriched in chemotaxis genes (Fig. 4b). The autoinducer enrichment was unexpected. These changes are consistent with the underlying commonalities of chemotaxis systems among diverse taxa and communities becoming increasingly structured over time38,39. We hypothesise that these enrichments will be universal to all systems. However, the enrichment of the various transporters is likely medium-specific. The capability of our meso-tube system to study functional enrichment will be useful when a library of transporter enrichments is created from many inoculants and compared to metagenomes directly from environments. Within the chemotaxis subsystem, cheA and MCP1 were extensively enriched (Fig. 4c). CheA is the primary control molecule for chemotaxis and has a default setting of smooth swimming with extended migration40.
MAGS-based species diversity. The 11 MAGs reveal the heterogeneity within the community and highlight the genomic divergence between L. adecarboxylata strains S1C52 and S1C53. This divergence has functional implications, as seen in the differences in relative abundance and gene content between these strains. Notably, S1C52 lacks key flagellar genes, including those for the flagellar basal body, biosynthesis, and hook proteins, which may suggest that S1C52 engages in 'cheating' by exploiting the motility of S1C53 through hitchhiking. The observed shifts in coverage patterns indicate that these two strains may occupy distinct ecological niches via niche differentiation41. For instance, the presence of flagellar genes in S1C53 may provide a selective advantage under the experimental conditions. Cheating behaviours add complexity and instability to microbial communities, preventing any single strain from dominating and thereby maintaining diversity and promoting the coexistence of multiple taxa42–44.
The recovery of high-quality MAGs for non-motile species supports their initial presence and aligns with the absence of motility and chemotaxis genes. The Tissierellales GPF-1 sp. and Peptostreptococcus sp. without flagellar genes require further investigation and may reflect the complexity of microbial community dynamics, where factors such as nutrient availability, interspecies interactions, and environmental stressors play critical roles in shaping community structure and function.
We extend classical microcapillary assays by increasing the assay's size, duration, and complexity from a single bacterial species to entire microbial communities with hundreds of species that co-migrate for metres. Using our continuous migration system, we observe nomadic communities composed of hundreds of bacterial and viral taxa migrating as chemotactic bands, leaving behind complex and ever-changing residual communities. The migrating bands and their residual communities enriched or depleted specific genes. Meanwhile, rare taxa maintained a near-constant abundance, while dominant taxa abundance fluctuated extensively. The invasion of tissue, soil, aquifers, and other nutrient sources by bacteria is key to biosphere existence. The experimental system presented is closer to natural systems and constitutes a powerful system to study mechanisms that sustain rare biospheres within different environments.
Sample Collection
Sewage samples were collected from the Glenelg SA Water Wastewater treatment facility, South Australia (34°57'49"S 138°30'38"E). Three, 250 mL samples of inlet sewage were collected on separate days (8/9/20, 13/04/21, 21/04/2021,15/06/2021). After collection, each sample was sealed and transported to the laboratory. A 10 mL aliquot was taken from each sample which was used for all further analysis. All samples were stored at 4°C before use.
Bacterial tube migration experiment
Three types of tube migration experiments were completed to investigate different qualities of migrating communities. We refer to these as community formation, long term migration and residual community tube experiments. The details of each type of tube experiment and the methods used to prepare the tubes are described below.
Tube preparation
For all three types of tube experiment, a length of transparent polyvinyl tubing with an internal diameter of 19 mm was sterilised by rinsing with 10% hydrochloric acid then rinsing with phosphate-buffered saline (65 mM NaCl, 1.2 mM KCl, 5.4 mM Na2HPO4, 0.52 mM KH2PO4) to neutralise the acid. This rinsing procedure was repeated three times for each tube. Once dry, each tube was sealed at one end with sterilised parafilm and attached using electrical tape. Tubes were then filled with media containing 0.25% agar and 0.5% tryptic soy broth and sealed at the opposite end again using sterile parafilm and electrical tape. Filled tubes were inoculated with 240 µL of sewage inoculum, 1 cm away from one end of the tube using a 25-gauge syringe. The inoculation site was then covered with clear adhesive tape to prevent leakage of media from the tube and any possible contamination from the external environment. Each experiment was conducted in a biosafety cabinet at ambient room temperature between 21–22°C.
Community formation tube experiments used tubes 35 cm in length, were inoculated with sewage and were left to migrate for 21 hours. Community formation experiments were repeated using three different sewage samples with five replicate experiments for each sewage sample. Methylene blue dye (3x10− 4%) was added to these tubes as a reduction-oxidation indicator to visualise oxidation due to bacterial metabolism changing the agar from blue to clear and representing the presence of bacteria within the system. Samples and observations of the tube were then taken at 3 hourly intervals for 21 hours until the community was seen to be moving as a band. Long-term migration tube experiments used the final sample taken from each community formation tube experiment as the inoculum for a 1m-long tube experiment (3 sewage replicates, 5 replicate experiments). After approximately 20 hours, a chemotactic band appeared in each of these tubes and migrated along the length of the tube. These bands were left to migrate for a week. If the migrating band reached the end of the band before the end of the week, the band was extracted using a sterile 25-gauge syringe and transferred by inoculating a fresh, 1m-long tube This transfer process was necessary to allow migrations over distances longer than 1 m due to the length of the biosafety cabinet. The residual community tube experiment was performed using a separate 1m-long tube experiment inoculated with a separate sewage sample, The chemotactic band in this tube was left to run to the end of the tube and samples were taken from the band at 2 cm intervals along the length of the tube to capture the microbial community lost during migration.
Samples extracted from the tube experiments were divided into two cryovials, one for sequencing and one for flow cytometry. For flow cytometry, 50 µL of the sample was added to a cryovial with 995 µL of Phosphate Buffered Saline (PBS) and 5 µL of glutaraldehyde (0.125% vol/vol concentration) to fix the sample. All remaining extract was reserved for sequencing. Both flow cytometry and sequencing samples were snap frozen in liquid nitrogen and stored at -80°C.
The volumes of samples collected from the experiments and used for extraction varied due to needing to maintain the integrity of the band in the community formation experiments. In community formation experiments, the volume of sample used in extraction and long-term migration experiments was 50 µL and 500 µL respectively. For long-term migration samples, the volume of samples extracted ranged from 500µL to 1mL depending on the width of the band. Samples from all experiments were extracted using the Qiagen DNeasy PowerSoil Pro Kit (Qiagen, Germany), as per manufacturer's instructions.
Speed Measurements
The speed of the migrating bands was quantified using speed measurements taken at regular intervals for the duration of the long-term migration experiment. These measurements were obtained by recording how far the band in each replicate tube experiment had travelled from the inoculation site at four-hour intervals using a tape measure. If the bands were not visible, a handheld light was used to improve visibility. Using these measurements, the average migration speed was calculated across replicates. The calculations were only performed over the first 1 m of migration due to the transfer process causing bands to disappear until the chemotactic band re-formed in a fresh tube.
Flow cytometry sample preparation
The flow cytometry samples were prepared as described by Paterson et al (2019). Glutaraldehyde fixed flow cytometry samples from the community formation experiment were prepared undiluted with 250 µL of sample and 6.3 µL of SYBR Green nucleic acid stain (1:20,000 final dilution; MolecularProbes). Samples from the residual community and long-term migration experiment were prepared as a 1:10 dilution with 450 µL of filtered 0.02 µm Tris EDTA (TE) buffer (10 mM Tris, 1 mM EDTA; pH 8), 50 µL of the sample, and 12.5 µL of SYBR Green. Samples were incubated at 80ºC in the dark for 10 minutes. One-micrometre fluorescent calibration beads (Molecular Probes) were added to all samples as a size and concentration standard 45. The prepared samples were examined using a Beckman Coulter CytoFLEX S flow cytometer to determine the number of populations and the abundance of microbial communities. Three negative control tubes, or noise samples, were prepared with 500 µL of filtered TE buffer and 12.5 µL of SYBR Green. The flow cytometry results were analysed using FlowJo v10.6.2 Software (BD Life Sciences). The noise events were subtracted from the final abundance counts 46. In this study, VLPs, as identified via flow cytometry, refer to events recorded with a small size and low nucleic acid content, typical characteristics of viruses.
DNA extraction
DNA was extracted from each tube sample using the DNeasy Powersoil Pro DNA Extraction Kit (QIAGEN, Germany). To concentrate the samples prior to extraction, samples were centrifuged for 15 minutes at 3056 x g to pellet the bacteria. The supernatant was then discarded, and the pellet was resuspended in 200 µL of PBS and then placed into the PowerBead Pro tube. All subsequent extraction steps were then performed as per the manufacturer’s instructions. The extracted DNA was quantified on a Quibit™ 2.0 fluorometer using the Qubit dsDNA HS Assay Kit (Invitrogen) following the manufacturer's protocol. Community formation samples with the highest DNA content were sent to the Australian Genome Research Facility (AGRF) for 16S sequencing, using the V3-V4 hypervariable region on the Illumina MiSeq platform.
DNA Sequencing
16S Sequencing: 16S Sequencing: QIIME2 software was used to process the 16S sequencing data as per their recommended workflow47. To remove background noise, remove primer sequences, filter chimeric sequences, and join paired reads, DADA 2 software was used 47. An alpha-rarefaction curve was used to determine the sequencing depth. All samples were rarefied to 29211 features per sample (Supplementary Fig. 9). Sequences were then classified using the Greengenes 13_8 database.
Metagenomic Sequencing
Long-term migration and residual community samples were sequenced using whole genome shotgun sequencing on separate sequencing runs. For both, Nextera XT libraries were prepared using the Nextera XT DNA library prep kit (Illlumina) with twelve rounds of amplification. The MGI DNBSEQ-G400 sequencer (DNBSEQ™ technology) was used to obtain 150 bp paired-end reads for each sample. Sequence datasets ranged from 6,297,882 bp (944,682,300 reads) to 7,914,048,000 bp (52,760,320 reads) with an average quality score (PHRED) of 32.53 and an average per quality base score above 30 for each base (Supplementary Fig. 10–11). Sequences can be accessed using the BioProject number PRJNA788639. Sample accession numbers range from SRR17326388-SRR17326409 and SRR29503425-SRR29503452.
Bioinformatic analysis
Whole genome sequencing reads were trimmed using fastp (version 0.23.4)48 to remove adapters with parameters specifying a required length of 100 base pairs and a limit of one ‘N’ per sequence (parameters n 1 -l 100). Further quality control was performed using PRINSEQ++ (PReprocessing and INformation of SEQuence data)(version 1.2.4)49 to remove sequencing tags, duplicates and reads containing N’s (parameters min_len 60 -min_qual_mean 25 -ns_max_n 1 -derep 1 -out_format 0 -trim_tail_left 5 -trim_tail_right 5 -ns_max_n 5 -trim_qual_type min -trim_qual_left 30 -trim_qual_right 30 -trim_qual_window 10). Taxonomic annotations were assigned by mapping reads to the standard kraken2 database (January 2024 update) using Kraken2 (version 2.1.3)50. The species-level taxonomy abundances and genus-level taxonomy abundances were estimated from the Kraken classification output using Bracken51, which is recommended to perform a Bayesian estimation of the taxonomy abundance after the use of kraken2. To stress test our results, confidence filter levels were used for running kraken between confidence levels of 0.1 and 0.5 at 0.05 increments and the relative abundance and richness at the species level were compared. An ‘elbow’ in species richness was identified at a confidence of 0.1 (Supplementary Fig. 12–15) and this cutoff was used for all further analysis. Rarefaction curves were constructed for each sample by randomly sampling the cleaned fastq read files with replacement using seqtk (https://github.com/lh3/seqtk) to extract between 10–90% of the sequence data at 10% increments. For the residual community, the number of ‘rare’ taxa was estimated in each sample using rank abundance curves. A break was observed in the rank abundance curves at a relative abundance of 10− 2.5. Species with a relative abundance below 10− 2.5 were considered rare and species with a relative abundance above this cutoff were considered dominant. Functional annotations were assigned by mapping reeds to the SEED functional classification system using SUPERFOCUS (version 1.6)52. Viral profiles were generated using the Hecatomb pipeline53 which uses the UniProt and UniClust amino acid databases to annotate each read to the closest viral taxon using mmseqs254. These read annotations were filtered to remove sequences lacking amino-acid alignment support to known viral proteins and applying an e value cutoff of 1×10− 20 to remove low quality alignments. The remaining hits were renormalised to the relative abundance of each viral taxa per sample. Cleaned, paired reads were assembled using MEGATHIT55 (version 1.2.9). Samples from residual community experiments were cross-assembled and samples arising from the same sewage replicate were cross-assembled. The resulting contigs were filtered to consider contigs longer than 2000bp and the coverage of each of the contigs across samples was computed using Koverage56. The resulting contigs were binned using VAMB (version 4.1.4)57 to generate metagenome-assembled genomes. The completeness and contamination of each genome was evaluated with Checkm258. Genomes were taxonomically annotated using GTDB-Tk (version 2.1.1)59 and annotated with SEED subsystems60 using BV-BRC61. Pyani was used to compute similarity between MAGs as the average nucleotide identity (%)62. To compute the coverage of each MAG in each sample, the coverage of each contig in each sample was calculated using Koverage56. The mean coverage and length of each contig belonging to each MAG was then used to compute the MAG coverage using Eq. 1.
Statistical analysis
Prokaryotic, viral and functional profiles were analysed using R. For abundance and coverage data, stacked area plots were generated using ggplot 63 and heatmaps were generated using the R pheatmap package64. Principal coordinates analysis based on square-root transformed, Bray-Curtis dissimilarity matrices were performed using the R ecodist and vegan packages and visualised using ggplot263,65,66. PERMANOVA tests were run with the factors of time, tube and sewage sample with 9999 permutations and unrestricted permutations of raw data. Differential analysis was performed in R with DESeq2 (version 1.38.3)67. Volcano plot was generated using EnhancedVolcano (https://github.com/kevinblighe/EnhancedVolcano). To identify significant differences in the abundance of bacterial particles and VLPs from flow cytometric data, two-way ANOVAs with Tukey’s post-hoc test were performed.
Data Availability
Data presented in this study can be accessed from the NCBI Sequence Read Archive, BioProject number PRJNA788639. Sample accession numbers range from SRR17326388-SRR17326409 and SRR29503425-SRR29503452. Intermediate files and output are available at https://figshare.com/s/953b050065ca18fae420.
Code Availability
Scripts used to analyse metagenomic sequencing data can be accessed at https://github.com/linsalrob/atavide_lite. Code used to generate figures and analyses from this project is available at https://github.com/susiegriggo/Comigration_TubeExperiment_Analysis.
Acknowledgements
This work was supported by Australian Research Council grants (www.arc.gov.au) DP150103018 (to J.G.M.) and DP220102915 (to R.A.E and J.G.M.). S.R.G was supported by a Playford Trust Honours and PhD scholarship. We acknowledge the Flow Cytometric Unit of the Flinders Medical Centre, South Australia and Dr Giles Best for the use of flow cytometry facilities. We acknowledge SA Water, the Glenelg Wastewater treatment facility and Ian Mackenzie, for providing wastewater samples and the Flinders Accelerator for Microbiome Exploration for funding sequencing. We acknowledge Adriana Gaeguta-Bierbaum for the use of PC2 Laboratory facilities at Flinders University and Bhavya Papudeshi, Dr Michael Doane and Dr Vijini Mallawarachchi from the Flinders Accelerator for Microbiome Exploration for assistance with bioinformatic analysis.
Contributions
J.G.M and S.K.G conceived the experiments. S.R.G, A.L.K.H, J.A.P.C.J and J.S.P performed flow cytometric analysis. S.R.G, A.L.K.H, S.K.G, A.N.A performed experiments. S.R.G, A.L.K.H, C.H and A.K.G performed DNA extractions. S.R.G, A.L.K.H and J.G.M wrote the paper. A.L.K.H and J.A.P.C.J analysed 16S data. S.R.G, M.J.R and R.A.E performed bioinformatic analysis. S.R.G performed statistical analysis.
J.G.M and P.G.S supervised the project. All authors read and approved the final manuscript.
Competing Interests
The authors declare no competing interests.
- Stubbendieck, R. M., Vargas-Bautista, C. & Straight, P. D. Bacterial Communities: Interactions to Scale. Front. Microbiol. 7, 1234 (2016).
- Mushegian, A. A., Arbore, R., Walser, J.-C. & Ebert, D. Environmental Sources of Bacteria and Genetic Variation in Behavior Influence Host-Associated Microbiota. Appl. Environ. Microbiol. 85, (2019).
- Hibbing, M. E., Fuqua, C., Parsek, M. R. & Peterson, S. B. Bacterial competition: surviving and thriving in the microbial jungle. Nat. Rev. Microbiol. 8, 15–25 (2010).
- Thanarajah, S. & Wang, H. Competition of motile and immotile bacterial strains in a petri dish. Math. Biosci. Eng. 10, 399–424 (2013).
- Ghoul, M. & Mitri, S. The Ecology and Evolution of Microbial Competition. Trends Microbiol. 24, 833–845 (2016).
- Gude, S. et al. Bacterial coexistence driven by motility and spatial competition. Nature 578, 588–592 (2020).
- Scheuerl, T. et al. Bacterial adaptation is constrained in complex communities. Nat. Commun. 11, 754 (2020).
- Ottemann, K. M. & Miller, J. F. Roles for motility in bacterial-host interactions. Mol. Microbiol. 24, 1109–1117 (1997).
- Xie, L. & Wu, X.-L. Bacterial motility patterns reveal importance of exploitation over exploration in marine microhabitats. Part I: theory. Biophys. J. 107, 1712–1720 (2014).
- Adler, J. Effect of amino acids and oxygen on chemotaxis in Escherichia coli. J. Bacteriol. 92, 121–129 (1966).
- Adler, J. & Templeton, B. The effect of environmental conditions on the motility of Escherichia coli. J. Gen. Microbiol. 46, 175–184 (1967).
- Adler, J. & Dahl, M. M. A method for measuring the motility of bacteria and for comparing random and non-random motility. J. Gen. Microbiol. 46, 161–173 (1967).
- Bakker, R. G., Li, C., Miller, M. R., Cunningham, C. & Charon, N. W. Identification of specific chemoattractants and genetic complementation of a Borrelia burgdorferi chemotaxis mutant: flow cytometry-based capillary tube chemotaxis assay. Appl. Environ. Microbiol. 73, 1180–1188 (2007).
- Bainer, R., Park, H. & Cluzel, P. A high-throughput capillary assay for bacterial chemotaxis. J. Microbiol. Methods 55, 315–319 (2003).
- Tena-Garitaonaindia, M., Llamas, I., Toral, L. & Sampedro, I. Chemotaxis of halophilic bacterium Halomonas anticariensis FP35 towards the environmental pollutants phenol and naphthalene. Sci. Total Environ. 669, 631–636 (2019).
- Cremer, J. et al. Chemotaxis as a navigation strategy to boost range expansion. Nature 575, 658–663 (2019).
- Dann, L. M., Clanahan, M., Paterson, J. S. & Mitchell, J. G. Distinct niche partitioning of marine and freshwater microbes during colonisation. FEMS Microbiol. Ecol. 95, (2019).
- Raina, J.-B. et al. Chemotaxis shapes the microscale organization of the ocean’s microbiome. Nature 605, 132–138 (2022).
- Ahmed, T., Shimizu, T. S. & Stocker, R. Bacterial chemotaxis in linear and nonlinear steady microfluidic gradients. Nano Lett. 10, 3379–3385 (2010).
- Parks, D. H. et al. A complete domain-to-species taxonomy for Bacteria and Archaea. Nat. Biotechnol. 38, 1079–1086 (2020).
- Adler, J. Chemotaxis in bacteria. Science 153, 708–716 (1966).
- Smriga, S., Fernandez, V. I., Mitchell, J. G. & Stocker, R. Chemotaxis toward phytoplankton drives organic matter partitioning among marine bacteria. Proc. Natl. Acad. Sci. U. S. A. 113, 1576–1581 (2016).
- Zhou, B., Szymanski, C. M. & Baylink, A. Bacterial chemotaxis in human diseases. Trends Microbiol. 31, 453–467 (2023).
- Seymour, J. R., Brumley, D. R., Stocker, R. & Raina, J.-B. Swimming towards each other: the role of chemotaxis in bacterial interactions. Trends Microbiol. (2024) doi:10.1016/j.tim.2023.12.008.
- Kiørboe, T., Grossart, H.-P., Ploug, H. & Tang, K. Mechanisms and rates of bacterial colonization of sinking aggregates. Appl. Environ. Microbiol. 68, 3996–4006 (2002).
- Stocker, R. & Seymour, J. R. Ecology and physics of bacterial chemotaxis in the ocean. Microbiol. Mol. Biol. Rev. 76, 792–812 (2012).
- Pelve, E. A., Fontanez, K. M. & DeLong, E. F. Bacterial Succession on Sinking Particles in the Ocean’s Interior. Front. Microbiol. 8, 2269 (2017).
- Bai, Y. et al. Spatial modulation of individual behaviors enables an ordered structure of diverse phenotypes during bacterial group migration. Elife 10, (2021).
- Fu, X. et al. Spatial self-organization resolves conflicts between individuality and collective migration. Nat. Commun. 9, 2177 (2018).
- Zeitz, M., Wolff, K. & Stark, H. Active Brownian particles moving in a random Lorentz gas. Eur. Phys. J. E Soft Matter 40, 23 (2017).
- Lynch, M. D. J. & Neufeld, J. D. Ecology and exploration of the rare biosphere. Nat. Rev. Microbiol. 13, 217–229 (2015).
- Clemmer, K. M., Bonomo, R. A. & Rather, P. N. Genetic analysis of surface motility in Acinetobacter baumannii. Microbiology 157, 2534–2544 (2011).
- Liu, Z. et al. The synergistic effects of plants and nitrogen on microbial hitchhiking. bioRxiv 2022.05.09.491057(2022) doi:10.1101/2022.05.09.491057.
- Finkelshtein, A., Roth, D., Ben-Jacob, E. & Ingham, C. J. Bacterial swarms recruit cargo bacteria to pave the way in toxic environments. MBio, 6 (3), e00074–15 (2015).
- Chang, C.-Y., Bajić, D., Vila, J. C. C., Estrela, S. & Sanchez, A. Emergent coexistence in multispecies microbial communities. Science 381, 343–348 (2023).
- Louca, S. et al. Function and functional redundancy in microbial systems. Nat Ecol Evol 2, 936–943 (2018).
- Madi, N., Chen, D., Wolff, R., Shapiro, B. J. & Garud, N. R. Community diversity is associated with intra-species genetic diversity and gene loss in the human gut microbiome. Elife 12, (2023).
- Gumerov, V. M., Andrianova, E. P. & Zhulin, I. B. Diversity of bacterial chemosensory systems. Curr. Opin. Microbiol. 61, 42–50 (2021).
- Keegstra, J. M., Carrara, F. & Stocker, R. The ecological roles of bacterial chemotaxis. Nat. Rev. Microbiol. 20, 491–504 (2022).
- Parkinson, J. S. cheA, cheB, and cheC genes of Escherichia coli and their role in chemotaxis. J. Bacteriol. 126, 758–770 (1976).
- Brochet, S. et al. Niche partitioning facilitates coexistence of closely related honey bee gut bacteria. Elife 10, (2021).
- Lenski, R. E. & Velicer, G. J. Games microbes play. Selection 1, 89–96 (2001).
- Luo, N. et al. The collapse of cooperation during range expansion of Pseudomonas aeruginosa. Nat. Microbiol. 9, 1220–1230 (2024).
- Leinweber, A., Fredrik Inglis, R. & Kümmerli, R. Cheating fosters species co-existence in well-mixed bacterial communities. ISME J. 11, 1179–1188 (2017).
- Gasol, J. M., Zweifel, U. L., Peters, F., Fuhrman, J. A. & Hagström, A. Significance of size and nucleic acid content heterogeneity as measured by flow cytometry in natural planktonic bacteria. Appl. Environ. Microbiol. 65, 4475–4483 (1999).
- Paterson, J. S. et al. A hydrocarbon-contaminated aquifer reveals a Piggyback-the-Persistent viral strategy. FEMS Microbiol. Ecol. 95, (2019).
- Estaki, M. et al. QIIME 2 Enables Comprehensive End-to-End Analysis of Diverse Microbiome Data and Comparative Studies with Publicly Available Data. Curr. Protoc. Bioinformatics 70, e100 (2020).
- Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890 (2018).
- Cantu, V. A., Sadural, J. & Edwards, R. PRINSEQ++, a multi-threaded tool for fast and efficient quality control and preprocessing of sequencing datasets. PeerJ (2019) doi:10.7287/peerj.preprints.27553v1.
- Wood, D. E., Lu, J. & Langmead, B. Improved metagenomic analysis with Kraken 2. Genome Biol. 20, 257 (2019).
- Lu, J., Breitwieser, F. P., Thielen, P. & Salzberg, S. L. Bracken: estimating species abundance in metagenomics data. PeerJ Comput. Sci. 3, e104 (2017).
- Silva, G. G. Z., Lopes, F. A. C. & Edwards, R. A. An Agile Functional Analysis of Metagenomic Data Using SUPER-FOCUS. Methods Mol. Biol. 1611, 35–44 (2017).
- Roach, M. et al. Hecatomb: An End-to-End Research Platform for Viral Metagenomics. bioRxiv 2022.05.15.492003(2022) doi:10.1101/2022.05.15.492003.
- Steinegger, M. & Söding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat. Biotechnol. 35, 1026–1028 (2017).
- Li, D., Liu, C.-M., Luo, R., Sadakane, K. & Lam, T.-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics 31, 1674–1676 (2015).
- Roach, M. J. et al. Koverage: Read-coverage analysis for massive (meta)genomics datasets. J. Open Source Softw. 9, 6235 (2024).
- Nissen, J. N. et al. Improved metagenome binning and assembly using deep variational autoencoders. Nat. Biotechnol. 39, 555–560 (2021).
- Chklovski, A., Parks, D. H., Woodcroft, B. J. & Tyson, G. W. CheckM2: a rapid, scalable and accurate tool for assessing microbial genome quality using machine learning. Nat. Methods 20, 1203–1212 (2023).
- Chaumeil, P.-A., Mussig, A. J., Hugenholtz, P. & Parks, D. H. GTDB-Tk v2: memory friendly classification with the genome taxonomy database. Bioinformatics 38, 5315–5316 (2022).
- Overbeek, R. et al. The SEED and the Rapid Annotation of microbial genomes using Subsystems Technology (RAST). Nucleic Acids Res. 42, D206–14 (2014).
- Olson, R. D. et al. Introducing the Bacterial and Viral Bioinformatics Resource Center (BV-BRC): a resource combining PATRIC, IRD and ViPR. Nucleic Acids Res. 51, D678–D689 (2023).
- Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T. & Aluru, S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9, 5114 (2018).
- Wickham, H. Ggplot2. Wiley Interdiscip. Rev. Comput. Stat. 3, 180–185 (2011).
- Kolde, R. Pheatmap: pretty heatmaps. R package version.
- Goslee, S. C. & Urban, D. L. The ecodist Package for Dissimilarity-based Analysis of Ecological Data. J. Stat. Softw. 22, 1–19 (2007).
- Dixon, P. VEGAN, a package of R functions for community ecology. J. Veg. Sci. 14, 927–930 (2003).
- Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).
- Flinders University. Deep Thought (HPC). 2021.
There is NO Competing Interest.