A mechanogenetic role for the actomyosin complex in branching morphogenesis of epithelial organs

The actin-myosin complex plays crucial roles in supporting various biological phenomena. The actomyosin network is first and most widely known as a functional unit for muscle contraction. In non-muscle cell types, this cytoskeletal assembly generates a mechanical force that primarily triggers diverse cellular functions, such as cell migration, mitosis and shape changes (Gorfinkiel and Blanchard, 2011). At the tissue and organ level, the actomyosin complex in the cell cortex mainly supports epithelial architectures, by balancing the forces generated by E-cadherin and the integrin family, and, thus, maintains tissue integrity and structural homeostasis in diverse epithelial organs (Guillot and Lecuit, 2013). The physical balance can be interrupted by a signaling input that induces morphogenesis, which is an exquisitely designed genetic program with complex spatiotemporal patterns. Under such circumstances, myosin acts as a force generator by interacting with cortical actin fibers; this contributes to the translation of the different types of forces into macroscopic changes in epithelial layers, resulting in organ deformation (Guillot and Lecuit, 2013; Maruthamuthu et al., 2010).

Besides the physical functions, the actomyosin complex is involved in mechanotransduction, in which extra/intercellular physical patterns are decoded into the perturbation of intracellular signaling pathways. In this process, the Yes-associated protein (YAP) and transcriptional co-activator with PDZ-binding motif (TAZ) open a new avenue for the interpretation of the mechanical-to-genetic changes. The stiffness of the extracellular matrix (ECM) and the shape changes in mammary epithelial cells regulate YAP and TAZ activity through a mechanism independent of the Hippo pathway, and this response determines the proliferation-differentiation switch of the cells (Dupont et al., 2011). Moreover, mechanically regulated YAP and TAZ activity was found to be pivotal for determining the differentiation fate of mesenchymal stem cells, ranging from osteocytes to adipocytes, according to the surrounding environment (Dupont et al., 2011; Engler et al., 2006). Upon both biological and physical stimulation, YAP and TAZ regulate their transcriptional activities through a nucleus-cytoplasm shuttling process. Nuclear YAP and TAZ proteins make a functional complex with diverse binding partners, especially with the TEA domain (TEAD) transcription factors, to establish core DNA-binding platforms (Vassilev et al., 2001). Such a versatile combination of complexes enables the dissemination of the genetic effect into broad biological contexts. The primordial cells in epithelial organs also show spatially patterned actomyosin foci, indicating that dynamic forces fluctuate throughout the developing structures (He et al., 2010; Rauzi et al., 2010). However, the exact role of the physical patterns and their relationship to the YAP and TAZ activity in the transcriptional changes during morphogenesis are not fully understood.

Branching morphogenesis is a crucial event for the morphological patterning that occurs in the initial process of epithelial organ development. The repetitive branching process of epithelial buds establishes the complex architectures that guarantee an efficient exchange of various molecules (e.g. nutrients and gases) with the mesenchymal components (Affolter et al., 2003).

The epithelial organs – the lung, kidney, pancreas, mammary gland and salivary gland – develop their own morphologies, although they share cognate extra/intracellular components for the genetic program (Lu and Werb, 2008). One plausible explanation for the morphological differences is the spatiotemporal differences of the components among the organs, which in turn guide the different epithelial architectures by regulating the cellular properties, such as the cell polarity, mitotic axis and location (Tang et al., 2011). In branching morphogenesis, actomyosin contractility is known as an important physical factor (Kim and Nelson, 2012). Distinct patterns have been reported for the actin cytoskeleton in developing epithelial buds according to the organ type (Ewald et al., 2008; Kim et al., 2015; Meyer et al., 2004), indicating that a locally concentrated force was established during the process. Moreover, perturbation of the actomyosin contractility results in impaired morphogenesis in diverse organs (Daley et al., 2009; Moore et al., 2005). Given evidence strongly suggests the importance of the mechanical factors in organ morphogenesis, an investigation of the physical effect on the developmental process could provide key information for our understanding of the core principle of organ shape patterning.

Here, we report the physical and genetic role of the actomyosin complex in the initial process of branching morphogenesis using an ex vivo epithelial organ explant culture system. Through the use of versatile techniques for the live imaging of developing organs, we found that the physical tension arranged by characteristic actin structures is essential for sculpturing the epithelial branching patterns. Moreover, the intercellular forces induce a YAP and TAZ-mediated mechanotransduction pathway, which maintains genes that are crucial for the expansion of epithelial progenitor cells.

To investigate the role of the actomyosin complex in branching morphogenesis, we first screened filamentous actin (F-actin) patterns in diverse organotypic cultures undergoing active branching. Phalloidin staining revealed locally concentrated patterns in the distal or peripheral area of epithelial buds, but the range and intensity were varied depending on bud shape (Fig. S1A). Broad and high signal patterns of F-actin were detected in organs showing lobular form (e.g. pancreas and salivary glands), whereas the high signal patterns were concentrated on the narrow area of epithelial curvature in tube-forming organs (e.g. kidney and lung) (Fig. S1A-C). There was no clear spatial difference of staining pattern in the organ with minimal architectural changes (e.g. tooth germ; Fig. S1A-C). This tendency highlights the potential importance of the spatial distribution of cytoskeleton in determining morphological characteristics. We then used an isolated culture model of the epithelial rudiment of mouse embryonic submandibular gland (epiSMG) (Nogawa and Takahashi, 1991), in which visualization and perturbation of the signaling activities during developmental processes is feasible. epiSMG cultures showed concentrated F-actin patterns on the cell cortices in the outermost epithelial layers when compared with the even signal distribution of cell membrane tracer (DiD) (Fig. 1A). The peripheral actin signals of epithelial buds were significantly higher compared with those in the central and pericleft areas (P=8.04×10⁻¹¹ and 1.34×10⁻¹⁰; Fig. 1B). The formation of peripheral actin patterns was apparent from embryonic day (E) 13, coincident with the initiation of cleft and bud formation (Fig. S1D). The concentrated actin patterns were more densely assembled and were maintained during the active morphogenic stage of epithelial buds (E14-15) (Fig. S1D). The functionality of the actin network was demonstrated by the colocalization of F-actin with phosphorylated myosin (pMyosin) (Fig. S1E,F). The peripherally concentrated F-actin and pMyosin patterns were additionally confirmed in the sectioned SMG samples (Fig. S1G) to avoid any biases from optical and chemical penetration problems. Moreover, we could not observe clear apico-basal difference in pMyosin signals in peripheral epithelial cells, which induces spontaneous curvature formation of epithelial organs (Fig. S1H). These data indicate that peripheral actomyosin complexes establish spatially distinct intercellular force distribution in developing epithelial buds.

