Transmural pressure signals through retinoic acid to regulate lung branching

The mammalian lung comprises a branched network of tubes that form a stereotyped tree (Fujii et al., 2020; Metzger et al., 2008). Two morphogenetic motifs build this tree: domain branching of one epithelial tube off the side of another and terminal bifurcation of the tip of an epithelial tube into two daughter branches (Metzger et al., 2008). Domain branching and bifurcation occur concomitantly with stereotyped differentiation of airway smooth muscle (ASM), a contractile tissue that sculpts the proliferating epithelium into new branches in lung explants (Goodwin et al., 2020; Goodwin et al., 2019; Kim et al., 2015). Branching morphogenesis and ASM differentiation are regulated by several soluble signals that form reciprocal feedback loops among the epithelium, mesenchyme and mesothelium, including fibroblast growth factor 10 (Fgf10), Fgf9, sonic hedgehog (Shh), Wnts and bone morphogenetic protein 4 (Bmp4) (Swarr and Morrisey, 2015; Yi et al., 2009). It remains unclear, however, what coordinates the relative rates of epithelial growth and ASM differentiation.

In addition to these biochemical signals, the embryonic lung is exposed to a complex mechanical microenvironment. The developing lung is filled with luminal fluid that is at a higher pressure than the fluid in the pleural cavity, resulting in a positive transmural pressure (ΔP) across the wall of the lung (Harding and Hooper, 1996). Peristaltic contractions of ASM push luminal fluid distally and, during later development, fetal breathing movements further increase ΔP (Harding and Hooper, 1996; Schittny et al., 2000). Clinically, loss of positive ΔP is associated with a decrease in the number of airway branches and pulmonary hypoplasia (Chandrasekharan et al., 2017; Shastry et al., 2012), but the underlying mechanism remains unclear. Using a microfluidic explant system, we recently found that increasing ΔP enhances ASM differentiation and the rate of branching morphogenesis (Nelson et al., 2017).

Here, we used a novel, unbiased RNA-sequencing (RNA-seq) analysis pipeline to identify signaling pathways that are activated by ΔP during lung development. We found that ΔP differentially affects the expression of genes that contain retinoic acid (RA)-response elements (RAREs) in their proximal enhancer regions. We therefore hypothesized that ΔP promotes branching morphogenesis at least in part through RA. Consistently, we found that elevating ΔP increases RA signaling and, through genetic ablation and pharmacological inhibition, that RA signaling is controlled by the mechanosensor Yap. Using predictions from a computational model of the developing airway, we found that RA signaling regulates branching morphogenesis by modulating the relative rates of ASM differentiation and epithelial proliferation. Our results indicate that ΔP signals in part through Yap and RA to regulate lung development.

Using microfluidic chest cavities to apply defined ΔP to embryonic mouse lungs (Fig. 1A), we previously found that airway epithelial branching and ASM differentiation both increase as ΔP increases (Nelson et al., 2017). Here, we carried out bioinformatics analysis on an RNA-seq dataset obtained from lungs held under low (20 Pa) or physiologic (high, 200 Pa) ΔP (Nelson et al., 2017; Olver et al., 2004; Vilos and Liggins, 1982) to identify transcription factors downstream of pressure. Combining the RNA-seq dataset with embryonic day (E)14 lung DNase-hypersensitivity data from the mouse ENCODE project (Yue et al., 2014) enabled us to identify the DNase-hypersensitive regions near the differentially expressed transcripts. By probing these regions for overrepresented motifs, we uncovered transcription factor-binding sites that are possibly associated with the changes in gene expression induced by ΔP (Fig. 1B; Table S1).

Although the DNase-hypersensitive regions that are closest to a transcription start site (TSS) are those most likely to contain cis-acting elements, important gene-regulatory regions can be located several kilobases away from a TSS (Cheng et al., 2014). We therefore investigated the length of each DNase-hypersensitive region and its distance from a TSS. We found that most DNase-hypersensitive regions were shorter than 4000 base pairs (bp) (Fig. 1C). We examined transcription factor-binding motifs that were significantly overrepresented in DNase-hypersensitive regions within different distances of the TSS in our dataset (Fig. 1D). At each distance cutoff, we ranked each transcription factor by its enrichment P-value (Heinz et al., 2010) and calculated the average difference between the rank of a transcription factor and its rank at the largest cutoff, which incorporates the most DNase-hypersensitive regions. As we increased the distance between the TSS and the DNase-hypersensitive regions that we considered, the transcription factors gradually approached their final ranks (Fig. 1E). We also carried out these analyses using DNase-hypersensitivity data from adult mouse lungs and found similar distributions of sizes, locations, and convergence upon their final ranks (Fig. S1A-C). We therefore decided to examine the transcription factor-binding motifs identified by including all DNase-hypersensitive regions that are associated with genes differentially expressed in response to ΔP, irrespective of distance from the TSS.

Using this approach, the top-identified motif was a putative binding site for Snai2 (also known as Slug) (Fig. 1F; Fig. S1E), which promotes epithelial-to-mesenchymal transition (EMT) and cell migration (Nieto, 2013). Although Snai2 has not been previously implicated in lung development, EMT plays a role in epithelial repair after lung injury (Vaughan and Chapman, 2013). We also found that Tead1-binding motifs were significantly overrepresented. Tead1 associates with the mechanosensitive transcription factors Yap and Taz (Han et al., 2018; Huraskin et al., 2016); Yap is crucial for proliferation of the lung epithelium during branching morphogenesis (Lin et al., 2017). In our analyses using DNase-hypersensitivity data from both embryonic and adult lungs, we found that putative binding sites for Nkx2-1, Nkx2-5, and retinoic acid receptors (RARs) were overrepresented (Fig. S1D,E). Nkx2-1 regulates airway epithelial progenitors during branching morphogenesis (Herriges and Morrisey, 2014; Wong et al., 2019), whereas Nkx2-5 regulates the differentiation of myofibroblasts and promotes heart development (Chen and Schwartz, 1997; Hu et al., 2010). Several pressure-responsive genes had at least one nearby motif for Snai2 (28.6% of genes), Nkx2-1 (41.9%), Tead1 (35.4%), Nkx2-5 (40.8%) and retinoid X receptor (RXR) DR1 (35.9%) (Fig. 1G). Most genes had one binding site for each transcription factor in their associated DNase-hypersensitive regions; however, for each transcription factor, several genes contained multiple motifs within their associated DNase-hypersensitive regions (Fig. 1H-J). Increasing ΔP therefore activates the expression of genes with binding sites for transcription factors known to be involved in lung development, indicating that pressure may signal through one or more of these pathways.

