Spatial-temporal dynamics of a microbial cooperative behavior resistant to cheating

Timeseries imaging tracks gene expression in spatial systems

Recent studies have shown it possible to identify the members of microbial consortia as well as their gene expression within spatially-structured systems^30,33,34. However, these methods capture data cross-sectionally and are unable to provide temporal insight into gene expression patterning as it emerges in these cell populations. To bridge this gap, we built a fluorescent imager inside an incubator (Supplementary Fig. 1). Our framework characterizes cellular growth and gene expression in spatially-structured environments with previously unattainable time-resolution and throughput. Fluorescently labeled cells are illuminated using LEDs connected to a custom-built control system (see methods). The images are background corrected and analyzed, tracking colony growth and gene expression information (Supplementary Figs. 2, 3) straight from the spatially-structured system.

In our experiments, we utilized a dual-labeled P. aeruginosa PA14 strain harboring PBad-DsRed(EC2)³⁵ driven by L-arabinose in the plate media, which cannot be metabolized by the cells³⁶, and PrhlAB-GFP^28,37. When grown in spatial structure, the constitutive expression of DsRed provided a measure of the local density of bacteria (Supplementary Fig. 4). In all our experiments, the dynamical expression of GFP, validated by RT-qPCR (Supplementary Fig. 5) (see methods), reported on the expression of rhlAB.

Using these data, we were able to characterize how the surroundings experienced by these microbes influence the dynamics of their cooperative behavior directly in a spatially-structured setting.

Rhamnolipid production differs in liquid and spatial environments

Rhamnolipids are necessary for cooperative swarming behavior in P. aeruginosa and for other traits related to virulence²⁶. Rhamnolipids can be produced in liquid culture^10,20,28,38, thus rhamnolipid production is often studied in detail there. Despite recent work indicating that gene expression related to quorum signaling systems in P. aeruginosa may differ in spatial structure²⁹, no studies assess how downstream genes, such as rhlAB, may be affected in spatially-structured colonies. Given the relevance of these diffusible inputs to the rhlAB system, we hypothesized that there could be differences between gene expression patterns in liquid and spatial environments.

We compared P. aeruginosa biomass growth and gene expression in the liquid and spatial environments (Fig. 1a). Liquid culture data was collected following prior methods²⁸. To interrogate the spatial system, we used the protocol from the classic Colony Forming Unit (CFU) assay. Cells were seeded with extreme dilution and we observed the behavior of the resultant colonies (cCFUs) across time and within the random configurations generated.

We observed differences in growth between cells grown in liquid culture (Fig. 1b) and spatial structure (Fig. 1c) with the same media composition. The growth pattern observed in liquid culture recapitulates previously reported data^22,28. In comparing WT growth (dark blue data in Fig. 1b, c) between environments, we observed that both achieve a period of exponential growth, followed by a period of slowed growth. This sub-exponential growth is prolonged and no period of biomass decay is observed in the spatially-structured environment during our observation window.

Quorum signal perturbation has long been an experimental tool to determine if a phenotype is responsive to social signaling^9,10. rhlAB gene expression in particular is known to be downstream of both the las and rhl quorum signal systems^39,40. However, it has previously been shown that liquid culture perturbation with additional C4-HSL and 3-oxo-C12-HSL, the rhl and las quorum signal system auto-inducers respectively, do not illicit significant change in growth or P_rhlAB dynamics in this strain of P. aeruginosa²². We replicated this liquid culture result (Fig. 1b, purple data). In the spatially-structured system, we performed this perturbation by including both quorum signal molecules in the plate media in the same concentration by volume as previously published²². This analysis was done using biological replicates with <70 colonies (Fig. 1c [Inset]). In comparing between colonies grown with or without quorum signals in the plate media, we observed that colonies perturbed by quorum signals may achieve a smaller final size after 48 h of growth (Supplementary Fig. 6a). We did not observe a difference between the specific growth rate of the colonies during the time interval when they came above detection (Supplementary Fig. 6b). However, we did observe that colonies given quorum signal perturbation show later colony detection (Supplementary Fig. 6c).