We next treated SMG cultures with blebbistatin (a myosin ATPase inhibitor) to disrupt actomyosin contractility in epithelial cells and monitored the morphological changes. Blebbistatin-treated SMGs and epiSMG cultures showed a remarkable architectural distortion of the epithelial buds, reminiscent of multiple bud formation (Fig. 1C). Through quantification of morphological parameters, we discovered that three points characterize the morphological distortion: shallower cleft depth [89.55±9.35 µm (control) versus 33.42±2.74 µm (blebbistatin); P=4.38×10⁻⁶], increased bud numbers (5.71±0.45 versus 37.00±1.92; P=5.85×10⁻⁷) and enlarged organ size (6.36±0.60×10⁴ µm² versus 8.35±0.42×10⁴ µm²; P=0.01) (Fig. 1D). The morphological alteration by blebbistatin showed dose-dependent patterns (Fig. 1E,F). Inhibiting actin polymerization with a low concentration of cytochalasin D treatment (100 nM) also resulted in similar morphological changes in the cultures (Fig. S1I).

Next, we investigated the epithelial morphology of developing buds at the single-cell level, by optically sectioning images of the midplane of each bud, showing maximum diameter (Fig. S2A). In the control group, epithelial cells in the peripheral sheet presented columnar features with highly polarized cells (aspect ratio=2.194±0.09; Fig. 2A-C and Fig. S2B). The relaxing of intercellular tension by blebbistatin transformed the cell morphology to a cuboidal shape (aspect ratio=1.65±0.09; Fig. 2C) and increased the cell area [37.19±3.99 µm² (control) versus 52.21±2.62 µm² (blebbistatin); Fig. 2D], but intrinsic cell polarity (marked by Golgi complex localization) was not affected (Fig. S2G). The typical columnar morphology of the peripheral cells is formed by an apico-basal contraction force and intercellular junctional integrity, which alter cell shape with a high aspect ratio (Hannezo et al., 2014). This background could provide a plausible explanation for the morphological changes: blebbistatin treatment impeded circumferential contractility in the cell cortex, but did not affect intercellular junctions (Fig. S2C), resulting in cuboidal-shaped cells. These results were confirmed through genetic perturbation of myosin function by expressing nonphosphorylatable mutants of the myosin regulatory light chain (MRLC-AA) (Vicente-Manzanares et al., 2008) in epiSMG cultures. MRLC-AA-expressing cells presented a lower aspect ratio (1.83±0.07) and increased cell area compared with controls (without expression) and EGFP-expressing cells (Fig. 2E-G). We additionally quantified the aspect ratio of the peripheral and central cells in both experimental setups for actomyosin perturbation [control (−BB) versus blebbistatin treatment (+BB) and MLC-AA− versus MLC-AA+ cells]. In the control groups (−BB and MLC-AA−), peripheral cells presented a significantly higher aspect ratio than central cells, whereas the discrepancies of the aspect ratio between peripheral and central cells were clearly reduced in actomyosin-perturbated groups (+BB and MLC-AA+) (Fig. S2D). Moreover, MLC-AA-expressing (MLC-AA+) cells and blebbistatin-treated cells (+BB) showed similar morphological changes compared with each control group – the significant decrease in the aspect ratio occurred in peripheral cells, but not in central cells (Fig. S2E). Moreover, the aspect ratio of peripheral cells was significantly higher compared with cells in the central and pericleft regions (Fig. S2F), which completely matched the F-actin density patterns (Fig. 1B). Collectively, the single-cell geometric data indicate that peripheral epithelial layers experience concentrated force generated by circumferential actomyosin contractility and that cleft formation results in de-stressing of the concentrated force. It is noteworthy that there is a limitation with respect to interpretation of the results, owing to absence of the data directly measuring intracellular forces; therefore, we could not conclude whether this force is applied actively or passively (from external factors). Meanwhile, low-concentration cytochalasin D also altered cell shape but did not induce a clear expansion of the cell area (Fig. 2C,D), which might be due to insufficient concentration for complete inhibition of the cortical actin assembly. Golgi complex polarity was sensitively disorganized by an equal concentration of cytochalasin D (Fig. S2H), implying the importance of actin structures for intracellular polarization as a binding platform of diverse polarity proteins (Li and Gundersen, 2008). Moreover, we investigated the pattern of fibronectin (FN), which is a crucial ECM component in cleft formation (Sakai et al., 2003). In SMG culture, overall FN signaling intensity around epithelial buds and clefts did not change upon blebbistatin treatment (Fig. S3A,B). Although blebbistatin-treated SMGs displayed the relaxed pattern of FN, the concentrated signals in the cleft region were not changed (Fig. S3A). In epiSMG cultures, there was no clear FN pattern around the outer epithelial layers in both groups (Fig. S3C). These data indicate that morphological changes of epithelial buds under actomyosin perturbation is mainly mediated by the alteration of internal stress in epithelial buds.