RA deficiency has been linked to congenital diaphragmatic hernia (CDH) (Keijzer et al., 2000; Major et al., 1998), wherein weaknesses in the developing diaphragm allow abdominal viscera to herniate into the pleural cavity and compress the lung, resulting in pulmonary hypoplasia at birth and respiratory complications throughout life (Montedonico et al., 2008a). RA signaling is thought to modulate ASM differentiation (Chen et al., 2014), which influences airway branching in cultured lung explants (Goodwin et al., 2019; Kim et al., 2015); consistently, ASM is disrupted in rodent models of CDH (Belik et al., 2003; Pederiva et al., 2012). RARγ translocates to the nucleus in response to mechanical signals (Swift et al., 2013). Furthermore, knocking out RARs disrupts cell-matrix adhesion and the actin cytoskeleton, suggesting an important role for this family of receptors in mechanotransduction (Al Tanoury et al., 2014). These observations led us to hypothesize that high ΔP increases the rate of branching morphogenesis in the embryonic lung at least in part by influencing RA signaling. Because our data suggested a previously unrecognized connection between the mechanical environment of the developing lung and RA signaling, we focused our subsequent analysis on this signaling pathway.

Several proteins are involved in the biosynthesis of RA (Fig. 2A). Retinol is transported through the bloodstream in a complex with retinol binding protein 4 (Rbp4; also known as Rbp) and is internalized when this complex binds to its cell-surface receptor, Stra6, or when free retinol, which is lipophilic, diffuses across the plasma membrane (Fernandes-Silva et al., 2020). In the cytoplasm, RA and RA precursors are bound by carrier proteins including cytoplasmic retinol-binding protein 1 (Crbp1; Rbp1) and cytoplasmic RA-binding proteins 1 and 2 (Crabp1/2). Retinol is oxidized to retinal, which is then further oxidized to RA by aldehyde dehydrogenase 1,2,3 (Aldha1/2/3; Raldh1/2/3). RA binds RAR/RXR dimers at RAREs to activate transcription of target genes. Retinal can also be stored through processing by lecithin retinol acyltransferase (Lrat) and RA is degraded by Cyp26b1.

The expression of Adh1, Rbp1 and Aldh1a1, genes that synthesize or bind to and stabilize RA intermediates, was significantly increased under high ΔP (Fig. 2B-D; Fig. S2A). Expression of Rdh14, which synthesizes retinal, and Crabp1, however, were both significantly decreased. Expression of Cyp26b1 was also significantly decreased. Overall, these data are consistent with increased production and decreased degradation of RA, which would lead to an increase in RA levels, concordant with our DNase-hypersensitivity analysis suggesting increased RA signaling under high ΔP.

To define the tissue compartments that express the RA-biosynthetic machinery, we combined in situ hybridization with single-cell RNA-sequencing (scRNA-seq) analysis of E11.5 lungs. As expected, both approaches revealed that Nkx2-1 is expressed in the airway epithelium (Fig. 3A,B; Fig. S3A). We found that Rara is expressed in all tissues of the lung, most prominently in the mesenchyme and mesothelium (Fig. 3C,D; Fig. S3B) and that Rarb is expressed in the mesenchyme and mesothelium but is most strongly expressed in the trachea and proximal airway epithelium (Fig. 3E,F; Fig. S3C). Our scRNA-seq data do not show epithelial expression of Rarb, likely because the tissue used for sequencing excluded the trachea and proximal bronchi where the Rarb probes most strongly hybridize. Consistently, scRNA-seq data from E15 lungs show that the levels of Rarb are higher in Sox2-expressing proximal epithelial cells than in Sox9-expressing distal epithelial cells (Fig. S3D). We found that Aldh1a2 is expressed primarily in the mesothelium (Fig. 3G,H; Fig. S3E-G). Rbp1 is expressed in all lung tissues (Fig. 3I). Stra6 is expressed in the mesenchyme and in a large proportion of ASM (Fig. 3J). The fact that these genes are expressed at this early stage of development is consistent with previous findings (Malpel et al., 2000) and with a possible role for RA downstream of ΔP in the developing lung.

To test directly whether pressure activates RA signaling, we isolated lungs from E11.5 mouse embryos that express lacZ under the control of the RARE from the Rarb promoter (RARE-lacZ), which reports RA signaling (Rossant et al., 1991), and cultured the explants in the presence of the adenylate cyclase agonist forskolin (Seamon et al., 1981). Forskolin promotes chloride and fluid secretion out of airway epithelial cells into the lumen of the lung (Dekkers et al., 2013), which increases ΔP and physically expands the airway (Fig. 4A). We found that treatment with forskolin increases β-galactosidase activity in RARE-lacZ lung explants (Fig. 4B; Fig. S4A), consistent with our bioinformatics analysis implicating RA downstream of high ΔP. Curiously, histological sections revealed that, although lung explants treated with vehicle control express β-galactosidase primarily in the mesenchyme, those cultured in the presence of forskolin show an increase in β-galactosidase specifically in the epithelium (Fig. 4C,D), without affecting epithelial proliferation (Fig. 4E,F). Increasing ΔP therefore enhances RA signaling in the epithelium of lung explants.

We previously found that high ΔP leads to increased ASM differentiation and epithelial bifurcations (Kim et al., 2015; Nelson et al., 2017). Consistently, immunofluorescence and qPCR analysis revealed that treatment with forskolin maintains the differentiation of ASM at sites of epithelial bifurcation (Fig. 4G; Fig. S4B). To quantify the spatial pattern of ASM differentiation, we correlated α-smooth muscle actin (αSMA) staining intensity with epithelial curvature around distal lung branches and found there was no difference between lungs treated with forskolin and controls (Fig. 4G,H). These data suggest that culture under high ΔP maintains ASM differentiation and allows continued sculpting of the epithelium into bifurcations.

