Root Hair Adhesion in Posidonia oceanica (L.) Delile Seedlings: A Numerical Modelling Approach

Animals and plants use adhesion to move, to anchor to a substrate, or to disperse seeds and fruits. Some plants developed a root pad as a common strategy to adhere to consolidated substrates. In the marine environment, the seagrass Posidonia oceanica attaches firmly to consolidated substrates via adhesive root hairs, forming a pad structure. We used novel morphological and ultrastructural data to develop a numerical model to study the dynamics of root hair adhesion during contact formation on rough consolidated substrates for this species. Morphological analysis, conducted using Scanning Electron Microscope, highlighted the role of root hair branching in pad formation. Transmission Electron Microscope microscopy allowed us to identify a glue-like substance at the pad/substrate interface. The numerical model highlighted the role played by the cell wall’s elasticity in pad formation and its importance in guaranteeing a firm adhesion. Furthermore, the effectiveness of these mechanisms was assessed at different simulated roughness levels. Increasing knowledge on the adhesion mechanism of seagrass to consolidated substrates could be pivotal in developing advanced seedling-based restoration protocols. The findings of this study could contribute to restoration activities planned to contrast seagrass regression. Transplanting initiatives using seedlings can now better address the search for suitable and low-impact ways to fix germinated plants to the substrate.


INTRODUCTION
Living organisms employ a great variety of attachment structures, such as hooks, suction cups, sticky pads, and glue-like substances (Scherge and Gorb, 2001), to give them the capability to climb (e.g., English Ivy), walk on smooth surfaces (e.g., geckos), be transported (e.g., zoochory), or be strongly fixed to consolidated substrates (e.g., mussels and barnacles). Plants use adhesion mainly to disperse fruits and seeds, but some species form adhesive structures to anchor to the substrate. These pad-like structures often consist of multi-layered multicellular formations that can originate from tendrils, modified leaves, or roots (Groot et al., 2003;Melzer et al., 2010;Steinbrecher et al., 2010;Steinbrecher et al., 2011;Melzer et al., 2012;Seidelmann et al., 2012;Bohn et al., 2015;Yang and Deng, 2017). The pad attachment occurs via mechanical interlocking which consists in replicating the surface of a given substrate, filling up the surface profile, and eventually becoming stiff by thickening the cell walls. Adhesion can act with or without the support of a glue. This mechanism has been shown to produce a firm attachment on a variety of surfaces in different plant species (Scherge and Gorb, 2001).
Seagrasses have thrived in the oceans for tens of millions of years (Larkum et al., 2018), but it is only in recent decades that they have begun to experience a severe regression due to human activities along coastal areas (Waycott et al., 2009;Telesca et al., 2015). Therefore, several attempts are underway to restore damaged meadows and to improve their ecological resilience (Unsworth et al., 2015). Seagrasses natural recovery usually occurs trough vegetative proliferation. In a species that lacks floating rhizomes, the recolonization of strongly damaged meadows may occur via sexual propagules. In this case, the success of seedling settlement and recruitment represents a bottleneck in the seagrass sexual reproduction strategy as its failure could impair the connectivity and genetic variability of these important habitat-forming species (Jahnke et al., 2016).
Most seagrass species have been reported to grow on soft bottoms (Green and Short, 2003;Balestri et al., 2015), but a few of them (i.e., Phyllospadix spp., Posidonia oceanica) possess specific morphological adaptations that enable them to colonize and grow on hard bottoms (Barnabas, 1994;Badalamenti et al., 2015). These species show a variety of morphological features, such as rudimental anchorage apparatuses or more sophisticated structures (e.g., adhesive root hairs) to adhere to consolidated substrates (Barnabas, 1994;Green and Short, 2003;Badalamenti et al., 2015).
The combination of knowledge on the specific morphology of the seagrass root system with the characteristics of the habitats in which each species settles and recruits is a prerequisite for developing projects aimed at advancing restoration protocols. In recent years, many transplanting projects have been planned to restore damaged meadows, and the use of sexual propagules could be pivotal in achieving restoration goals. In this context, an increase in the current knowledge on the mechanisms of root adhesion, which provide the plantlet with a strong anchorage and persistence against drag forces on the sea floor, is of fundamental importance .
In the Mediterranean Sea, seedlings of the endemic seagrass Posidonia oceanica (L.) Delile successfully colonize hard bottoms and their settlement occurs preferentially on consolidated and firm substrates via adhesive root hairs . Such a strategy represents a mechanism for early settlement on vegetated and unvegetated rocky shores, which favors plantlet persistence on consolidated substrates compared to unconsolidated ones Badalamenti et al., 2015). Root hairs appear soon after germination at the hypocotyl and both at the principal and adventitious roots. When not in contact with the substrate, root hair tips remain long and straight, but they quickly begin to branch after contact formation Badalamenti et al., 2015;Zenone et al., 2020). P. oceanica may thrive in habitats characterized by high hydrodynamism (Montefalcone et al., 2016) hence its anchorage system has to be highly efficient, especially in the early life stages, when seeds reach the recruitment sites and settlement occurs.
Posidonia oceanica roots form a pad starting from the branching of a single cell, the root hair, which determines the anchoring via mechanical interlocking (Zenone et al., 2020). Root hairs adapt themselves to the characteristics of the substrate on which the attachment process occurs. The pad morphology depends on the substrate micro-roughness and optimal adhesion is achieved at a micro-roughness of 3-52 µm (asperities diameter), while the maximum adhesion force occurs at a micro-roughness of 12 µm, a value similar to the average diameter of the root hair (Zenone et al., 2020).
In recent years, numerical models have been elaborated to integrate a variety of important biological and physical processes involving plant growth (Godin and Sinoquet, 2005;Dupuy et al., 2010). This approach is useful to address biological questions, as models may show hidden properties and behaviors that can be validated with experimental observations. In the present paper, we report on a numerical model that is inspired by the morphological characteristics of the P. oceanica root hairs to investigate the dynamics of P. oceanica root hair adhesion during contact formation on rough consolidated substrates.
The model was constructed to answer the following questions: 1) Does the root hair growth change vary the cell's elastic properties when substrate micro-roughness is kept constant? 2) Does the contact formation change, varying the levels of substrate roughness, when the cell's elastic properties are kept constant? and 3) Can attachment properties of the root hair be explained exclusively by mechanical interlocking? In this work, we hypothesize that: 1) root hair softness plays a fundamental role in the adhesion process; 2) the substrate micro-roughness at to the scale of the root hair diameter improves proper contact area formation and additionally contributes to the mechanical interlocking; and 3) a glue-like substance may be additionally involved in the adhesion mechanism, as hypothesized by Zenone et al. (2020).