How is a highly stressed epithelial sheet involved in patterning an epithelial architecture? To address this question, we first monitored the growth patterns of epithelial buds through time-lapse live imaging. While control epiSMGs showed restricted expansion of epithelial buds, rapid outgrowth was observed in blebbistatin-treated epiSMGs [outgrowth velocity: 8.33 µm/h (control) versus 25.13 µm/h (blebbistatin)] (Fig. 3A,B; Movies 1 and 2). The peripheral epithelial margin of the latter group was radially expanded without directional tendency, resulting in the bleb-like irregular shape (Fig. 3A,B; Movie 2). The irregular outgrowth pattern was more apparent in a salisphere culture model using an SMG acinar cell line (SMG-C6) [Fig. S4A-D; circularity: 0.82±0.01 (control) versus 0.44±0.02 (blebbistatin)]. Cell division profoundly affects the organization and size control of epithelial tissues (Gibson et al., 2006), but we observed no significant difference in proliferation rate between the groups (Fig. S4E,F). Based on these results, we found that actomyosin contractility in peripheral layers plays a key role in the restriction of outward growth; furthermore, we hypothesized that the stress generated by increased cell density might be resolved by transmitting the force to another direction. To confirm this conjecture, the spatial force distribution around the surface of epithelial buds was measured via a three-dimensional (3D) force distribution analysis using fluorescent microbeads (Fig. 3C). We quantitatively mapped the time-lapse 3D displacement field in the vicinity of the epithelial layer for 8 h (Fig. 3D-F; see Materials and Methods). First, we tracked the individual microbeads throughout the whole experimental duration. Interestingly, lateral flow of the microbeads alongside the epithelial layer was detected in a bidirectional manner, and turbulence-like flow was generated towards the inside of the epithelial bud at the point with converging displacement vectors (Fig. 3D; Movie 3). The points of concentrated inward forces were perfectly matched to cleft-propagating positions (Fig. 3D and Fig. S4G). In contrast, blebbistatin-treated epiSMGs showed an even force distribution around the epithelial layer (Fig. 3D; Movie 4). Next, we quantified the full spatiotemporal dynamics of the displacement fields in 3D (see Materials and Methods). Visualizing the 3D distribution of the tangential traction revealed a highly organized and segregated spatial pattern (Fig. 3E). Strikingly, the temporal dynamics of this force distribution predicted cleft formation. After a pair of local and transient tractions are applied inwards, a gradual and prolonged surface retraction is induced, eventually forming a cleft (Fig. 3E). Meanwhile, blebbistatin treatment completely abolished the spatiotemporal organization of the force distribution and resulted in diminished cleft formation (Fig. 3F). Decomposing the measured displacement vectors identified converging tangential traction alongside the epithelial layer at the cleft locations (Fig. 3G). The tracking patterns demonstrated that actomyosin complexes efficiently transmit the internal stress in the epithelial buds to tangential forces, and continuous cell proliferation acts as a major driving force of cleft propagation (Fig. 3H). This effect could provide compelling evidence for a buckling-infolding model of cleft formation (Kim et al., 2018).

Cleft formation and propagation are crucial events for the compartmentalization of the functional units of salivary glands, i.e. the acinus and duct (Harunaga et al., 2011). Together with the halted cleft propagation, we found that salivary glands (SMGs and sublingual glands) treated with blebbistatin presented an immature and undistinguishable ductal structure (Fig. 4A). Notably, tubulogenesis was blocked by blebbistatin treatment, as illustrated by impaired midline condensation of zonula occludens-1 (ZO-1), the key step of ductal lumen formation (Nedvetsky et al., 2014) (Fig. 4A). Compared with the control group, blebbistatin-treated SMG cultures featured a more obtuse transition angle between the acinar and ductal borders [116.40±5.22° (control) versus 140.41±3.12° (blebbistatin)] with a wider ductal diameter [86.35±9.86 µm (control) versus 148.76±25.46 µm (blebbistatin)] (Fig. 4B). We applied an additional model, in which actomyosin contractility is modulated by adjusting ECM concentration around epiSMG cultures (Paszek et al., 2005) (see Materials and Methods). Low-concentration matrix reduced myosin phosphorylation levels (Fig. 4C,D) and elicited immature morphogenesis, undefined cleft formation and impaired ductal structures (Fig. 4E), similar to blebbistatin-treated SMG cultures (Fig. 4A). We verified the differentiation state of SMG cultures by investigating the pattern changes of proximal (K19⁺)/distal (K14⁺) progenitor cells upon actomyosin perturbation. Strikingly, blebbistatin-treated epithelial buds showed less condensed patterns of K14⁺ progenitor cells in the distal region of the acinus (Fig. 4F-H). Moreover, K19⁺ progenitor cells showed diffused patterns from proximal to distal region, resulting in a mixed population of K14⁺ and K19⁺ cells (Fig. 4F-H).