We analyzed the promoter activity^28,41 of P. aeruginosa grown both in liquid culture (Fig. 1d) and spatial structure (Fig. 1e). In liquid culture, we found rhlAB promoter activity to be low during periods of high specific growth rate as seen previously^22,28. Promoter activity increased as the specific growth rate decreased and below a threshold growth rate promoter activity dropped as expected during prolonged stationary phase²⁸. Unexpectedly, in spatial structure, we observed a strong positive correlation between specific growth rate and promoter activity (R² = 0.96) (Fig. 1e).

Previous work done in liquid culture captured no significant change in rhamnolipid production in WT bacteria grown with quorum signals added to the media²² and our data agree (Fig. 1d, purple data). Conversely, we found that in spatial structure, WT colonies expressed even higher levels of rhlAB during periods of high growth rate when quorum signals were added to the same plate media recipe (two-sided rank sum test p-value < 1e−4, see methods). This presents in our data as a steeper positive slope in the association between specific growth rate and promoter activity (R² = 0.98). We conclude that not only are there phenotypic differences between rhlAB gene expression in the liquid and spatial systems, but that there is a phenotypic difference under quorum signal perturbation that is specific to the spatially-structured system.

Cellular response to diffusive quorum signals is distance-dependent

Next, we sought to understand whether the differences we observed in liquid culture and spatially-structured gene expression could be driven by diffusible quorum signals. Previous work in liquid culture revealed that rhlAB expression can integrate nutrient and quorum signal information from at least three diffusible small molecules: a growth-limiting nutrient^22,28 and the hierarchical quorum sensing pathway involving the autoinducer molecules 3-oxo-C12-HSL and C4-HSL^{40,42,43,44,45,46}. When bound to their cognate receptors, LasR and RhlR respectively, these complexes may act as transcription factors, instigating systemic gene expression change^47,48 (Fig. 2a). All three diffusive inputs are of similar molecular size and thus may act on similar length and time-scales. The ratio of the diffusion coefficients and decay rates for these molecules in bacterial growth media³⁸ indicate that the quorum signals could reach biomass that is multiple millimeters away, though whether their physiologic concentrations could influence biomass at that distance was unknown.

P. aeruginosa has been shown capable of micrometer length-scale communication in constrained microfluidic experiments³¹. However, it is difficult to extrapolate these results to the full spatial-temporal system. With the knowledge that P. aeruginosa cCFUs respond to systemic quorum signal perturbation, we asked: over what macro- spatial-temporal scales are these cells capable of responding to quorum signal perturbation? To address this, we utilized a signal-mute mutant, PA14 ΔlasIΔrhlI, that cannot produce the 3-oxo-C12-HSL and C4-HSL molecules, but is able to respond when these signals are exogenously provided (Fig. 2a). In these experiments, this strain was double-labeled in the same way as the WT PA14. We focused on response to C4-HSL. We added the upstream quorum signal (3-oxo-C12-HSL) directly to the plate media and loaded 4 µL of 5 µM C4-HSL on a filter paper in the center of the plate (Fig. 2b).

We tracked the growth and rhlAB expression in colonies seeded around the filter paper (Fig. 2c). In this experimental configuration, observing P_rhlAB activity in a colony indicated that it had encountered both quorum signals at concentrations high enough to trigger a rhlAB response. We found that the response in the signal-mute mutants varied with the distance between the colony and the center of the filter paper (Fig. 2c, d).

We found that maximum colony promoter activity was inversely proportional to the colony’s distance to the filter paper (Fig. 2d) (R² = 0.41). The highest maximal promoter activities we observed occurred in colonies <2.5 cm from the quorum signal source. In colonies between 2.5 and 4 cm away, the maximal promoter activity scaled with the distance to the quorum signal source more strongly (R² = 0.63). We investigated the presence of rotational biases in our data by comparing the distributions of maximal promoter activity within 45° increments with the two-sided rank sum test (Supplementary Fig. 7). The strongest bias our investigation revealed was between 90 and 135°. However, in our dataset, this region had fewer samples (Supplementary Fig. 7b) and all colonies were within 2.75 cm from the quorum signal source (Supplementary Fig. 7c), a region of high variability across all our data.

