Mechanical instability and interfacial energy drive biofilm morphogenesis

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.


Introduction
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.

A mechanical instability model for biofilm morphogenesis
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., eLife digest 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. 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 Figure 1. Mechanical instability drives V. cholerae biofilm morphogenesis. (A) Bright-field images of biofilms grown for 2 days on the designated percentages of agar. (B) Schematic of the wrinkling and delamination processes that occur during biofilm expansion. Red with a black outline denotes the biofilm. Gray denotes the substrate, agar in this case. (C) Three-dimensional (3D) profile of two colliding biofilms, initially inoculated 9 mm apart, grown on a 0.6% agar plate for 36 hr. (D) Transmission image of a V. cholerae biofilm grown for 35 hr (top) and 48 hr (bottom) on a 1.0% agar plate. (E) Transmission image of a V. cholerae biofilm inoculated as a line and grown for 30 hr on a 0.5% agar plate. In panels (D) and (E), blue arrows denote the expansion directions, and black arrows denote the tangential directions along which compressive stress accumulates. All scale bars are 5 mm. DOI: https://doi.org/10.7554/eLife.43920.003 The following source data and figure supplements are available for figure 1: 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/m 2 ), which is the same order of magnitude as the deformation energy of the soft substrate . 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). Video 1. Part 1: A V. cholerae biofilm grown for 24 hr on 0.6% agar medium was peeled off of the substrate by the capillary method using LB medium as the liquid starting from the bottom left. The movie is played in real time. Part 2: The peeled biofilm from Part 1 grew at the air-liquid interface over time. Imaging began immediately after peeling and its total duration is 6 hr with 5-min time steps. The field of view is 73.0 mm Â 48.3 mm. DOI: https://doi.org/10.7554/eLife.43920.007 A trilayer mechanical model predicts the biofilm wrinkling wavelength 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 2figure supplement 1). The slope has a geometrical origin: N = (2p/l)r in which l 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 l 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 l 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 Video 2. Part 1: Collision of two V. cholerae biofilms grown on medium containing 0.6% agar. Imaging began 5 hr after inoculation and has a total duration of 75 hr with 15 min time steps. Biofilms were separated by 9 mm at the time of inoculation. At t = 20 hr, the biofilms begin to contact one another. The additional compressive stress present at the collision front leads to the formation of a large blister in the middle. The field of view is 41.5 mm Â 27.7 mm. Part 2: Growth of a V. cholerae biofilm on medium containing 0.5% agar after cells were inoculated in a line. Imaging began 5 hr after inoculation and has a total duration of 72 hr with 15 min time steps. The field of view is 50.2 mm Â 33.3 mm. DOI: https://doi.org/10.7554/eLife.43920.008 Video 3. Growth of a V. cholerae biofilm on medium containing 0.4% agar. Imaging began 5 hr after inoculation and has a total duration of 75 hr with 15 min time steps. The field of view is 41.5 mm Â 27.7 mm. DOI: https://doi.org/10.7554/eLife.43920.019 Video 4. Growth of a V. cholerae biofilm on medium containing 0.7% agar. Imaging began 5 hr after inoculation and has a total duration of 75 hr with 15 min time steps. The field of view is 41.5 mm Â 27.7 mm. DOI: https://doi.org/10.7554/eLife.43920.020 predicts that l, normalized by the film thickness h, should be equal to 2p(G f /3G s ) 1/3 , in which G f and G s 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 l, h, G s , and G f for all growth conditions. G f varies minimally over a wide range of agar concentrations, whereas G s varies by almost three orders of magnitude for agar concentrations from 0.4% to 3% (Supplementary file 1 Table S1). Plotting l/h versus G f /G s 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 G f /G s < 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 G f /G s 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 G f /G s 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) . 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 h r ( Figure 2C). Our preliminary analysis suggests that this layer consists primarily of matrix polysaccharide (Figure 2figure supplements 2 and 3). Thus, the corrected biofilm thickness h f was obtained as h f = h -h r . We replotted our data using h f ( 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, G r . In our theoretical model, we use a residual layer thickness h r = 0.3h f , which was obtained from our experimental measurements, and we left G r /G f 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 (G s > G f ). 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: l/h f scales according to the bilayer model as 2p(G f /G s ) 1/3 for large G f /G s ratios, but increasingly deviates from the 1/3 scaling law for smaller G f /G s values. In the low G f /G s regime, wrinkling is increasingly controlled by the soft intermediate layer. An intermediate layer stiffness of G r = 0.1G f allows the trilayer model to best fit our experimental data over all conditions.
The biofilm wrinkling-to-delamination transition is controlled by interfacial energy and substrate stiffness 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, G , and the elastic energy in the substrate. A dimensionless term G * , defined as G/(h f G' s ) in which G' s is the effective substrate modulus taking into account the residual layer (Lejeune et al., 2016a; see also Materials and  The corrected biofilm thickness h f was obtained by subtracting the residual thickness h r from the total thickness h. The solid portion of the black line corresponds to the prediction from the bilayer model, which applies only to x coordinates greater than 4.75 (Wang and Zhao, 2015). The dashed portion of the black line is an extrapolation to zero from the bilayer prediction provided as a guide to the eye. The red line is the fitted data from the trilayer model in which the stiffness contrast between the residual and biofilm layers G r /G f is treated as a fitting parameter while holding h r /h f = 0.3. Inset: finite-element simulation of the trilayer model undergoing wrinkling instability. Red denotes the biofilm. Gray denotes the substrate. Blue denotes the residual layer. Simulation parameters were chosen to mimic the growth condition on 1.0% agar (black arrow). Data are represented as mean ± std with n = 3. DOI: https://doi.org/10.7554/eLife.43920.009 The following source data and figure supplements are available for figure 2: Source data 1. Experimental measuremants of biofilm residual layer thicknesses and wavelengths and predictions from trilayer wrinkling theory. DOI: https://doi.org/10.7554/eLife.43920.010 Figure supplement 1. Analysis of intrinsic wavelengths in the morphologies of biofilms. Figure 2 continued on next page 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 G~5-10 mJ/m 2 . Hence, G * 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, G * . On stiff substrates, G * is small, so delamination is favored over wrinkling. Blisters form extensively but they are small because they share the overall compression. On soft substrates, G * 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).

Interfacial energy controls blister development dynamics and interactions between blisters
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   Figure 3D. (B) Developing profile of a single blister, extracted from side view images at successive time points after delamination. Profiles are shown at 2.5 hr (gray line), 10 hr (gray dotted line) and 17.5 hr (black dashed line) after the onset of delamination. The distance between the red arrows corresponds to W, which, over time, approaches twice the biofilm thickness (2h f ). Regions near the blister become flatter as cell mass is pulled into the blister. Agar concentration: 0.4%. (C) Representative merging of adjacent blisters (white arrows) at specified times (top). Cross-section image from a biofilm producing fluorescent mKate2 reveals blister peak-topeak contact (bottom; designated by the white arrow). Agar concentration: 0.7%. Scale bars: 1 mm (top) and 0.5 mm (bottom). (D) Interfacial energy of the biofilm-air interface g fa , biofilm-liquid interface g fl , and the adhesion energy between the biofilm and the substrate G for WT V. cholerae biofilms. Data are represented as mean ± std with n = 3. Inset: schematic of different interfaces. (E) Schematic of blister development in a WT V. cholerae biofilm. White stars and dashed black lines denote interface annihilation events. For panels (D) and (E), the color code is the same as that in Figure 3D. 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 . They are: biofilm blister-liquid underneath, g fl~4 9 mJ/m 2 ; biofilm blister-air above, g fa~3 0 mJ/m 2 ; and the energy needed to separate the biofilm from the residual layer underneath, G~5 mJ/m 2 ( 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 Dbap1DrbmC 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 Dbap1DrbmC blister dynamics, we measured the relevant interfacial energies ( Figure 5C). The adhesion energy G between the Dbap1DrbmC biofilm and the substrate is below the detection limit, meaning that delamination occurs more easily in the Dbap1DrbmC biofilm than in the WT biofilm. Indeed, Dbap1DrbmC 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: g fl approaches zero whereas g fa is large, consistent with the hydrophilicity of the Dbap1DrbmC 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 Dbap1DrbmC 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.

Mechanical instability and biofilm expansion feed back onto one another
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 a = P 2 /4pA, in which P is the perimeter of the biofilm and A is the area  (Asally et al., 2012). a = 1 for a perfect circle. For a biofilm growing on soft agar (0.4%, Figure 6A), there is a sharp increase in a at t c , 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 a for mutant biofilms lacking the matrix structural polysaccharide (DvpsL) (Figure 4-figure supplement 3; Hammer and Bassler, 2003) or lacking matrix structural proteins (DrbmADbap1DrbmC) (Figure 2-figure supplement 3C; Berk et al., 2012;Yan et al., 2017). In both cases, the biofilm has no surface features and a 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 k and expansion velocity V f along the biofilm periphery ( Figure 6B,C, Figure 6-figure supplements 1 and 2). Both k and V f 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 V f 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 V f , thus establishing the correlation between blister locations and negative local curvature. After this transient difference, V f 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 nonhomogeneous 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.

Discussion
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 (Dphz) 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 Dphz mutant, all of the cells overproduce extracellular polysaccharides (Madsen et al., 2015), so we speculate that the Dphz 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 . 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 slowgrowing 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.

Bacterial strains
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 l-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 P tac 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).

Biofilm growth
Biofilm growth on agar plates LB medium solidified with different percentages of agar was used as the solid support to grow biofilms. V. cholerae strains were streaked onto LB plates containing 1.5% agar and grown at 37˚C overnight. Individual biofilms were selected and inoculated into 3 mL of LB liquid medium containing~10 glass beads (MP Biomedicals Roll and Grow Plating Beads, 4 mm diameter) and the cultures were grown with shaking at 37˚C to mid-exponential phase (5-6 hr). Subsequently, the cultures were mixed by vortex to break clusters into individual cells, OD 600 was measured, and the cultures were back-diluted to an OD 600 of 0.5. 1 mL of these preparations were spotted onto prewarmed agar plates. Subsequently, the plates were incubated at 37˚C. Typically, four biofilms were grown per agar plate. For time-lapse imaging, one or two biofilms were grown on each plate.

Biofilm growth on substrates with defined defects
On prewarmed agar plates, syringe needles (BD PrecisionGlide needles, 0.6 mm Â 2.5 mm) were used to punch holes at eight locations, equally separated by 45˚around a circle. Marks were made on the bottoms of the Petri dishes to guide our eyes for placement of holes in the agar surface. 1 mL of V. cholerae cultures at OD 600 = 0.5, prepared as described in the preceding paragraph, were spotted at the center of the circle. The diameter of the circle was~1 cm for biofilms grown on 0.6% agar and~0.6 cm for biofilms grown on 1.0% agar. Different circle diameters were used to accommodate the differently sized biofilms that form on soft and stiff agar, and to guarantee that, in both cases, when biofilms expanded to cover the pre-defined defects, they remained flat. Following biofilm growth, the positions of these defects were inferred from the marks drawn on the bottoms of the Petri dishes.

Biofilm growth in a line geometry
A V. cholerae culture at OD 600 = 0.5 was prepared as described above. A sterile razor blade was carefully dipped into this culture and dried in air for 1 min. The razor blade was gently touched to the surface of a prewarmed agar plate to initiate biofilm growth.

Biofilm growth at the liquid-air interface
First, a biofilm was grown for 24 hr following the procedure describe above. Subsequently, 25 mL of LB medium was gently added from the edge of the agar plate. When the liquid reached the biofilm, the liquid lifted the biofilm off the substrate by capillary force.

Biofilm imaging Bright-field imaging
Biofilms were imaged with a Leica stereoscope in the reflective (bright field) mode. For biofilms larger than the field of view, multiple overlapping images were acquired manually (3 by 3 or 3 by 2) at different locations in the biofilm. Images from multiple locations in biofilms were stitched together with the Image Composite Editor software from Microsoft to yield the full images of the biofilms while preserving the original resolution. Raw images from the stereoscope contain iridescence as the result of reflections from agar, which were removed by setting the color saturation to zero (i.e. converting to black-and-white images).

Transmission imaging
A custom transmission imaging setup was built in a 37˚C environmental room to follow biofilm growth. Briefly, an agar plate containing the inoculum was placed on an LED illumination pad (Huion L4S Light Box) and imaged with a Nikon D3300 SLR camera equipped with a Sigma 105 mm F2.8 Macro Lens. The entire setup was covered to exclude light. The camera was controlled using Digi-CamControl software. Imaging was started 5 hr after inoculation, at which time the camera was capable of focusing on the growing biofilm. Imaging was performed automatically every 15 min for 3 d. The growth of the biofilm floating at the air-liquid interface was monitored with images acquired at 5 min intervals.

Side view imaging
A similar setup to the one described in the preceding paragraph was used to image biofilms from the side, with the following changes. First, the LED illumination pad was placed on the side so that the camera received scattered light from the biofilm surface. Second, an additional camera (Nikon D3300 SLR equipped with DX Zoom-Nikkor 18-55 mm lens) was also placed on the side of the biofilm, at an~90˚angle with respect to the first optical path. To remove the optical obstruction from the wall of the agar plate, an imaging window (~1 cm Â 1 cm) was created using a hot razor blade. Imaging started immediately before the onset of the wrinkling-to-delamination transition, and the time interval between images was 5 min. From time to time, the focus in the side view was adjusted manually.

3D optical profiling
Biofilms were imaged with a Keyence VR-3200 optical profiler using a telecentric multi-triangulation algorithm. Subsequent analyses related to obtaining the 3D profiles of biofilms were performed with the Keyence Analyzer software. In brief, noise was first removed from the raw data using the built-in function in the Keyence Analyzer software to give smooth, continuous surface profiles. Surfaces corresponding to agar were excluded by setting upper and lower height thresholds. 3D views of biofilms were rendered with a built-in function in the software. The corresponding line profiles were extracted along an arc centered at the center of the biofilm.

Cross-sectioning of biofilms
Biofilms of V. cholerae strains expressing mKate2 were grown on agar plates as described above.
Where indicated, 0.5 mM 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.

Rheological measurements Shear rheology of biofilms
All rheological measurements were performed with a stress-controlled shear rheometer (Anton Paar Physica MCR 301) at 37˚C. For each measurement, 100-960 biofilms were collected with a pipette tip or a razor blade and transferred onto the lower plate of the rheometer. After sandwiching the biofilm cells between the upper and lower plates with a gap size of 0.5 mm, silicone oil (5 cSt at 25˚C, Sigma Aldrich) was applied to surround the biofilm. Sandblasted surfaces were used for both the upper and lower plates to avoid slippage at the boundary. Oscillatory shear tests were performed with increasing amplitudes of the oscillatory strain "' from 0.01 to 2000% at a fixed frequency of 6.28 rad/s. The storage modulus G' was extracted with the RheoPlus software as a function of "'.
To extract the plateau shear moduli of biofilms, segmented linear fittings were applied to G'-"' curves on a log-log scale. G' varies minimally in the plateau region. We used the fitted G' value at "' = 1% as the modulus of the biofilm G f . All rheological properties of the biofilm remained roughly constant for at least 48 hr.

Shear rheology of agar
LB medium containing different agar concentrations was freshly prepared in 100 mL bottles. The semi-solid medium was heated in a microwave, cooled to~55˚C, and added (2 mL) to the lower plate of the rheometer preheated to 60˚C. The heated agar solution was subsequently sandwiched between the two rheometer plates with a gap size of 0.5 mm and sealed with silicone oil. The preparation was cooled to 22˚C using a cooling rate of 1˚C/min. Subsequently, the solid agar was heated to 37˚C for measurement. This procedure mimics the sequence of events that agar plates were exposed to during preparation and biofilm growth. Smooth surfaces with TrueGap technology were used. Oscillatory shear tests were performed in the linear elastic region at a fixed frequency of 6.28 rad/s. For data obtained with agar, we averaged 10-20 points in the plateau region of the G'("') curve to give G s .

Poisson ratio measurement
The Poisson ratio n of the biofilm was estimated by compressing the biofilm in the vertical direction and measuring its bulk modulus. Briefly, a home-built hollow cylinder made of polytetrafluoroethylene with a diameter of 25.5 mm was placed between two parallel plates of a rheometer. The biofilm was loaded into the cylinder to fill its volume. The upper plate of the rheometer (with a diameter of 25 mm) was subsequently lowered with a constant velocity (of between 8 mm/s and 12 mm/s). During this measurement, the shaft does not rotate, but rather acts as a piston to measure the normal force. Using the relationship between normal force and shaft displacement, we calculated the bulk modulus K of the biofilm to be~130 kPa; much larger than the shear modulus G'. From these data, we could calculate the Poisson's ratio n = (3K -2G') / 2(3K + G') » 0.495, close to the incompressible limit (n = 0.5).

Biofilm thickness measurements
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 mm. 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 h r 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 h r . The three measured values were averaged to give the value for one biological replicate. The biofilm thickness h f was obtained by h f = h À h r .

SEM sample preparation and imaging
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).