Different modes of force (e.g. compressive, tensile and shear) have been known to elicit transcriptional changes in diverse genes, especially in cell types exposed to dynamic stress fields (Diamond et al., 1990; Hatton et al., 2003; Rysa et al., 2018). Although there is a possibility that blebbstatin treatment physically perturbed the ductal formation, we hypothesized that actomyosin contractility plays a significant role in the regulation of genetic programs for acinus-duct differentiation of salivary glands. Therefore, we profiled transcriptional changes of epiSMGs upon actomyosin perturbation (blebbistatin treatment) by conducting a whole-transcriptome analysis via RNA-sequencing (RNA-seq) (Fig. 5A). RNA-seq revealed that blebbistatin treatment significantly altered the expression of 150 genes (greater than 1.5-fold change, P<0.05), including 53 up- and 97 downregulated genes. Among the 150 differentially expressed genes (DEGs), multiple transcripts were related to developmental processes (Tgfb1, Tgfb3, Wnt4 and Wnt7b) (Cadigan and Nusse, 1997; Nishita et al., 2000) and ECM-cell interaction (Col2a1, Col9a2 and Itgb3) (Table S1). Gene ontology (GO) analysis of biological processes showed a number of annotations containing the categories of morphogenesis, differentiation and development with high significance (P<0.001) (Fig. S5A). The GO analysis and DEG results support the fact that diverse genes were dynamically regulated upon perturbation of actomyosin function. Interestingly, volume plots of the expression values revealed a striking decrease in Kit expression (Fig. 5B), which is a crucial mediator of epithelial stem cell differentiation and morphogenesis (Leong et al., 2008). Kit was also listed in the multiple gene set of highly enriched GO annotations (e.g. anatomical structure morphogenesis) and in the phosphoinositide 3-kinase pathway in the KEGG database (Fig. 5C-E, Fig. S6). The decreased Kit levels induced by blebbistatin were confirmed by quantitative PCR (Fig. 5F), and decreased signal with disrupted localization of KIT protein was also observed via immunostaining (Fig. 5G). Surprisingly, different epithelial organ cultures in the initial morphogenic stage (SMG, E14; lung, E12; kidney, E13; pancreas, E13) commonly showed a clear decrease in KIT expression upon blebbistatin treatment [SMG, −36.6%; lung, −14.3%; kidney, −16.4%; pancreas, −29.3%, P<0.005] (Fig. S7), implying the general involvement of actomyosin contractility in KIT expression in diverse organs.

During SMG development, KIT is known to transmit signals for the expansion of distal progenitor cells (K14⁺) in epithelial buds, which is crucial for determining cell fate for acinus and duct formation (Lombaert et al., 2013). The mixed pattern of K14⁺/K19⁺ progenitor cells produced by blebbistatin treatment that we observed (Fig. 4F) was highly similar to the pattern of epithelial buds seen in response to inhibition of KIT activity and in KIT knockout models that have been previously reported (Lombaert et al., 2013). In addition to these results, we performed structural evaluation of developing ducts under the condition of KIT inhibition by investigating ZO-1 midline patterns. As expected, ISCK03 (a KIT inhibitor) treatment elicited the diffused patterns of ZO-1 without midline condensation at the site of ductal formation (Fig. 5H), which is analogous to the case of blebbistatin treatment (Fig. 4A). These data indicate that actomyosin is required for a normal level of KIT expression through transcriptional regulation, and that this cascade can be a major reason for immature acinus and duct formation of SMGs.

Next, we dissected the signaling mechanism underlying this genetic event induced by mechanical perturbation. We identified Cyr61, a representative target gene of the YAP/TAZ-TEAD complex, among the top five upregulated genes [Tnc, Cpe, Aldh3a1, Nsnl1 and Cyr61 (Ccn1) in descending order; Fig. 5B]. We then verified the localization pattern of YAP proteins and confirmed the increased nuclear localization of YAP in SMG epithelial cells after blebbistatin treatment [nucleus/cytoplasm=0.84 (control) versus 1.45 (+blebbistatin), P=4.72×10⁻¹¹] (Fig. 6A,B). The same pattern of YAP translocation was also observed in lung and kidney cultures (Fig. S8A-C). Moreover, inhibition of ROCK activity induced the nuclear translocation of YAP (Fig. S8D,E). The localization pattern of YAP sensitively responded to the alteration in ECM stiffness (Fig. 6C,D). The epithelial buds in low-stiffness Matrigel had a high proportion of nuclear YAP proteins, confirming the results from the blebbistatin treatment. Transcriptional activity of YAP was verified by the increased expression of the YAP and TAZ target genes Cyr61 and Ctgf (Ccn2; Fig. 6E). The enhanced YAP and TAZ activity was not wired to conventional Hippo signaling, as demonstrated by the absence of detectable changes in large tumor suppressor kinase (LATS) and YAP phosphorylation upon blebbistatin treatment (Fig. 6F,G). To validate the regulatory role of YAP and TAZ in KIT expression, we suppressed YAP and TAZ activity by knocking down both genes in epiSMG cultures (Fig. S9A,B) and detected increased Kit expression (Fig. 6H). In addition, verteporfin (an inhibitor of the YAP/TAZ-TEAD interaction) treatment completely rescued the suppression of Kit expression by blebbistatin treatment (Fig. 6I). K14 expression levels followed the trend of Kit expression (Fig. 6I). Collectively, these results confirmed the following scenario: actomyosin perturbation by blebbistatin induces YAP and TAZ translocation, and the activated complex negatively regulates KIT expression.