These experiments carried out with the signal-mute mutant confirm that P. aeruginosa can respond to diffusible quorum-signal perturbation on a centimeter length-scale. They also illustrate that P. aeruginosa is capable of a concentration-dependent dose-response to diffusible quorum signals. As our experimental protocol uses physiologically relevant quorum signal concentrations and time-scales²², these data indicated that these cells may be capable of configuration-dependent behavior.

Inference of the cellular spatial environment by model selection

Given our results suggesting that spatially distinct multicellular aggregates may be capable of communication over macro-scale distances, we next looked to test whether similar centimeter length-scale interactions could be observed in the WT. The diffusion coefficients for the quorum signals we have investigated here, C4-HSL (MW 171.9 g/mol) and 3-oxo-C-12-HSL (297.37 g/mol), are on the order of ~7 × 10⁻⁶ cm²/sec, slightly slower for the larger 3-oxo-C-12-HSL⁴⁹. This means that over the course of 24–48 h, these signals are capable of traveling 1–2 cm away from their source. Based on this, we predicted that P. aeruginosa colonies may be able to detect and respond to each other within small macro-scale distances within the viewing timeframe of our experiments. However, the integration of these signals in liquid is known to be non-trivial^20,22,28. This infrastructure poses a system that may be highly sensitive to fluctuations in the diffusive environment. Therefore, a data-driven and unsupervised approach was required to provide unbiased insight into the spatial-temporal scale over which spatially segregated bacterial communities may influence one another.

To test these predictions, we performed CFU experiments with between 60 and 150 WT colonies seeded throughout the plate. Here, the variability in colony location inherent to extreme dilution seeding provided an experimental advantage. These plates explore a variety of colony configurations. Each colony had a unique location with respect to every other colony on the plate. Each biological replicate had colonies spread over a similar total plate area. As a result, this approach generated a large amount of variation in growth pattern and gene expression. This self-generated variation allowed us to leverage unbiased and data-driven methodology to uncover, 1—whether colony-colony interactions could explain the variation we observed in colony growth and gene expression and if so, 2—how did these interactions develop with time and over what spatial scale? We define colony-colony interactions here by statistical association, examining whether focal colony promoter activity can be explained by the growth patterns of non-focal colonies.

To answer these questions, we applied a spatial kernel approach—a method long used in biophysical systems to model the patterns generated by the interactions between spatially segregated or dispersing individuals^50,51,52,53. As an example (Fig. 3a), arid landscapes often show a patchy pattern of vegetation due to the presence and removal of limiting resources such as water and soil nutrients⁵³. A single shrub in an arid environment can preserve some water and nutrients in the soil surrounding their roots. Other shrubs in close proximity impact the focal individual positively. By making the local root system dense, the subsequent preservation of nutrients facilitates the survival of all nearby plants. However, shrubs farther away have a competitive (negative) impact, drawing nutrients away from the localized hub^54,55. This methodology can capture behaviors occurring simultaneously across multiple length-scales while allowing the flexibility to encapsulate a wide range of shapes⁵¹ besides the Gaussian shape typical of diffusional processes. We apply this same idea to describe the positive and negative effects that the spatial configuration of biomass may have on the cooperative behavior of a focal colony. However, where previous work has used the data to fit the parameter values for spatial kernels of a specified shape^50,51,56, we use a data-driven approach to fit the shape of the spatial kernel itself. This is the value of the heterogeneity generated by our application of the CFU assay protocol. The variety in spatial configuration allows us to sample a wide range of colony arrangements, giving us the statistical opportunity to infer the spatial-temporal length-scales at play directly from the data.

The spatial kernel can be modeled as a collection of concentric annuli of fixed radius emanating from each colony (Fig. 3b). The promoter activity of a focal colony (Fig. 3b, black arrow) was investigated with respect to the surrounding colony configuration. We discretized the kernel with a 1 mm distance between the outer and inner radii of each annulus. We fit the following linear model to the colony promoter activity:

where ({P}_{{C}_{F}}) is the promoter activity of a focal colony, ({mu }_{C_{F}}) is the specific growth rate of the focal colony and ({B}_{{C}_{F}}) is the amount of biomass in the focal colony at time (t). ({B}_{{C}_{i,j}}) is the amount of biomass in colony (j) in annulus (i) where there are (A) total annuli and ({n}_{i}) colonies in the annulus of interest. In this formulation, focal colony promoter activity may scale with the colony’s growth rate as well as the number of cells present in the focal colony at time (t). The total biomass in each annulus is used as a series of features in the model to explain the variation in promoter activity between colonies of similar size and growth rate. All annulus features were normalized for annulus area before fitting to ensure that all annuli contributed equally to the fit. We assumed that all cells that founded a CFU landed on the plate at the same time. Diffusion processes that impacted a colony were then assumed to occur across the same time-scale for all colonies simultaneously. It stands to reason that a focal colony may experience a time-delay in the impact of distant colonies; we found our implementation to be a decent approximation.

We independently fit models for data taken every 30 min between 20 and 48 h. At each timepoint, a series of models were fitted where each new model included one annulus farther from the focal colony than the previous model. We compared models using the Aikake Information Criterion (AIC) to assess the trade-off between model simplicity (fewer features) and the quality of fit. In accordance with standard practice, we chose the best model as the one with the lowest AIC. We interpret the distance of the largest annulus included in that optimal model as the longest length-scale of colony-colony interaction at that timepoint (Fig. 3c).

Our AIC-selected spatial kernel approach fit our data well (Supplementary Fig. 8) and revealed the surprisingly clear result that the colony-colony interaction length-scale lengthens with time (R² = 0.54) (Fig. 3d). Further, focal colony feature coefficients showed internal consistency across our independently fit models. (Fig. 3e). Model selection identified interactions between colonies up to 1–1.5 cm apart early in the timeseries (22–26 h) and up to 1.8–2.3 cm apart later in the timeseries (40–44 h). This data-driven approach leads us to conclude that WT PA14 are capable of centimeter length-scale colony-colony interactions within a 48 h window.

Finally, we reviewed the spatial kernels predicted by our model to see how biomass localized in each annulus was predicted to influence focal colony promoter activity at various timepoints. Earlier in the timeseries, at 24 h, we observed that biomass more than 0.5 cm away from the focal colony negatively impacted promoter activity. However, this relationship shifted with time. By 44 h, all colonies within 1.5 cm of the focal colony had a negative impact, while colonies more than 1.5 cm away may have positively impacted focal colony promoter activity (Fig. 3f). We do not claim that these interactions are due only to the diffusion of C4-HSL and 3-oxo-C12-HSL, though these results do match the length-scales of interaction predicted by our investigation of quorum signal gradients (Fig. 2). All together, these results characterize a general length-scale of gene expression association between colonies on the order of 1–2 cm that lengthens and changes shape with time.

Swarm tendrils achieve exponential growth despite constant velocity

Uncovering associations between gene expression patterns in spatially-distinct biomass aggregates led us to ask whether these findings could extend to motile P. aeruginosa swarms. This swarming behavior has long been of general interest due to its example as a cooperative behavior that is not invadable by non-cooperating strains in competitive assays²². Given our success in interrogating gene expression directly in spatially-structured systems, we looked next to extend our investigations to the motile swarming system.

Specifically, we wanted to know whether the WT motile swarms would show growth and rhlAB promoter activity patterns more similar to classic well-mixed liquid culture or the new dynamics found in the immotile spatial system. To do this, we fluorescently imaged swarms (Supplementary Fig. 9) and isolated cross-sectional biomass (DsRed) and P_rhlAB (GFP) measurements along the length of three tendrils in each of four independent swarms. We first investigated growth in swarming tendrils (Fig. 4a, b[top]). In a striking departure from both our liquid and cCFU results, we found that despite attaining an average constant velocity of 3.56 mm/h (±0.65 mm/h), these tendrils were capable of achieving and sustaining periods of exponential growth (Fig. 4a, dashed line—linear growth trajectory). This phenomenon may be related to cell motility as seeding cells in a tendril configuration on motility-preventing agar showed growth dynamics similar to cCFUs (Supplementary Fig. 10).

Swarming tendril P_rhlAB activity matches cCFUs

