Bi-directional crosstalk between cells and extracellular matrix leads to network morphogenesis in multi-layered tissues

Cell-generated mechanical forces drive many cellular and tissue-level movements and rearrangements required for the tissue or organ to develop its shape1, 2, 3, 4, 5. The prevalent view of tissue morphogenesis relies on epithelial folding resulting in compressed epithelial monolayers, overlooking the involvement of stroma in morphogenesis1, 4, 6, 7. Here, we report a giant web-like network formation of stromal cells in the epithelium-stroma interface, resulting from a multi-scale mechano-reciprocity between migrating cells and their extracellular environment. In multi-layered tissues, surface wrinkles form by a stromal cell-mediated tensional force exerted at the basement membrane. The topographical cue is transmitted to the stromal cell, directing its protrusion and migration along the wrinkles. This inductive movement of the cells conveys traction forces to its surrounding extracellular matrix, remodeling the local architectures of the stroma. In this manner, stromal cells and wrinkles communicate recursively to generate the cellular network. Our observation provides a rational mechanism for network formation in living tissues and a new understanding of the role of cellular-level tensional force in morphogenesis.

Cell-generated mechanical forces drive many cellular and tissue-level movements and 20 rearrangements required for the tissue or organ to develop its shape 1, 2, 3, 4, 5 . The prevalent view 21 of tissue morphogenesis relies on epithelial folding resulting in compressed epithelial 22 monolayers, overlooking the involvement of stroma in morphogenesis 1,4,6,7 . Here, we report 23 a giant web-like network formation of stromal cells in the epithelium-stroma interface, 24 resulting from a multi-scale mechano-reciprocity between migrating cells and their 25 extracellular environment. In multi-layered tissues, surface wrinkles form by a stromal cell-26 mediated tensional force exerted at the basement membrane. The topographical cue is 27 transmitted to the stromal cell, directing its protrusion and migration along the wrinkles. This 28 inductive movement of the cells conveys traction forces to its surrounding extracellular matrix, 29 remodeling the local architectures of the stroma. In this manner, stromal cells and wrinkles 30 communicate recursively to generate the cellular network. Our observation provides a rational 31 mechanism for network formation in living tissues and a new understanding of the role of 32 cellular-level tensional force in morphogenesis. 33 34 Fibroblasts are the primary constituents of connective tissues, accounting for most stromal 36 cells embedded in the extracellular matrix. The fibroblast-derived contractile force is sufficient 37 to generate the deformation of soft tissue 8 . Such a contractile force-mediated deformation has 38 been documented in several systems, primarily in single-layered tissues 3,9,10,11 , but the 39 morphogenesis in a multi-layered biological system remains elusive. In this study, we 40 investigate how stromal cells lead to tissue morphogenesis forming a cellular network at the 41 epithelium-stroma interfaces (ESIs). 42 First, we visualized the morphology and cellularity of the ESIs in various tissues after whole-43 tissue clearing, including the intestine, dermis, esophagus, and lung tissue (Extended Data Fig. 44 1a, and Supplementary Video 1). We observed an aligned cellular network in a three-45 dimensional en-face view. The stromal cells connected and formed a web-like structure aligned 46 along the direction of grooves (white dotted arrow) (Extended Data Fig. 1a). Unlike previous 47 observations of stromal condensate and follicle patterns at the ESI 4, 7 , the stromal cells were 48 not restricted to the condensate, instead forming a three-dimensional cellular network at the 49 interfaces. 50 Indeed, the contribution of the stroma is distinguishable during tissue morphogenesis. For 51 example, at gestational week 8 12 , the multi-layered tissue surface emerges with fibrous type I 52 collagen stroma development 13 (Extended Data Fig. 1b). The critical elements of the dermal-53 epidermal junction (DEJ) are known to form by 8-10 weeks of gestation 12 , and matured 54 wrinkles in the DEJ are shown in histologies at late gestation compared to that of early stages 14, 55 15, 16 ; similar changes are observed in villi, formed at 9 -10 weeks of gestation 17 . Thus, the 56 emergence of multilayered tissue is a crucial determinant of the surface morphogenesis of ESI. 57 To explore the principles of the stromal cell-mediated morphogenesis, we developed a 58 tissue-equivalent model consisting of the epithelium (Epi), basement membrane (BM; reconstructed by Collagen Vitrigel Membrane (CVM), and stroma (Str; fibroblasts embedded 60 grooves that served as the main path in the cellular network b). We therefore sought to elucidate such a new class of morphogenesis driven by stromal 86

cells. 87
To exclude the contribution of the Epi layer, we conducted bilayered tissue composed of the 88 stromal cell-embedded collagen matrix covered with a thin CVM to act as the BM (BM-Str 89 model). Mouse fibroblasts (NIH 3T3, CellTracker; red) were embedded into the collagen 90 matrix (gray) beneath the CVM (white or green). To observe the onset of wrinkle formation 91 over time, we visualized cells and collagen fibers using fluorescence confocal microscopy and 92 reflection imaging (Fig. 2a,Supplementary Video 6). The cells near the CVM-Str interface 93 extended pseudopods, forming wrinkles in the early period (1.5h) (Fig. 2a). At the initial stage 94 of wrinkle onset, the cells made a proteolytic path while migrating, thereby generating 95 irreversible deformation in the 3D matrix ( Fig. 2a inset, Supplementary Video 7). The 96 neighboring cells then moved around the proteolytic path (Extended Data Fig. 4). 97 At the wrinkle onset, the stromal cells showed polarization of F-actin and phosphorylated 98 myosin II light chain 2 (pMLC2) in the lamellipodia or tips (Fig. 2b). F-actin was concentrated 99 in each tip (white star), whereas MLC2 was highly phosphorylated between the tips (white 100 arrowhead). The cells are attached directly to the CVM interface by focal adhesion molecules, 101 as revealed by paxillin staining. When the myosin inhibitor and blebbistatin (Bleb) were used 102 to block cell contraction before the formation of wrinkles, no wrinkles were produced. By 103 contrast, the contraction inhibition following the development of wrinkles did not affect the 104 pre-formed wrinkles due to the plastic nature of the CVM (Extended Data Fig. 5a). It indicates 105 that the actomyosin-mediated contractility of stromal cells is a key driving factor to form the 106 wrinkle (Fig. 2c). 107 Using a fibronectin-Forster resonance energy transfer (Fn-FRET) mechanical strain sensor, 108 we examined the strain with high spatiotemporal resolution (Extended Data Fig. 5b,c). The FRET ratio (Ia/Id) decreased (shown by the green pseudocolor) in the regions with high strain. 110 The wrinkle propagated biaxially around the entire interface, leading to changes in the FRET 111 ratio. (Extended Data Fig. 5d, e). This observation was validated with an in-silico 3D Agent-112 Based Model (ABM) comprising contractile cells, fiber matrix, and membrane on the top (Fig  113   2g, Methods). Membrane deformation in silico was quantified and showed noticeable wrinkle 114 formation across two cells (Extended Data Fig. 5f). Particle image velocimetry (PIV) technique 115 was adopted to visualize the surface strain during wrinkling (Extended Data Fig. 5g, h). The 116 displacement map presents the strain field of the surface around the wrinkles. The network 117 pattern defined by the grooves (bulging part in en-face view) exhibited relatively high strain. 118 Given the fact that cell attachment modifies the local strain and propagates the contractile 119 force across the matrix, cell-to-cell distances influence the wrinkle onset (Fig. 2d,120 Supplementary Video 8,9). Cells with a radial distance, r, < 70 μm, formed a condensate 121 instantly. However, cells within a range between 70 and 90 μm migrated toward each other, 122 forming condensates. With r > 270 μm, the cells did not interact directly but instead formed a 123 long groove that connects two cells (Fig. 2e, f). The principle is validated in silico by the ABM 124 ( Fig. 2g). As cells were initially located closer, they migrated toward each other while rested 125 stagnant or even migrated away when placed further away (Fig. 2h). Our results suggest that 126 nearby cells are mechanically coupled in a distance-dependent manner and can promote the 127 directional movement of cells along their axis of interaction, forming in topographically 128 confined condensates. 129 Active cell migration was initiated from the confined condensates within grooves. To 130 migrate, a cell extends protrusions, such as lamellipodia and filopodia, along the wrinkles in 131 response to the topography ( Fig.3a. 8h, Extended Data Fig. 6a). The distant condensates were 132 spread and connected as cell migration continued along the wrinkles (Fig. 3a. 24-32h. white 133 star, Supplementary Video 10). Consequently, most cells were interconnected and aligned in parallel to the wrinkles (Fig. 3b). The aligned cells (lower circularity) within the wrinkle had 135 polarized focal adhesions (Fig. 3c), implying that they generated more traction force than non-136 stretched cells far from the wrinkles (higher circularity) 18 . 137 During cell migration, the initial wrinkles were reorganized and evolved into a network. (Fig.  138 3d, Extended Data Fig. 6b). The two separated wrinkles were joined (black and white star; 139 joining). At the junction, a new perpendicular wrinkle formed (black and white arrowhead; 140 branching). The phenomena were also validated in silico with cell migration (Fig. 3e). In 141 joining, the front cells in two neighboring wrinkles ( Fig. 3f. i and ii) continued to connect. 142 Combining the two stresses resulted in a wrinkling rearrangement (black dotted arrow) with an 143 acute angle (Fig. 3f,Supplementary Video 11). In branching, the migrated cells tugged the 144 membrane in parallel to the wrinkle, forming wrinkles perpendicular to the initial wrinkle ( Fig.  145 3g, Supplementary Video 12). These observations indicate that stromal cells form an interstitial 146 cellular network through repetitive interactions between cell migration and wrinkling. 147 It is noteworthy that the observed network morphogenesis resembles a hierarchical network 148 of folds in elastic membranes subjected to biaxial compressive stress 19 in terms of visual 149 appearance and dynamics. In the elastic membrane system, each fold reorganizes the local 150 stress field over time, resulting in a repetitive wrinkle-to-fold transition. Likewise, in our 151 biological model, individual cells generate a force that modifies the stress distribution in spatial 152 and temporal domains, forming a cellular network. It suggests that spatiotemporal force 153 localization and stress relaxation through wrinkling may play a universal role in the formation 154 of the network in layered biological and material systems. 155 Next, we examined how the cells respond to the dynamic remodeling of the interstitial space. 156 We found that the collagen fibers, filopodia (F-actin) of cells, and wrinkles of CVM were 157 aligned in the same orientation (Fig. 4a). This observation accorded with the results of in-vivo 158 analyses. In the dermis, for example, collagen fibers were aligned through the wrinkles under the epithelium interface (Extended Data Fig. 7a). To understand the causality among the 160 aligned components (wrinkle, pseudopod, collagen fiber), we visualized the dynamics of 161 matrix remodeling during wrinkle propagation.  The stromal cell was mixed into col I gel premix in 2x10 5 /ml (for human cells, 1x10 5 was 290 enough). The col I gel premix was made by 2mg/ml rat-tail col I, 43x 1M NaOH to col I, 10x 291 PBS, and the cell containing full DMEM. 292 The excessive TR2 was removed by pipetting, then the hollow cylinder-shaped PDMS mold 293 coated with Pluronic F-127 (P2443, 0.5% in 70% EtOH) was fitted into the inner hole of the 294 donut-shaped CVM mold. 85μl of the col I gel premix with the stromal cell was poured into 295 the gel mold (※ should be done as fast as possible to prevent evaporation). The mold was 296 incubated for 2h in an incubator for gelation and adherence. After incubation, the gel mold 297 was mildly pressed by a tweezer and the outer CVM mold was removed by striping out. 298 Pouring DMEM would mildly detach the CVM+gel tissue (in gel mold) from the parafilm 299 surface. 300 To make the full Epi-BM-Str model, the epithelial cell was seeded in 1x10 5 /well directly on 301 the CVM of BM-Str model and incubated at 37℃ for 15 min to adhere. After checking the 302 adherence, the whole tissue was transferred into 35pi filled with full media (for the 303 developmental dermis model, add 10ng/ml TGFβ2 (B266458)). The transferred tissue was 304 inverted and agitated smoothly to remove the mold. The free-floating tissue was incubated for 305 the desired time to make self-assembled wrinkles. Mice used in this study were maintained in a specific pathogen-free facility of KAIST 331 Laboratory Animal Resource Center. C57BL/6J Mice, ages 6-8 weeks, were used in this 332 study. After anesthetization with isoflurane, mice were transcardially washed with ice-cold 1× phosphate-buffered saline (PBS) supplemented with heparin with a concentration of 10 334 units/mL, followed by perfusion with ice-cold 4% paraformaldehyde (PFA) in 1× PBS. 335 Dermal tissue was harvested and post-fixed with 4% PFA in 1× PBS at 4°C for 6-12 hours. 336 Dermis, intestine, and ovary tissues were treated with SHIELD processing as described 337 previously 27 . Briefly, the heart of a mouse was exposed to perfuse PBS. Then, the SHIELD 338 perfusion solution (20ml for 4min) is applied. After perfusion, each target organ was 339 extracted and incubated in SHIELD perfusion solution for 2 days at 4℃. The organs were 340 incubated in SHELD-OFF solution for 1 day at 4℃. The organs were then incubated in 341 SHELD-ON solution for 1 day at 37℃ (epoxy polymerization). After the SHIELD process, 342 the organs were incubated in an SDS clearing buffer for 3~5 days at 47℃. The SDS was 343 washed out by 1% PBSTN for 1 day at 37℃. Finally, The PBSTN was washed out by PBST 344 for 12h at 37℃. The whole organs were treated with an optical clearing solution for 12~24h 345 for transparency. 346 For expansion microscopy, tissues were first blocked and permeabilized with a blocking 347 buffer (5% NDS (w/w) 0.2% Triton X-100, 1× PBS) for 3 hours. Then the tissues were 348 incubated with the primary antibody in blocking buffer for 2 weeks at 4 °C, and washed in 349 0.1% PBST 3 times for 30 min each time. Then they were incubated with secondary antibody 350 in blocking buffer for 1 week at 4 °C, and washed in 0.1% PBST 3 times for 30 min each 351 time. AcX was then treated in 1× PBS with a concentration of 0.1mg/mL 6-12 hours at 4°C. 352 The samples were washed in 0.1% PBST 3 times for 30 min each time with gentle shaking. 353 Then the tissues were subjected to a gelation process. In this process, (2.5% (w/w) 354 acrylamide, 7.5% (w/w) sodium acrylate, 0.15% (w/w) N,N′-methylenebis ImageJ. Z-stack projection is done with summation or max intensity projection depending on 380 the data requirement. The raw image was 16-bit, but the processed images were stored in 8-381 bit or RGB for compatibility with documentation. 382

Calculation of the Alignment Percentage 384
The alignment percentage of subjects A and B were acquired by the difference of orientation 385 between them: DA-B. Then the alignment percentage was: 1 180 ⁄ 100 %. The fabricated FN-DA was simply mixed with collagen pre-mixture in 50μg/ml -1 (5μg/ml -1 of 422 FN-DA and 45μg/ml -1 of pure FN) to form FN-FRET CVM. The mixture was incubated for 2 423 hours at 37℃ for gelation, then 1 week at 18℃ for self-assembly (vitrification). 424 To measure the relative FRET analysis, two image sets were acquired in confocal 425 microscopy. Two photomultiplier tubes (PMTs) were used to acquire the donor and acceptor 426 signals. The first set of images was donor to donor (DD). The sample was excited with the 427 488 laser and PMT received 515-525nm filtered wavelength for the donor channel. The 428 second set of images was donor to acceptor (DA). The sample was excited with the 546 laser 429 and PMT received 567-577nm filtered wavelength for the acceptor channel. where , and , are the length and diameter of segments located between endpoints in 482 equilibrium, respectively, and is the viscosity of a surrounding medium. The positions of all 483 points evolve every time step using the Euler integration scheme: 484 includes extensional forces to regulate equilibrium lengths, bending forces to maintain 486 equilibrium bending angles, and repulsive forces to account for volume-exclusion effects. 487 Extensional, bending, and repulsive forces are determined by the following potential equations: 488 where and are extensional and bending stiffness, respectively. and are instantaneous 492 length and bending angle, respectively, and the subscript 0 indicates their equilibrium values. 493 represents the strength of repulsive force, is the shortest distance between two 494 neighboring elements, and , is a threshold distance below which the repulsive force starts 495 acting. 496 497

Fibers and cross-linkers in the matrix 498
Each matrix fiber consists of cylindrical segments connected in series. The length of 499 each segment and an angle formed by two interconnected segments are regulated by Eqs. 5, 6 500 with extensional stiffness ( , = 0.277 N/m) and bending stiffness ( , = 8.12 × 10 -20 Nꞏm), 501 and their equilibrium values are , = 1.0 µm and , = 0 rad. Each cross-linker is comprised of two segments. The length of each segment and a bending angle between the two segments 503 are regulated by Eqs. 5, 6 with the extensional stiffness ( , = 2.0 × 10 -3 N/m) and bending 504 stiffness ( , = 1.04 × 10 -19 Nꞏm), and their equilibrium values are , = 230 nm and , = 0 505 rad. The repulsive force acts only between fiber segments. The shortest distance between two 506 neighboring cylindrical segments is calculated as , and , is equal to the diameter of the 507 segments, , = 60 nm. The repulsive force calculated via Eq. 7 with , = 1.69×10 -4 N/m is 508 distributed to two endpoints of each segment as in our previous study 34 . 509 510

Membrane and its interaction with the matrix 511
The length of chains between node points in the triangulated mesh is regulated by Eq. 512 5 with the extensional stiffness ( , = 0.001 N/m), and its equilibrium value is , = 1.0 µm. 513 A dihedral angle formed by adjacent triangles is regulated by Eq. 6 with the bending stiffness 514 ( , = 1.0 × 10 -19 Nꞏm), and its equilibrium value is , = 0 rad. The node points can be 515 permanently linked to the endpoints of fiber segments if they are located within 380 nm. The 516 repulsive force acts between the matrix and the membrane by considering the shortest distance 517 between fiber segments and triangular elements in the mesh, . , is equal to the sum of 518 the radius of the fiber segments, , /2, and the half of membrane thickness, , /2. The 519 repulsive force calculated via Eq. 7 and , = 4×10 -4 N/m is distributed to two endpoints of 520 the fiber segment and three vertices of the triangular element. 521 522

Cells and their interaction with the matrix 523
Cells are modeled as spheres, and they can co-exist in the same position with matrix 524 elements (i.e., phantom cells without physical volume) (Fig. 2d, left). Thus, the repulsive force 525 between cells and matrix fibers is not considered. Each cell can mechanically interact with matrix fibers if a distance between the cell center and the endpoint of a fiber segment falls 527 within a range between and , where = 2 μm is a cell radius and dI = 4 μm is an 528 interaction distance. All fiber endpoints located within the interaction range can be physically 529 linked to a cell at a constant rate, k+,FC = 0.03 s -1 . The cell exerts a constant force (100 pN) to 530 fiber endpoints linked to it toward the cell center, and the cell also experiences a reaction force 531 with the same magnitude in an opposite direction. The physical links can be dissociated in a 532 force-dependent manner 35 where , is a force-dependent dissociation rate, , is a zero-force dissociation rate 535 constant, , is a spring force vector on the physically-linked fiber endpoint, and , is force 536 sensitivity. Note that it is assumed that physical links behave as a catch bond, so physical links 537 formed on fiber endpoints with high tensile forces last for longer time. The sum of reaction 538 forces acting on each cell is calculated, and the position of its center point is updated using Eq. 539 1. This leads to cell migration toward high-force or stiffer regions. 540 541

Domain and its interactions with elements 542
A rectangular domain (120 × 120 × 60 µm) is employed for all simulations. The 543 boundaries of the domain keep elements within the domain by exerting repulsive forces to 544 elements that attempt to cross the boundaries. It is assumed that the node points of of the 545 membrane on edges and the endpoints of fiber segments close to the boundaries of the domain 546 can be permanently linked to the boundaries. On each boundary, the sum of forces acting on 547 the membrane node points and the fiber endpoints is calculated, then divided by the area of the 548 boundary to calculate stress. By assuming that each boundary represents an elastic material 549 with imposed Young's modulus, the position of boundaries is updated over time; strain is calculated by dividing stress by Young's modulus, and then the displacement of the boundary 551 by multiplying the thickness of the elastic material to the calculated strain. The boundary 552 gradually moves toward the displacement level. A rationale for considering the boundaries an 553 elastic material is that a system in our computational model represents a small portion of the in 554 vitro experimental system, so the surrounding elastic material in simulations represents the rest 555 of the real large system. 556 557

Simulation procedures 558
For the first 100 s, a simulation system is created via self-assembly of fibers and cross-559 linkers and the interaction between the matrix and the membrane. The membrane is initially 560 flat and located at z = 30 µm, and the matrix is assembled between z = 0 µm and the membrane. 561 Cells are allocated in predetermined positions in a symmetric manner relative to the center of 562 the domain; their y and z positions are identical, whereas their x position are 60 µm -dC/2 and 563 60 µm + dC/2, where dC is an initial distance between two cells. As the value of dC, either 20, 564 30, 40, 50, or 60 µm is used. After creating the system, cells were allowed to interact with 565 surrounding matrix fibers and migrate for 200 sec. The rest of parameters and theirs values are 566 summarized in Extended Data Table 1. 567 568

Quantification of fiber alignment 569
To assess the alignment of fiber segments, we first measured alignment for every grid, 570 defined by 〈cos 〉 , where is the angle between cell-cell orientation and each fiber's 571 direction and angle bracket averages all angle deviations within a particular grid. As fiber 572 segments are aligned toward cell-cell orientation, alignment gets 1. Furthermore, regardless of 573 cell-cell orientation, order parameter was considered, defined as 〈2cos 1 〉, where 574 indicates the angle between the mean fiber orientation (of a particular grid) and each fiber's 575 direction. Order parameter approaches 1 as fibers themselves are aligned in any direction.