Incoherent Feedforward Regulation via Sox9 and ERK Underpins Mouse Tracheal Cartilage Development

Tracheal cartilage provides architectural integrity to the respiratory airway, and defects in this structure during embryonic development cause severe congenital anomalies. Previous genetic studies have revealed genes that are critical for the development of tracheal cartilage. However, it is still unclear how crosstalk between these proteins regulates tracheal cartilage formation. Here we show a core regulatory network underlying murine tracheal chondrogenesis from embryonic day (E) 12.5 to E15.5, by combining volumetric imaging of fluorescence reporters, inhibitor assays, and mathematical modeling. We focused on SRY-box transcription factor 9 (Sox9) and extracellular signal-regulated kinase (ERK) in the tracheal mesenchyme, and observed a synchronous, inverted U-shaped temporal change in both Sox9 expression and ERK activity with a peak at E14.5, whereas the expression level of downstream cartilage matrix genes, such as collagen II alpha 1 (Col2a1) and aggrecan (Agc1), monotonically increased. Inhibitor assays revealed that the ERK signaling pathway functions as an inhibitory regulator of tracheal cartilage differentiation during this period. These results suggest that expression of the cartilage matrix genes is controlled by an incoherent feedforward loop via Sox9 and ERK, which is supported by a mathematical model. Furthermore, the modeling analysis suggests that a Sox9-ERK incoherent feedforward regulation augments the robustness against the variation of upstream factors. The present study provides a better understanding of the regulatory network underlying the tracheal development and will be helpful for efficient induction of tracheal organoids.


INTRODUCTION
The mammalian trachea is a tubular organ of the respiratory system and is composed of several tissues from different origins, such as endoderm-derived epithelium and mesodermderived cartilage (Cardoso and Lü, 2006). The tracheal cartilage, also known as the tracheal ring, exhibits a C-shaped semi-ring architecture that surrounds the epithelial airway on the ventral side and provides structural support. Abnormal formation of the tracheal rings can collapse the airways and obstruct breathing, which leads to congenital defects, including tracheomalacia and tracheal stenosis (Arooj Sher and Liu, 2016). Thus, a fundamental understanding of the processes underpinning tracheal ring development is essential, yet it is still incomplete.
The development of tracheal rings has been investigated using mouse genetics, which has revealed the importance of multiple transcription factors and signaling pathways. Among these, SRYbox transcription factor 9 (Sox9) is known to be a master regulator that plays a critical role in cartilage differentiation by inducing gene expression of key cartilage matrix molecules, such as collagen II alpha 1 (Col2a1), and aggrecan (Agc1) (Bi et al., 1999;Han and Lefebvre, 2008). It has been demonstrated that SOX9 functions in each successive step of the cartilage differentiation processes, including mesenchymal condensation, commitment to the chondroprogenitor, and maintenance of proliferating chondrocytes (Bi et al., 1999;Akiyama et al., 2002). The importance of Sox9 in tracheal development was revealed by reports which described a complete absence of the tracheal rings in mesenchymal Sox9 knockout mice (Hines et al., 2013;Turcatel et al., 2013). In addition, the haploinsufficiency of Sox9 caused hypoplastic cartilage formation, indicating that SOX9 dosage is critical to tracheal ring formation (Bi et al., 2001). In sonic hedgehog (Shh) knockout mice, Sox9 mRNA expression was lost in the developing tracheae at a later stage of tracheal ring development, leading to the failure of tracheal ring formation, suggesting that Shh plays an important role in controlling Sox9 expression (Park et al., 2010).
Another important signaling pathway is the fibroblast growth factor (Fgf)-extracellular signal-regulated kinase (ERK) signaling axis. Previous studies showed that severe malformations of the tracheal rings were observed when Fgf10 or its main receptor Fgfr2b were either ubiquitously knocked out or overexpressed in the tracheal mesenchyme (Tiozzo et al., 2009;Sala et al., 2011). Abnormal tracheal cartilage formation was also reported in the FGF18 overexpressing mice (Elluru et al., 2009). These results suggest that FGF signaling level must be within a specific range to ensure correct formation of the tracheal rings. The most downstream kinase of the signaling cascade, ERK, is considered to be essential for tracheal development because the genetic deletion of both Mek1 and Mek2 (upstream kinase of ERK) in the mesenchyme resulted in defective tracheal rings (Boucherat et al., 2014). These studies make a strong case for the necessity of Fgf-ERK signaling in the normal development of the tracheal rings; however, the mechanisms through which ERK activation regulates cartilage differentiation are still unknown. It is noteworthy that there are several studies which have presented conflicting results; on one hand, the ERK signaling enhances the expression of cartilage matrix molecules, shown using mouse primary chondrocytes (Murakami et al., 2000), while on the other, the ERK signaling suppresses it, demonstrated using chicken embryonic limb buds (Oh et al., 2000;Bobick and Kulyk, 2004;Zákány et al., 2005).
In this study, we explore the crosstalk between the Shh-Sox9 and Fgf-ERK signaling pathways and the regulatory network that contributes to tracheal ring development. We first show a synchronized temporal profile of Sox9 expression and ERK activity using volumetric imaging of fluorescence reporters. Combined with inhibitor assays, we then show that cartilage matrix genes are positively regulated by Sox9, and in parallel, negatively regulated by ERK activity, in an incoherent feedforward manner. Finally, a mathematical model demonstrates that an incoherent feedforward loop via Sox9 and ERK can explain the dynamics of cartilage matrix gene expression during tracheal cartilage formation.