We next looked to compare the relationship between promoter activity and growth rate within a swarming tendril. To do this, we calculated these metrics along our tendril cross-sections and examined the data for spatial localization (Fig. 4b). We separated our cross-sections into three segments: the swarm center, the swarm edge, and the mid-tendril region between them (Fig. 4c inset). The edge of the tendril was a region that showed biomass localization, located near the tip of each tendril (Fig. 4b, Supplementary Figs. 11b and 12). It was typically between 2.75 and 4.5 mm in length. This analysis again revealed a positive correlation between growth rate and promoter activity R² = 0.79 in the tendril tip, the front-most 0.86 mm of the tendril edge (Fig. 4c). This finding continues to be highly counter-intuitive given previously published work^22,28 as well as our own liquid culture data (Fig. 1b, d). However, these trends fit with our new spatially-driven expectations (Fig. 1c, e). We performed our analysis conservatively, only examining a pixel after biomass had been present in that location for three or more timepoints (15 min) to prevent artifacts. This approach did not account for biomass flux into or out of any pixel. We assumed that the mid-tendril region is seeded by cells left behind as the edge proceeds away from the swarm center. As there will be non-negligible flux out of the tendril edge, the growth rate calculated for the edge of a tendril may be an underestimate. By contrast, the mid-tendril and the swarm center exhibited a much narrower range of growth rate and corresponding promoter activity (Fig. 4d, e).

Quorum signal perturbation reveals swarm biomass redistribution

Finally, we wanted to know how perturbation with quorum signals impacted swarming behavior. In light of our previous results, we hypothesized that in swarms provided with exogenous quorum signals we would see higher promoter activity related to rhamnolipid production, resulting in faster spreading and earlier tendril formation in P. aeruginosa swarms.

We extracted data from three tendrils in each of three independent swarms grown with quorum signals in the plate media and examined the data for quorum signal-induced phenotypic differences. We found that, similarly to the non-perturbed swarms, these swarms were able to achieve a sustained period of exponential growth along the length of each tendril and reach the same total biomass at the end of our 24 h timeseries (Supplementary Fig. 11).

Surprisingly, we found that localization of biomass within these tendrils differed significantly from the original WT swarming tendrils (Supplementary Figs. 11–12). We found that under quorum signal perturbation, biomass increasingly localized to the swarm center and the tendril edge as the tendrils lengthened (Fig. 5a). To our knowledge, a socially-driven spatial segregation phenotype such as this has not been previously identified in P. aeruginosa, with the closest comparisons being the social regulation involved in facilitating fruiting body formation in Myxococcus xanthus⁵⁷ or Bacillus subtilis⁵⁸.

We found that, indeed, swarms provided with quorum signals form tendrils sooner than non-perturbed swarms, p-value < 1e−8 by Kolmogorov–Smirnov test (Fig. 5b), and moved with a faster average velocity (4.58 mm/h ± 1.00 mm/h) than swarms not given quorum signals (3.56 mm/h ± 0.65 mm/h). Cells in the tendril tips of these swarms achieved higher growth rates, p-value < 1e−10 (Fig. 5c) and higher promoter activities (Supplementary Fig. 13), scaling linearly with growth rate, R² = 0.89. We did not detect a change in the slope of the relationship between promoter activity and growth rate in P. aeruginosa swarms perturbed with quorum signals.

Quorum signal perturbation does not facilitate invasion by defectors

Lastly, we investigated swarming competition in the presence of quorum signals. Rhamnolipid production has long been posited as a possible competitive weak point in P. aeruginosa cooperative swarming as it represents a large resource investment. Once secreted, the rhamnolipids can be utilized by other cells in the vicinity²². However, competitions conducted between the WT and a rhamnolipid defector (ΔrhlA) have never shown the WT to definitively lose. In light of our quorum signal perturbation data in swarms (Supplementary Fig. 13), it was unclear whether the increased P_rhlAB activity could make these cells susceptible to invasion by this defector strain in a competitive setting.

Swarming competition experiments were performed as previously reported²², now with plate media containing quorum signals. The exact initial mix of WT to defector was calculated to account for variability due to mixing and dilution (left(frac{{{WT}}_{i}}{{{Total; Cells}}_{i}}right)). After the competition, the final ratio (left(frac{{{WT}}_{f}}{{{Total; Cells}}_{f}}right)) was calculated. We found that, in our hands, the WT continued to resist invasion by the defector strain despite quorum signal perturbation (Fig. 5d).

Source: Ecology - nature.com