Throughout this work, we reported the mechanical and genetic effects of actomyosin contractility in the initial phase of branching morphogenesis through optical, genetic and biochemical assays (Fig. 6J). We identified the different actin patterns that characterize various developing epithelial organs and dissected the mechanical role of actomyosin contraction in driving cleft formation by using epithelial cultures of salivary glands. We introduced a buckling mechanism, a key physical process in complex morphology of diverse organs (Nelson, 2016), in the context of salivary gland branching morphogenesis. The core mechanism is actomyosin contractility of the outer epithelial cells efficiently transmitting the forces from mechanical instability to an inward buckling force (Fig. 3H). When the actomyosin contractility was perturbed, striking morphological irregularity of epithelial buds was observed (Fig. 1C,E). Based on our observational data and background knowledge, we can anticipate the underlying processes in the morphological distortion of epithelial buds. First, decreased actomyosin contractility allows epithelial outgrowth with disrupted integrity of epithelial layers. Moreover, perturbated contractility decreases epithelial tension, then lateral compressive stress, the triggering factor of the buckling mechanism, potentially passively increasing in epithelial layers. In addition to compressive stress, decreased stiffness of epithelial cells by blebbistatin treatment can alter buckling patterns more sporadically and transiently (Martens and Radmacher, 2008; Nelson, 2016). We speculate that the aforementioned processes collectively elicit the morphological changes of epithelial buds.

Genetically, the physical balance in the epithelial buds intricately sustains KIT expression, a critical mediator of progenitor cell differentiation for acinus/duct formation, by regulating YAP and TAZ transcriptional activity. In branching morphogenesis, the role of YAP and TAZ activity has been recently highlighted in lung and kidney culture models (Lin et al., 2017; Reginensi et al., 2016). Notably, a recent article reported the role of YAP in ductal progenitor patterning in mouse SMGs (Szymaniak et al., 2017). The article provided RNA-seq results evaluating differences in the expression profiles between wild-type and YAP-conditional knockout (Yap-cnull) mice. This study reported an inverse relationship between YAP and KIT expression, which is consistent with our results. However, they detected a decrease in K14 expression in SMGs of Yap-cnull mice, which is different from our result. The authors demonstrated that the Hippo pathway is a key upstream signal of YAP-mediated epiregulin expression, and the paracrine effect of epiregulin regulates K14 expression through EGFR signaling pathway. In this study, we suggest a different mechanism underlying K14 regulation mediated by the YAP-KIT signaling cascade. In terms of the regulatory mechanism of YAP activity, a previous study suggested the canonical Hippo pathway, while our investigation suggested the Hippo pathway-independent mechanotransduction process as an upstream signal. Based on the background, we anticipate that a differential upstream signaling input has a multimodal effect on YAP activity that causes differences in K14 expression.

The concept of YAP and TAZ-mediated mechanotransduction has largely contributed to provide a plausible explanation for various biological events regulated by physical factors. Previous reports were generally based on 2-dimensional in vitro culture systems, and have found that well-defined actin stress fibers engaged to adhesion clutches are a key factor for mechanosensing and YAP and TAZ transcriptional activity. Through the experimental models, diverse regulators involved in the mechanotransduction process have been introduced: YAP and TAZ activity was enhanced when cells were on a stiff substrate and suppressed on soft conditions or upon actomyosin perturbation (Case and Waterman, 2015; Dupont et al., 2011). The mechanical checkpoint modulated by ECM-stress fibers is known to be regulated by F-actin capping and severing proteins (e.g. CapZ and cofilin 1) and by focal adhesion complexes, including adhesion kinases (Aragona et al., 2013; Elosegui-Artola et al., 2016). Moreover, geometric factors, including adhesion area and cell shape, also affect the mechanogenetic processes (Dupont et al., 2011).

However, the experimental model in this study presents a distinct biological context in which developing epithelial clusters interact with ECM (or mesenchymal tissues) in a 3D manner. At a subcellular level, the cells presented the following morphological characteristics: first, F-actin structures consisting of a circumferential cortex-dominant actin pattern; and, second, well-established intercellular adhesions and a highly polarized morphology with low ECM adhesion area (Fig. 2A). According to these differences, YAP and TAZ also showed a unique response to physical changes (Fig. 6A-D), indicating that this system employs a distinct operating mechanism distant from previously reported mechanotransduction modules. Notably, recent work also introduced different modes of YAP and TAZ regulation in the epithelial monolayer (Furukawa et al., 2017; Totaro et al., 2018). The authors reported highly organized cortical actin (circumferential actin belt) in cells cultured at high densities, and a low basal activity of YAP and TAZ. In addition, blebbistatin treatment of the monolayer induced YAP and TAZ nuclear translocation, which is consistent with our results. These results support our findings, and the authors suggest that the structure of actin patterning is a crucial determining factor of this phenomenon. Taken together, we can conclude that YAP and TAZ are differentially regulated, depending on their biological context.

The reciprocal interaction of epithelial cells and ECM is a key event to sculpture the organ morphologies (Dzamba and DeSimone, 2018). In salivary glands, ECM remodeling and related proteins have been extensively studied for cleft formation (Harunaga et al., 2011). Among ECM components, fibronectin has been proposed as a key inducer of cleft formation (Sakai et al., 2003). Fibronectin was concentrated on cleft-forming sites, and regulates the expression of Btbd7 to reduce E-cadherin integrity that is an essential event for cleft formation (Onodera et al., 2010). Fibronectin as well as other ECM proteins also act as a binding partner of integrin families, leading to outside-in signals in epithelial cells. This signaling event is essential to establish the cortical actin patterns of epithelial layers, which is a core structural mediator of our findings. This evidence suggests the necessity of in-depth study of the interplay between epithelial cells and ECM proteins in future research. It will enable researchers to address one of the fundamental questions in the branching processes of epithelial organs: how are clefts localized, initiated and propagated at specific sites?

Including this work, various models and hypotheses explaining mechanosensing and transduction process are now being introduced in broad fields of biology (Ray et al., 2013). The mechanical aspects of diverse biological processes are beginning to be understood, although a number of questions remain, especially regarding the underlying molecular mechanisms. Recently introduced mechanistic studies greatly help to unravel the mechanisms underlying mechanical processes (Furukawa et al., 2017; Meng et al., 2018). We expect that our model of branching morphogenesis will also contribute to elucidating previously unknown aspects of mechanobiology.