MATERIALS AND METHODS
Posidonia oceanica seeds were collected along the Sicilian coast (NW Sicily) between May and June 2017 and left growing ( Figure 1A) for 6 months on 3 cm × 3 cm tiles of epoxy resin replicas (Spurr, 1969) of sand paper (12 µm asperities size, Starcke Ersta Abrasives, Germany). Resin replicas were created using a two-step molding method (Gorb, 2007).
Different microscopy techniques were used to investigate the root hairs morphology and their principal characteristics. All the morphological features were measured on images from a set of 30 seedling samples (see Table 1 for details) by means of ImageJ software (Schneider et al., 2012).
Mean curvature and mean pad radius were calculated according to Eqs. 10-14 (see below for details) using Matlab (ver. R2019b, The Math Works, Inc., (2018) United States). For this, a batch of homogeneously distributed points (N ≥ 100) were selected along the pad perimeter. The coordinates of the points were used for calculation of the mean pad radius and the mean curvature.

Root Hair Tips Morphology (Scanning Electron Microscope)
Root hair morphology was analyzed using a Scanning Electron Microscope (SEM, Hitachi TM3000; Hitachi High-Technologies, Corp., Tokyo, Japan). Portions of resin tiles with root hairs adhered on their surfaces were dehydrated in hexamethyldisilane (HMDS/1,1,1,3,3,3-Hexamethyldisilazan, Carl Roth GmbH & Co. KG, Karlsruhe, Germany). For HMDS drying, tiles were cut into pieces of 1 mm 3 , fixed in a solution of 2.5% Glutardialdehyde and 2% Formaldehyde in 0.1 M cacodylate buffer overnight, and washed two times in 0.1 M cacodylate buffer, each for 20 min. They were then stained in a solution of 1% OsO4 in a 0.1 M cacodylate buffer for 1 h, then washed two times again in the 0.1 M cacodylate buffer for 20 min. Samples were then dehydrated by a graded series of ethanol (30 and 50%, 15 min each) and transferred to a 70% ethanol solution. Samples were then air-dried in a desiccator overnight. Finally, samples were mounted on SEM aluminum stubs, coated with gold palladium (10 nm thickness; Leica Bal-TEC SCD500, Leica Microsystems GmbH, Wetzlar, Germany), and immediately visualized in the SEM.

Root Hair Ultrastructure (Cryo-Scanning Electron Microscope and Transmission Electron Microscope)
Fresh portions of roots were visualized at high resolution using a cryo-SEM Hitachi S-4800 (Hitachi High-Technologies Corp., Tokyo, Japan) equipped with a Gatan ALTO2500 cryopreparation system (Gatan, Inc., Abingdon, United Kingdom). Fresh root portions were fixed to the holder and immediately placed in the cryo-chamber (−140°C). The frozen samples were sputter-coated with gold-palladium (layer thickness 10-20 nm) and transferred to the SEM.
The root hair tips at the interface with the substrate were examined using a Transmission Electron Microscope (TEM) (Tecnai Spirit BioTWIN; FEI Company, Eindhoven, Netherlands). Samples were cut into small pieces of about 1 mm 3 , fixed in 2.5% glutardialdehyde and 1% OsO4, dehydrated by a graded alcohol series to 100%, and gradually transferred to Spurr resin which was then cured at 70°C for 12 h. Resin blocs were trimmed and cut into semithin (1-3 µm) and ultrathin (50-60 nm) sections with a Leica EM UC7 ultramicrotome (Leica Microsystems GmbH; Wetzlar, Germany). Ultrathin sections were post-stained with gadolinium (III) acetate tetrahydrate and lead citrate and examined in TEM.