We next sought to define the mechanism by which pressure promotes RA signaling in the developing lung. RARs could be directly mechanosensitive, for example by translocating to the nucleus in response to force. Alternatively, RA signaling could be downstream of another mechanosensor, such as Yap, which associates with Teads to induce transcription in response to mechanical force (Liu-Chittenden et al., 2012). Yap has been implicated in RA signaling in the adult lung epithelium (Ng-Blichfeldt et al., 2018) and in other systems (Mezquita et al., 2018; Xiao et al., 2018). Based on our findings that Tead1-binding motifs are over-represented in cis-regulatory elements of pressure-responsive genes (Fig. 1F), we hypothesized that ΔP signals through Yap to increase the expression of the RA-biosynthetic machinery. Consistent with this hypothesis, we found that increasing ΔP by culturing lungs in the presence of forskolin induces nuclear localization of Yap in the epithelium (Fig. 4I,J). To determine whether Yap is upstream of the RA-biosynthetic machinery, we evaluated a published RNA-seq dataset acquired from Shh-Cre;Yap^fl/fl (epithelial Yap KO) lungs and littermate controls (Lin et al., 2017). This dataset revealed that epithelial knockout of Yap leads to a significant decrease in the expression of Rbp1, Adh1, Aldh1a1 and Crabp1, among other genes in the RA-biosynthetic pathway (Fig. 2E,F; Fig. S2B), consistent with our hypothesis that pressure signals through Yap to induce RA signaling.

To investigate directly whether Yap signals to RA, we cultured RARE-lacZ lung explants in the presence of verteporfin, which blocks Yap-Tead interactions and Yap-dependent gene expression (Liu-Chittenden et al., 2012). As predicted, we found that verteporfin reduces nuclear localization of Yap (Fig. S4C,D) and decreases X-gal staining in RARE-lacZ lung explants (Fig. 4K,L; Fig. S4E) without increasing apoptosis (Fig. S4C). To disrupt Yap specifically in the epithelium, we generated RARE-lacZ embryos in which we conditionally deleted Yap from the airway epithelium (Shh-Cre;Yap^fl/fl;RARE-lacZ). Consistent with our verteporfin experiments, Shh-Cre;Yap^fl/fl;RARE-lacZ lungs showed significantly decreased X-gal stain within the epithelium compared with littermate controls, both in histological sections (Fig. 4M,N; Fig. S4F) and in epithelial rudiments denuded of mesenchyme (Fig. S4G). RA signaling is therefore downstream of Yap in the developing lung.

Given that high ΔP increases both RA signaling and epithelial branching, we hypothesized that RA promotes morphogenesis of the embryonic lung. First, we confirmed that retinoids in the culture medium would not interfere with our ability to detect changes in branching in lung explants (Fig. S5A,B). Next, we inhibited RA signaling by culturing E11.5 lung explants in the presence of the pan-RAR inverse agonist BMS493 (Chen et al., 2014), which decreased β-galactosidase activity in RARE-lacZ reporter lungs (Fig. S5C-H) and reduced expression of the RA pathway target Rarb (Fig. 5A), but did not affect Yap nuclear localization in the epithelium (Fig. S5I,J). We found that explants treated with BMS493 develop dilated airways (Fig. 5B) and form fewer branches (Fig. 5C) than those treated with vehicle control, indicating that RA signaling promotes airway branching. Consistently, treatment with BMS493 also decreases branching in lungs cultured under high ΔP (Fig. S5K,L). Aberrant RA signaling has been linked to CDH (Beurskens et al., 2007) and the nitrofen rodent model of CDH acts, at least in part, by disrupting RA biosynthesis (Chinoy et al., 2001; Mey et al., 2003; Nakazawa et al., 2007). The underdeveloped lungs of the nitrofen model are thought to result from a combination of direct and indirect effects on the lung, the latter through malformation of the diaphragm (Keijzer et al., 2000). Consistently, we found that explants cultured in the presence of nitrofen form fewer, more dilated branches than those cultured with vehicle control (Fig. 5D,E), confirming that nitrofen can directly affect lung branching independently of its effects on the diaphragm (Chinoy et al., 2001). RA signaling therefore promotes branching in lung explants.

We previously found that altering the rate of epithelial proliferation significantly alters the pattern of branching in mesenchyme-free lung explants (Varner et al., 2015). We therefore asked whether RA signaling affects branching by tuning epithelial proliferation. We found that exposure to BMS493 does not significantly affect 5-ethynyl-2′-deoxyuridine (EdU) incorporation in the epithelium in lung explants (Fig. 5F,G), suggesting that RA signaling promotes the formation of branches without altering epithelial proliferation.

Branching is influenced by patterned differentiation of ASM (Goodwin et al., 2019; Kim et al., 2015) and deposition of fibronectin (Sakai et al., 2003) around the airway epithelium in cultured lung explants. We therefore sought to determine whether RA signaling is necessary for ASM differentiation or fibronectin localization around developing epithelial branches. Immunofluorescence staining showed that treatment with BMS493 led to a reduction in ASM localizing to sites of epithelial bifurcation (Fig. 5H,I) but did not affect the distribution of fibronectin around epithelial branches (Fig. S5M,N). Instead, inhibiting RA signaling induced the formation of dilated and buckled epithelium lacking ASM. We therefore hypothesized that inhibiting RA signaling reduces epithelial branching in part by decreasing ASM differentiation. Consistent with this hypothesis, qRT-PCR analysis of E11.5 lung explants cultured in the presence of BMS493 revealed a significant reduction in the expression of αSMA (Acta2), which is a marker of contractile smooth muscle cells (Jaslove and Nelson, 2018), transgelin (Tagln; SM22α) (Li et al., 1996; Liu et al., 2017) and calponin (Cnn1). In contrast, the level of serum response factor (Srf), the main regulator of ASM differentiation, is not significantly affected by treatment with BMS493 (Fig. 5A), suggesting that RA signaling is required for the differentiation of ASM independently of Srf.

Physical stretch induces embryonic lung mesenchymal cells to elongate and differentiate into ASM in culture (Yang et al., 2000). We therefore cultured primary embryonic lung mesenchymal cells on microcontact-printed squares or lines of fibronectin to vary cell elongation and thus mimic the morphologies observed before and after application of mechanical strain (George et al., 2015; Yang et al., 2000). We found that elongated cells develop αSMA-containing fibers along their long axes (Fig. S6A,B), which is characteristic of contractile smooth muscle (Jaslove and Nelson, 2018), and express higher levels of αSMA protein than rounded cells (Fig. S6C,D), consistent with previous reports (Yang et al., 1999). Surprisingly, we found that blocking RA signaling does not affect the expression of αSMA (Fig. S6E). These data suggest that, in the intact lung, RA signaling indirectly promotes ASM differentiation in a non-cell-autonomous manner. This conclusion is consistent with our observation that increasing ΔP enhances RA signaling specifically within the epithelium.