Measurement of colony-forming units
Biofilms grown for 2 days were peeled off of agar substrates using a phosphate-buffered saline solution (PBS) as described previously . The floating biofilms were collected with clean pipette tips and the corresponding residual layers were removed from the agar using a sterile razor blade. All samples were transferred to 1.5 mL microcentrifuge tubes containing 1 mL PBS and~0.2 mL small glass beads (acid-washed, 425-600 mm diameter, Sigma), vigorously mixed by vortex for 15 min at 37˚C to break apart aggregates, serially diluted in PBS, and plated onto LB plates. The LB plates were incubated overnight at 37˚C and subsequently assessed for colony forming units (CFU). Four biological replicates were measured, each with two technical replicates. Raw CFU values were normalized by the volume of each biofilm and residual layer, calculated from the radius and thickness of each biofilm and residual layer, respectively.

India Ink staining
Biofilms grown for 2 days were peeled off of agar substrates with PBS as described above. 1 mL of Higgins Black India ink solution (10% in PBS) was added to the agar to cover the area containing an intact biofilm or a residual layer, and the preparation was air-dried at room temperature for 30 min. The stained residual layer was subsequently imaged with a Leica stereoscope in the bright-field mode.

Antibiotic killing assay
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 mM SytoX Green stain with or without 50 mg/mL tetracycline.