Numerical Modelling
In order to understand interactions between individual root hairs and rigid substrates of different roughness levels, we developed a  Frontiers in Mechanical Engineering | www.frontiersin.org November 2020 | Volume 6 | Article 590894 simplified two-dimensional numerical model aimed at solving the complex original 3-dimensional problem. We substituted the root hair material by its boundary with an effective "stiffness" reminiscent of the surface tension. This approach is based on an assumption that young root hair tips, being thin-walled bladders with a fluid inside, might behave similarly to a fluid-like system in contact with a solid substrate. Further, we defined the boundary conditions as a movable and growing "elastic chain." It allows us to naturally describe a multivalued function of the boundary and its evolution over time. The rigidity of the chain corresponds to the surface tension of real material and allows simulating the system under consideration.
The model includes the following elements. The boundary of the root hair is modelled by an elastic fiber whose elasticity varies along its length, depends on time, and can change from initially soft to much stiffer or almost rigid (but still with some degree of flexibility). The fiber is constructed as an array of nodes connected by the segments r → j {x j , y j }. Longitudinal F → jj' and transversal F →⊥ j stiffness of each segment is simulated by the following interactions between the segments: and The coordinates at the beginning and end of the segment represent the two nearest neighbors j; j ' j ± 1 and j 1, . . . , N s . Here, N s is the total number of segments which varies with time according to the rules which will be described below. We estimated the P. oceanica root hair elastic constant K ⊥ j as 225 N/m ( Table 1).
The longitudinal force, F → jj ' is described by a two-valley potential which tends to keep the distance between the nodes r → j and r → j±1 close to the equilibrium length of the segment dr. The transversal elastic force, F →⊥ j , holds r → j close to the position between its nearest neighbours ( r → j+1 + r → j−1 )/2, and therefore tends to keep the angle between the neighbouring segments close to 180°.
It is important to note that the "chain" simulates a boundary of the root hair pad, so its beginning and its end must be connected. It means that the last element of the array r → Ns {x Ns , y N s } must be initially placed at a distance from the first one r → 1 {x 1 , y 1 } equal to the segment length dr and it must elastically interact with it, according to the same interactions given by Eqs. 1 and 2, as all other segments.
A further demand from the model is that the root should grow with time. To reproduce root hair elastic properties, the length of each segment must be kept more or less fixed. So, an increase of root hair length will be simulated by changing the number of the FIGURE 3 | Typical configurations of the boundary for three different elastic constants, K ⊥ 01 , K ⊥ 02 , and K ⊥ 03 , are shown in the subplots (A), (B), and (C), respectively. Different colors corresponding to these three cases (blue, green, and red, respectively) are maintained in the subsequent plots. The substrate structure is shown by gray scale map, where white and black correspond to the highest and lowest asperities height, respectively.
FIGURE 4 | Time dependencies of the curvature radius (A) and mean radius of the root hair (B) calculated for the same initial elastic constants: K ⊥ 01 , K ⊥ 02 , and K ⊥ 03 as above. The meaning of the colors is the same as in Figure 2.
Frontiers in Mechanical Engineering | www.frontiersin.org November 2020 | Volume 6 | Article 590894 5 segments N s . In other words, the size of the nodes array increases with time. Moreover, according to the experimental observation, the elasticity of each segment starting from its appearance changes with the time too. In other words, both elastic coefficients K jj ' and K ⊥ j are time-dependent. The resulting distribution of the elasticity along the root at every particular moment should reflect the gradients of the elasticity between softer newly generated segments and harder ones in the older parts of the boundary. In this case, we considered a sort of aging of the root hair in the model, intended as a hardening with a characteristic time of about 40 min ( Table 1).
Counting on the particular nature of the system, we will insert a new node r → k {x k , y k } between two already existing nodes r → k and r → k+1 after time interval Δt. We will do this randomly along the whole chain. If we suppose that the boundary grows isotropically (i.e., branching occurs at different portions of the root hair pad), the random number j should be uniformly distributed in interval j [1: N s ]. However, the boundary of real root has preferable centers of growth independently on the substrate asperities. To simulate this, we defined some segments, where insertion of the new segments occurs with a higher probability. To be definite, we defined five symmetrically placed orientations (in originally isotropic boundary), where new segments can be inserted with 10% higher probability than into other ones.
The ratio between Δt and discrete time steps dt of the calculation, Δt/dr regulates the rate of the root growth. From a mathematical point of view, the only limitation is Δt/dr ≫ 1, but biologically the rate must also be adjusted to the total time of the process (real time of root growth: 7,000 min, Table 1).
In general, root hairs are very soft at the growing tip. However, their hardiness increases along the shaft (Grierson et al., 2014). In the model, the elastic constant K ⊥ k of each newly created segment is very small and increases with time. Mathematically, it means that every segment has its life counter τ j which starts from τ j 0, when the segment is "just born," and increases synchronously with common time steps: τ j (t + dt) τ j (t) + dt. In other words, each value of K ⊥ j starts from min(K ⊥ j ) K ⊥ k at τ j 0 and goes asymptotically at τ j → ∞ to the maximal rigidity max(K ⊥ j ) λK ⊥ k , where coefficient λ > 1. For a real system, the coefficient λ can be very large, but for the transparency of final plots of the simulation results, we used λ 3. Of course, resulting structures in all the numerical experiments are defined by a competition between the The contact surface can be modelled in different ways depending on the particular problem. A widely accepted approach is based on the supposition of scale-invariant fractal distribution of the surface heights. In this case, one can construct the surface using Fourier transform with a large number of corresponding waves. Depending on the number of waves and on both limits, the surface will be, more or less, close to the real fractal one and include or not include irregularities with small sizes (Filippov and Popov, 2007).
Like many other natural surfaces, here, the surface was physically created by a set of deposited particles (spheres with radius R, 12.75 ± 2.42 µm, Table 1). A natural way to produce such a surface mathematically is to apply the depositing of spheres procedure. In the purely numerical approach, different roughness of the surfaces can be achieved by depositing mathematically defined spheres with radius R on an originally flat surface. In this case, one takes an array of spheres with a radius close to R and puts them successively on randomly chosen positions {x n , y n }, with n 1, 2 . . . , N. Total number N of the spheres should be chosen to cover the surface with a density close to the incrementally used one. Each sphere is added virtually to the corresponding segment of the surface: δz n (x; x n ; y; y n ) R 2 − (x − x n ) 2 − (y − y n ) 2 . Due to the random deposition, the spheres can fall either on top of already existing coverage or onto the flat base surface. As a result of such a deposition, the total surface gradually accumulates all the spheres z(x, y) N n 1 δz n (x; x n ; y; y n ) which are, generally speaking, placed on top of each other. It essentially complicates the procedure, which should be repeated many times during the numerical simulations. FIGURE 6 | Equilibrium z-positions of the elastic boundaries plotted along closed curve as the functions of the index j of nodes normalized to the number N c (thin lines). The subplots (A-C) correspond to the three different roughness levels (red -171 µm, green -12 µm, blue -6 µm asperities size, respectively). Random realizations of the rigid surfaces are plotted by the bold lines.
Frontiers in Mechanical Engineering | www.frontiersin.org November 2020 | Volume 6 | Article 590894 One can simplify the procedure when considering that, for our aim, the specific form of the deposed particles is not very important. In such a case, one can treat the random deposition procedure as an abstract way to construct the desired surface. One of the simplest ways to do this is to model the rough surface Z(x, y) by a random deposition of Gaussians: The advantage of this approach is its mathematical transparency, since the simple manipulation of the position, amplitude, and width of the Gaussian functions will suffice to generate the desired kind of surface irregularity. In particular, the typical distance between the hills and valleys of the randomly generated surface Z(x, y) n G n can be regulated by the number of Gaussian curves. An additional convenience of this method is the fact that the amplitude of the asperities is not controlled during accumulation of the sum n G n . Instead, the amplitude of the roughness after accumulation is regulated by the normalization: where the desired amplitude A can be chosen from the limit of completely flat surface A 0 to values comparable with the characteristic heights of the system under consideration. This approach has been implemented into the model here.
Two basic mechanisms only will be involved in our minimalistic model: mechanical (potential) interaction with the surface and energy dissipation. These mechanisms restrict the scale of the objects which we include into the model from working in dimensionless units. The irregularities, which one can incorporate into the model, have to correspond to the size of the typical curvature, which can be described by the chain formed by the discrete segments. In this context, the surfaces with irregularities much smaller or larger than the size of the discrete segments must be treated as almost flat. That is why we will use the segment length approximately 10 times shorter than the numerically generated surface irregularities. Formal model, based on this assumption, can well explain the results of experimental studies on real surface profiles, depending on the effective "surface tension" determined by the "chain" stiffness.
In contrast to our previous studies (Filippov and Popov, 2007;Gorb and Filippov, 2014) , where the elastic curve of the connected segments represented the real fiber of an attachment pad contacting with a surface by the adhesive force of its terminal segments, now the "chain" has a meaning of the soft boundary separating internal material of the root from the external space. So, its extension, caused by the insert of the additional segments, does not mean the simple elongation of the "fiber" but rather the addition of the 2-dimentional internal area of a root hair cell enveloped by the cell wall. Further, we will return to the question and show how to implement it into the model.
In any case, expanding the boundary mechanically interacts with the relief, repulsing from its bulges and attracting to the valleys. This interaction is described by the force caused by continual potential: Here, the derivatives are supposed to be calculated for each node at its discrete position: x x j , y y j . However, the "continuous" surface is defined on the mesh. Technically, operation Eq. 5 means that we find a cell of the mesh closest to the position x x j , y y j and calculate the derivative Eq. 5 for this cell.
Moving different segments according to the potential interactions, in principle, can result in their overlapping, what is unphysical for the boundary, containing practically incompressible material inside, with its own internal pressure. The repulsion between nodes, which could avoid overlapping between FIGURE 7 | Relaxation of the fractions in contact F r F r (t z ) to the equilibrium values for three different roughness (see Figure 6 for details). Meaning of the colors is the same as above.
Frontiers in Mechanical Engineering | www.frontiersin.org November 2020 | Volume 6 | Article 590894 segments, should be incorporated into the of interactions. For isotropic roots, the easiest way to simulate internal pressure is to use repulsion of all the nodes from the common center. However, the form of the root boundary is very complex. It expands to the plenty of random potential valleys and produces many "local centers," the positions of which are impossible to predict. It is preferable to use a mutual short-range repulsion between all the nodes of the chain: This interaction should be weak and short-range, just enough to prevent the segments from overlapping. For isotropic roots, it automatically reduces to the repulsion from the common center, but for many protrusions, it causes several independent repulsions from the local centers.
Counting on the microscopic nature of the problem, we will use the so-called "over-damped" equations of motion without inertial terms. Accumulating all the forces of the problem, one can now write the equation of motion as follows: Here, index j runs from 1 to current number N s of the segments j [1: N s ]. A formal solution of these equations must be accompanied by the creation of new nodes in random r → k , according to the above-described digital procedure, and to the aging of the already existing ones, which leads to an increase of their rigidity K ⊥ j K ⊥ j (τ j ). One can determine a mean rigidity of the root boundary from the beginning of its growth to the current time moment. To do this, one has to integrate each elastic constant K ⊥ j (τ j ) over time starting from the moment of the appearance of particular segments. It gives an array of the time average elasticity of all the segments: Now, the mean rigidity of the root can be found as the statistical average over the ensemble of the segments: However, both these averaging Eqs. 8 and 9 are linear operations. Their results are proportional to the starting elastic constant K ⊥ j (τ j 0). The fluctuations of the initial elastic constant are small, and one can ignore the difference between K ⊥ j (τ j 0) values of different segments and of K ⊥ j (τ j 0) K ⊥ 0 ≈ const. Therefore, one can characterize different variants of the system by the sole value of K ⊥ .