To determine whether a decrease in ASM wrapping can explain the dilated epithelial morphologies we observe when RA is inhibited, we used the finite element method (FEM) to solve a mechanical model of the developing airway. We modeled the lung as three layers of tissue: the epithelium, ASM, and undifferentiated mesenchyme (Fig. 6A,B; Fig. S7A-C), using measurements of epithelial and mesenchymal thicknesses from confocal slices of E11.5 lungs and published estimates of the relative stiffnesses of these tissues (Alcaraz et al., 2003; Chevalier et al., 2016; Smith et al., 2005) (Fig. 6C). We accounted for tonic ASM contractions at this developmental stage (Jaslove and Nelson, 2018) by modeling this tissue as the stiffest of the three layers. As a first approximation, we assumed that the epithelium grows isotropically in the tangential plane to the tissue at every point. Viscoelastic relaxation of embryonic tissue occurs over minutes (Chevalier et al., 2016), whereas epithelial branches form over 12-24 h (Fig. 4A). Therefore, the rate of tissue growth is slow relative to the rate of mechanical relaxation, and we assumed that the tissue is in quasi-static equilibrium at each time point (Goodwin et al., 2019). Finally, to make the model computationally tractable, we used this quasi-static assumption to approximate all three tissue layers as hyperelastic materials.

In the absence of ASM, our simulations revealed that the epithelium dilates and eventually buckles (Fig. 6D,E; Movie 1). Given that the thickness of the airway epithelium relative to its surrounding mesenchyme can vary (Fig. 6C), we explored the parameter space of our models by changing the relative thickness of the epithelial layer. We found that epithelia of all thicknesses dilate in the absence of ASM, but that decreasing the relative thickness of the epithelium (t_ep/R_out) leads to more buckling (Fig. 6F). We next used the computational model to define the effects of ASM wrapping on the geometry of the epithelium and found that the presence of ASM is sufficient to induce bifurcations (Fig. 6G,H; Movie 2), regardless of epithelial thickness (Fig. 6I). The dilated and buckled geometries of models without ASM (Fig. 6E,F) are reminiscent of the epithelia that we observed in explants cultured in the presence of RA inhibitors (Fig. 5H). In contrast, the bifurcated geometries we observed in models with ASM (Fig. 6H,I) are similar to the bifurcated branches in lungs cultured in control conditions (Fig. 5H), supporting our conclusion that the dilated, buckled epithelial morphology observed in explants treated with RA inhibitors results from losses of ASM wrapping. Loss of Yap decreases actomyosin contractility and mechanical tension in the developing lung epithelium (Lin et al., 2017). To determine whether changes in epithelial mechanics could account for the observed effects on lung branching, we varied the stiffness of the epithelial layer in our model and found that the overall geometry of the branches is defined by smooth muscle wrapping regardless of epithelial stiffness (Fig. S7D-I). The similarities between our computational and experimental results suggest that our approximation of the three tissue layers as hyperelastic materials captures the most salient aspects of branch initiation.

An alternative (but not mutually exclusive) hypothesis is that localized increases in epithelial proliferation lead to epithelial outgrowth and terminal bifurcation (Varner and Nelson, 2014). Using our computational model, we found that patterned increases in epithelial proliferation lead to flattening of the epithelium, reminiscent of the early stages of bifurcation (Kim et al., 2015), but are not sufficient to generate bifurcations in the absence of ASM (Fig. 6J,K; Movie 3). In addition, this flattened morphology depends on the relative thickness of the epithelium: when epithelial thickness is set near the low end of our experimentally measured range, the entire structure buckles (Fig. 6L; Movie 4). Our simulations therefore suggest that robust bifurcations require an external mechanical constraint, such as that from ASM differentiation.

The results of our simulations led us to hypothesize that the rate of ASM differentiation relative to that of epithelial growth (proliferation) is a key parameter for determining whether the lung branches or dilates. We experimentally tested this hypothesis by performing EdU analysis on lungs treated with forskolin or BMS493 and comparing the amounts of ASM localizing to sites of bifurcation to the rates of epithelial proliferation. We found that inhibiting RA signaling decreases ASM wrapping relative to epithelial proliferation (Fig. 5J) and leads to dilated epithelial branches. In contrast, increasing ΔP by culturing with forskolin maintains the balance between ASM wrapping and epithelial proliferation (Fig. 4O), resulting in normal branch morphology.

We therefore investigated the effect of changing these relative rates in the computational model. Our results revealed that decreasing the rate of ASM differentiation relative to epithelial proliferation promotes the formation of a dilated tissue; conversely, increasing the rate of ASM differentiation relative to epithelial proliferation promotes bifurcation (Fig. 7A-C). These simulations predicted that inducing proliferation specifically in the epithelium would also lead to a dilated morphology. To test this prediction experimentally, we cultured E11.5 lung explants in the presence of exogenous Fgf1 to promote epithelial proliferation (Varner et al., 2015; Wang et al., 2005) (Fig. 7D,E), which caused the lungs to form fewer branches (Fig. 7F,G) that lack ASM wrapping (Fig. 7H,I). As predicted, inducing the epithelium to proliferate faster than the ASM differentiates (Fig. 7J) led to dilated epithelial branches.

By E12, Shh-Cre;Yap^fl/fl lungs form dilated epithelial cysts instead of branches (Fig. 7K) (Isago et al., 2020; Lin et al., 2017). We found that epithelial deletion of Yap results in decreased RA signaling in the epithelium and decreased expression of ASM markers and epithelial signals that promote ASM differentiation, including Shh, Tgfb2 and Tgfb3 (Fig. S8A,B). We therefore hypothesized that the failure of these lungs to form branches was due, at least in part, to a decrease in the rate of ASM differentiation relative to epithelial proliferation. EdU analysis revealed that epithelial deletion of Yap had no effect on epithelial proliferation (Fig. 7L,M), but decreased the correlation between epithelial curvature and αSMA staining intensity (Fig. 7N) relative to controls. Overall, ASM wrapping relative to epithelial proliferation is greatly decreased in Shh-Cre;Yap^fl/fl lungs (Fig. 7O) and is not rescued by treating explants with forskolin to increase ΔP (Fig. 7P; Fig. S8C). The severe branching defects observed in these lungs is consistent with decreased mechanical constraint on the epithelium from ASM downstream of reduced RA signaling. Altogether, these data reveal that ΔP signals through Yap to activate RA in the epithelium, which increases ASM differentiation and promotes epithelial branching.

