article: Surface-attached bacterial communities called biofilms display a diversity of morphologies. Although structural and regulatory components required for biofilm formation are known, it is not understood how these essential constituents promote biofilm surface morphology. Here, using Vibrio cholerae as our model system, we combine mechanical measurements, theory and simulation, quantitative image analyses, surface energy characterizations, and mutagenesis to show that mechanical instabilities, including wrinkling and delamination, underlie the morphogenesis program of growing biofilms. We also identify interfacial energy as a key driving force for mechanomorphogenesis because it dictates the generation of new and the annihilation of existing interfaces. Finally, we discover feedback between mechanomorphogenesis and biofilm expansion, which shapes the overall biofilm contour. The morphogenesis principles that we discover in bacterial biofilms, which rely on mechanical instabilities and interfacial energies, should be generally applicable to morphogenesis processes in tissues in higher organisms.
  Many of the stunning morphologies that distinguish living entities do not arise exclusively from gene expression programs, but rather from overarching contributions from mechanical forces (Heisenberg and Bellaïche, 2013; Thompson, 1992; Yamada and Cukierman, 2007). Such morphomechanical processes include the formation of ripple-shaped leaves (Liang and Mahadevan, 2009), tendrils and flowers (Gerbode et al. , 2012; Liang and Mahadevan, 2011), as well as the dorsal closure and apical constriction-mediated epithelial folding processes that take place during Drosophila embryonic development (He et al. , 2014; Solon et al. , 2009). One key feature is common to many of these morphogenic transformations: two or more layers of biomaterials are attached to one another but each grows at a different rate (Wang and Zhao, 2015). Inevitably, such growth mismatches generate mechanical stresses, and corresponding shape instabilities, which depend on the mechanical and other material properties of the biological constituents, as well as their geometries. Some examples include villi formation during the development of the human gut and formation of gyri and sulci during cerebrum development (Shyer et al. , 2013; Budday et al. , 2015; Tallinen et al. , 2016). Though ancient in their evolutionary origin, bacterial cells can also display intricate developmental patterns, particularly when they exist in the community lifestyle known as biofilms (Hobley et al. , 2015; Humphries et al. , 2017; Persat et al. , 2015). Biofilms are surface-associated bacterial communities that are embedded in a polymer matrix (O' Toole et al. , 2000; Thongsomboon et al. , 2018) and are a predominant growth mode for bacteria in nature (Hall-Stoodley et al. , 2004; Humphries et al. , 2017). Biofilms can be beneficial, for example in waste-water treatment (Nerenberg, 2016), but they also cause significant problems in health and industry (Costerton et al. , 1999; Drescher et al. , 2013) because they are resistant to physical perturbations and to antibiotics (Kovach et al. , 2017; Meylan et al. , 2018). Biofilms on surfaces undergo morphogenic transformations, beginning as smooth colonies and, over time, developing complex morphological features (Beyhan and Yildiz, 2007). Genes specifying matrix components that enable polysaccharide production, cell-surface adhesion, and cell–cell adhesion are required for the morphological transition (Hobley et al. , 2015). However, the underlying mechanisms that dictate how these biofilm matrix components direct overall morphology are not well-understood. One model focuses on the differential spatial regulation of genes encoding matrix components as the key driver of biofilm morphogenesis (Okegbe et al. , 2014). Another model suggests that localized cell death serves as an outlet for mechanical stresses and thus determines biofilm morphology (Asally et al. , 2012). Most recently, theory has been put forward to suggest the possibility that global mechanical instabilities are involved in the development of biofilm morphology (Zhang et al. , 2016; Zhang et al. , 2017). Here, by combining quantitative imaging, biomaterial characterization, mutant analyses, and mechanical theory, we show that the mismatch between the growing biofilm layer and the non-growing substrate causes mechanical instabilities that enable the biofilm to transition from a flat to a wrinkled film, and subsequently to a partially detached film containing delaminated blisters. The sequential instabilities that the film undergoes, coupled with the generation and annihilation of interfaces, drive the evolution of biofilm topography. Our results demonstrate that bacterial biofilms provide a uniquely tractable system for the quantitative investigation of mechanomorphogenesis.
  Our central hypothesis is that biofilm morphogenesis is driven by mechanical instabilities that arise from the growth mismatch between an expanding biofilm and the non-growing substrate to which it adheres. To garner evidence for this idea, we grew biofilms on agar plates, which enabled us to control the mechanical properties of the substrate by changing the agar concentration (Nayar et al. , 2012). We employed a commonly used Vibrio cholerae strain that lacks motility and constitutively produces biofilms (Beyhan and Yildiz, 2007; Yan et al. , 2017). This strain (denoted WT in the present work) produces biofilms that have disordered cores decorated with radial features extending to the rims (Figure 1A). Indeed, biofilm surface morphology changes with increasing agar concentration: the spacing between the peripheral, radial features is reduced and their amplitudes become more homogeneous (Figure 1—figure supplement 1). Encouraged by the observations described above and inspired by models developed to describe mechanical instabilities in abiotic materials systems (Li et al. , 2012), here we propose a mechanomorphogenesis model for biofilms (Figure 1B). The biofilm originates as a flat film. Its volume increases over time due to cell proliferation and matrix production. If the biofilm were not attached to a substrate, it would grow into a stress-free state to cover a large area (Figure 1B, top, ‘virtual state’). However, the non-expanding agar substrate constrains biofilm expansion. Thus, biofilms are always subject to compressive stress (Figure 1B, middle right), which we hypothesize drives the surface morphology. Indeed, a biofilm growing at an air–liquid interface, not limited or compressed by a substrate, exhibits no surface features (Video 1). According to mechanical instability theories, surface-adhered films under compression have several pathways to release compressive stress (Wang and Zhao, 2015). For example, the film can buckle out of the growth plane and deform together with the substrate into a periodically wrinkled pattern (Figure 1B, bottom left). In this mode, the compressive stress is released by film bending and substrate deformation. Alternatively, the film can directly delaminate from the substrate to form ‘blisters’ (Figure 1B, bottom right) (Vella et al. , 2009), leaving the substrate essentially undeformed. An extra interfacial energy penalty is paid for delamination because new interfaces are generated, so direct delamination occurs in systems with film–substrate adhesion energies that are much smaller than their elastic deformation energies. Biofilms possess finite adhesion strength (~ 5 mJ/m2), which is the same order of magnitude as the deformation energy of the soft substrate (Yan et al. , 2018). Thus, we suggest that biofilms could first wrinkle, and subsequently delaminate as growth gradually builds up compressive stress (Figure 1—figure supplement 2). According to this mechanomorphogenesis model, we should be able to change the biofilm topography by changing the spatial distribution of the mechanical stress. To this end, we inoculated two V. cholerae biofilms onto the same agar plate and allowed them to collide. Indeed, a large localized blister formed at the collision front where mechanical stress is most concentrated (Figure 1C; Video 2). Our mechanomorphogenesis model provides an intuitive explanation for the commonly observed biofilm surface pattern of a disordered core surrounded by radial features at the edge (DePas et al. , 2013; Okegbe et al. , 2014; Wilking et al. , 2013). Soon after the initial expansion of the biofilm, growth occurs primarily at the edge of the biofilm because of nutrient limitation at the center of the biofilm (Liu et al. , 2015; Yan et al. , 2017; and Figure 1—figure supplement 2). At the biofilm center, cell death has been shown to drive pattern formation (Asally et al. , 2012). However, in the biofilm periphery, which is the region of focus of the current study, wrinkling and delamination occur with no preceding localized cell death (Figure 1—figure supplement 2). In this outer region, mechanical instabilities dominate the pattern formation and its wavelength. Directionality at the edge stems from the asymmetry between radial and tangential compressive stresses on the expanding front (Figure 1D). During cell proliferation, radial compressive stress is partially relieved by new biomass extending the biofilm boundary (Zhang et al. , 2016). By contrast, in the tangential direction, compressive stress becomes concentrated because there is no analogous relaxation mechanism. Therefore, starting from a flat film, a growing biofilm will undergo mechanical instabilities predominantly in the tangential direction, leading to radial wrinkling, and later, to delamination patterns (Figure 1D). By contrast, in the interior region of a biofilm, compressive stress occurs in both the radial and tangential directions, giving rise to a network containing both radially and tangentially oriented features (Figure 1A, D). To demonstrate that pattern directionality is determined by expansion anisotropy, we changed the biofilm growth geometry by inoculating cells starting from a line so that the biofilm would extend quasi-unidirectionally (Video 2). In this geometry, compressive stress along the inoculation line is higher than that perpendicular to the line (the expanding direction). Therefore, wrinkles or blisters occur perpendicular to the biofilm line (Figure 1E). Mechanical instability theory predicts that, for a film–substrate system that is subject to compressive stress, the wrinkling wavelength is determined exclusively by the thickness and mechanical properties of the relevant materials (Huang et al. , 2005). If so, we would expect the wrinkle wavelength to change with the mechanical properties of the biofilm and substrate but to be independent of the growth stage and geometry of the biofilm. To extract the wrinkle wavelength, we imaged the biofilm morphogenesis process over 72 hr and quantified the periodicity of radial stripes (Figure 2—figure supplement 1; Videos 3–5). We note that blisters emerge from wrinkles and that they inherit the wavelength of wrinkles, so we do not distinguish between the two in this analysis. We quantified the number of wrinkles or blisters N as a function of radial distance r from the biofilm center at different times. We found a linear relationship between N and r (Figure 2A, Figure 2—figure supplement 1). The slope has a geometrical origin: N = (2π/λ) r in which λ is the inherent wavelength of the system irrespective of the time in the developmental process or the location in the overall pattern (except at the biofilm core). A constant wavelength λ also means that radial wrinkles or blisters must bifurcate to maintain constant spacing as r increases, and indeed, we observed this to be the case (Figure 2A, inset). We also confirmed that the same λ was maintained when cells were inoculated in the line geometry and grew quasi-unidirectionally (Figure 2—figure supplement 1). We conclude that the wavelength of wrinkles or blisters reflects an intrinsic physical property of the biomechanical system. Mechanical instability theory also predicts how the wavelength varies with the stiffness contrast between the biofilm and the substrate. Classical linear stability analysis for bilayer film–substrate systems predicts that λ, normalized by the film thickness h, should be equal to 2π (Gf/3Gs) 1/3, in which Gf and Gs are the shear modulus of the film and the substrate, respectively (Chen and Hutchinson, 2004; Huang et al. , 2005). The 1/3 power law is a result of the competition between the energy cost to deform the film and that to deform the substrate. To test whether this relationship applies to biofilms, we measured λ, h, Gs, and Gf for all growth conditions. Gf varies minimally over a wide range of agar concentrations, whereas Gs varies by almost three orders of magnitude for agar concentrations from 0. 4% to 3% (Supplementary file 1 Table S1). Plotting λ/h versus Gf/Gs on a log-log scale (Figure 2B) reveals the characteristic scaling power law of 1/3, indicating the applicability of mechanical instability theory to biofilm morphogenesis. One key discrepancy exists between the experimental measurements and the bilayer model. Bilayer theory predicts that, if Gf/Gs < 1. 3, the substrate is too stiff for the flat-to-wrinkling transition to occur (Wang and Zhao, 2015). However, wrinkling occurs in our experiments for Gf/Gs well below 1. 3, corresponding to agar concentrations ≥ 0. 7%. To reconcile this discrepancy, we considered that a third soft, intermediate layer could exist between the growing biofilm and the non-growing substrate, which has been shown to allow wrinkling behavior even at low Gf/Gs ratios (Lejeune et al. , 2016a). To acquire evidence for an intermediate layer, we employed a capillary peeling method in which biofilms are gently dipped into water and the capillary force peels the biofilm off the substrate without destroying the biofilm or the underlying surface (Figure 2—figure supplement 2) (Yan et al. , 2018). Prior to peeling, using reflective confocal microscopy, the total biofilm thickness h was measured. After peeling, a residual layer remained on the substrate with a thickness hr (Figure 2C). Our preliminary analysis suggests that this layer consists primarily of matrix polysaccharide (Figure 2—figure supplements 2 and 3). Thus, the corrected biofilm thickness hf was obtained as hf = h – hr. We replotted our data using hf (Figure 2D, Figure 2—figure supplement 4). To rationalize the replotted curve, we took advantage of recent modeling efforts concerning multi-layer wrinkling phenomena (Lejeune et al. , 2016a). The only unknown parameter in our work is the shear modulus of the residual layer, Gr. In our theoretical model, we use a residual layer thickness hr = 0. 3hf, which was obtained from our experimental measurements, and we left Gr/Gf as a fitting parameter (Figure 2—figure supplement 4). The trilayer model qualitatively and quantitatively captures our experimental observations. Qualitatively, with a soft intermediate layer, the wrinkling pattern persists even when the substrate becomes stiffer than the biofilm (Gs > Gf). Unlike the bilayer model, in which the substrate is deformed by the wrinkling film, in the trilayer model, the soft interfacial layer assumes the major share of the deformation, effectively reducing the substrate stiffness (Figure 2D, Figure 2—figure supplement 4) (Lejeune et al. , 2016a). Quantitatively, predictions from the trilayer model recapitulate the prominent features of the revised plot: λ/hf scales according to the bilayer model as 2π (Gf/Gs) 1/3 for large Gf/Gs ratios, but increasingly deviates from the 1/3 scaling law for smaller Gf/Gs values. In the low Gf/Gs regime, wrinkling is increasingly controlled by the soft intermediate layer. An intermediate layer stiffness of Gr = 0. 1Gf allows the trilayer model to best fit our experimental data over all conditions. We next investigated the second transition predicted by our mechanomorphogenesis model: wrinkling-to-delamination. Whether and when a film–substrate system undergoes delamination is controlled by a competition between the adhesion energy between layers, Γ, and the elastic energy in the substrate. A dimensionless term Γ*, defined as Γ/ (hfG' s) in which G' s is the effective substrate modulus taking into account the residual layer (Lejeune et al. , 2016a; see also Materials and methods), was used previously to quantify the relative importance of the two energies (Wang and Zhao, 2015). We recently measured the biofilm–agar interfacial adhesion energy Γ ~ 5–10 mJ/m2 (Yan et al. , 2018). Hence, Γ* is in the order of 0. 01–1 in the current system, making delamination highly likely to occur during biofilm growth. In the context of the trilayer model, delamination takes place at the weakest interface, which is between the biofilm and the residual layer. To access the wrinkling-to-delamination transition experimentally, we simultaneously imaged the growing biofilm from the top and the side (Figure 3A, Figure 2—figure supplement 1). Radial wrinkles developed into blisters when growth proceeded beyond ~ 36 hr. At low agar concentrations, large amplitude blisters emerged among small amplitude wrinkles (Figure 3A). At higher agar concentrations, additional wrinkles developed into blisters, although with amplitudes smaller than those on low concentration agar substrates. We verified these findings using optical profiling to capture the full three-dimensional (3D) height information of the entire biofilm (Figure 3B). To peer inside blisters, we imaged cross-sectioned biofilms grown from cells producing fluorescence from mKate2 (Figure 3C). At low agar concentration (i. e. , 0. 6%), only a small fraction of wrinkles were detached from the substrate in the form of blisters (Figure 3—figure supplement 1). By contrast, at high agar concentration (i. e. 1. 0%), nearly all wrinkles had developed into blisters. In the cross-sectional images, voids were clearly present underneath the blisters, which were presumably filled with liquid (Wilking et al. , 2013). Figure 3D quantifies the positive correlation between the percentage of wrinkles that converted to blisters at the biofilm edge and the substrate agar concentration. To rationalize the dependence of the delamination pattern on agar concentration, it is useful to recall the notion of normalized adhesion energy, Γ*. On stiff substrates, Γ* is small, so delamination is favored over wrinkling. Blisters form extensively but they are small because they share the overall compression. On soft substrates, Γ* is large, so blisters form only infrequently while the majority of the biofilm remains attached to the substrate. In this case, the isolated blisters concentrate the compressive strain and become larger than those on a stiff substrate. We hypothesized that the locations of blisters on soft substrates are defined by surface defects that trigger local delamination. This hypothesis is consistent with the observed heterogeneous sizes of blisters in biofilms grown on soft substrates. Specifically, we argue that blisters emerge at different times and at different locations in growing biofilms depending on when a surface defect is encountered during biofilm expansion. The different ages of blisters naturally lead to their heterogeneous heights. To test this possibility, we made surface imperfections in the soft agar substrate at defined positions (Figure 3E). Indeed, these imperfections dictated the exact locations at which blisters formed as the biofilm expanded. By contrast, on stiff substrates, delamination occurred along the entire biofilm rim, irrespective of the predefined surface imperfections (Figure 3E). In conventional materials systems, a blister initially assumes a sinusoidal profile and then continues to grow in both width and height as the strain mismatch between the film and substrate increases (Vella et al. , 2009). We wondered how blister width and height would develop in a living biofilm as the biofilm expands and accumulates strain mismatch. To examine this, we tracked isolated blisters by imaging the rim of the expanding biofilm (Figure 2—figure supplement 1). The width of each biofilm blister decreased while its height increased over time until the final width of the blister reached twice that of the thickness of the biofilm (Figure 4A, B). This final value for the blister width indicates that the two sides of the blister come into contact with one another. Subsequently, blisters continue to develop only in height. Moreover, large blisters suppress nearby wrinkles from delaminating (Figure 4B), presumably because the biofilm and the substrate can slide relative to one another such that a blister captures nearby compressed biofilm material, and in so doing, releases compressive stress in the vicinity. Neighboring blisters tend to merge during late stages of biofilm development (>48 hr), forming single dark features in the transmission images (Figure 4C (top) and Figure 4—figure supplement 1). Indeed, cross-sectional images reveal that head-to-head contact occurred (Figure 4C (bottom) ). The sequential biofilm blister dynamics described above involve the generation or annihilation of new or existing interfaces, which have energy penalties or payoffs. To understand the order of these events, we measured their interfacial energies in WT V. cholerae biofilms (Yan et al. , 2018). They are: biofilm blister–liquid underneath, γfl ~ 49 mJ/m2; biofilm blister–air above, γfa ~ 30 mJ/m2; and the energy needed to separate the biofilm from the residual layer underneath, Γ ~ 5 mJ/m2 (Figure 4D). This energy hierarchy determines the sequence through which interfaces are generated or annihilated (Figure 4E). First, compressive stress leads to delamination of the biofilm from the residual layer, forming a local blister. This step generates an additional high-energy interface between the blister and the liquid underneath it. To eliminate this high-energy interface, the two sides of the inner face of the blister come into contact with each other as the blister grows. Indeed, electron microscopy imaging of the cross-section of a blister shows this to be the case (Figure 4—figure supplement 2). After internal contact occurs, the blister can only develop in the vertical direction. However, blister growth enlarges the interface between the biofilm and the air. Subsequent merging of neighboring blisters (Figure 4C) eliminates biofilm–air interfaces, and in so doing, lowers the free energy of the entire system. An added benefit to the bacteria stems from these blister dynamics: cells in blisters are less susceptible to the lethal effects of antibiotics that diffuse in from the substrate than are cells residing in the base of the biofilm, presumably because cells in blisters are located further away from the antibiotic source (Figure 4—figure supplement 3). If the above interpretations concerning the involvement of interfacial energy in blister development are correct, changing the relative magnitudes of the three interfacial energies should modulate blister dynamics, and, in turn, the global biofilm morphogenesis process. To test this idea, we deleted bap1 and rbmC, which encode proteins that are responsible for cell-surface interactions and biofilm hydrophobicity (Fong and Yildiz, 2007; Berk et al. , 2012; Hollenbeck et al. , 2014). Rather than forming isolated blisters, when formed on soft agar substrates, the Δbap1ΔrbmC biofilm exhibits a star-shaped morphology with flat regions between the facets of the stars (Figure 5A (top) ) (Yan et al. , 2017). The cross-section of a single facet shows that it consists of a group of congregated blisters (Figure 5A (bottom) ). Curiously, in contrast to the WT blisters, in the mutant, only the external surfaces of neighboring blisters are in contact with one another, leaving the internal spaces under each blister intact. Indeed, transmission images show that multiple stripes exist within one facet, corresponding to multiple blisters (Figure 5B, Figure 5—figure supplement 1). To rationalize the Δbap1ΔrbmC blister dynamics, we measured the relevant interfacial energies (Figure 5C). The adhesion energy Γ between the Δbap1ΔrbmC biofilm and the substrate is below the detection limit, meaning that delamination occurs more easily in the Δbap1ΔrbmC biofilm than in the WT biofilm. Indeed, Δbap1ΔrbmC biofilm blisters emerge directly from the expanding flat film, skipping the wrinkling state (Video 6). Second, the relative order of interfacial energies changes in the mutant: γfl approaches zero whereas γfa is large, consistent with the hydrophilicity of the Δbap1ΔrbmC biofilm (Hollenbeck et al. , 2014). These alterations in interfacial energies have profound consequences for blister dynamics (Figure 5D). Instead of annihilating biofilm–liquid interfaces inside of the blisters, in the mutant, neighboring blisters prefer to collapse against each other, which eliminates the high-energy interface between the biofilm and the air. Indeed, during the development of the mutant biofilm, newly emergent blisters move towards, and ultimately join, existing blister groups (Figure 5E; Video 6). The triangular shape of each facet in the Δbap1ΔrbmC biofilm is therefore a consequence of the merging of multiple blisters, whose ages and radial lengths decrease from the center to the edge of the aggregate. We wondered whether the emergence of the 3D biofilm surface topography affected biofilm expansion in the growing plane. One common morphological feature of bacterial biofilms is their irregular petal-shaped 2D contours (Videos 3 and 4). We hypothesized that the evolution of contours could also be a consequence of blister formation. To quantify the contour undulation, we define the acircularity parameter α = P2/4πA, in which P is the perimeter of the biofilm and A is the area (Asally et al. , 2012). α = 1 for a perfect circle. For a biofilm growing on soft agar (0. 4%, Figure 6A), there is a sharp increase in α at tc, the time at which the 3D surface morphology forms at the edge (Figure 6—figure supplement 1). To show that blisters are required for contour undulations, we tracked α for mutant biofilms lacking the matrix structural polysaccharide (ΔvpsL) (Figure 4—figure supplement 3; Hammer and Bassler, 2003) or lacking matrix structural proteins (ΔrbmAΔbap1ΔrbmC) (Figure 2—figure supplement 3C; Berk et al. , 2012; Yan et al. , 2017). In both cases, the biofilm has no surface features and α remains close to 1 (Figure 6A). To investigate the coupling between contour undulations and biofilm morphogenesis in the z direction, we followed the time evolution of growing biofilm borders in different geometries (Figure 6B, Figure 6—figure supplement 2). Visually, the indentations along the contours always correspond to the locations of large blisters. To quantify this finding, we measured the local curvature κ and expansion velocity Vf along the biofilm periphery (Figure 6B, C, Figure 6—figure supplements 1 and 2). Both κ and Vf are negatively correlated with the positions of blisters. Monitoring the evolution of a single blister and a nearby flat region shows a transient large difference in Vf when the blister initially forms at the edge (~ 45 hr in this case; Figure 6D), which triggers the local contour indentation. The emergence of a blister creates an extra dimension into which newly produced biomass can be distributed, which causes local slowing in Vf, thus establishing the correlation between blister locations and negative local curvature. After this transient difference, Vf becomes comparable for boundaries with and without blisters, and the local curvature reaches a steady value, provided that there is no nearby blister (Figure 4—figure supplement 1). In this steady state, the petal-like contour propagates radially without changing the overall shape of the contour. This explanation for the formation of the biofilm petal shapes suggests that contour undulations require non-homogeneous blister distribution along the biofilm rim and indeed, WT biofilms that are grown on stiff agar (>1. 0%) remain nearly circular because they possess regularly and closely spaced blisters (Figure 6A, blue line). As additional evidence for the connection between blister formation and boundary undulation, we show that we can control the number and positions of the petals by specifying the positions of the blisters using patterned substrates (Figure 3, Figure 6—figure supplement 2). We conclude that the 3D surface topography that arises owing to mechanical instabilities caused by biofilm expansion feeds back to slow down expansion and drive contour evolution.
  We show here that mechanical instabilities, including wrinkling and delamination, underlie biofilm morphogenesis. Moreover, differences in interfacial energies drive mechanomorphogenesis by dictating the creation or annihilation of new or existing interfaces. Finally, feedback between mechanomorphogenesis and biofilm expansion shapes the overall biofilm contour. Collectively, our findings concerning the connections between a biofilm’s surface morphology and its mechanical and material properties suggest that new genes and/or new compounds that alter biofilm morphology by altering mechanics could be discovered and investigated to address biofilm-related problems. Morphological patterns can certainly involve gene regulation programs. Nonetheless, we expect our mechanical instability findings in V. cholerae biofilms to apply to other systems — from bacteria to humans — because they reveal links between the specific material properties of the biological components and morphological transitions. Regarding bacterial systems, we have already commented on how localized cell death underpins pattern formation at the core of Bacillus subtilis biofilms (Asally et al. , 2012). In fact, in light of our mechanomorphogenesis model, localized cell death can be viewed as a source of surface defects that functions to trigger delaminations, similar to the defined surface imperfections that drive delaminations shown in Figure 3E. Another example concerns biofilms of Pseudomonas aeruginosa, an opportunistic pathogen (Costerton et al. , 1999). WT P. aeruginosa develops biofilms with a labyrinthine inner pattern surrounded by flat rims (Madsen et al. , 2015). By contrast, P. aeruginosa mutants that are incapable of phenazine production (Δphz) form biofilm topography similar to those that we examined here for V. cholerae with disordered cores surrounded by radial features (Dietrich et al. , 2013). We suggest that the mechanical principles uncovered here could also drive the morphological transitions in P. aeruginosa biofilms. The WT P. aeruginosa biofilm pattern occurs because cells at the biofilm center display upregulated matrix production (Madsen et al. , 2015), whereas cells located at the periphery are downregulated for matrix production. In the case of the Δphz mutant, all of the cells overproduce extracellular polysaccharides (Madsen et al. , 2015), so we speculate that the Δphz P. aeruginosa mutant forms peripheral radial wrinkles and subsequently delaminations because of the same mechanical instability described here in V. cholerae. These examples illustrate how gene regulation and spatially differentiated cell physiology can be coupled to mechanical instability to promote biofilm surface morphologies. Recent theoretical work on bacterial biofilms has considered mechanical instabilities. Zhang et al. (2016) used simulations to suggest that anisotropic growth coupled with wrinkling instability could generate the surface topography observed in bacterial biofilms, and most recently they considered the possibility of delamination (Zhang et al. , 2017). Wang and Zhao (2015) introduced competition between adhesive and elastic energies and computed a phase diagram of the different modes of instability for a film–substrate system. These inspiring theories will be made more valuable by the inclusion of measured biophysical parameters and additional observations generated through experiments. For example, the thin intermediate residual layer that we discovered here is not accurately considered in biofilm simulations, but is required to explain the wrinkling instability in biofilms (Figure 2D). In addition, interfacial energies play a predominant role in driving the morphologies of biological materials that possess soft layers, whereas their roles are minor in classical mechanical systems (Qi et al. , 2011). To date, contributions from interfacial energies have been suggested in contexts such as cell sorting in tissues (Brodland, 2002; Foty and Steinberg, 2005), but we are not aware of any work incorporating interfacial energies into mechanical instability models for morphogenesis. Future theoretical analyses can now incorporate measured parameters to understand the rich hierarchical dynamics and the history dependence of mechanomorphogenesis, taking into account biofilm viscoelasticity, interfacial energies, and the consequences of sliding and friction between the biofilm and the substrate (Beroz et al. , 2018; Peterson et al. , 2015). Though more sophisticated, eukaryotic organisms often employ similar mechanical instability principles to generate fascinating morphologies. Thus, our findings for biofilms are potentially generalizable and relevant for studies of development in higher organisms (Kim et al. , 2015). A close analogy is presented by cerebellum development. The cerebellum possesses a thin, soft layer of Purkinje cells that is sandwiched between the rapidly growing external granular layer and the slow-growing internal granular layer (Lejeune et al. , 2016b). Through wrinkling instabilities, the cerebellum develops finely spaced parallel grooves called folia. This hard-soft-hard geometry and the associated wrinkling instabilities directly mirror the configuration that we discovered in V. cholerae biofilms. Hence, our work suggests that exploiting mechanical principles to drive key morphogenic events is ancient: it occurs in bacteria, and evolution, as is often the case, has reused prokaryotic processes and principles in eukaryotes. In summary, biofilms represent an intriguing and highly tractable model system to investigate the general role of mechanical forces in morphogenesis, and they provide a convenient system for morpho-engineering.
  All of the V. cholerae strains used in this study are derivatives of V. cholerae O1 biovar El Tor strain C6706str2 (Thelin and Taylor, 1996), harboring a missense mutation in the vpvC gene (VpvC W240R) (Beyhan and Yildiz, 2007). Bacterial cultures were grown at 37°C under constant shaking in standard lysogeny broth (LB) medium. Genetic engineering of V. cholerae was performed using allelic exchange with pKAS32 (Skorupski and Taylor, 1996). All plasmids used in the current study have been reported previously (see' Key resources table' ). pKAS32-derived plasmids were introduced into V. cholerae by conjugation with Escherichia coli S17 λ-pir (de Lorenzo and Timmis, 1994), selection on plates containing ampicillin (100 mg/L) and polymyxin B (6 mg/L), and subsequent counterselection on plates containing streptomycin (500 mg/L). Deletions were verified by PCR and phenotypic analysis. The constitutive mKate2 gene (Shcherbo et al. , 2007) is driven by Ptac and was inserted into the V. cholerae chromosome at the lacZ locus (as previously described) with X-Gal (50 mg/L) present in the counterselection step (Nadell and Bassler, 2011). Biofilms of V. cholerae strains expressing mKate2 were grown on agar plates as described above. Where indicated, 0. 5 μM SytoX Green Nucleic Acid Stain (ThermoFisher) was added to the agar to stain dead cells. The region of the agar substrate containing a biofilm (~ 2. 5 cm × 2. 5 cm) was removed and transferred to an empty petri dish. O. C. T. agent (Tissue-Tek, Sakura) was applied to the surface of the biofilms, and the entire Petri dish was rapidly dipped into a dry ice–ethanol mixture to solidify the O. C. T. agent together with the biofilm. Razor blades were used to cut through the solidified samples. Samples with exposed cross-sections were immediately transferred to a homemade T-shaped sample holder and kept frozen in a dry ice–ethanol mixture. These samples were transferred to a Leica stereoscope and imaged in bright-field mode or in fluorescent mode with an mCherry or GFP filter set. The surface profiles of biofilms grown for 48 hr were analyzed with a Leica DCM 3D Micro-optical System. A 10× objective was used to image a 3 mm x 3 mm region covering roughly one quarter of the biofilm, with a z step size of 2 μm. To measure the thickness of the residual layer, agar plates containing biofilms were slowly vertically lowered into water to peel the biofilms from the substrate. The entire agar plate was allowed to air dry for 5–10 min to remove liquid remaining from the peeling step. After drying, the above analysis procedure was performed to measure the thickness of the residual layer. The total thickness of the biofilm h and the thickness of the residual layer hr were measured using Leica Map software. A three-point flattening procedure was first performed on the agar surface to level the image. Next, line profiles were generated at three different locations spanning the agar surface to the surface of the biofilm or the residual layer. An automatic step-size detection procedure was performed with a built-in function in the software to extract h or hr. The three measured values were averaged to give the value for one biological replicate. The biofilm thickness hf was obtained by hf = h − hr. Biofilms were grown on 0. 6% agar plates for 2 days as described above. The region of the agar containing a biofilm (~ 2 cm × 2 cm) was separated from the remainder of the agar plate, transferred to a piece of glass, and placed horizontally in a 50 mL conical tube and frozen at −80°C overnight followed by overnight lyophilization (Millrock Technology, BT85A-A). The biofilm samples were sliced with a razor blade to expose blisters, sputter-coated with a 5 nm layer of Pd (VCR IBS/TM200S ion beam sputterer), adhered to an upright SEM stub with conductive tape, and imaged with a scanning electron microscope (FEI XL30 FEG-SEM). Biofilms of V. cholerae strains constitutively expressing mKate2 were inoculated onto semipermeable membranes (EMD Millipore VSWP02500) that had been placed on top of 0. 6% agar. The plates were incubated at 37°C for 2 days. The semipermeable membranes were gently removed from the agar surface using tweezers, and subsequently floated at room temperature overnight on top of 3 mL LB medium containing 1. 7 μM SytoX Green stain with or without 50 μg/mL tetracycline. Image analyses were performed with custom codes written in MATLAB and with ImageJ software. Raw transmitted light image data were first converted into intensity images. From the pixel intensity distributions, we identified the peak with the highest intensity Ib and used it as background. We set the minimum intensity Imin = 0 and the average background intensity Ib = 0. 9 to standardize the contrast of the images. Images were then smoothed with a median filter. From the intensity distribution, we also identified the intensity value IV of the valley immediately adjacent to the background peak and used it as the thresholding value to binarize the image (using a built-in thresholding function in MATLAB). We separated the biofilm object F from the background. We used the image of each biofilm at t = 12 hr after inoculation to define the center OF for all time points. When mutations affecting biofilm morphology arose, they were manually excluded from the image analysis. To quantify variations in the amplitudes of biofilm morphological features, we extracted the intensity profiles IE (θ) along a circle near the biofilm edge. We use a built-in function in MATLAB to identify the positions and the prominence ΔIp of the peaks in –IE (θ). We set the minimum peak prominence to be 0. 02 to eliminate noise. To extract the periodicity of the wrinkling or delamination pattern, we tracked the time evolution of these patterns from images. For wavelength analysis, we applied fast Fourier transformation (FFT) to intensity functions Ir (θ) in a ring at time t and radial coordinate r, and identified N (r, t) from the peak frequency in the power spectrum. We also verified the values by autocorrelation and manual counting. We plotted all data from different time points and fitted them with a linear function N (r) =2πr/λ to obtain the intrinsic wavelength λ. The radial coordinate at which N decreases to zero was defined as Rp. For images of biofilms grown in a line geometry, several values of N were extracted from multiple lines at different distances from the central line, averaged, and subsequently used to extract λ. For contour analyses, we first obtained the biofilm object F from the binarized image. From the binarized object F, we extracted the perimeter P and the area A of region F. At each time point, we calculated the acircularity α as α = P2/4πA. To define the radii for biofilms that were not strictly circular, we used <Rf> = <|ri–rO|>i, averaged over all the points ri on the circumference ∂F. <Rf> was then calculated over time to give <Rf (t) > versus t. Segmented linear regression with two segments was used to quantify the expansion velocity of the biofilm <Vf> before and after the critical time tc and to define the critical time itself. To capture local curvature κ and expansion velocity Vf, the smoothed boundary ∂F was locally approximated by quadratic polynomials ri, 2 (t) at ri. The parametrized curve xi, 2 (t) and yi, 2 (t) allowed us to calculate the analytical curvature κi and normal ni locally using the weighted central difference. Coarse-grained contours at time points t and t+Δt were then connected by joining ri (t) to its nearest neighbor ri (t+Δt) in ∂Ft+Δt, yielding local velocities Vf, i = |ri (t+Δt) – ri (t) |/Δt. To analyze the side-views of blisters, blister contours were manually extracted with ImageJ software and then smoothed. The baseline of the blister was obtained by averaging the z coordinate of the left and right bottom region of the blister. The blister height H was calculated as the distance between the peak of the blister to the baseline. The width of the blister W was measured at half of the blister height. We adapted a trilayer model from previous work (Lejeune et al. , 2016b), and modeled the biofilm system with the following three elastic components: the biofilm (top), the residual layer (middle), and the agar substrate (bottom) denoted with subscripts f, r, and s, respectively. V. cholerae biofilms harbor an active growing top cell layer and a dead cell layer underneath (Figure 4—figure supplement 2). The live and dead cell layers are connected to each other, and they were removed together for our mechanical measurements; so in the model, we do not distinguish between the two and we treat them as a single biofilm layer. Biofilm and residual layers were modeled as thin elastic sheets with thickness hf and hr, whereas the agar substrate was modeled as an elastic body with a thickness hs, much larger than that of the other two layers. The relevant scale for the continuum model is about the thickness of the film (>50 μm). Therefore, we could neglect potential structural and materials heterogeneities in the biofilm, which exist on a much smaller scale (~ 5 μm, see Yan et al. , 2018). The shear modulus and Poisson’s ratio of the materials are denoted by G and ν, respectively. For theoretical calculations, we treated all three layers as incompressible materials and hence, ν = 0. 5 (see the above experimental measurements). In the simulation, the residual layer grows at the same rate as the biofilm layer, while the substrate does not grow (as confirmed by comparing the locations of the edge of a biofilm and the residual layer; see Figure 3—figure supplement 1). This growth difference induces a strain mismatch ε between the biofilm/residual layer and the substrate. Following previous studies (Lejeune et al. , 2016a), we applied the Föppl-von Kárman equation to the biofilm model. Assuming a sinusoidal profile of the surface undulations, we can write the longitudinal stress S in the film as: S (n) = Gfhf2n23+K~hfn2, where n is the wave number and K~ is the combined stiffness of the residual layer and the substrate layer: K~= 4Gsnnhr (Gs/Gr−1) +2, and the effective substrate modulus of the composite substrate can be calculated by Gs′=K~hs′ in which h' s is the total depth of the strained region (see Lejeune et al. , 2016a for details). By numerically solving the nonlinear equation dS/dn = 0, we determined the minimal critical value of S for mechanical instability and the corresponding n gives the critical wavenumber ncr. The wavelength at the onset of wrinkling was then calculated as λcr=2π/ncr. The critical stress and strain were obtained by Scr = S (ncr) and εcr = Scr/3Gf, respectively. Theoretical predictions from the bilayer model can simply be calculated by setting Gs = Gr. The model described above, despite assuming only small strains, accurately predicted the wavelength and critical stress/strain for finite strains (Lejeune et al. , 2016a). We verified that the analytical predictions were in reasonable agreement with results obtained from finite element simulations. The only unknown parameter in the model is the shear modulus of the residual layer Gr, which is difficult to probe experimentally. Therefore, we treated Gr as the only fitting parameter. We used hr/hf = 0. 3 as an average value from the relevant experimental data and fit the model against the experimental data for wavelength versus stiffness contrast between the biofilm and the agar substrate. Fitting was carried out by minimizing the least-square error between the theoretically predicted and the experimentally measured wavelengths. A bisection method was employed that converged in fewer than 10 iterations. A plane-strain computational model was developed to take into account growth, large deformations, and the nonlinear elasticity of the system. We considered the same planar three-layer structure as above. According to finite strain theory, we define the deformation gradient tensor as Fij=∂xi/∂Xj, where xi and Xi denote the coordinates in the deformed and undeformed configurations, respectively (Ogden, 1997). To incorporate the effect of growth, we further introduced the decomposition of the deformation tensor F = FeFg as the product of the growth deformation Fg and the elastic deformation Fe (Figure 2—figure supplement 4) (Rodriguez et al. , 1994). We used Fg= (1+g001) for the biofilm and residual layers to describe their 1D growth (g > 0) in the X1 direction, and we set Fg to be the identity matrix I for the non-growing agar substrate. The growth-induced compressive strain is thus ε=g/ (1+g). To account for the nonlinear stress-strain behavior of materials undergoing large deformations, all three layers were modeled as neo-Hookean materials. The strain energy density of each layer is given by Ogden (1997): Ψ (Fe) =μe2IC-2-2ln⁡J+λe2ln⁡J2, where μe and λe are the Lamé parameters, and they are related to the shear modulus G and Poisson’s ratio ν byμe=G, λe=2Gυ1-2υ. IC = tr (FeTFe) is the first invariant of the right Cauchy-Green deformation tensor C = FeTFe, and J = det (Fe). The total elastic energy of the system can thus be calculated byΠ=∫Ωf ΨFe, fJg, fdX+∫Ωr ΨFe, rJg, rdX+∫Ωs ΨFe, sJg, sdX, where Ωf/r/s denotes the volume occupied by biofilm/residual/substrate in the initial undeformed reference configuration, and Jg = det (Fg) specifies the volume element change following growth. We assumed that the present instability pattern always seeks the lowest potential energy among all possible configurations at any time during biofilm growth, neglecting the viscoelasticity and plasticity of the biomaterials that could potentially lead to hysteresis in mechanical instability. For the computational model, we considered a rectangular domain Ω = Ωf∪Ωr∪Ωs = [0, L]×[0, hf+hr+hs] composed of three layers, where L denotes the size of the system. We use subscripts 1 and 2 to denote the horizontal and vertical components, respectively. Numerically, the task is to calculate the displacement field ui = xi - Xi that minimizes the total potential energy, that is u=argminu∈Vu⁡Π, where Vu is the function space that satisfies the boundary conditions on u. Without loss of generality, we considered a scenario in which the biofilm and residual layers grow together but are confined by the left and right walls of the bottom fixed domain Ω, that is, the boundary conditions were set by u1|X1=0=u1|X1=L=u2|X2=0=0 (Figure 2—figure supplement 4). The nonlinear constrained minimization problem was implemented in the open-source computing platform FEniCS (Alnæs et al. , 2015). The computational model was discretized by first-order triangular elements generated by Gmsh (Geuzaine and Remacle, 2009), and the accuracy of the results was verified by mesh refinements. A growth increment of ∆g = 0. 002 was employed in the simulations, up to a maximum of 1. For each step, we computed the equilibrium configuration x and the Green-Lagrange strain tensor e = 0. 5 (FeTFe – I) of the system. The critical condition for wrinkling instability was identified as a vertical displacement of the biofilm that surpassed the threshold value (0. 01hf). We further calculated the deviatoric strain tensor e´ij = eij – 0. 5δijekk and the von Mises equivalent strain εvM = (2e´ij e´ij /3) 1/2 (Jones, 2009) to visualize the strain distribution among the three layers. All results were visualized by Paraview software (Ahrens et al. , 2005). For the model parameters, we set hr/hf = 0. 3 based on the measured thickness values from experiments, and hs/hf = 10 to represent the thick substrate. The stiffness contrast Gr/Gf = 0. 1 was used according to the optimal fitting value from theoretical curves, and we varied Gf/Gs from 0. 02 to 10 to correspond to the experimental conditions. In all simulations, L was set to be larger than 10 times the wavelength to minimize the finite size effect, and the Poisson’s ratios of all three layers were set to be 0. 45 to ensure convergence of the algorithm. Error bars correspond to standard deviations of the means. Standard t-tests were used to compare treatment groups and are indicated in each figure legend. Tests were always two-tailed and unpaired/paired as demanded by the details of the experimental design. All statistical analyses were performed using GraphPad Prism software. The custom-written MATLAB scripts and simulation codes used in this study are available at https: //github. com/f-chenyi/biofilm-morphogenesis (Fei, 2019; copy archived at https: //github. com/elifesciences-publications/biofilm-morphogenesis).
summary: Engineers have long studied how mechanical instabilities cause patterns to form in inanimate materials, and recently more attention has been given to how such forces affect biological systems. For example, stresses can build up within a tissue if one layer grows faster than an adjacent layer. The tissue can release this stress by wrinkling, folding or creasing. Though ancient and single-celled, bacteria can also develop spectacular patterns when they exist in the lifestyle known as a biofilm: a community of cells adhered to a surface. But do mechanical instabilities drive the patterns seen in biofilms? To investigate, Yan, Fei, Mao et al. grew biofilms of the bacterium called Vibrio cholerae – which causes the disease cholera – on solid, non-growing ‘substrates’. This work revealed that as the biofilms grow, their expansion is constrained by the substrate, and this situation generates mechanical stresses. To release the stresses, the biofilm initially folds to form wrinkles. Later, as the biofilm expands further, small parts of it detach from the substrate to form blisters. The same forces that keep water droplets spherical (known as interfacial forces) dictate how the blisters evolve, interact, and eventually shape the expanding biofilm. Using these principles, Yan et al. could engineer the biofilm into desired shapes. Collectively, the results presented by Yan et al. connect the shape of the biofilm surface with its material properties, in particular its stiffness. Understanding this relationship could help researchers to develop new ways to remove harmful biofilms, such as those that cause disease or that damage underwater structures. The stiffness of biofilms is already known to affect how well bacteria can resist antibiotics. Future studies could look for new genes or compounds that change the material properties of a biofilm, thereby altering the biofilm surface.
keywords: physics of living systems
title: Mechanical instability and interfacial energy drive biofilm morphogenesis