Root Hair Tip Morphology (Scanning Electron Microscope)
Posidonia oceanica seedlings produced root hairs with a total length ranging between 200 and 1,000 µm and a mean width of 12.31 ± 1.94 µm (mean ± SD, Table 1). The portion of root hairs in contact with the substrate always showed a strongly branched and dactyliform pad-shaped structure (insert of Figures 1B and  2A). These branches were very thin (6.11 ± 0.74 µm), with a width equal on average to about half that of the root hair shank. The mean pad radius was 28.95 ± 4.06 µm.
The contour curvature of the root hair pads attached to the substrate (resin tiles replica of a polishing paper with 12 µm particle size) was 0.1317 ± 0.0317 µm −1 ( Table 1, mean ± SD).
A thin layer of a solidified substance was discovered in SEM images of the root hair pad base and of the resin tile where some of the root hairs had been detached ( Figure 2B). Such an observation was confirmed by TEM analysis (Figures 2C,D).

Root Hair Ultrastructure (Cryo-Scanning Electron Microscope and Transmission Electron Microscope)
The average cell wall thickness, measured along the shaft in fractured root hairs using Cryo-SEM, was 1.50 ± 0.24 µm. In the apical region, the cell wall thickness was lower, at 0.86 ± 0.14 µm. Cross sections of root hairs adhering to the substrate were visualized in TEM, and results confirmed the presence of a glue-like substance filling the empty spaces at the root hair/ substrate interface. The observed substance created a layer of different thickness in all samples tested ( Figures 2C,D). The layer was thin when the root hair was in close contact with the substrate, and thicker when interstices or small gaps were present and the root hair could not penetrate into the gaps between the asperities. Moreover, TEM images allowed an estimating of the actual portion of root hair in real contact with the substrate.

