Stiffness transitions in new walls post-cell division differ between Marchantia polymorpha gemmae and Arabidopsis thaliana leaves

Significance During morphogenesis, plant cells divide and undergo significant shape changes. The mechanical properties of the cell wall during these two processes are important for plant morphogenesis. We introduce a systematic method to map cell wall age and growth to its bulk elasticity. We demonstrate that the stiffness of the wall correlates with their growth and that the new walls in Marchantia polymorpha gemmae become transiently stiffer and slower-growing compared to the older walls, a phenomenon not seen in Arabidopsis thaliana leaves. Using computational modeling, we were able to link the properties of the new walls to a local cell shape change in the two species—junction angle between the new and the old walls.

Table S5.Computational model results.The pinch-in angles, the new division walls' lengths and areas from the inflation study of a patch of hexagonal cells showing the consequence of different wall properties compared to a control after a hexagon has been divided in two.All results were taken when the cells of interest grow and exceed the same set area.Initial conditions used as the starting point for the simulations, mesh was inflated with no growth and uniform parameters of a Young's modulus= 100 MPa, turgor pressure= 0.2 MPa that will be the value used for all the simulations.

Methods
Immunostaining and image analysis of Marchantia Gemmae.Immunofluorescence of cell wall epitopes was performed as previously described (Alonso-Serra et al., 2020) with slight modifications.Four 48h old gemmae were fixed in a solution consisting of 4% formaldehyde (freshly prepared from paraformaldehyde powder, Sigma) and 0,5% glutaraldehyde (Sigma) in a 0.1M phosphate buffer at pH 7. Fixation, dehydration and resin infiltration steps were all micro-wave (MW) assisted using a PELCO BioWave Pro (Ted Pella, Redding, CA).Fixation was realized at 150W, under vacuum (20Hg) (5x 1').Samples were left in the fixative overnight at 4°C and then washed 3 times in PBS.Samples were then processed through increasing dehydration steps (25%, 50%, 70%, 90%, 96%, 3x 100% Ethanol, vacuum 20Hg, MW 150W 5').Resin infiltration (LR White medium grade, Agar scientific) was then realized through increasing resin concentration: 33% Resin in ethanol 100%, 66% Resin in ethanol 100%, and 3 times 100% Resin (20Hg, MW 200W 5').Samples were left overnight in 100% resin for effective resin penetration in the samples.Samples were then aligned parallelly prior to resin polymerization, which was realized at 60°C during 17h. 1 µm thick sections were cut using an EM UC7 ultramicrotome (Leica Microsystems), mounted on Menzel Ultrafrost Plus microscopy slides on water droplets and air dried.Mounted samples were blocked with blocking solution (2% BSA in PBS supplemented with 0.05% Tween 20) for 30 min at room temperature.Primary antibodies LM19, LM20, LM12, LM15, LM25, LM23 and LM21 (Plant Probes, Leeds, UK) were applied as 1/100 dilutions in blocking solution over night at 4°28 C and slides washed 5 times in blocking solution.Alexa Fluor 488-conjugated Goat anti-Rat IgM (ThermoFisher Scientific, A-21212) secondary antibodies were as 1/200 dilutions in blocking solution for 2 h at room temperature and washed 5 times in blocking solution.For callose immunolocalization, callose primary antibody (Biosupplies) was applied in 500 µl as 1/400 dilution and secondary antibody Alexa Fluor 555-conjugated Goat anti-Mouse (ThermoFisher Scientific, A21425) secondary antibodies were applied as 1/200 dilutions in blocking solution for 2 at room temperature and washed 5 times in blocking solution.Slides were finally mounted in a 1:1 solution of AF1 antifadent (Citifluor) with PBS containing 0.1% (w/v) Calcofluor White (Merck/Sigma, F3543) as a cell wall counterstaining and then imaged by confocal laser scanning microscopy (Zeiss LSM700) and slides were sealed with nail varnish.Sections were subsequently imaged by confocal laser scanning microscopy as z-stacks (Zeiss LSM700).Fluorescence analysis of callose immunolocalizations was performed using a custom-made macro on ImageJ software.The sequential steps of the macro are as follows: images, acquired as z-stacks were first processed to maximal projections; a region of interest restricted to the cell division area is first hand drawn by the user; within the defined cell division region of interest, a binary mask is applied with a user defined threshold to the cell wall counterstaining channel (here calcofluor) in order to create a region of interest (ROI) restricted to cell walls in the cell division area (the creation of a region of interest restricted to the cell wall avoid to underestimate the levels of fluorescence measured in different conditions, as the studied epitopes are exclusively located to the cell wall); a second region of interest representing the older cell walls is defined as the inverse of the cell division area, and restricted to the cell walls using again a binary mask to restrict the ROI to the cell wall counterstaining; defined ROIs are finally applied to the immunolocalized channels, to measure the average fluorescence levels.The detailed macro code can be consulted in the github repository (3).

Statistical analysis of immunostaining experiments.
Data analysis used the Tidyverse R package collection and the ggplot2 package for boxplots and barplots.Statistical analysis was performed using the R functions for ANOVA and Tukey HSD tests for multiple comparisons after having determined that parametric tests were applicable using R functions for Bartlett and Shapiro test.Significance values for P < 0.05 were grouped using the agricolae package (with alpha = 0.05).In the case parametric tests could not be applied, statistical analysis was performed using the R functions for Kruskal-Wallis one-way Modelling.Hexagonal shaped cells were generated by defining the nodes' positions and connections.The 3D meshes of the generated hexagonal outlines were obtained using a custom code (4) a Delaunay-based unstructured mesh generation algorithm from Darren Engwirda (5).All meshes were created so that the triangle edges on the cell boundaries were approximately 0.5 µm long, where they smoothly increase in size as they move away from the boundary to reduce the degrees of freedom in the system.This mesh was then converted into a compatible form to inflate in Tissue, (1,6,7) software developed by the Jönsson group.
Tissue uses a finite element method approach to approximate solutions to the hyperelastic continuous mechanical equations over the discretized mesh surface, which are being stretched by a normal (outward from each cell) pressure force applied to each mesh triangle.It uses the 'Triangular Bi-Quadratic Springs' (TRSB) method, which uses a St Venant Kirchoff formulation and approximates it using biquadratic springs, which resist triangle edge and inner angle deformation (8) and allows us to simulate cell walls efficiently.In more detail, the St Venant Kirchoff formulation stress-strain relation for an isotropic material is, where W is the total strain energy in a domain Ω, E the Green-Lagrange strain tensor and λ and µ are the Lame coefficients and defined in plane elasticity as λ = Eν 1−ν 2 and µ = E(1−ν) 1−ν 2 with E being the Young's modulus and ν the Poisson coefficient.
Using TRBS, W is approximated over a triangle T as where k T i and c T k are the tensile and angular stiffness of the biquadratic springs and Li and Li,0 are the current and resting lengths of edge i in the triangle.These stiffnesses are defined as where αi is the angle opposite edge i in the resting triangle configuration and AT the resting triangle area.A force on each node can then be calculated from this strain energy as a result of the stretching triangles.Additionally, deriving the stress tensor for each triangle allows us to find the maximum stress magnitude and direction from the maximum corresponding eigenvector and eigenvalue.
The pressure force is applied to all the triangles proportional to their area in the outward normal direction.A solution to the simulation is found when the system has reached mechanical equilibrium, i.e. when the pressure forces on each node acting on them from their surrounding triangles are balanced with their surrounding stretched and strained triangles, such that the total force is now 0. We find this mechanical equilibrium using the Newton-Raphson method (9) and setting the tolerance to 1e − 6.
The hexagons' growth simulations were performed by first inflating to mechanical equilibrium and saving their state.The cross wall of interest was then given a higher stiffness and lower growth rate to replicate the effect of a higher apparent young modulus from the AFM data.The system was then allowed to grow using a spring with a variable resting length, a formulation previously introduced to model plant growth (10,11).The resting length L t i,0 of the element i at time t + δt with time step (δt) is given by the resting length in the previous time step L t i,0 , plus the increment that is proportional to their strain and their current length.This formulation takes the form (1, 12) where gi is the growth rate of that extensibility rate edge and Li is the current edge length.
All simulations were inflated with a turgor pressure of 0.2 MPa, Young's modulus of 100 MPa and a Poisson's ratio of 0.2 and the growing hexagons had a extensibility rate of 0.2 hr −1 .The values chosen for the turgor pressure and the walls Young's modulus are in line with previous studies (13,14).The stiffer cell wall has a extensibility rate of 0.1 hr −1 and a Young's modulus of 140 MPa when we vary both the extensibility rate and stiffness with the third contrapositive simulation having the same parameters but just reversed.In the second simulation where we vary just the stiffness, the stiffer cell wall has a Young's modulus of 500 MPa (Table S4).
The simulations depicted in Figure 4A-D were stopped when the cells of interest (the one with the different cell wall properties); area, A satisfies where A0 is the area of the cell at the beginning of the simulation to match with when the angles in Figure 4(G) were measured.
The black lines depicting the cell boundaries in Figure 4A-D were found using the indicator labels stored for each triangle in our mesh letting us know whether it starts the simulation on top of the cell or not (on the periclinal surface or not).From this we can then find the edge elements which are between one triangle on the top (the periclinal wall) and another on the side (the anticlinal wall).
Fig. S2.(A) Comparison of cell wall epitopes immunostaining signal intensity.Visual comparison of immunostaining signals from several cell wall epitopes on 48h old gemmae cross-sections through the meristematic area.Images are represented in Fire LUT with a 8 bit colour coding of original grey levels (256 levels) represented on the figure top left.An example of cell wall counterstaining (calcofluor) is also represented.Arrowhead indicates the position of the meristematic area.Scale bar: 100 µm.(B) Immunolabeling quantification in early gemmae (2 days old) shows that LM20 fails to detect methyl-esterified pectin in M. polymorpha gemmae.(C) Bright field image of the sectioned gemma for immunostaining where the blue line indicates the position of the section and the black contour is the gemma edge (scale bar 100 µm).

Fig. S3 .BFig. S4 .Fig. S5 .
Fig. S3.LM19 immunostaining cross sections on 48 HAG M. polymorpha gemmae.Images are represented in Fire LUT with a 8 bit colour coding of original grey levels (256 levels).A dashed-white box has been placed to indicatively highlight the dividing zone.Bright field images of the sectioned gemma for immunostaining, where the black curves indicate the contour of the gemmae and the blue thick lines indicate the location of the sections for each sample.

Fig. S7 .MPaFig. S8 .
Fig. S7.Cell division in A. thaliana first true leaf follows geometrical clues.(A) Heat map of the ratio of cell area over the 12 hour time period shown in days after stratification, scale bar 100 µm.White lines show the major axis of expansion.(B) Heatmap of the number of cell divisions in the period of time shown, scale bar 100 µm.(C) Sequential imaging of cortical microtubules (MT) for extraction of cell shapes and MT orientation of an active cell division area, scale bar 50 µm.Maximum projection of cortical MT signal on extracted cell shapes shows the main MT orientation (yellow lines on left picture), while comparison of cell shapes between the two different time points allows the calculation of the maximum growth direction (green lines on left picture).Analysis of the cell shapes allows calculation of the shortest path connecting two non-consecutive cell walls (blue lines on left picture).Comparison of the MT main orientation and shortest path direction with the actual plane of division for the dividing cells (see same letters between top and bottom picture) shows that shortest path orientation predicts the orientation of the division plane.(D) Quantification of angles between the MT main orientation and the actual division orientation (yellow distribution), and of the angles between the shortest path and the actual division orientation (blue distribution), and of the angles between the Principle Growth Direction and the actual division orientation (pink distribution) for the cells for which the MT orientation and shortest path do not coincide with the division plane (cells for which both MT orientation and shortest path possess an angle of divergence to division plane < 15 have been removed), total number of cells 86.The yellow, blue and pink solid lines represent respectively the cumulative distribution for the three predictions-MT orientation, shortest path and perpendicular direction to PGD.A comparison of the cumulative distributions for the MT orientation and the shortest path shows that cell division in A. thaliana better follows the shortest path (Anderson-Darling test, p-value <0.0001).

Fig. S9 .
Fig. S9.Modelling details. (A) Quantification of the pinch-in angle 24 HAD in M. polymorpha gemmae and quantification of the cell growth area of the dividing cells 24 HAD.(B) Mesh convergence study performed on a regular hexagonal cell with increasing mesh density.

Table S4 . Computational model parameters. Properties of the new wall (stiffness and growth rate) and properties of the surrounding walls (stiffness and growth rate).
0019.All p-values are computed as Wilcoxon Rank Sum Test).Note that the cells undergoing AFM testing are plasmolysed.Hence, this analysis illustrates a trend, and the specific quantitative values obtained are not directly comparable to physiological conditions.
Alessandra Bonfanti, Euan Thomas Smithers, Matthieu Bourdon, Alex Guyon, Philip Carella, Ross Carter, Raymond Wightman, Sebastian Schornack, Henrik Jönsson, Sarah Robinson A Fig. S6.(A) No correlation between height and Young's modulus is observed in AFM maps -4 independent patches of cells from M. polymorpha gemmae.The plots refer to anticlinal walls only.These plots suggest that the stiffer new cell walls in M. polymorpha are not correlated to the fact that new cells might have a lower curvature.(B) Correlation between Apparent Young's modulus and growth for the M. polymorpha samples that underwent two AFM tests (AFM time course), for M. polymorpha gemmae, and for A. thaliana leaves.