The chemical reagents and plasmids used in this study are as follows: 50 µM blebbistatin (Sigma-Aldrich, B0560), 100 nM cytochalasin D (Sigma-Aldrich, C8273), 10 µM verteporfin (Sigma-Aldrich, SML0534), 50 µM Y27632 (Sigma-Aldrich, I2764), Texas Red-X Phalloidin (Thermo Fisher Scientific, T7471), 10 µM ISCK03 (Sigma-Aldrich, I6410) and DiD (Thermo Fisher Scientific, D7757). pEGFP-MRLC T18A, S19A was a gift from Tom Egelhoff (Beach et al., 2011) (Addgene plasmid 35681). AAV-CAG-GFP was a gift from Karel Svoboda (Mao et al., 2011) (Addgene plasmid 28014). pHelper (340202) and pAAV-RC2 plasmids (VPK-422) were purchased from Cell Biolabs. AAV-CAG-MRLC T18A/S19A was cloned by inserting the MRLC PCR amplicon into the AAV-CAG-GFP vector with flanking BamHI/AgeI restriction sites.

ICR strain mice were used for embryonic epithelial organ culture. Animal studies were approved by the CHA University Institutional Animal Care and Use Committee (approval number IACUC180069). All experiments were conducted according to the relevant guidelines and regulations. Submandibular glands (SMGs) were explanted from E13 embryos under a stereo light microscope (Leica, KL300 LED). The explanted organs were plated on 0.1 μm pore-size polycarbonate membranes (Whatman, 110405) and cultured in 50 mm microwell dishes (MatTek Corporation, P50G-1.5-14-F) in Dulbecco's modified Eagle's medium/nutrient mixture F-12 (DMEM/F-12, Gibco, 21041-025) containing 100 U penicillin-streptomycin (Gibco, 15140122), 150 μg/ml ascorbic acid (Sigma-Aldrich, A5960) and 50 μg/ml transferrin (Sigma-Aldrich, T8158). Epithelial rudiment culture was conducted based on previous methods (Rebustini and Hoffman, 2009). Briefly, explanted SMGs were immersed in dispase I (0.5 U/ml; Life Technologies, 17105-041) for 20 min and washed with DMEM/F12 containing 5% bovine serum albumin (Sigma-Aldrich, A7906). The mesenchymal regions of SMGs were then dissected under a stereo microscope, and the isolated epithelial rudiments were embedded in growth factor-reduced Matrigel (Corning, 356231) diluted with DMEM/F12 culture media supplemented with 100 ng/ml Fgf7 (R&D Systems, 251-KG) and 10 ng/ml EGF (R&D Systems, 236-EG) at a 1:1 ratio. A polycarbonate membrane was placed on the gel matrix, and MEM/F12 culture medium was added at 15 min intervals for complete gelation. To modulate the elasticity of Matrigel (Figs 4C and 6D), we decreased the Matrigel concentration from 100% to 10% by increasing the amount of culture media. Other epithelial organs from mouse embryo were cultured according to previously reported methods (Barak and Boyle, 2011; Del Moral and Warburton, 2010; Gaete et al., 2015; Greggio et al., 2013; Voutilainen et al., 2013).

SMG cultures extracted from E13 mouse embryos were embedded in paraffin blocks and sliced into 4 µm sections. Morphology of organ cultures was captured on an inverted light microscope (Nikon, Ti) equipped with a digital camera (Nikon, DS-Ri2) and a CFI Plan Fluor 4× objective (Nikon). Time-lapse images of morphological changes were monitored using JuLI Br live-cell movie analyzer (NanoEnTek). Optical sectioning images were acquired on a confocal laser scanning microscope (Carl Zeiss, LSM700) equipped with Plan-Apochromat 10×, Plan-Apochromat 20× and C-Apochromat 40× objectives (Carl Zeiss), and with 405, 488, 555 and 647 nm excitation lasers. All images of developing buds at the single-cell level were taken by optically sectioning images of the midplane of each bud, showing maximum bud diameter (Fig. S2A). Images were analyzed using Fiji software (Schindelin et al., 2012). Signal intensity, angle and geometric parameters (bud number, organ area, aspect ratio and circularity) were analyzed using the measure and ROI manager plug-ins based on phase contrast and fluorescence images. F-actin signals were quantified by manual detection of cortical actin in a single-cell region (Fig. S1A). Line scan analyses were performed using the Multi Kymograph plugin for quantifying the patterns of pMyosin/F-actin (Fig. S1F) and K14/K19 (Fig. 4H). Nucleus- or membrane-located fluorescent signals were automatically detected using a wand tracing tool (supported by Fiji software) based on the organelle marker signals (DAPI and Did). Cleft depth was measured as the distance between the contact point of two inner bud surfaces and the tangential line of the bud outer surface based on phase contrast images (Fig. 1D, upper panel). For measuring KIT signal intensities (Fig. S7), the regions of each single bud were manually selected and analyzed by ROI manager plug-ins (Fiji software).