Numerical Simulation
In order to model the behavior of the system with different elasticity levels, we used three representative values of the initial elastic constant: small K ⊥ 01 21.7 N/m, medium K ⊥ 02 76 N/m, and large K ⊥ 03 217 N/m. These values are comparable with the experimental expectations, but the differences between them are one order of magnitude. The numerical modelling approach permitted us to study the possible behavior of the system, even if the values are not observed in real plants.
Typical configurations of the boundary for these three elastic constants -K ⊥ 01 , K ⊥ 02 , and K ⊥ 03 -are shown in Figure 3. It is important to note that each new repetition of the numerical run produces new particular realizations of the surface and resulting structures. Numerous generations of these structures confirm a general qualitative difference between the systems with different surface tensions (corresponds to K ⊥ 0 ), which is already seen from Frontiers in Mechanical Engineering | www.frontiersin.org November 2020 | Volume 6 | Article 590894 Figure 3: The curvature of the boundary increases with an increasing value of K ⊥ 0 . A typical dynamic scenario of the expansion of the root boundary on the substrate is reproduced also in Supplementary Movie S1. It is directly seen how the boundary dynamically reacts on the surface structure, repulsing from the asperities and filling the valleys between them.
The visual difference between the smoothness of the boundaries plotted in the subplots (A-C) must be characterized quantitatively. For this goal one can calculate the curvature of each curve. Local curvature in arbitrary points χ(l) (x(l), y(l)) of parametrically defined curves is defined by the radius R c of the circle touching the curve in this point. An analytical formula for the curvature can be written as follows: where derivatives are calculated along the curve: 3 2 in Eq. 10 can give both positive and negative values of the curvature. For our goal, we need its absolute value. The curvature was calculated as discrete array Curv j .
The curvature varies with time, because the root starts to grow from a smooth circular boundary and develops into a more and more complex structure. We have performed this calculation for all three cases: K ⊥ 01 , K ⊥ 02 , and K ⊥ 03 . As a timedependent integral characteristic, the mean value of the curvature could be: It is convenient for us to represent the value of the curvature radius as R c (t) 1/Curv(t). The resulting curvature radiuses for all three cases, K ⊥ 01 , K ⊥ 02 , and K ⊥ 03 , are shown in subplot (A) of Figure 4 by the blue, green, and red curves, which correspond to the colours used in the subplots (A-C) in Figure 3, respectively. Another important characteristic is the space to which the root hair expands in the current moment. It can be estimated from the mean radius R Area of the boundary calculated starting from the "center of mass" (x c , y c ) of all its segments, where: Individual components of this radius can be calculated as follows: And its mean value is equal to: For the expanding root hair, this radius also depends on time R Area R Area (t). But its increase has a meaning other than the time dependence of curvature Curv(t). The curvature increases because the boundary of the root becomes more complex, but the radius R Area (t) increases simply because the area covered by the root is expanding.
The mean radius (R Area ) of the hair root increases linearly with time for the three modelled typologies of the hair stiffness; this is shown in the subplot (B) of Figure 4. Soft root hairs show the steepest slope, followed by root hairs with intermediate stiffness and then by hard root hairs. The soft and intermediate soft root hairs show a similar trend, with a rapid growth at the beginning of the development and a more constant growth in the sequent phases. At large times, the softest root hair reaches the largest pad size due to greater elasticity. The simulated root hair with intermediate stiffness reaches values similar to the one measured in the real root hair pad observed at SEM.
As we already discussed, the elastic constants evolve independently for different segments of the boundary depending on the individual time elapsed from the moment of their appearance K ⊥ j K ⊥ j (τ j ). Moreover, the number of the segments N s N s (t) increases with time and the array of the individual times τ j becomes longer as well. It makes the visualization of the time evolution of all individual values K ⊥ j K ⊥ j (τ j ) in one static plot difficult. However, one can use the following approach. First of all, we plot, as a set of dots, all the values K ⊥ j K ⊥ j (τ j ) vs. index j. Each new member of the array appears with an arbitrary index inside the array, but it automatically finds its own place in the plot expanding the array. As a result, we get a cloud of points independently filling the range from the minimal value K ⊥ to the maximal value λK ⊥ k . An instant plot of this configuration is shown by colored dots for all three variants in subplot (A) of Figure 5. One can see that the positions of the points are not sorted along the index j. Obviously, they are different for different realizations. There is a well pronounced compaction of the points near λK ⊥ k , where all the points concentrate at the final stage of their motion for each of the three cases: K ⊥ 01 , K ⊥ 02 , and K ⊥ 03 . This process in dynamics is reproduced in Supplementary Movie S2.
To provide a statistical description of the points distribution, one can calculate the histogram of this distribution along the vertical axis. The number of segments N s N s (t) increases with time and the histogram becomes smooth enough to find out the regularity of the process. Histograms for all three cases are plotted in subplot (B) of Figure 5. Perfect correlation between the direct observation of the compaction of the points near final stiffness and calculated maximum of each histogram is obvious. The only difference between the cases is the final value λK ⊥ k , which is proportional to the initial one (K ⊥ ).
What is more important here, is that all three curves have generally the same structure, which confirms the universality of the process independently on the K ⊥ . One can expect that this universality takes place for the whole time interval from t 0 to current t 700 min. Having an array of the histograms of K ⊥ (t) for every time moment t, one can collect them into a twodimensional array and plot them together as a color-map. This Frontiers in Mechanical Engineering | www.frontiersin.org November 2020 | Volume 6 | Article 590894 procedure has been performed by us for all three cases and corresponding maps are presented in the subplots (C-E) of Figure 5. Zero value of the probability is represented here by a deep blue color, while maximum value is represented by deep red. In fact, these maps show in static form (despite the changeable length of the array!) the complete process of evolution of all the values K ⊥ j (t) together. As expected, for all the cases, universality is obvious. The only difference is in the starting value K ⊥ .
However, an important question remains for the validation of the model. As we stressed from the very beginning, the original problem is a complex 3-dimensional problem, and the present study is just an attempt to model it using a 2-dimensional approach. However, we still can advance the model for such a compromised approach. Let us stop the expansion of the root at some (final or intermediate) stage. Its boundary, at this moment, has a certain shape determined by a competition of all the forces involved in the problem. In particular, it is a competition between the repulsion from the surface asperities Z(x, y) n G n and the bending stiffness. As a result, the height Z(x j , y j ), corresponding to the different nodes r → j (x j , y j ), is different. However, the ability of the boundary to bend in Z-direction is also limited. One can suppose that transversal elastic constant in this direction is practically the same as it is for the horizontal plane: K z ≈ K ⊥ . It causes vertical elastic force: In Z-direction, the hair root boundary chain also tends to find a compromise between this force and interaction with the rigid surface in the same manner as it was done in (x, y)-plane. We can calculate this compromise using the relaxation approach. For definiteness, we apply constant vertical attraction force f →z j f →z const < 0 and exponential repulsion from the surface F →z j A (−(zj−Z)/az . Now, by considering all the vertical forces, one can solve the relaxation equation: and find an equilibrium configuration. Here, adaptation time t z runs from t z 0 and stops when the derivative zz j /zt z becomes negligibly small zz j /zt z → 0. We have solved this equation for three different roughness levels (6, 12, and 171 µm asperities size) using the average K ⊥ constant corresponding to the real situation. The results are plotted using three different colors for these cases: red, green, and blue corresponding to high, medium, and low roughness, respectively ( Figure 6). As expected, the result is different for different roughness levels. Because the boundary is a closed curve, one cannot plot corresponding functions along any of the horizontal coordinates (x, y). The only possibility is to plot them along the index of the nodes: j 1, . . . , N s . Counting that, the number of the segments N s varies with the time and the evolution of the root, since the process can be stopped at any arbitrary moment. It is convenient to normalize the horizontal axis in the plots of Figure 6 to this number j/N s and plot in the interval [0, 1]. The substrate surface in all the subplots is shown by the bold line.
In principle, the difference between the cases with different roughness levels is seen directly from Figure 6. However, it is necessary to extract this information quantitatively; it is also essential that this extraction should be done for any adaptation time t z moment, but not for the only fixed one, which was used in Figure 6 to illustrate the final results of the solution of the relaxation equation. We decided to use the fraction of the segments in contact with the substrate as a measure of the root hair boundary configuration difference. In fact, it means that we have to formalize a procedure of finding a fraction of the segments in contact, solve Eq. 16 at every step of t z , and calculate the fractions using this formalized procedure.
Natural formalization is to apply the following condition: if distance between the boundary and surface is less than the threshold (z j − Z) < Δz tresh , the segment is in contact. In the opposite case, it is not. We can consider the Δz tresh as the thickness of the solidified substance. Formally, this condition can be written in the form: Using this condition, the fraction of segments in contact can be easily calculated as the following sum: This procedure is quite efficient and it allows for performing the calculation of the fraction depending on time Fr Fr(t z ). The results for the three different roughness levels (6, 12, and 171 µm asperity size) used for the simulations are presented in Figure 7. Root hairs, due to their specific elasticity and their morphological characteristics, cannot establish good contact with surfaces with too high or too low roughness levels. Optimal contact was produced by the root hair on the average simulated roughness.

DISCUSSION
Root hairs originate from trichoblasts by a highly polarized growth process that occurs at the apex of the cell, referred to as tip growth, and the concomitant repression of proximal growth along the shank (Hirano et al., 2018). Cytoarchitecture of the root hair is highly polar-organized during growth but it disappears in mature cells, when growth usually ceases (Ryan et al., 2001;Ketelaar et al., 2002;Emons and Ketelaar, 2009;Balcerowicz et al., 2015). The root hair's cell wall is a complex network of polysaccharides. The major components are cellulose microfibrils, the polysaccharides matrix made up of hemicelluloses and pectins, and various cell wall proteins (Somssich et al., 2016). A primary cell wall, which is a thin layer that continually expands outwards, is deposited at the root hair tip, while a strong multi-layered secondary cell wall develops inside the primary cell wall, 20-30 μm from the tip, and continues along the shank (Dumais et al., 2004). Once tip growth ends, hair tips also become covered by a secondary cell wall (Galway et al., 1999).
Cellulose microfibrils are randomly distributed in the primary cell wall of root hairs, while in the nonexpanding tubular portion they present an oriented distribution (Galway, 2006). Hemicelluloses and pectins serve as cementing material and any change in their chemical nature influences the plasticity and elasticity of the cell wall. Both the microfibrils arrangement and the chemical properties of polysaccharides in the cell wall create an anisotropic gradient along the cell, which is a matter of considerable importance in the root hair adhesion process.
Our analysis of the morphology and ultrastructure of P. oceanica seedlings showed a unique pad formation process characterized by the root hair tip branching behavior and by the presence of a solidifying (filler) substance at the substrate/ pad interface. These two features, reported here for the first time, together with the presented numerical model simulations allow for better understanding of the adhesion mechanism of P. oceanica described in previous works Badalamenti et al., 2015;Tomasello et al., 2018;Zenone et al., 2020). The numerical model simulated the contact formation of the P. oceanica seedling root hair on hard rough substrates, clarifying what role the physical properties of the root hairs play in pad formation and in the facilitation of the attachment process.
In contrast with most terrestrial plants, where root hairs increase the root surface area to enhance mineral and water absorption, in P. oceanica, and possibly also in other seagrass species, the main role played by root hairs is to provide anchorage to seedlings. However, similar anchorage systems have been described in some terrestrial climbing plants [eg, Ficus pumila (Groot et al., 2003), Parthenocissus tricuspidata (Scherge and Gorb, 2001;Steinbrecher et al., 2010;Steinbrecher et al., 2011), Hedera helix (Melzer et al., 2010, Amphilophium crucigerum (Seidelmann et al., 2012), Passiflora discophora (Bohn et al., 2015), and Syngonium podophyllum (Yang and Deng, 2017)]. To the best of our knowledge, no similar detailed studies exist for most seagrass species, even though the morphology and ultrastructure of root hair have been described for some of them (Cooper and McRoy, 1988;Roberts, 1993;Tomasello et al., 2018).
The terrestrial species P. tricuspidata, P. discophora, and A. crucigerum, build a multi-layered and multicellular pad at the base of tendrils to increase the contact area and produce a firm adhesion to climb on a given surface. Other terrestrial species (e.g., S. podophyllum, H. helix, and F. pumila) climb by means of clusters of adventitious roots which are covered by root hairs that emerge from internodes of the plant stem and form an adhesive pad according to a multistep process (Groot et al., 2003;Melzer et al., 2010;Melzer et al., 2012;Yang and Deng, 2017). In both cases, the pad formed is a complex structure and its adhesion occurs at different spatial scales and in different moments. For these species, root/root hair adhesion occurs according to different steps. Generally, the first step starts with one of the plant adhesive structures that grows towards the substrate. In the second step, the adhesive structure establishes contact with the substrate. In this phase, adhesion is very weak. In the third step, adhesion become stronger as a consequence of mechanical interlocking, and the presence of a glue-like substance may facilitate adhesion strength in some species. The fourth and last step is the hardening of the glue, a fundamental process for bonding and strengthening the attachment. In F. pumila, S. podophyllum, and H. helix, root hairs adhere to the substrate by mechanical interlocking, swelling, or flattening the root hair tip to increase the contact area.
In P. oceanica, the pad structure is produced by the root hair tip branching. According to our results, root hair branching increases the real contact area about 50 times compared to the contact area of a single tubular root hair. The branching of P. oceanica root hairs is a pivotal mechanism for pad formation, and supports the hypothesis that adhesive root hairs may represent an adaptive trait to life in marine rocky habitats characterized by strong turbulent water movements Montefalcone et al., 2016). Interestingly, in terrestrial plants, root hair branching is reported to happen only under stress conditions (Yang et al., 2011;Bobrownyzky, 2015) or in mutant phenotypes (Engstrom et al., 2011).
The role played by the pad in anchoring the root hair is most likely improved by the presence of a filler substance at the interface between the pad base and the substrate. The "glue" thickness observed at the root hair/substrate interface in TEM in this study was not homogeneous. It was thinner where the root hair pad was in close contact with the substrate, and thicker where the contact was looser due to the presence of gaps or interstices. The glue can be seen as a form filler whose role is related to the improvement of root hair adhesion to the substrate by the enhancement of the gap closure at the micro-scale. On rough substrates, root hairs grow between the asperities or inside gaps and small crevices. Under these conditions, root hairs were observed filling up the empty spaces either by the branches, by the compliancy of the cell wall, or by the glue. Such a strategy improves the plant surface-replicating mechanisms. Without the glue presence, attachment would solely be based on mechanical interlocking (between the root hair pad and substrate). It can be hypothesized that the glue may facilitate the first adhesion of the root hair to the substrate and to some extent may reduce the friction of the root hair tip during growth (until the glue is solidified). A low friction environment at the root hair tip region could represent an advantage for the plant to penetrate small interstices and to follow the topography of the substrate, which we clearly see in model simulations.
The elasticity of the root hair tip has strongly affected the effectiveness of formation of the pad in our model simulations performed on a constant micro-roughness (i.e., asperities diameter of 12 µm). When the root hair tip is too soft (low value of the initial elastic constant K) it produces a large and complex pad (i.e., characterized by several branches) but, due to its softness, strong mechanical interlocking can be inefficient. On the contrary, a root hair tip that is too hard (high value of the initial elastic constant) improves mechanical interlocking but, due to its high stiffness, it produces a small and simple pad, ultimately reducing the area in contact with the substrate. The curvature radius (R c , i.e., the inverse of curvature) of the root hair pad increases with increasing elastic constant K ⊥ 0 . Our model showed how root hair radius increases with time, reaching the greatest pad apparent contact surface at lowest stiffness and the smallest at highest stiffness.
In the model, the medium-level elasticity (initial constant elasticity, K ⊥ 02 corresponding to the estimated root hair natural elasticity) proved to be optimal for developing a pad capable of increasing the amount of surface in contact with the substrate and maximizing the mechanical interlocking at the roughness provided. Substrate roughness is a key factor affecting root hair adhesion strength. No general rules define the relation between the roughness and the adhesion strength, as it is highly dependent on the observed system (Budhe et al., 2015). Zenone et al. (2020) reported that the highest adhesion strength in P. oceanica was achieved by seedling root hairs grown at a roughness of 12 µm (asperities diameter), a value roughly corresponding to the root hair diameter. Roughness lower than 3 µm and greater than 52 µm negatively affected the adhesion strength. No adhesion was recorded on a smooth glass replica (supposed no roughness), highlighting the important role played by substrate roughness in the adhesion process of P. oceanica root hairs.
Compared to a smooth substrate, roughness allows root hairs to grow between the interstices, offering paths of minimum growth resistance and inducing an increase in the curvature of the pad which, in turn, contributes to generating a stronger adhesion to the substrate. In the model simulation performed under the estimated natural elasticity constant (K ⊥ 02 ), the highest adhesion (i.e., the number of simulated root hair segments in contact with the substrate surface) was reached at the intermediate surface microroughness level. At this roughness, the fraction of segments in contact with the substrate surface was higher than that observed at lower and higher roughness levels. This high real contact surface is very likely responsible for a stronger adhesion and corroborates the results reported in Zenone et al. (2020). At the high roughness level, the low number of contact segments recorded by the model may limit the mechanical interlocking and, consequently, the adhesion strength. Under this condition, curvature does not cause adhesion enhancement.
The shape and size of the simulated root hair pads properly reflected the real system, shedding light on the role played by pads in the attachment process of P. oceanica seedlings' root hairs to the substrate. Root hair specific elasticity and pad formation seem to be the best combination to produce the optimal contact to a substrate with a roughness level at the scale of the root hair diameter (i.e., 12 µm). It is important to highlight that Δz tresh introduced in the model corresponds to the layer of the filler substance reported in this work. It is very likely that, during the course of evolution, P. oceanica has found a perfect fit between the size of the root hair, the elasticity gradient of the cell wall, the secretion of a glue, and, presumably, the natural roughness of the substrates, to determine through the pad formation a strong anchorage to consolidated substrates and to allow seedlings to colonize and grow in very high energy environments.
Several restoration attempts have been made in the Mediterranean Sea to contrast the general regression caused by human impact on P. oceanica meadows, such as transplanting plugs, sods, or bare-root rhizomes collected from donor plants (Meinesz et al., 1993;Molenaar and Meinesz, 1995;Piazzi et al., 1998;Carannante, 2011;Calvo et al., 2014;Alagna et al., 2019a). These initiatives have raised concerns in relation to their environmental sustainability (Balestri et al., 2019) and also because their large-scale application has rarely achieved high success rates (van Katwijk et al., 2016). To reduce the environmental impact on donor meadows, an approach known as "nursery seagrass" has been proposed which involves the propagation of plants from seeds reared in aquaculture facilities (Balestri and Lardicci, 2012;Balestri and Lardicci, 2014;Balestri et al., 2015). Research has recently focused on finding a suitable and low-impact way to fix germinated plants to the substrate (Alagna et al., 2019b). The adhesion mechanism of P. oceanica seedlings and the model we described here may give a boost to the use of particular consolidated substrates with specific roughness levels for new restoration initiatives. Identifying the right substrate where the seedlings adhere before transplantation could represent an effective approach to avoid the limits highlighted by other experiences (Alagna et al., 2019a). The model developed in this study could be useful in identifying the best characteristics of natural substrates, where the seedlings adhere and obtain a first stable anchorage before releasing them at sea.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
AZ, SG, and FB conceived and designed the research. AZ conducted the experiments and analysis. AF and AK developed the model. SG contributed analytical tools. AZ and AK analyzed data. AZ and AF wrote the manuscript. SG, FB, and AK revised the manuscript. All authors read and approved the final manuscript.

FUNDING
This paper has been partially supported by the project Marine Hazard, PON03PE_00203_1, Italian Ministry of Education, University and Research (MIUR). We acknowledge financial support by DFG within the funding program Open Access Publizieren.

ACKNOWLEDGMENTS
The Authors would like to thank Esther Appel for her collaboration in laboratory work at the Functional Frontiers in Mechanical Engineering | www.frontiersin.org November 2020 | Volume 6 | Article 590894 Morphology and Biomechanics lab in Kiel and Giuseppe Di Stefano for the technical support at IAS-CNR lab in Castellammare del Golfo. This study was developed during a short-term scientific mission granted to AZ by Cost Action ENBA (CA15216).