Biofilm image analyses
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 I b and used it as background. We set the minimum intensity I min = 0 and the average background intensity I b = 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 I V 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 O F 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 I E () along a circle near the biofilm edge. We use a built-in function in MATLAB to identify the positions and the prominence DI p of the peaks in -I E (). 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 I r () 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) =2pr/l to obtain the intrinsic wavelength l. The radial coordinate at which N decreases to zero was defined as R p . 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 l.
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 a as a = P 2 /4pA. To define the radii for biofilms that were not strictly circular, we used <R f > = <|r i -r O |> i , averaged over all the points r i on the circumference qF. <R f > was then calculated over time to give <R f (t)> versus t. Segmented linear regression with two segments was used to quantify the expansion velocity of the biofilm <V f > before and after the critical time t c and to define the critical time itself.
To capture local curvature k and expansion velocity V f , the smoothed boundary qF was locally approximated by quadratic polynomials r i , 2 (t) at r i . The parametrized curve x i , 2 (t) and y i , 2 (t) allowed us to calculate the analytical curvature k i and normal n i locally using the weighted central difference. Coarse-grained contours at time points t and t+Dt were then connected by joining r i (t) to its nearest neighbor r i (t+Dt) in qF t+Dt , yielding local velocities V f,i = |r i (t+Dt) -r i (t)|/Dt.
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.