Total RNA was extracted from SMG cultures using the RNeasy Mini Kit (Qiagen, 74140). Total RNA (1 μg) was applied for complementary DNA (cDNA) synthesis using the SuperScript III First-Strand Synthesis System (Thermo Fisher Scientific, 18080) with an oligo-dT primer. Real-time PCR was performed using SYBR PCR master mix (Applied Biosystems, 4309155) and a real-time PCR machine (Applied Biosystems, 7200). Nested PCR was performed by following previous protocols (Won et al., 2017). Briefly, the collected SMG cultures were used directly for cDNA synthesis using the Superscript III system. The cDNA was used for two-step PCR procedures with Platinum Taq DNA Polymerase (Thermo Fisher Scientific, 10966) and nested primers. The primer sequences were as follows (5′ to 3′): Kit-outer-F, CAAAGGAAATGCACGACTGC; Kit-outer-R, TTGACCCTCACGGAATGGTC; Kit-inner-F, GGTCGACTCCAAGTTCTACAAG; Kit-inner-R, CTTTTCAAGGGGTCAGCGTC; Kit-F, TCCTCTGGGAGCTCTTCTCCTT; Kit-R, GTTGGACAACCTGCTTGAATGTT; Gapdh-outer-F, ACTTGAAGGGTGGAGCCAAA; Gapdh-outer-R, TTCAGCTCTGGGATGACCTT; Gapdh-inner-F, TCCTGCACCACCAACTGCTT; Gapdh-inner-R, TGGCAGTGATGGCATGGAC; Ankrd1-F, AAACGGACGGCACTCCACCG; Ankrd1-R, CGCTGTGCTGAGAAGCTTGTCTCT; Ctgf-F, GCTTGGCGATTTTAGGTGTC; Ctgf-R, CAGACTGGAGAAGCAGAGCC; Cyr61-F, GCCGTGGGCTGCATTCCTCT; Cyr61-R, GCGGTTCGGTGCCAAAGACAGG; Yap-outer-F, TCCTCCTTTGAGATCCCTGA; Yap-outer-R, TGCTCCAGTGTAGGCAACTG; Yap-inner-F, TTTGCCATGAACCAGAGGAT; Yap-inner-R, TTCAACCGCAGTCTCTCCTT; Taz-outer-F, ATGAATCACCAACACCAGCA; Taz-outer-R, CAAAGTCCCGAGGTCAACAT; Taz-inner-F, TCCCCACAACTCCAGAAGAC; Taz-inner-R, GAAGGCAGTCCAGGAAATCA; K14-outer-F, TGGAGACGTCAATGTGGAGA; K14-outer-R, GGATGACTGAGAGCCAGAGG; K14-inner-F, CAAGGATGCTGAGGAATGGT; K14-inner-R, AGGGATGCTTTCATGCTGAG.

The isolated epithelial rudiments were embedded in a Matrigel-DMEM/F-12 mixture containing 5% carboxylate-modified microspheres (FluoSpheres, 2.0 µm diameter, Invitrogen, F8826). Live imaging was performed on a confocal microscope (Carl Zeiss) with a customized live-cell chamber (Live Cell Instruments). The images were acquired by excitation of 591 nm laser with transmitted light detector (T-PMT). The image analysis software IMARIS was used to track bead movement using SPOT analysis (intensity threshold, 170; size threshold, 5 µm; distance threshold, 3 µm). The tracking results were analyzed using custom Microsoft Excel Visual Basic and MATLAB scripts. The full 3D image area was subdivided into 10×10×10 µm voxels; 3D displacement of each microbead per frame was determined by simple Euclidian distance calculation (the square-root of (dx^2+dy^2+dz^2), where dx, dy, dz indicate displacement along the x-, y- and z-axis), and the average displacement within a voxel was calculated. For 2D projections, the average 3D displacement and direction were calculated for voxels on the same z-axis.

For spatiotemporal analysis of 3-D traction, a custom MATLAB script was written and used. First the imaged volume was resized into 2×2×2 μm voxel resolution. The surface of the epithelial rudiments was determined by thresholding with frame-wise thresholds determined by Otsu's method with time-smoothed (over four frames) and heavily Gaussian-filtered (kernel standard deviation 20 μm) volume (Otsu, 1979). The frame-wise 3D optical flow vectors were calculated from the Lucas-Kanade optical flow equation (Bruce and Lucas, 1981), with a neighborhood radius of 6 μm, and then filtered with the same parameters. The flow vectors at the surface were decomposed into tangential and normal components, with respect to the surface of the epithelial rudiments segmented earlier, and the tangential components were used as the color codes for the surface. The edge protrusion/retraction dynamics at a specific depth was quantified using the previously determined surface boundaries, by applying a marker tracking scheme similar to a previous report (Yang et al., 2016), and then was compared spatiotemporally to the in-plane tangential optical flow dynamics.

Prof. Guang-Yan Yu (Peking University, China) kindly shared SMG-C6 cells. The cell line is authenticated and has been tested for contamination. SMG-C6 cells were cultured in 5% CO₂ at 37°C with DMEM/F12 (Sigma-Aldrich, D8900) containing 2.5% FBS, 5 μg/ml transferrin, 1.1 μM hydrocortisone, 100 nM retinoic acid, 2 nM T3, 5 μg/ml insulin, 80 ng/ml EGF, 5 mM glutamine, 50 μg/ml gentamicin sulfate and 1% penicillin-streptomycin. For salisphere cultures, 10 µl of Matrigel was injected into the microwells of µ-Plate Angiogenesis 96-well plate (Ibidi, 89646) and 1.0∼3.0×10⁴ cells in culture media were plated on the gel after setting.

To produce adeno-associated virus (AAV), human embryonic kidney (HEK)-293 cells were grown to ∼70-80% confluence in DMEM (WelGene, LM-001-05) containing 10% fetal bovine serum (WelGene, S001). Transfection of the plasmids containing AAV constructs was conducted using Polyethylenimine (PEI, Polyscience, 23966) at a 1:2.5 DNA:PEI ratio. AAV expression plasmid, pAAV-RC2 and pHelper were transfected at a 1:1:2 ratio. After 48 h, the AAV products were purified using a pAAVpro Purification Kit (Takara, 6232) according to the manufacturer's instructions. Isolated epithelial rudiments were immersed in DMEM/F-12 media containing AAV particles (∼1:5-1:50) for 1 h at room temperature, washed three times with DMEM/F-12 and incubated in Matrigel.