RESULTS
Mesenchymal Condensation Begins Between E12.5 and E13.5 To morphologically characterize the developing murine trachea, we dissected tracheae from the embryos and processed the tissues using whole mount immunohistochemistry to examine the staining of the nuclei and cell membrane of the tracheal epithelium. We focused on the ventral region of the tracheal mesenchyme, where the tracheal cartilage rings are formed, and we show the frontal and sagittal planes of the ventral region ( Figure 1A). We also show the flattening of nuclei, which were used to evaluate how flat the best fitting ellipse should be in comparison to a circle for nuclei evaluations (Figure 1B).
At E12.5, the mesenchymal nuclei were almost round and homogeneously distributed, in both the frontal and the sagittal planes ( Figure 1A). At E13.5, some mesenchymal nuclei were more elongated and surrounded the cells close to the epithelium, forming a template of the chondrogenic nodule (Figures 1A,C). This process is known as mesenchymal condensation, the initial step of chondrogenesis in general (Goldring et al., 2006). These condensations were observed between the ridges of wavy tracheal epithelium ( Figure 1A'). At E14.5, the peripheral nuclei of the mesenchymal condensations were further elongated and the azimuthal polarization pattern was observed (Figure 1A-C). This phenomenon is consistent with the formation of the perichondrium, the next step of mesenchymal condensation during chondrogenesis (Goldring et al., 2006). These observations indicate that mesenchymal cell differentiation into chondrocytes of the tracheal cartilage begins between E12.5 and E13.5.

Sox9 Expression Increases up to E14.5, and Decreases Thereafter
We next examined the expression of Sox9, the early differentiation marker of chondrocytes, in the developing murine tracheae from E12.5 to E15.5. For this purpose, we employed Sox9-EGFP knock-in mice, in which the expression level of EGFP has been shown to correlate with the endogenous expression of Sox9 (Nel-Themaat et al., 2009;Nakamura et al., 2011). EGFP fluorescence in the ventral side of the epithelium was observed using two-photon microscopy.
Observation by 3D imaging showed the appearance of distinct cell clusters, defined by EGFP intensity, from E13.5 onward (Figure 2A). These high-EGFP clusters are located between epithelial ridges, corresponding to the mesenchymal condensations ( Figure 1A and Supplementary Figure 1A). The condensed mesenchymal cells with high EGFP intensity also exhibited a C-shaped semi-ring structure, indicating formation of the chondrogenic nodule (Supplementary Figure 1B). We then quantified the EGFP intensity in each mesenchymal cell. Although EGFP expression was uniform along the anteroposterior axis in the median line at E12.5, the periodic pattern of EGFP intensity was confirmed at E13.5 and the amplitude of this periodic intensity became larger at E14.5 ( Figure 2B). The median EGFP intensity in the mesenchymal condensations increased 1.18-fold from E12.5 to E13.5, and 1.65-fold from E13.5 to E14.5 ( Figure 2C), which linearly correlates with Sox9 protein levels (Supplementary Figures 1C,D). However, the EGFP intensity in the mesenchymal condensations decreased 0.56-fold from E14.5 to E15.5 (Figures 2A-C), indicating that Sox9 expression is suppressed at E14.5.
We also quantified the expression levels of various cartilage matrix genes (the major determinants of chondrogenesis), including Col2a1 and Agc1, as well as Sox9, by RT-qPCR analysis. Sox9 expression in the whole trachea showed an increasing profile from E12.5 to E14.5 and a slight decreasing profile from E14.5 to E15.5 ( Figure 2D). Note that the RT-qPCR measures cellensemble of gene expressions in the whole trachea. Because the growth of the chondrogenic nodule was evident from E14.5 to E15.5 (Figures 2A,B), the value obtained by RT-qPCR at E15.5 was most likely to be larger than the net expression level in the chondrogenic nodule. Thus, Sox9 expression would exhibit a non-monotonic profile with a peak at E14.5, similar to that obtained from our imaging analyses of Sox9-EGFP ( Figure 2D). However, Col2a1 and Agc1 exhibited monotonic increases during progression through the developmental stages ( Figure 2D). These results suggest that other signaling inputs, together with Sox9, primarily regulate the expression of the cartilage matrix genes.
We also performed time-lapse imaging of the dissected tracheae using the FRET biosensor-expressing mice to continuously monitor the dynamics of the mesenchymal cells and the changes in ERK activity. In tracheae cultured from E12.5 embryos, the mesenchymal cells in the ventral epithelium formed some condensations and ERK activity gradually increased during this process ( Figures 3C,C').

Single cell quantification indicated that averaged ERK activity
Frontiers in Cell and Developmental Biology | www.frontiersin.org increased to a maximum of 1.31-fold when cultured from E12.5 for 40 h ( Figure 3D). The increase in averaged ERK activity was also observed when E13.5 tracheae were cultured for 24 h (1.17-fold, Figure 3D). In contrast, ERK activity in the chondrogenic nodules decreased 0.78-fold when cultured from E14.5 (Figures 3D,E,E'). These observations confirmed that the gradual activation of ERK switches at E14.5 to inactivation.

ERK Activation Is Dispensable for Mesenchymal Condensation
We then explored the role of ERK activation in the development of tracheal cartilage. For these experiments, we treated tracheae dissected from E12.5 embryos with a MEK inhibitor under ex vivo culture conditions for 1 day. Since how the nuclei change shape in a way that is characteristic of mesenchymal condensation in the initial differentiation step occurring between E12.5 and E13.5 (Figures 1A-C), we focused on assessing the shape and distribution of the mesenchymal nuclei, as well as Sox9 levels in the condensations. We observed elongated mesenchymal nuclei, that were distributed concentrically in tracheae treated with the MEK inhibitor, similar to the control samples ( Figures 4A,B,D) and to the E13.5 samples ( Figure 1A,A ′ ). Furthermore, the increase in Sox9 levels during mesenchymal condensation was confirmed in PD0325901treated tracheae as well as in the control explants (Figures 4A,B).
These results indicate that the activation of ERK between E12.5 and E13.5 is dispensable for mesenchymal condensation, which is consistent with a previous report (Oh et al., 2000).
We also examined the effect of mesenchymal condensation on the suppression of cell-generated contractile forces, since it has been shown that cell contraction drives mesenchymal condensation in the development of other organs (Mammoto et al., 2011;Shyer et al., 2017). Tracheae dissected at E12.5 were treated with blebbistatin, an inhibitor of non-muscle myosin II, and cultured for 1 day. Most mesenchymal nuclei remained  round, and the tissue did not exhibit clear condensations or increased Sox9 expression near the epithelium (Figures 4C,D), suggesting that cell-generated contractile forces are required for mesenchymal condensations in tracheal development, but that ERK activation is dispensable for the regulation of contractile forces in this process.

ERK Inactivation Promotes Expression of Cartilage Matrix Genes
To investigate the role of ERK activation in tracheal cartilage formation from E13.5 onwards, we measured the expression levels of Sox9 and several cartilage matrix genes, which show a marked elevation in expression at this stage, in response to treatment with the MEK inhibitor ( Figure 2D). The tracheae were dissected at E13.5 and cultured for 2 days in the presence of PD0325901 under ex vivo conditions. The expression of Spry2, an indicator of Ras-ERK cascade activity (Mason et al., 2006), was decreased by PD0325901 treatment, confirming the inhibition of ERK activation in this assay ( Figure 5A). The expression of Sox9 was unaffected, as anticipated from the data shown in Figure 4. We also confirmed that ERK inactivation did not affect Sox9 expression using imaging measurements of Sox9-EGFP tracheae (Supplementary Figure 2). These data regarding ERK signaling are similar to those from a previous study which reported that inactivation of Fgf10 did not affect Sox9 expression (Sala et al., 2011). Interestingly, ERK inhibition significantly upregulated the expression of Col2a1 and Agc1 (Figure 5A). These findings suggest that ERK inactivation promotes gene expression of these cartilage matrix proteins in a manner independent of Sox9 expression.
We next examined the impact of ERK inactivation on the tracheal cartilage phenotype. To this end, PD0325901 was administered to pregnant mice from E13.5 periodically for 2 days by oral gavage. The morphology of tracheal cartilage rings from dissected embryos was assessed by Alcian blue staining (Figure 5B), and the tracheae were evaluated in terms of anteroposterior (AP) and left-right (LR) lengths ( Figure 5C). In the anterior region of trachea, close to larynx, ERK inactivation significantly increased the AP length (1.33-fold in median, Figure 5D), while it did not affect the LR length (0.98fold in median, Figure 5E), indicating that cartilage matrix accumulation was enhanced due to ERK inactivation. Relative variance in each sample was <5% in either treatment (averaged in-sample coefficient of variance: 0.023 in the DMSO and 0.042 in the PD0325901). In the posterior region of trachea, close to the carina, the AP length was slightly increased by the ERK inactivation (1.07-fold in median) with large variance (Supplementary Figure 3). Together, our results suggest that active ERK represses the transcriptional activity of the cartilage matrix genes, thereby influencing the AP length of the cartilage rings, especially in the anterior region.

Incoherent Feedforward Loop Can Explain the Regulation of Sox9-ERK-Cartilage Matrix Genes
Our results so far present two pathways which have antagonistic effects on the transcriptional regulation of the cartilage matrix genes, namely activation via Sox9 and repression via ERK. Furthermore, since Sox9 expression and ERK activity exhibited similar non-monotonic temporal profiles (Figures 2C, 3B), it suggests that there may be a common upstream regulator X (Figure 6A, left) or positive regulation of ERK by Sox9 (Figure 6A, right). These network topologies are examples of a so-called incoherent feedforward loop (IFFL), known as a network motif, which represent recurrent patterns in transcriptional regulations (Shen-Orr et al., 2002;Mangan and Alon, 2003). We therefore considered the question of how these regulatory networks might function in tracheal development.
We modeled the possible core regulatory networks ( Figure 6A) with a linear relation between X, Sox9, and ERK (see Mathematical analysis in Materials and Methods) because of the similar temporal variations, and found that, in either case, the expression level of the cartilage matrix genes in the steady stateĈ converged to the following form: (1) The first term represents the activation by Sox9 with the coefficient K ′ S , and the second term represents the repression via ERK with the coefficient K ′ E , each of which is modeled by the Hill function with coefficient h S or h E , respectively, with input level X as the simplest form. The function of Equation 1 shows a biphasic inverse U-shape response, in which the output levelĈ increases to a high level, and then decreases, as input level X increases, meaning that the output level reaches at its maximum at an intermediate input level (Figure 6B and Supplementary Figure 4). From Equation 1, the biphasic response modeled by the IFFL can explain the mapping of the non-monotonic temporal profiles of both Sox9 and ERK to the monotonic increases of the cartilage matrix genes ( Figure 2E). We next considered the role of ERK as a negative regulator of the pathways involved in chondrogenesis and queried the benefit of an IFFL involving ERK as a core regulatory design, compared with a single pathway regulated only by Sox9. To answer these questions, we analyzed the robustness of the steady state level of cartilage matrix gene expressionĈ with respect to the input level X in the IFFL compared to the hypothetical single pathway (Figures 6A,C). To do this, we calculated the parameter sensitivity coefficient (Goldbeter and Koshland, 1981) ofĈ with respect to X in both the IFFL (S 1 ) and in the single pathway (S 0 ), and obtained a condition of parameters which satisfies S 1 > S 0 . That is, IFFL is more robust than the single pathway for input variations in the following inequality: Here we introduced a fraction of parameter sets that satisfy Eq. 2, designated Robust index of IFFL. When this value is more than 0.5, the IFFL is a more robust design compared with the single pathway, and vice versa. We then numerically investigated the dependency of the parameters X and K ′ S , both of which are common to the IFFL and the single pathway, without fixing other parameters, and found that the parameter space where the index is more than 0.5 was larger than the area representing an index below 0.5 ( Figure 6D). Collectively, our model analyses show that adopting ERK as a negative regulator of the regulatory chondrogenic pathways makes the system more robust against variations of the upstream regulator of Sox9 during chondrogenesis.

DISCUSSION
In this study, we demonstrate that there is an inverted U-shaped temporal change in Sox9 expression and ERK activity, while the expression of target cartilage matrix genes monotonically increases, during murine tracheal development. Using assays with a MEK inhibitor, we found that inhibiting ERK activity significantly promotes the expression of cartilage matrix genes. The role played by ERK during cartilage differentiation has been controversial in several in vitro studies (Murakami et al., 2000;Oh et al., 2000;Bobick and Kulyk, 2004;Zákány et al., 2005).
Our results show that the ERK signaling pathway functions as an inhibitory regulator of tracheal cartilage differentiation, which has also been demonstrated in embryonic chicken limb buds (Oh et al., 2000;Bobick and Kulyk, 2004;Zákány et al., 2005). Our findings led us to propose a model in which cartilage differentiation is controlled by an IFFL via positive and negative contributions of Sox9 and ERK activity, respectively. To delineate the core regulatory pathways, we tested two networks, one of which included a common upstream regulator without interactions between Sox9 and ERK, while the other network included a positive regulation of ERK by Sox9 (Figure 6A), although both of them resulted in the same class of output in the steady state (Equation 1). It has previously been demonstrated that the expression level of Fgf10 does not change significantly even in Sox9 knockout mice (Turcatel et al., 2013). Furthermore, there is evidence to imply that SHH, an upstream regulator of Sox9, could be a potential activator of Fgf10 in tracheal development (Sala et al., 2011), and it raises the possibility that SHH could also be the common upstream regulator of Sox9 and ERK. It is therefore likely that the former regulatory network is adopted during tracheal ring formation.
The proposed core regulatory networks would be controlled by different types of stimulus, such as mechanical loading (O'Conor et al., 2013), oxidative stress (Zuscik et al., 2008), and electrical signals (Atsuta et al., 2019), through complex signaling pathways (Kozhemyakina et al., 2015). In chondroprogenitor cells derived from chicken limb buds, uniaxial cyclic compression increased the expression level of the cartilage matrix genes and elevated the activity of cyclic AMP-dependent protein kinase A (PKA) (O'Conor et al., 2013;Juhász et al., 2014), known to phosphorylate SOX9 (Huang et al., 2000). In addition, the mechanical loading stimuli activated the pituitary adenylate cyclase activating polypeptide (PACAP) pathway, of which the classical downstream targets include PKA and MAPK signaling cascades (Juhász et al., 2015a,b), suggesting the common input to the Sox9 and the ERK activity. PACAP activation also inhibited the hedgehog signaling activity in the chondroprogenitor cell cultures. Thus, there would be an another IFFL at the upstream of Sox9-ERK signaling layer. Speaking of oxidative stress, it was shown that the hydrogen peroxide (H2O2) inhibited the chondrogenesis using the chicken limb bud cell cultures; H2O2 led to a concentration-dependent decrease of Sox9 expression and an increase of ERK activity in a calcineurin-dependent manner (Zákány et al., 2005). This result suggests that the H2O2 would be a common input regulator to the Sox9 and the ERK activity, which brings a complex signaling crosstalk in the system.
In this study, we primarily focused on the expression of the cartilage matrix genes as an output, but the spatiotemporal dynamics of Sox9 expression and ERK activity are likely to play roles in other aspects of cartilage differentiation. From E14.5 to E15.5, we observed a decrease in Sox9 expression ( Figure 2C). This may reflect the fact that the downregulation of SOX9 has multiple roles in the progression of cartilage differentiation, such as the transition of proliferating chondrocytes to a hypertrophic state, and cartilage vascularization (Akiyama et al., 2002;Hattori et al., 2010;Lefebvre et al., 2019). We also observed that endogenous ERK activation was maintained at a relatively low level from E12.5 to E13.5 ( Figure 3B) and did not affect mesenchymal condensation, an initial step in tracheal ring patterning (Figure 4B). Since Fgf10 overexpression in the mesenchyme between E11.5 and E13.5 was shown to cause disruption of the tracheal ring patterns (Sala et al., 2011), it is possible that ERK activity at below-moderate levels may be a key factor in normal tracheal development. Future studies should clarify the physiological importance of temporal ERK activity dynamics in this process. In addition, further systematic investigations of other signaling pathways, such as Wnt and bone morphogenetic protein (Bmp) signaling (Snowball et al., 2015;Kishimoto et al., 2018), as well as physical interactions with epithelia and smooth muscles (Hines et al., 2013;Yin et al., 2018), would improve our understanding of tracheal development and in vitro organoid systems (Conway et al., 2020;Kishimoto et al., 2020).

Animals
For FRET imaging, we used transgenic mice that ubiquitously express an ERK biosensor with a long flexible linker, which has been described elsewhere (Harvey et al., 2008;Komatsu et al., 2011Komatsu et al., , 2018. Sox9-EGFP mice were provided from RIKEN BRC through the National BioResource Project of the MEXT/AMED, Japan (RBRC05651). Otherwise, we used ICR mice purchased from Japan SLC, Inc. The midnight preceding observation of a plug was designated as embryonic day 0.0 (E0.0), and all mice were sacrificed by cervical dislocation to minimize suffering. All the animal experiments were approved by the local ethical committee for animal experimentation (MedKyo 19090 and 20081) and were performed in compliance with the guide for the care and use of laboratory animals at Kyoto University.

Whole-Tissue Fluorescence Staining and Imaging
Staining and optical clearing of dissected tracheae were performed as described in a previous study (Hirashima and Adachi, 2015). Briefly, the samples were fixed with 4% PFA in PBS overnight at 4 • C. For anti-SOX9 staining, the samples were incubated in 25 mg/mL hyaluronidase (Nacalai Tesque, #18240-36) for 1 h at 37 • C, to digest hyaluronic acid. The samples were then blocked in 10% normal goat serum (Abcam, #ab156046) diluted in 0.1% Triton X-100/PBS (PBT) for 3 h at 37 • C. The samples were treated with primary antibodies overnight at 4 • C, washed in 0.1% PBT, and subsequently incubated in secondary antibodies conjugated to either Alexa Fluor 546 or Alexa Fluor 647 overnight at 4 • C. DAPI was used for nuclear counterstaining (Dojindo Molecular Technologies, #D523-10, 1:200 dilution). The samples were mounted with 10 µL of 1% agarose gel onto a glass dish (Greiner Bio-One, #627871) for stable imaging. Then, the samples were immersed in CUBIC-R+ (Tokyo Chemical Industry Co., # T3741) solution for optical clearing. Images were acquired using the confocal laser scanning platform Leica TCS SP8 equipped with the hybrid detector Leica HyD, using a ×40 objective lens (NA = 1.3, WD = 240 µm, HC PL APO CS2, Leica) and the Olympus FluoView FV1000 with a ×30 objective lens (NA = 1.05, WD = 0.8 mm, UPLSAPO30XS, Olympus).

Alcian Blue Staining
The dissected tracheae were fixed in 4% PFA in PBS overnight at 4 • C, and stained in Alcian Blue Solution (FUJIFILM Wako Pure Chemical Corporation, # 013-13801) for 60 min at 23 • C. The samples were then washed with 20% acetic acid in PBS overnight at 23 • C, and finally clarified in 50% glycerol in PBS for 2 hours at 37 • C. The samples were visualized by the stereo microscopy (SZX16, Olympus).

Drug Administration
For the drug administration under ex vivo culture conditions, (-)-Blebbistatin (FUJIFILM Wako Pure Chemical Corporation, #021-17041) and PD0325901 (FUJIFILM Wako Pure Chemical Corporation, #162-25291) were mixed in the culture medium. Equivalent amounts of DMSO were used as a vehicle control for each drug. For administration by oral gavage, PD0325901 (ChemieTek, #CT-PD03) in 30% PEG400, 0.5% Tween 80 in PBS was administered to pregnant mice at a dose of 25 mg/kg body weight twice a day (at 8:00 a.m and 6:00 p.m) for 2 days from E13.5, i.e., 4 times in total.

Target Region for Analysis
As the mesenchymal condensation simultaneously occurs in the whole trachea, we especially focused on anterior-ventral region, where the future 3rd to 7th tracheal cartilages below the cricoid ring are formed, for analysis. These regions were also convenient for imaging because the frontal plane is parallel to the optical plane.

FRET Image Analysis
The median filter of a 3×3 window was processed to remove shot noises, and the background signal was subtracted each in FRET and CFP channel. Then, the ratio of FRET intensity to the CFP intensity was calculated using a custom-designed MATLAB (MathWorks) script.

Quantification of Nuclear Shape
Since the nuclear configuration exhibited no significant differences throughout the trachea, we randomly chose one representative mesenchymal condensation from the ventral region where the future 3rd to 7th tracheal cartilages below the cricoid ring are formed. The DAPI staining images were smoothened using the median and gaussian filters, and the nuclei were then manually extracted. Flattening or ellipticity is defined as 1-b/a, where a and b are the major and the minor axis length of the best fitting ellipse, respectively. The value is 0 for a circle, and it approaches 1 as it is compressed. All processing was done using ImageJ.

Single Cell Measurement
For the EGFP and FRET measurement at single cell resolution, we first manually extracted each cell or nucleus and measured the average intensity within the extracted region. Then, the background signal was subtracted from averaged signal. For the EGFP signal, the obtained intensity was normalized by the maximum value of the bit depth. All image processing was done using ImageJ.

Software
For digital image processing, MATLAB (MathWorks) and ImageJ (National Institute of Health) were used. For graphics, MATLAB (MathWorks), Imaris (Bitplane) and ImageJ (National Institute of Health) were used. MATLAB (MathWorks) was used for statistical analysis and Mathematica (Wolfram Research) for mathematical analysis.

Graph
For the boxplot, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using the "+" symbol. All of the graphs were prepared in MATLAB.

Modeling
For model construction, we let X, S, E, C be input level, Sox9 concentration, ERK activity, and the expression level of cartilage matrix genes, respectively. Since Sox9 expression and ERK activity showed similar temporal profiles, they would be regulated linearly by input level X. Thus, where α and γ denote production rate and decay rate, respectively. Regarding the alternative regulation, i.e., X indirectly, but not directly, regulates ERK activity through Sox9, the regulation of ERK activity can be represented, instead of Equation S2, asĖ As Sox9 and active ERK antagonistically regulate the expression level of cartilage matrix genes, the dynamics of C can be represented in the simplest form as follows: where K S is the activation coefficient via Sox9, K E is the repression coefficient via ERK, and h denotes the Hill coefficient. Both combinations of Equations S1, S2, and S4, and Equations S1, S3, and S4 led to the following function of C in the steady state: where parameters with prime represent integrated parameters, and the case with α ′ C =1 is shown as the Equation 1 in the main text.

Parameter Ranges in Numerical Investigation
Owing to the importance of relativity for X, K ′ S , and K ′ E , we set those ranges as 0.01 to 1.0. For h S and h E , we set the range from 1 to 5 to consider non-linearity of reactions. α ′ C was arbitrarily set to 1.

Parameter Dependence in the Steady State
Assuming h S = h E = 1 for simplicity and feasibility of analysis, we found that the input value at the output peak X * was determined by the two parameters K ′ S and K ′ E as follows which corresponded to the numerical results (Supplementary Figure 4B'). Also, the function ofĈ was convex upward at the peak as ∂ 2Ĉ Moreover, the maximum level ofĈ is represented aŝ This clearly indicates that the output peak level decreases with increases in K ′ S and it increases with increasing K ′ E , corresponding to the numerical results (Supplementary Figure S4B').

Sensitivity Analysis for Robustness
The parameter sensitivity coefficient with respect to the input value X, denoted S(Ĉ, X), is defined as the relative change inĈ for a given small relative change in X: From Equations S5 and S9, we obtained the explicit form S in the IFFL regulation.
As an alternative regulation, we assumed a hypothetical single pathway regulation, instead of Equation S4 C and from Equations S9 and S11, we obtained the parameter sensitivity coefficient for a single pathway as follows: We then evaluated the robustness of the IFFL compared to the single pathway by the inequality of the sensitivity coefficients: With Equations S10 and S12, Equation S13 reached the inequality condition shown in Eq. 2 in the main text.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by all the animal experiments were approved by the local ethical committee for animal experimentation (MedKyo 19090 and 20081) and were performed in compliance with the guide for the care and use of laboratory animals at Kyoto University.

AUTHOR CONTRIBUTIONS
TY and TH: conceptualization, methodology, validation, formal analysis, investigation, data curation, and writing-original draft. TH: software, visualization, and project administration. MM and TH: resources, writing-review & editing, supervision, and funding acquisition. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the JSPS KAKENHI 17KT0107 and 19H00993, 16H06280, by the JST PRESTO JPMJPR1949, by the CREST JPMJCR1654, by the AMED 19gm5010003h0003, and by the Medical Research Support Center of Kyoto University.

ACKNOWLEDGMENTS
We would like to thank Yu Kurata and Akane Kusumi for technical assistances, and Keisho Hirota and Mitsuru Morimoto for fruitful discussion.

SUPPLEMENTARY MATERIAL
The