Theoretical modeling procedures
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 h f and h r , whereas the agar substrate was modeled as an elastic body with a thickness h s , much larger than that of the other two layers. The relevant scale for the continuum model is about the thickness of the film (>50 mm). Therefore, we could neglect potential structural and materials heterogeneities in the biofilm, which exist on a much smaller scale (~5 mm, see Yan et al., 2018). The shear modulus and Poisson's ratio of the materials are denoted by G and n, respectively. For theoretical calculations, we treated all three layers as incompressible materials and hence, n = 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: where n is the wave number andK is the combined stiffness of the residual layer and the substrate layer:K and the effective substrate modulus of the composite substrate can be calculated by G 0 s ¼Kh 0 s 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 n cr . The wavelength at the onset of wrinkling was then calculated as l cr ¼ 2p=n cr . The critical stress and strain were obtained by S cr = S(n cr ) and " cr = S cr /3G f , respectively. Theoretical predictions from the bilayer model can simply be calculated by setting G s = G r .
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 G r , which is difficult to probe experimentally. Therefore, we treated G r as the only fitting parameter. We used h r / h f = 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.

Computational modeling procedures
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 F ij ¼ qx i =qX j , where x i and X i 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 = F e F g as the product of the growth deformation F g and the elastic deformation F e (Figure 2-figure supplement 4) (Rodriguez et al., 1994). We used for the biofilm and residual layers to describe their 1D growth (g > 0) in the X 1 direction, and we set F g to be the identity matrix I for the non-growing agar substrate. The growthinduced 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): where e and l e are the Lamé parameters, and they are related to the shear modulus G and Poisson's ratio n by e ¼ G; l e ¼ 2G 1 À 2 : I C = tr(F e T F e ) is the first invariant of the right Cauchy-Green deformation tensor C = F e T F e , and J = det(F e ). The total elastic energy of the system can thus be calculated by where W f/r/s denotes the volume occupied by biofilm/residual/substrate in the initial undeformed reference configuration, and J g = det(F g ) 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.