Here, we describe an unbiased approach that combines publicly available information about genome-wide regulatory regions with RNA-seq analysis to uncover novel transcriptional regulators of morphogenesis. We used data on chromatin accessibility as a proxy for enhancer regions. As the ENCODE project approaches its goal of acquiring ChIP-seq data for every known transcription factor and provides other measures of chromatin accessibility such as ATAC-seq in a wide array of tissues and at more developmental time points (Encode Project Consortium et al., 2020; He et al., 2020), investigations based on our approach will be able to uncover signaling pathways involved in specific developmental events with more precision and fewer false-positives.

Consensus motifs for Nkx2-1 are overrepresented in DNase-hypersensitive regions near genes increased under high ΔP in the lung, which is expected given that Nkx2-1 is expressed throughout the epithelium during lung development (Herriges and Morrisey, 2014) and that high ΔP increases expression of Nxk2-1 (Nelson et al., 2017). Although we also identified consensus motifs for Nkx2-5, this is likely a false-positive result due to the similarity between its motif and that of Nkx2-1 and the lack of a known role for Nkx2-5 in lung development. Our top identified motif was for Snai2. Snai2 is a transcription factor expressed by migratory cells during gastrulation and development of the neural crest (Mayor et al., 1995; Nieto et al., 1994), and was previously found to be expressed in the mesenchyme of the developing mouse lung (Savagner et al., 1998). Although Snai2 has not been previously implicated in lung development, we speculate that this transcription factor is involved in mesenchymal motility or remodeling during branching morphogenesis. Unexpectedly, we found that binding sites for RARs and Tead1 are also overrepresented in DNase-hypersensitive regions near pressure-responsive genes. The results of our analysis support a model in which high ΔP activates Yap to increase RA signaling in the developing lung (Fig. 8A,B). Our data show that pressure enhances RA signaling in the epithelium and tunes the relative rates of ASM differentiation and epithelial proliferation to promote branching (Fig. 8C). Conditions that disrupt this balance lead to decreased branching and a dilated or buckled epithelial morphology (Fig. 8D).

Dysregulated RA signaling is thought to contribute to weakening of the diaphragm during development in humans (Loo et al., 2018), rats (Montedonico et al., 2006, 2008b; Nakazawa et al., 2007) and mice (Clugston et al., 2010; Nakamura et al., 2020). Diaphragmatic weakening permits herniation of the abdominal viscera and compression of the lung. Previous investigations have also suggested that decreased RA signaling directly decreases airway branching (Chinoy et al., 2001; Keijzer et al., 2000). Our data show that high ΔP enhances the expression of the RA-biosynthetic machinery and induces RA signaling in the developing lung. These data suggest that diaphragmatic herniation and consequent compression of the lungs cause decreased RA biosynthesis within the developing lung tissues themselves, potentially in a positive-feedback loop that worsens defects in pulmonary morphogenesis in the context of CDH. Future studies are needed to define the mechanical signaling pathways, including those involving Yap, that are altered during lung hypoplasia secondary to CDH.

RA serves several functions in the lung – inducing the primary lung bud (Chen et al., 2010, 2007), stimulating proliferation of the mesenchyme (Schuger et al., 1993), shifting alveolar epithelial differentiation from type I towards type II cells (Wongtrakool et al., 2003) and maintaining ASM homeostasis later in life (Chen et al., 2018). Nonetheless, the precise role of RA signaling during airway branching is controversial: some authors report that treatment with exogenous RA up to a maximum dose increases branching in lung explants (Fernandes-Silva et al., 2020; Pereira-Terra et al., 2015; Schuger et al., 1993), whereas others using higher concentrations of RA or longer durations of exposure have concluded that RA decreases branching (Cardoso et al., 1995; Chazaud et al., 2003; Malpel et al., 2000; Mollard et al., 2000). Our finding that RA signaling balances the relative rates of ASM differentiation and epithelial proliferation provides a possible explanation for these apparently contradictory results: either increasing or decreasing RA signaling has the potential to disrupt this balance and lead to a dilated and buckled morphology instead of a properly branched epithelium.

Branching morphogenesis of the mouse lung thus appears to require growth of the airway epithelium and patterned stiffening of the surrounding mesenchyme, which serves as an external constraint to direct the growing epithelium into the bifurcations necessary to form a tree, akin to the mechanism proposed for branch formation in the human lung (Danopoulos et al., 2018; Warburton, 2021). The computational and experimental analyses presented here and in previous work (Goodwin et al., 2019; Kim et al., 2015) reveal that this stiffening results from patterned differentiation of ASM. Curiously, a recent study reported that conditional knockout of the smooth muscle-associated transcription factor myocardin (Myocd) does not affect the ability of the lung to form branches (Young et al., 2020); based on this phenotype, the authors of that study concluded that the airway epithelium branches independently of ASM but did not propose an alternative mechanism for airway branching. What sculpts the branches in the Myocd knockout? It is possible that a non-cellular component of the mesenchyme, such as patterned deposition of ECM, provides sufficient constraint to direct airway epithelial branching in the absence of Myocd. However, our analysis here shows that the morphology of the epithelium is more closely correlated with the pattern of ASM than with that of fibronectin. Furthermore, when we replace the ASM layer in our computational model with a layer that has the thickness and stiffness of ECM, the epithelium fails to bifurcate (Fig. S7J-M), suggesting that matrix deposition is insufficient to act in this capacity.

It is also possible that the subepithelial mesenchymal cells in Myocd-knockout lungs retain their ability to sculpt epithelia by modifying their cytoskeletons and ECM interactions independently of Myocd. Although Myocd knockout decreases the expression of Acta2 and Myh11, grouping genes by their functions and their biochemical pathways reveals that most genes involved in smooth muscle differentiation, cytoskeletal organization, contraction, and ECM adhesion are unaffected by loss of Myocd (Fig. S9). This observation is also consistent with our previous findings from fluorescent reporter mice, which revealed that the epithelium begins to bifurcate at the earliest stages of ASM differentiation, after the Acta2 gene promoter has been activated but before αSMA protein can be detected in lungs imaged at low magnification (Kim et al., 2015). Finally, bioinformatics analysis of the ASM differentiation program reveals that nascent immature smooth muscle cells activate the expression of cytoskeletal components associated with cell stiffening before they begin expressing Myocd (Goodwin et al., 2020). Given the coupled dynamics of epithelial growth, ASM differentiation and branching morphogenesis, we postulate that the shape of the growing epithelium is influenced by this presumed gradual stiffening that occurs at the earliest stages of the ASM differentiation program; we predict that these stiffer, immature smooth muscle cells direct airway epithelial branching and that their presence is unaffected by knockout of Myocd (Fig. 8E). Therefore, patterned ASM differentiation directs the branching of the growing airway epithelium but Myocd expression does not.