Epithelial cultures were fixed with 4% paraformaldehyde for 20 min at room temperature and washed with phosphate-buffered saline containing 1% Tween-20 (PBST). The cultures were permeabilized with PBS containing Triton X-100 (PBSX) for 20 min at room temperature and washed again with PBST. Blocking was performed by treatment with PBS containing 0.1% BSA and 10% FBS. After 1 h, the blocking solution was changed to PBS containing primary antibody (1:200) and incubated on a laboratory shaker at 4°C overnight. The cultures were incubated with secondary antibody-PBST solution (1:500) for at least 6 h at room temperature. The following antibodies were used in this study: rabbit polyclonal phospho-myosin light chain 2 (Ser19) (Cell Signaling Technology, 3671), rabbit polyclonal ZO-1 (Thermo Fisher Scientific, 40-2200), mouse monoclonal E cadherin (Abcam, ab76055), mouse monoclonal GM130 (BD Bioscience, 610822), goat polyclonal KIT (Santa Cruz Biotechnology, sc-1494), rabbit polyclonal cytokeratin 14 (BioLegend, 905301), rat monoclonal cytokeratin 19 (Developmental Studies Hybridoma Bank, Troma111), rabbit polyclonal YAP (Cell Signaling Technology, 4912), rabbit monoclonal phospho-YAP (Ser127) (Cell Signaling Technology, 13008), rabbit polyclonal phospho-LATS1 (Ser909) (Cell Signaling Technology, 9157), mouse monoclonal phospho-Histone H3 (Cell Signaling Technology, 9706) and rabbit polyclonal fibronectin antibody (Abcam, ab2413). The following secondary antibodies were used in this study: Alexa Fluor 488 goat anti-mouse IgG (H+L) (Invitrogen, A11001), Alexa Fluor 594 goat anti-mouse IgG (H+L) (Invitrogen, A11005), Alexa Fluor 488 donkey anti-goat IgG (H+L) (Invitrogen, A11055), Alexa Fluor 488 goat anti-rabbit IgG (H+L) (Invitrogen, A11008) and Alexa Fluor 647 donkey anti-rat IgG (H&L) (Abcam, ab150155).

Submandibular salivary glands (SMGs) epithelial rudiments were immersed in 200 µl of Opti-MEM (Gibco, 31985) containing 5 µl of Oligofectamine Transfection Reagent (Invitrogen, 12252) and 500 nM siRNA products for 4 h at 37°C and 5% CO₂. The transfected cultures were embedded in Matrigel plated on glass-bottomed microwell dishes (MatTek Corporation). siYAP, siTAZ and control siRNA were purchased from Dharmacon (ON-TARGETplus siRNA, L-046247-01, L-041057-01 and D-001810-01), and Block-iT Alexa Red fluorescent oligo was purchased from Invitrogen (14750).

In order to construct cDNA libraries with the TruSeq RNA Sample Prep kit, 1 μg of total RNA was used. The protocol consisted of polyA-selected RNA extraction, RNA fragmentation, random hexamer primed reverse transcription and 100 nucleotide paired-end sequencing by Illumina HiSeq4000. The libraries were quantified using qPCR according to the qPCR Quantification Protocol Guide and qualified using an Agilent Technologies 2100 Bioanalyzer. We processed reads from the sequencer and aligned them to the Mus musculus (mm10, http://genome.uscs.edu) using Tophat v2.0.13 incorporating the Bowtie v2.2.3 algorithm to perform the alignment. Transcript assembly and abundance were estimated using Cufflinks. Cufflinks v2.2.1 was used to assemble aligned reads into transcripts and to estimate their abundance. To correct sequence expression count bias, ‘--max-bundle-frags 50000000’ options were used. The transcript counts in isoform level were calculated, and the relative transcript abundances were measured in FPKM (fragments per kilobase of exon per million fragments mapped) from Cufflinks. Gene level expression values were also calculated from the transcript counts. We performed the statistical analysis to find differentially expressed genes. Genes with one more than zeroed FPKM values in the samples were excluded. To facilitate log2 transformation, 1 was added to each FPKM value of filtered genes. Filtered data were log2 transformed and subjected to quantile normalization. Statistical significance of the differential expression data was determined using an independent unpaired two-tailed t-test and fold change in which the null hypothesis was that no difference exists among groups. False discovery rate (FDR) was controlled by adjusting P using the Benjamini-Hochberg algorithm. For the DEG set, hierarchical clustering analysis also was performed using complete linkage and Euclidean distance as a measure of similarity to display the expression patterns of differentially expressed transcripts that are satisfied with fold change≥1.5 and an independent unpaired two-tailed t-test raw value of P<0.05. Gene-enrichment and functional annotation analysis of the significant gene list was performed using Gene Ontology (www.geneontology.org/).

All experiments were performed at least three times for independent tests. Differences between two groups were determined using unpaired two-tailed t-tests. Multiple groups were compared using one-way ANOVA, followed by post-hoc analysis via an unpaired two-tailed t-test with Bonferroni correction. Statistical outliers were excluded using a 2×s.d. cut-off. The Analysis ToolPak-Visual Basic for Application (VBA) supported in Excel was used for all statistical analyses.

We gratefully thank Dr Nury Kim (Institute of Basic Science) for helping to analyze and visualize the trajectories from 3D force distribution analysis.

The peer review history is available online at https://dev.biologists.org/lookup/doi/10.1242/dev.190785.reviewer-comments.pdf