Finite element simulations
For the computational model, we considered a rectangular domain W = W f S W r S W s = [0, L]Â[0, h f +h r +h s ] 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 u i = x i -X i that minimizes the total potential energy, that is u ¼ arg u2Vu min P, where V u 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 W, that is, the boundary conditions were set by u 1 j X1¼0 ¼ u 1 j X1¼L ¼ u 2 j X2¼0 ¼ 0 (Figure 2-figure supplement 4). The nonlinear constrained minimization problem was implemented in the open-source computing platform FEniCS (Alnaes 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 Dg = 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(F e T F e -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.01h f ). We further calculated the deviatoric strain tensor e´i j = e ij -0.5d ij e kk and the von Mises equivalent strain " vM = (2e´i j e´i j /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 h r /h f = 0.3 based on the measured thickness values from experiments, and h s /h f = 10 to represent the thick substrate. The stiffness contrast G r /G f = 0.1 was used according to the optimal fitting value from theoretical curves, and we varied G f /G s 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.

Statistical methods
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.
. Supplementary file 1. Summary of biomaterial parameters for V. cholerae biofilms. Table S1 reports the measured biomaterial parameters for V. cholerae biofilms grown on different concentrations of agar substrates. These measurements include the shear modulus of the substrate and the biofilm, the thickness of the biofilm and the residual layer, and the wavelength of the biofilm surface pattern. Data availability All data generated or analyzed during this study are included in the manuscript and supporting files. Source data files have been provided for all figures, tables and figure supplements.