It remains unclear how the spatial pattern of ASM is specified in the developing lung (Goodwin et al., 2019). In general, smooth muscle differentiation is controlled by both mechanical stretch (Jakkaraju et al., 2005) and soluble signals (Kim et al., 2015). Surprisingly, one study found that inhibiting RA led to modest increases in the expression of ASM markers (Chen et al., 2014). However, this previous work cultured explants in the presence of 10% serum, which can significantly alter the concentration of retinol (Engedal et al., 2006) and would also be expected to alter cell proliferation, which can affect RA responsiveness in some tissues (Marcelo et al., 2013; Qiu et al., 2020). Here, we measured the distribution of αSMA protein, which forms the core of the smooth muscle cytoskeletal machinery (Jaslove and Nelson, 2018), and found that inhibiting RA decreases αSMA-positive cells at the tips of epithelial branches. Our observations support the conclusion that inhibiting RA decreases the rate of ASM differentiation.

Although the lung clearly develops in the presence of mechanical forces, our understanding of its mechanosensory and mechanotransduction machinery is still in its infancy. Our data reveal that ΔP promotes the nuclear translocation of Yap in the developing lung epithelium and activates the expression of genes with Tead-binding sites in their nearby regulatory regions. A central role for this mechanosensory pathway is supported by the phenotype of Yap-knockout animals, which fail to form lungs or initiate branching (Isago et al., 2020; Lin et al., 2017). It is likely that some of these defects are due to dysregulation of RA signaling. Consistent with our results, recent work has uncovered links between Yap nuclear localization and the expression of RA receptors and biosynthetic machinery in colorectal cancer cells (Bauzone et al., 2021).

Our findings identify RA signaling as one mechanism by which mechanical forces regulate morphogenesis in the developing lung. RA contributes to several developmental processes, from inner ear development to body-axis elongation (Cunningham and Duester, 2015; Ono et al., 2020), and is crucial for the formation of organs in which mechanical forces are also required for tissue growth, including the eye (Molotkov et al., 2006). However, few studies have investigated whether mechanical forces influence RA signaling or vice versa (Guo et al., 2016; Xiao et al., 2018). Our finding of a role for RA signaling in mechanotransduction downstream of Yap in the developing lung suggests that future studies might reveal a similar interplay between RA and Yap in the development of other tissues.

Differential gene-expression analysis was carried out as previously described for an RNA-seq dataset of lungs cultured under high or low ΔP, available from the Gene Expression Omnibus with the GEO accession GSE90148 (Nelson et al., 2017). A custom Python script was created to identify sequences of DNase-hypersensitive regions near genes that were significantly differentially expressed under high versus low ΔP. Using the locations of DNase-hypersensitive regions from the mouse ENCODE lung DNase-hypersensitivity dataset (Yue et al., 2014), we generated a list of the locations of DNase-hypersensitive regions in the mm10 mouse genome build downloaded from the University of California, Santa Cruz (UCSC) genome browser (Kent et al., 2002). We then used the list of all NCBI Refseq mouse genes from the UCSC genome browser to locate the TSS of each gene. We compared the TSS of each gene to the start and end points of each DNase-hypersensitive region on the same chromosome and assigned each DNase-hypersensitive region as associated to the gene with the nearest TSS. We then extracted the sequences of all DNase-hypersensitive regions assigned to genes that were significantly increased under high ΔP as a FASTA file. This process was repeated on a list of random mouse genes generated by the Regulatory Sequence Analysis Tools web interface (Nguyen et al., 2018) to obtain a FASTA file of background sequences. The high ΔP FASTA file was assayed for motifs that were overrepresented relative to the background file using the HOMER findMotifs.pl FASTA motif analysis script (Heinz et al., 2010). Gene identifiers were converted between Refseq and official gene-symbol formats using DAVID (Huang et al., 2009a,b). To determine whether genes involved in RA signaling depend on the expression of Yap in the airway epithelium, we analyzed a previously published bulk RNA-seq dataset (GSE93339) from E12.5 lungs of Shh-Cre⁺;Yap^fl/fl embryos and littermate controls (Lin et al., 2017). For each gene of interest, we retrieved the log₂ fold changes and adjusted P-values provided in this GEO entry. We also obtained published bulk RNA-seq data from E13.5 lungs (five controls and five mutants) in which Myocd had been deleted from the mesenchyme (GSE143394; Young et al., 2020). Data were processed using the DESeq2 package (Love et al., 2014), which was used to import data and estimate size factors and dispersions. DESeq2 was then used to calculate differential expression of genes and their significance using the Wald test. We focused on genes from GO biological process terms and KEGG signaling pathway terms involved in cell contractility and adhesion. We then used volcano plots to visualize adjusted P-values and log₂ fold changes for each of these gene sets.

Left lobes from E11.5 mouse lungs were treated with dispase for 10 min, dissociated with tungsten needles, and resuspended in DMEM without HEPES supplemented with 5% fetal bovine serum (FBS). The resulting solution was passed through a 40-μm-diameter pore filter. Samples were then processed at the Princeton Genomics Core Facility on a 10x Genomics Chromium Controller and with an Illumina Nextera library prep kit to barcode and reverse-transcribe RNA from single cells, and the cDNA was sequenced using an Illumina HiSeq 2500 Rapid Flow Cell to generate paired-end reads. FASTQ base calling was performed using 10x CellRanger software. The resulting sequence files were processed using the R Bioconductor package Seurat 3.0 to scale the data, find variable features and perform cell-stage correction and cell clustering based on expression of known markers of cell types in the developing lung (Butler et al., 2018; Stuart et al., 2019). Blood cells were removed at this stage to focus the analysis on native lung cell types. scRNA-seq data are available from the Gene Expression Omnibus with the GEO accession GSE153069 (Goodwin et al., 2020). Transcriptional differences between proximal and distal epithelial cell identities were investigated using a previously published dataset (GSE144170) of scRNA-seq data from airway epithelial cells isolated at E15.5 (Gerner-Mauro et al., 2020) also processed with Seurat. Only cells expressing Cdh1 were included. To focus on proximal and distal cell types, cells were classified based on whether they expressed Sox2, Sox9 or neither.

