Mechanical instability and interfacial energy drive biofilm morphogenesis
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 ) .
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 .
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 .