CD1 mice were obtained from Charles River. RARE-hsp68lacZ RA receptor activation reporter (RARE-lacZ; JAX 008477) (Rossant et al., 1991), Shh-Cre-EGFP (Shh-Cre; JAX 005622) and Yap-flox/flox (Yap^fl/fl; JAX 027929) mice were obtained from the Jackson Laboratory. Breeding pairs were set up locally for timed pregnancy. Yap^fl/fl;RARE-lacZ females were generated by breeding Yap^fl/fl mice with RARE-lacZ mice and Shh-Cre⁺;Yap^fl/+ males were generated by breeding Yap^fl/fl mice with Shh-Cre⁺ mice. Yap^fl/fl;RARE-lacZ females were mated with Shh-Cre⁺;Yap^fl/+ males to generate timed pregnancies of Shh-Cre;Yap^fl/fl;RARE-lacZ (lung epithelial Yap knockout) and littermate control RARE reporter embryos. Pups and embryos were genotyped by PCR and gel electrophoresis on DNA extracted from pup tail-snips or embryo hindlimbs and tails. All primers used for genotyping are listed in Table S2. Embryos were isolated in accordance with institutional guidelines following the National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by the Princeton University Institutional Animal Care and Use Committee.

Lungs were dissected and cultured as explants as previously described (Kim et al., 2015). Briefly, lungs were dissected from embryos in cold phosphate buffered saline (PBS) supplemented with 1% penicillin/streptomycin (Gibco, Thermo Fisher Scientific). Explants were cultured on microetched nucleopore filters (8 μm pore size, 25 mm diameter; Whatman) in DMEM/F12 medium without HEPES buffer (HyClone, GE) supplemented with 5% filtered, heat-inactivated FBS and 1% penicillin/streptomycin (Gibco). Explants were treated with BMS493 (1 µM, Sigma-Aldrich), nitrofen (5 μM, Sigma-Aldrich), forskolin (5 µM, Sigma-Aldrich), verteporfin (2 µM, Sigma-Aldrich), or dimethyl sulfoxide (DMSO) vehicle control (ATCC), or Fgf1 (1 μg, R&D Systems) or PBS vehicle control. Medium was replaced every 24 h in all explant culture experiments. Cell proliferation was measured by EdU incorporation using a Click-iT EdU Imaging Kit (Invitrogen) by incubating lungs at 37°C with 20 μM EdU for 30 min. In culture experiments, EdU was added directly to the culture medium; in Yap epithelial knockout experiments, lungs were dissected and cultured on nucleopore filters for 1 h before adding EdU to the culture medium.

Primary mesenchymal cells were harvested from E11.5 CD1 or RARE-lacZ mouse lungs using a differential plating technique (Schuger et al., 1993). Briefly, lungs were digested with 0.25% trypsin and 0.05% EDTA (Gibco) while pipetting vigorously. Homogenized mixtures of embryonic lung cells were allowed to attach directly onto microcontact-printed islands of fibronectin for 1 h before unattached cells were washed off with PBS and new medium was added, allowing only enough time for mesenchymal cells to attach. Cells were cultured in Eagle's Minimal Essential Medium (EMEM; ATCC) supplemented with 10% FBS, 1% penicillin/streptomycin, and 50 µg/ml gentamycin (Gibco). Cultures were treated with BMS493 (1 µM) or DMSO vehicle control where indicated.

Regularly spaced 30-μm-wide squares or lines of fibronectin were microcontact printed onto UVO-treated polydimethylsiloxane (PDMS)-coated glass coverslips using sterile stamps of PDMS that were adsorbed with fibronectin (25 μg/ml, Corning) in PBS at 4°C overnight and dried under compressed nitrogen (Gomez and Nelson, 2011). Coverslips were then passivated using a solution of 1% synperonic F108 (Fluka) for 20 min and washed twice with PBS.

Total RNA was extracted using TRIzol^® reagent according to the manufacturer's instructions, followed by cDNA synthesis using the Verso cDNA synthesis kit with poly(dT) primers (Thermo Fisher Scientific). Transcript levels were measured by qRT-PCR using an Applied Biosystems StepOne Plus instrument and iTaq Universal SYBR Green Supermix (Bio-Rad). Amplification was followed by melting-curve analysis to verify the presence of a single PCR product. Primers specific for each target are listed in Table S3. The relative expression of each transcript was quantified using the ΔΔC_T method with the 18S ribosomal subunit as the calibrator in each sample (Schmittgen and Livak, 2008). The efficiency of each primer was between 90 and 105% (Schmittgen and Livak, 2008) according to a standard curve of serially diluted template.

Cells were fixed with 4% paraformaldehyde (PFA; Fisher) for 15 min, permeabilized in 0.5% Triton X-100 (PBS-T), blocked in PBS-T supplemented with 10% FBS, incubated with primary antibody for 1 h at room temperature, washed, blocked, incubated with secondary antibody for 1 h at room temperature, and mounted with Fluoromount-G (Southern Biotech). Lungs were fixed, permeabilized in 0.1% PBS-T, blocked in PBS-T supplemented with 5% FBS and 0.1% bovine serum albumin (BSA; Sigma-Aldrich), incubated overnight with primary antibody at 4°C, washed, and incubated overnight with secondary antibody at 4°C. Cells and lungs were counterstained with Hoechst 33342 (Invitrogen). The antibody sources, targets and dilutions are listed in Table S4. Lungs were dehydrated progressively in methanol and cleared in 1:2 benzyl alcohol:benzyl benzoate (Murray's clear).

Immunostained samples were imaged on a Nikon Ti-U inverted microscope with a Hamamatsu camera and a spinning disk confocal attachment (X-Light) using a Plan Apo 4× 0.2 NA, S Plan Fluor 20× 0.45 NA or Plan Fluor 40× oil 1.3 NA objective. Tiled images were stitched together using the stitching plugin (Preibisch et al., 2009) for the ImageJ FIJI distribution (Schneider et al., 2012). Branching morphogenesis was quantified by counting the number of terminal branches in confocal stacks using the FIJI cell counter plugin. Brightfield stacks taken on a Leica M205FA stereomicroscope were processed with the FIJI linear stack alignment with SIFT and extended depth-of-field plugins (Forster et al., 2004). Displayed images of sections were white-balanced in FIJI. X-gal staining intensity was quantified with a custom MATLAB script by outlining the whole lung, epithelium or mesenchyme, adding together the average red, green and blue (RGB) values of the selected area and subtracting the sum of the RGB values of a random background area from the bright-field images. Smooth muscle wrapping around epithelial branches was measured by finding the Pearson correlation between epithelial curvature (concave versus convex regions) and αSMA staining intensity along the outline of each epithelial branch. αSMA staining intensity was determined by outlining the epithelium in a confocal slice through the widest portion of a branch and taking a maximum projection of the αSMA staining at each point along the outline. The magnitude of curvature was calculated by constructing a triangle using the points on the path 5 pixels before and after each point and finding the radius of the circle circumscribed around the triangle. After skeletonizing the lung based on the epithelial outline, a positive curvature was assigned if the vector between the skeleton and the point was in the same direction as the vector between the point and center of curvature, and negative otherwise. EdU incorporation was quantified by selecting a confocal slice through the airway epithelium, counting the number of EdU⁺ nuclei per length of epithelium, and averaging over one to three separate branches on each lung using a custom MATLAB script. Intensity of αSMA staining in micropatterned cells was measured by sum-projecting confocal stacks and quantifying cell shape and the total intensity of staining using a custom MATLAB script. Yap1 quantifications were also performed using a custom MATLAB script. First, images were background-subtracted, then nuclear masks were generated based on Hoechst staining. Yap1 pixel intensities within or outside of nuclei in the airway epithelium were then measured in five randomly selected regions of interest per confocal stack. We then took the ratio of nuclear to cytoplasmic intensity as a readout of Yap activity.

LacZ-expressing tissues were stained with X-gal by fixing for 30 min in 1% formaldehyde, 0.2% glutaraldehyde, 2 mM MgCl₂, 5 mM EGTA, 0.02% tergitol (NP-40), permeabilizing with 0.02% tergitol, staining at 4°C or room temperature from 4 h to overnight in a solution of 5 mM K₃Fe(CN)₆, 2 mM MgCl₂, 0.01% sodium deoxycholate, 0.02% tergitol and 1 mg/ml X-gal, and post-fixing with 4% PFA. Stained tissues were either cleared in 75% glycerol and imaged in wholemount on a Leica M205FA stereomicroscope or embedded in optimal cutting temperature compound (OCT), cut into 10-μm-thick sections on a Leica CM3050S cryostat, mounted with Fluoromount-G and imaged on a Nikon Eclipse TS100 with an attached Nikon CoolPix camera or Leica M205FA stereomicroscope. To examine X-gal stain in epithelial rudiments denuded of mesenchyme, lungs from RARE-lacZ reporter mice were explanted and incubated at room temperature for 5 min in dispase (Corning), which was subsequently inactivated in DMEM:F12 media supplemented with 5% FBS and 0.1% BSA. Tungsten needles (Fine Science Tools) were then used to remove the mesenchyme and the resulting epithelial rudiments were fixed and stained with X-gal as described above.

Probe sequences were generated by PCR (primers listed in Table S5) from embryonic mouse lung cDNA and ligated into a dual promoter pCRII-TOPO vector (Invitrogen) that was transformed into Max Efficiency DH5α chemically competent cells (Invitrogen) for amplification. Plasmids were harvested using a Qiagen Midi Prep kit and digested using the appropriate restriction enzymes for each promoter. Antisense and sense digoxigenin-uridine-5′-triphosphate (UTP)-labeled riboprobes were synthesized using a DIG RNA Labeling Kit (SP6/T7) (Roche) and purified using quick spin columns for radiolabeled RNA purification (Roche). Explanted lungs were fixed in 4% paraformaldehyde overnight, dehydrated in methanol, rehydrated and digested with Proteinase K (Ambion) before incubating with probe overnight at hybridization temperature in a Bambino Hybridization Oven (Boekel Scientific). Hybridized tissues were probed with AP-anti-DIG antibody Fab fragment (Roche) and color was developed using nitro blue tetrazolium chloride (NBT; Roche) and 5-bromo-4-chloro-3-indolyl-phosphate (BCIP; Roche) solutions. Stained tissues were cleared in 75% glycerol and imaged in wholemount on a Leica M205FA stereomicroscope.

An idealized 3D model and finite-element mesh of an epithelial branch surrounded by smooth muscle and undifferentiated mesenchyme were created using the computer-aided design functionality of Gmsh (Geuzaine and Remacle, 2009). As described previously (Goodwin et al., 2019), the epithelium, mesenchyme, and smooth muscle compartments were modeled as connected hyperelastic solids, each with different elastic moduli. Epithelial growth was modeled as an isotropic expansion of the size of the epithelial layer in a plane tangential to the epithelium at each point. The FEniCS Python module (Logg et al., 2012) was used to numerically minimize the potential energy at each timestep as the epithelium grew. The results were visualized and plotted using ParaView (Squillacote and Ahrens, 2006) and the Matplotlib Python package. Additional details on the mathematical formulation of the problem are available in supplementary Materials and Methods.

Statistical analyses were carried out using the Python SciPy Stats and statsmodels packages or Microsoft Excel. Unpaired, two-tailed Student's t-tests were used for pairwise comparisons and ANOVA with Tukey's post-hoc test was used to compare more than two groups. Statistical significance in figures is denoted as *P<0.05, **P<0.01 and ***P<0.001 and exact values for significant results are noted in the figure legends. In box-and-whisker plots, the midline, box, and whiskers represent the median, upper and lower quartiles, and furthest points, respectively, with the whiskers excluding outliers defined as points further than 1.5× interquartile range outside of the box limits. Individual data points on graphs represent individual lungs from three independent experiments, except where noted. The standard error of the mean (s.e.m.) of ASM wrapping to epithelial proliferation ratios was calculated by error propagation of the s.e.m. from the correlation of ASM to curvature and epithelial proliferation. Permutation tests to compare ASM wrapping to epithelial proliferation were carried out by defining a test statistic ⁠, then randomly reassigning the ASM wrapping and proliferation data points to treatment or control groups and calculating the test statistic 10,000 times. The P-value was calculated by finding the proportion of permuted T that had a larger absolute value than the experimentally determined T.

We are grateful to Sarah Paramore and other members of the Tissue Morphodynamics Group for helpful discussions and experimental assistance, and Lance Parsons of the Princeton Genomics Core for assistance in processing bulk RNA-seq data.