- Research article
- Open Access
Salivary gland branching morphogenesis: a quantitative systems analysis of the Eda/Edar/NFκB paradigm
BMC Developmental Biology volume 9, Article number: 32 (2009)
Ectodysplasin-A appears to be a critical component of branching morphogenesis. Mutations in mouse Eda or human EDA are associated with absent or hypoplastic sweat glands, sebaceous glands, lacrimal glands, salivary glands (SMGs), mammary glands and/or nipples, and mucous glands of the bronchial, esophageal and colonic mucosa. In this study, we utilized EdaTa(Tabby) mutant mice to investigate how a marked reduction in functional Eda propagates with time through a defined genetic subcircuit and to test the proposition that canonical NFκB signaling is sufficient to account for the differential expression of developmentally regulated genes in the context of Eda polymorphism.
The quantitative systems analyses do not support the stated hypothesis. For most NFκB-regulated genes, the observed time course of gene expression is nearly unchanged in Tabby (EdaTa) as compared to wildtype mice, as is NFκB itself. Importantly, a subset of genes is dramatically differentially expressed in Tabby (Edar, Fgf8, Shh, Egf, Tgfa, Egfr), strongly suggesting the existence of an alternative Eda-mediated transcriptional pathway pivotal for SMG ontogeny. Experimental and in silico investigations have identified C/EBPα as a promising candidate.
In Tabby SMGs, upregulation of the Egf/Tgfα/Egfr pathway appears to mitigate the potentially severe abnormal phenotype predicted by the downregulation of Fgf8 and Shh. Others have suggested that the buffering of the phenotypic outcome that is coincident with variant Eda signaling could be a common mechanism that permits viable and diverse phenotypes, normal and abnormal. Our results support this proposition. Further, if branching epithelia use variations of a canonical developmental program, our results are likely applicable to understanding the phenotypes of other branching organs affected by Eda (EDA) mutation.
Branching morphogenesis is a common mechanism of mammalian development (salivary glands, lungs, mammary glands, pancreas, kidney, etc.), and has been a classic topic of study for generations of developmental biologists . Based on recent findings regarding signal transduction pathways and transcriptional control, it has reasonably been proposed that all branching systems use variations of a canonical developmental program . Ectodysplasin-A, a protein required for epithelial differentiation, appears to be an important constituent of such a program.
Ectodysplasin-A (Eda in mouse, EDA in human) is mapped to the X-chromosome [3, 4]. Eda is a glycosylated, oligomeric type II membrane protein with three collagenous repeat domains and a TNF homology domain [3–6]. Eda is shed from the cell membrane and binds as a trimer to its trimerized cognate receptor (Edar) [7, 8]. Like TNF/TNFR signaling, Eda/Edar signaling is thought to be primarily through the canonical NFκB pathway [9–13]. The initial appearance of Eda and Edar proteins in Late Pseudoglandular/Canalicular Stage (E15) submandibular salivary glands (SMGs) in vivo indicate that they participate in late branching morphogenesis and histodifferentiation . In vitro study of the Eda/Edar pathway in embryonic SMG development indicates that this pathway is important for epithelial cell proliferation, lumina formation, and histodifferentiation; exogenous Eda delivered to SMG explants upregulate NFκB activation and nuclear localization .
Mutations in mouse Eda or its human homologue EDA result in hypohydrotic ectodermal dysplasia (HED), a syndrome variably characterized by absent or hypoplastic teeth, hair, sweat glands, sebaceous glands, lacrimal glands, salivary glands, mammary glands and/or nipples, and mucous glands of the bronchial, esophageal and colonic mucosa [14–22].
The mouse EdaTa(Tabby) allele is characterized by a ~2 kb deletion [3, 4]. Specifically, genomic DNA hybridized with an exon 1 probe shows a deletion including the coding region, and primers for DNA flanking exon 1 fail to amplify in a PCR assay (Jackson Laboratories; http://www.informatics.jax.org). Importantly, RT-PCR assays of embryonic EdaTa(Tabby) skin reveals that EdaTatranscript levels are about 10–20% of that seen in the wildtype [23, 24]. This may reflect the fact that DNA sequences that regulate gene transcription occupy no fixed position relative to coding DNA regions and are often diffuse and widely dispersed, including secondary ("shadow") enhancers [25, 26]. Secondary enhancers map far from the target gene and mediate activities overlapping the primary enhancer, accounting for why deletions of well-defined enhancers are sometimes associated with weak or no phenotypic abnormality .
In this study, we utilized EdaTa(Tabby) mutant mice to investigate the in vivo relationship between Eda/Edar signaling and progressive (E13-NB) SMG morphogenesis and histodifferentiation. Our investigation reveals that, from the outset (E13), embryonic Tabby SMGs are smaller, exhibit fewer branches, and are developmentally delayed compared to wildtype (WT) glands. By E18, Tabby glands remain smaller with fewer presumptive acini than WT, though both display a similar degree of terminal differentiation (presumptive functional maturation), as evidenced by mucin (MucCAM/Muc10) protein expression. The key question is how altered Eda function is related to abnormal SMG ontogeny.
Developmental biology has progressed from studying one or two molecules at a time to studying scores of molecules. Such studies require complex experiments and mathematical analyses of their results. While the preferred approach for data analysis has been statistical, its value is limited for mechanistic understanding of signal transduction, mostly because the approach is correlational. To understand how a specific ligand/receptor binding (or lack thereof) produces a change in cell/tissue behavior, one is compelled to utilize a dynamic model of a larger relevant genetic subcircuit. Such dynamic models allow one to test different sets of unbiased assumptions and determine if the predicted behavior matches the actual.
Thus, the wider goal of the present study was to understand how a marked reduction in functional Eda propagates with time through a defined subcircuit that includes 5 signaling pathways that are both critical for SMG ontogeny and share post-activation downstream targets. More specifically, it was the central objective of our extensive quantitative study to demonstrate the way in which a mechanistic pathway diagram can be programmatically transformed to a corresponding system of ordinary differential equations in order to quantitatively test the proposition that the canonical NFκB signaling cascade is sufficient to account for the differential quantitative expression of developmentally-regulated genes in the context of Eda polymorphism (Eda v. EdaTa).
Quantitative gene expression analysis and mechanistic kinetic modeling reveal that the EdaTaallele is associated with significantly downregulated Edar, Shh and Fgf8 expression, and vastly upregulated Egf/Tgfα/Egfr expression. Further, the results revealed that Eda/Edar signaling is not a major determinant of NFκB signaling in normal and mutant Tabby SMG development; rather, TNF is. Moreover, the evidence strongly points to the existence of an important, alternative Eda-initiated transcriptional control pivotal to SMG development. These results are likely to be applicable to other branching organs affected by Eda (EDA) mutation.
The SMG in the adult Tabby mouse is smaller than the wildtype (WT) gland and is characterized by decreased granular convoluted ducts and acini; there is a significant decrease in mucin protein in the hypoplastic Tabby SMGs compared to WT, but not more than expected given the gland size differences [14, 15]. To elucidate the natural history of the pathogenesis, we compared the developmental phenotypes of embryonic (E) days 13–18 Tabby and WT SMGs and relate these to multiple gene expression and putative functional integration of related pathways.
Progressive pathogenesis of Tabby SMGs
From the outset, embryonic Tabby SMGs are smaller, exhibit fewer branches and are developmentally delayed compared to WT glands (Fig. 1). At E13, the Tabby gland is characterized by a solid epithelial stalk ending in a bulb and achieving the Initial Bud stage (Fig. 1A), whereas the WT gland exhibits deep clefts and the formation of 3–6 epithelial buds (branches) (i.e. Late Initial Bud stage) (Fig. 1B). Moreover, although the Pseudoglandular stage, composed of a network of epithelial branches and terminal buds, is achieved in E14 WT SMGs (Fig. 1D), the Tabby gland has only progressed to Late Initial Bud stage on E14 (Fig. 1C) and does not achieve the Pseudoglandular stage until E15 (Fig. 1E). This developmental delay of approximately 1 day persists in Tabby glands, with the Canalicular stage being seen in E15 WT and E16 Tabby glands (compare Fig. 1F to 1G) and the Early Terminal Bud stage being seen in E16 WT and E17 Tabby glands (compare Fig. 1H to 1I, 2A to 2B). The initial detection of immunolocalized mucin protein in Tabby E17 glands, as well as the similarity of its distribution pattern to E16 WT glands, indicates that Tabby proacinar maturation is also delayed, but not precluded (compare Fig. 2A to 2B). At E18, the Tabby glands remain smaller and exhibit fewer presumptive acini than E18 WT glands (compare Fig. 1K to 1L), a hypoplastic phenotype. Nevertheless, both E18 WT and Tabby SMGs display distinct proacinar lumina bounded by a single layer of cuboidal epithelium and filled with mucin protein, as well as continuity between ductal and proacinar lumina (compare Fig. 1K to 1L, and 2C to 2D). This similarity in E18 Tabby and WT phenotypes indicates that terminal differentiation (presumptive functional maturation) is less affected by the EdaTamutation.
A key principle of signaling and transcriptional regulation of development is that a wide array of ontogenic effects can emerge from a relatively small number of individual molecular components . This nonlinear process is dependent on the functional integration of information transmitted by diverse pathways. Many such pathways in SMG ontogeny have been identified (see reviews, [28–31]). The putative functional relationships within and between pathways have been modeled for mouse SMGs as the estimative equivalent of a "café napkin" sketch, namely a diagrammatic genetic network [28, 32].
Mapping genes to their function is called the "genotype-to-phenotype problem," where phenotype is whatever is changed in the organism when a gene's function is altered . Because developmental regulation is both robust and degenerate, it is limiting to simply investigate the average effects of single network genes across samples, large or small . Rather, we are compelled to map genotype to phenotype within the context of the underlying complexity of the networks that regulate cellular functions .
Here we investigated a subcircuit (Fig. 3) that includes 31 probative genes (Tables 1, 2 and 3) in 5 signaling pathways that are both critical to SMG ontogeny and share post-activation downstream targets: Eda, Tnf, Il6, Egf, Fgf [14, 28–31, 35–42]. This approach has proven very useful, not least because it is experimentally constrained and computationally accessible . Below we present novel emergent properties that would not be otherwise evident.
WT gene expression
We determined the sequential time-dependent changes in the expression of 31 genes from the earliest time of immunodetectable Eda and Edar to newborn by quantitative RT-PCR. The presentation of wildtype (WT) gene expression data for embryonic day 16 (E16) to newborn (NB) is found in Table 1. The relative expression ratio (R) is the mean increase or decrease in gene expression in WT glands compared to a single standard, WT E15 glands, the earliest time of Eda and Edar protein expression. These data are particularly heuristic in our attempt to more fully expand upon what is already known about the control of embryonic SMG branching from prior genotype to phenotype analyses [14, 28, 30–32, 37, 38, 41, 42, 44].
Normal Terminal Bud Stage (E16-NB) maturation is critical to subsequent acini and ductal terminal differentiation, postnatally. The steady state gene expression of Fgf8 and Fgfr2, and the more dramatic 20+ fold upregulation of Fgf10, suggests that the Fgf pathways are important to Terminal Bud Stage maturation in the same way they are prior to this stage. This also may be said of the functionally related Shh gene expression. Progression from the Pseudoglandular Stage to the Canalicular Stage, and continued epithelial branching, also includes Eda/Edar, Tnf and Il6 signaling through the canonical NFκB pathway. The importance of these pathways for Terminal Bud Stage maturation appears considerably diminished: Eda, Edar, Edaradd, Tnf, Il6, and Nfkb1 gene expression are significantly (P < 0.01) downregulated. Finally, the Egf/Tgfα/Egfr pathway regulates the rate of branching and histodifferentiation in preparation for progression from the Canalicular Stage to the Terminal Bud Stage. The steady state of Egf expression and the significant upregulation of Tgfα and Egfr suggest the importance of this signaling pathways through Early Terminal Bud Stage (E16); however, beyond this point there is a highly significant (P < 0.01) decline in Egf, Tgfα, and Egfr expression, suggesting a greatly diminished role in SMG maturation.
Tabby gene expression
To determine if the sequential expression of these important 31 genes differs in Tabby SMGs from that seen in WT SMGs, we also compared Tabby E16 to NB to the same standard (WT E15). The presentation of Tabby SMG gene expression data for E16 to NB, and its comparison to WT is found in Tables 2 and 3. It is readily apparent that embryonic EdaTaSMGs have many highly significant (P < 0.01) gene expression differences as compared to WT. Such differences have always to be viewed in the context of progressive SMG pathogenesis (Figs. 1, 2).
Though Eda gene expression in embryonic EdaTaSMGs are 3–19% that of WT, in the Early Terminal Bud Stage there is a 4-fold increase in Nfkb1 gene expression relative to WT (Table 3). This may largely be explained by the concomitant 11-fold increase in Tnf expression, Tnf being a powerful inducer of NFκB activation in SMGs . By Mid to Late Terminal Bud stage, this response is replaced by near normal Nfkb1 and Tnf expression (Table 3), and dramatically upregulated Egf/Tgfα/Egfr expression: Egf by a relative 10–100 fold; Tgfα by a relative 10–30 fold, Egfr by a relative 2.5–5 fold (Table 3). This unexpected upregulation of Egf/Tgfα/Egfr message is reflected in the expression of immunodetectable proteins (Fig. 4), with a substantial increase in Egf, Tgfα, and Egfr proteins being seen in Tabby SMGs as compared to WT glands. At the same time, diminished Eda expression is associated with significantly (P < 0.01) downregulated Shh and Fgf8 expression, an outcome previously shown to be correlated with severe SMG pathology [38, 44]. Shh and Fgf8 are known to positively and reciprocally regulate one another in SMGs . Still, while Shh and Fgf8 expression in WT SMGs is highly correlated (P < 0.01) and the variation of one accounts for nearly all of the variation in the other, in EdaTaSMGs they are not at all correlated (P > 0.90), one accounting for less than 4% of the variation in the other.
In summary, then, the pathology of EdaTaSMGs, relative to the circumscribed signaling subcircuit (Fig. 3), is associated with downregulated Shh and Fgf8 gene expression. Perhaps, as a mark of this subcircuit's robustness, the potential severity of the SMG pathology may be mitigated by vastly upregulated Egf/Tgfα/Egfr gene expression (Table 3) and protein expression (Fig. 4). These results naturally raise questions regarding the underlying mechanism that links the loss-of-function EdaTamutation to the positive and negative regulation of the 5 genes just noted (Shh, Fgf8, Egf, Tgfα, Egfr), and the overall putative relationship to NFκB-mediated gene regulation.
The expression of an individual gene in a developing organ is not a soliloquy; rather, it acts in a chorus of quantitative functional relations appropriately termed a genetic circuit, network or connections map. It is the task of systems biology to quantitatively define and analyze the parts (subcircuits) of the whole, the goal being to put it all together in the future [45–49]. To do this, two effective strategies have emerged. One seeks to infer biologic pathways from large data sets and increasingly powerful software tools for managing and searching literature and pathway databases . The second seeks to construct and test increasingly complex mechanistic models of biologic systems . Here we have combined these strategies to determine the way in which a mechanistic pathway diagram (derived from our prior studies, the literature, and pathway databases) can be programmatically transformed to a corresponding system of ordinary differential equations in order to quantitatively test the proposition that the canonical NFκB signaling cascade is sufficient to account for the differential quantitative expression of developmentally regulated genes in the context of Eda polymorphism (Eda v. EdaTa).
Pathway and literature databases were mined for processes that interrelate a quantitatively measured panel of genes which are critical to SMG development (Tables 1, 2 and 3; Fig. 3). Mined processes were assembled into a mechanistic kinetic framework and the model was tested for consistency with the WT and Tabby quantitative RT-PCR derived expression data. The resulting best-fit mechanistic system diagram contains 138 states (a state being a molecule or a complex in a physical place) in 5 cellular locations (nucleus, cytoplasm, secretory pathway, plasma membrane and extracellular space), and 217 processes (transport, chemical transformation or binding). The system diagram is too large to be displayed legibly on a journal page, but is provided as a scalable PDF file in Additional file 1, along with access to all the equations and parameters for both WT and Tabby data sets, as well as the fits of the quantitative experimental data (see Additional file 2; Additional file 3; Additional file 4; Additional file 5).
Previous systems analyses have utilized computational methods that allow ready deduction of genetic network connectivity and functional properties solely from gene expression data, temporal or not . Still, recent studies in yeast and E. coli compel the caveat that the ratios of protein levels between mutant and WT may not have a one-to-one correlation with those of the corresponding mRNAs [53, 54]. Further, since the mechanistic diagram inevitably includes a large number of proteins as well as expressed mRNAs, we concede at the outset that gene expression data alone are not sufficient to allow highly detailed parameter identification of the downstream protein networks.
Fortunately, the kinetic characteristics of signaling circuits provide general constraints that support quantitative analysis of gene expression data in isolation. First, it is generally agreed that protein transport from the site of synthesis in the ER to the cell surface requires no more than 30 minutes . Second, it is widely acknowledged that signaling pathways require on the order of only two minutes to convey information from the plasma membrane to the nucleus. Since the developmental time courses we analyze include transients on a time scale of 4 to 5 days, we adopted the working hypothesis that signal transduction events (mediated by protein transport in the secretory pathway, receptor binding, and binding to response elements) are fast relative to the dynamics of transcription, mRNA turnover, translation and protein turnover. This simplification allows a full mechanistic test of a complex hypothetical model even in the absence of experimental protein kinetic data on secreted cytokines and growth factors, activated enzymes, and transcription factors (see Methods).
It is readily apparent that Eda is not a major determinant of NFκB signaling. This conclusion is inescapable given the large number of genes in the model (Fig. 3; see Additional file 1) whose transcription is reported to be significantly regulated by canonical NFκB. Given our working assumption that all differences between WT and Tabby developmental mRNA profiles are secondary to differential Eda expression, it must follow that if Eda is the dominant determinant of Eda/Edar signaling, then the Tabby state of many-fold less Eda must propagate in this kinetic model (see Additional file 1) into much reduced expression of genes whose transcription is significantly regulated by NFκB. To the contrary, only 4 (Edar, Fgf8, Shh, Egf) of the 17 genes thought to be NFκB response genes (Fig. 3) have dramatically different time courses in WT and Tabby SMGs (Fig. 5; see Additional file 2). This strongly suggests that Eda is not a major determinant of NFκB activation in SMG branching morphogenesis, and presages the existence of an additional Eda-initiated transcriptional control of the 4 genes: Edar, Fgf8, Shh, Egf.
Concomitantly, it is informative to consider Nfkb1 gene expression itself. The mechanistic model (see Additional file 1) accounts for the observed significant, but transient, decrease in WT and Tabby SMGs (Tables 1 and 2); Nfkb1 expression has essentially the same time course in WT and Tabby SMGs (Fig. 6). Interestingly, the mechanistic model (see Additional file 1; Additional file 3; Additional file 4; Additional file 5) demonstrates that the resulting decrease in protein synthesis is largely buffered by the pre-existing cytoplasmic pool of NFκB1, and thus the downregulation of Nfkb1 has little effect on the propagation of Eda, Fgf10, and Tnf induced signals to the nucleus. Further, since prior SMG studies demonstrated that a significant portion of the activation and nuclear translocation of NFκB can be accounted for by the variation in Tnf signaling , we tested and corroborated the hypothesis that Tnf, not Eda, is the dominant controller of IKK activation in developing WT and Tabby SMGs (see Additional file 1; Additional file 3; Additional file 4; Additional file 5).
To address this question further, we employed an in vitro loss-of-function strategy to determine if loss of canonical NFκB function abrogates Eda-enhanced SMG branching and maturation. We utilized SN50, a cell permeable inhibitor of NFκB translocation into the nucleus , at a concentration (100 μg/ml) previously shown to entirely preclude Tnf- and Eda-induced NFκB nuclear translocation and to significantly inhibit embryonic SMG branching ([32, 57], unpublished). To enhance SMG branching, we used 250 ng/ml Eda-A1 in a manner previously reported . Thus, E14 explants were cultured for 7 days (E14 + 7) in the presence of 250 ng/ml Eda-A1, 100 μg/ml SN50, or 250 ng/ml Eda-A1 + 100 μg/ml SN50; controls consisted of E14 explants cultured in control medium alone.
Eda supplementation induces a substantial increase (Jaskoll et al., 2003; data not shown) and SN50 induces a marked decrease (compare Fig. 7C, D to 7A, B) in embryonic SMG branching morphogenesis and terminal differentiation, as indicated by mucin protein expression. The presence of exogenous Eda supplementation in SN50-treated explants "rescued" the SN50-induced abnormal phenotype and restored it toward that seen in controls (compare Fig. 7E, F to 7A, B). Note that Eda + SN50-treated glands exhibit a marked increase in epithelial terminal buds, lumina formation, and immunodetectable mucin protein compared to glands treated with SN50 alone (compare Fig. 7E, F to 7C, D). Taken together, our results would indicate that a fully functional NFκB pathway is not the sine qua non for Eda signaling, and that Eda signaling likely uses additional, more critical and yet unidentified, pathway(s).
These conclusions were supported by gene expression analysis (Fig. 8). SN50-treated explants show a decline in gene expression of 9 of 10 NFκB1 response genes (Fig. 8A). An optimized (neural network) gene expression model was derived, resulting in a molecular signature that was able to distinguish between SN50-treated and Eda-treated explants with 100% sensitivity and specificity (Fig. 8B). More importantly, when we used the derived model to classify the gene expression data of Eda + SN50-treated explants, we found that all Eda + SN50-treated explants (n = 9) were classified as "Eda-treated" and none were classified as "SN50-treated." This is consistent with the histologic "rescue" (Fig. 7E) and again suggests Eda signaling uses pathways in addition to NFκB.
Within our gene expression datasets (Tables 1, 2 and 3), there are four growth factors and two growth factor receptors whose quantitative expression time courses are dramatically different in Tabby SMGs relative to WT: Egf, Tgfα, Shh, Fgf8, Edar, Egfr (Figs. 5, 9). The mechanistic model (see Additional file 1; Additional file 3; Additional file 4; Additional file 5) is internally consistent with some of these but totally incapable of accounting for others. For Shh and Fgf8, the small differences in nuclear NFκB1 caused by differential Tnf signaling plus the purported feedback of these genes on one another are sufficient to explain the temporal phenotypes of Tabby gene expression reasonably well. For Egf and Edar, the necessity for a second (non-NFκB) Eda-induced transcriptional regulator is clearly indicated by the discrepancies between model solution and mRNA expression data. For Tgfα and Egfr, the model solutions are internally consistent but the link to Eda remains unknown. That is, the model can correctly propagate the different temporal phenotypes of WT and Tabby gene expression, but the model cannot explain the causal chain of events that links these genes to Eda. These discrepancies and missing links compelled us to consider an alternative mechanistic explanation and refine the model.
To wit, we hypothesized and tested the alternative model that a second pathway is activated by Eda/Edar signaling, and that its activated transcription factor ("TFx") regulates the six differently expressed genes in unique ways (Fig. 10). With this addition, all six quantitative gene expression time courses are now internally consistent with and accounted for by the mechanistic model. Further, it is possible to deduce whether the added transcriptional regulation is positive or negative. If gene expression is decreased in Tabby SMGs, TFx must be a positive regulator of that gene. Conversely, if gene expression is increased in Tabby SMGs, then TFx must be a negative regulator of transcription. Examination of the quantitative gene expression data (Tables 1 and 2) reveals that there are three of each: Edar, Shh and Fgf8 are positively regulated by the Eda-activated TFx; Egf, Tgfα, and Egfr are negatively regulated (Fig. 10). The identity of TFx is presently uncertain. However, further in silico, in vivo and in vitro investigations suggest a promising candidate.
Searching for TFx
Regulatory trans-acting transcription factors (TF) activate or repress transcription by physical interaction with genomic cis-regulating DNA elements that may be found in promoters or at some distance from the target gene's start site(s) (see review, ). Putative interactions between TFs and their target DNA sequences can be identified by web-based tools for searching TF binding sites in DNA sequences. Using AliBaba 2.1, TRANSFAC 12.1, MATCH 11.2, and P-MATCH 1.0 (see Methods), we determined that the most probable identity of TFx is CCAAT/enhancer binding protein alpha (C/EBPα). C/EBPα is part of a family of leucine zipper proteins that control the differentiation of many cell types, as well as regulate cell proliferation . Gene and protein analyses are consistent with our in silico result and model prediction (Fig. 10). Quantitative RT-PCR reveals a 60% decline in Cebpa message in Tabby glands relative to WT at E16, and a 25% decline at E17. In addition, there is a dramatic reduction in activated, nuclear-localized C/EBPα protein in Tabby SMGs relative to WT SMGs at E17 (compare Fig. 11B to 11A). To further delineate the relationship between Eda/Edar signaling and C/EBPα, we compared the spatial distribution of C/EBPα protein in Eda-treated and control E14 + 7 explants (compare Fig. 11D to 11C). Enhanced Eda/Edar signaling in vitro induces a notable increase in activated, nuclear-localized C/EBPα protein. Taken together, our in vivo and in vitro data suggest that C/EBPα is an Eda responsive gene, thus a promising TFx candidate.
The ability of C/EBPα to regulate differentiation and proliferation in a context-specific manner often depends on the presence of specific collaborating TFs [59, 60]. Composite cis-regulating elements are combinations of two or more TF binding sites with synergistic regulatory action . Our model predicts that 4 key genes (Shh, Fgf8, Edar, Egf) would contain NFκB/C/EBPα composite elements (Fig. 10). Scanning all known human 5'-flanking sequences, Shelest et al.  found that the most abundant composite elements were of the NFκB/C/EBPα type, including one in the Egf sequence and another in the sequence of a TNF receptor superfamily member. Applying the model for composite elements of the NFκB/C/EBPα type , we determined that such putative composite elements were also present in the murine sequences of Shh, Fgf8, Edar, and Egf.
That C/EBPα is the identity of postulated TFx in the model (Fig. 10) is far from certain. Still, it is an informed hypothesis which can be tested further with protein-protein and protein-DNA interaction assays of high sensitivity and specificity, among other strategies.
In 1875, Charles Darwin presented a meticulous description of what we now know as X-linked Hypohydrotic Ectodermal Dysplasia (XLHED): "I may give an analogous case, communicated to me by Mr. W. Wedderburn of a Hindoo [sic] family in Scinde, in which ten men, in the course of four generations, were furnished, in both jaws taken together, with only four small and weak incisor teeth and with eight posterior molars. The men thus affected have very little hair on their body, and become bald early in life. They also suffer much during hot weather from excessive dryness of their skin. It is remarkable that no instance has occurred of a daughter being affected...though they transmit the tendency to their sons; and no case has occurred of a son transmitting it to his sons" .
Patients with X-linked HED have since been shown to also have SMG hypoplasia and variably reduced salivary secretion [19–21, 63]. The reduced saliva flow results in dryness of the oral mucosa and predisposes XLHED patients to dental caries and Candida albicans infections. Tabby (EdaTa) is a mouse homologue of human XLHED, a model of the human syndrome that displays near identity. Adult EdaTaSMGs are hypoplastic with decreased granular convoluted ducts and acini, and there is a notable decrease in immunodetectable mucin protein [14, 15]. Prior studies of in vitro SMG ontogeny suggested that Eda/Edar signaling effects epithelial cell proliferation, lumen formation, and histodifferentiation via the canonical NFκB pathway . Of great interest, although EdaTaresults in a severely diminished functional Eda ligand [3, 4], hemizygosity (males) and homozygosity (females) for the mutant allele result in variable (mild to moderate) adult (mouse and human) SMG hypoplasia, not aplasia [14, 19–21, 63].
The reproductive success of mus musculus, homo sapiens and other thriving organisms is in no small measure due to their robustness against perturbations, including gene mutation. This robustness may be seen at all levels of biologic organization from gene expression to terminal organ differentiation and function. What might be the basis of "robustian" rescue (partial or full) from perturbations? Several recent lines of experimental evidence suggest that cells, and their genetic regulatory networks, are dynamically critical [64–68]. Such critical dynamical systems, poised between dynamical order and chaos, maximize the correlated behavior of variables in systems of many variables and maximize the diversity of what they can do as they become larger . In the present study, we sought to detect and measure the degree of SMG developmental robustness, and to elucidate its underlying genetic mechanism.
We began by clarifying the in vivo ontogeny of EdaTaassociated SMG pathology (Figs. 1, 2). From the Initial Bud stage (E13) on, Tabby SMGs are smaller, exhibit fewer branches and are developmentally delayed compared to WT glands. This is consistent with the findings that in immunodetected in the Late Pseudoglandular/Early Canalicular stage (~E15) . Of note, though Eda protein is detectable through the Late Terminal Bud Stage (E18-19), comparison of E18 Tabby and WT proacinar phenotypes indicates that terminal differentiation (presumptive functional maturation) is less affected by Eda loss-of-function.
Our prior in vitro study  clearly showed that control glands express very little activated NFκB, but enhanced Eda/Edar signaling induces a very significant increase in activated NFκB. What we did not show was whether or not canonical NFκB activation was necessary or sufficient for in vivo Eda/Edar signaling. We had our doubts because Baltimore's group  had shown that p50-/- mice have normal development, except for some postnatal (6 weeks) defects in immune responses. Further, we investigated the IkBα transgenic mouse created by Schmidt-Ullrich's group [12, 71] and the SMGs are normal (unpublished).
Branching morphogenesis is not a simple dichotomous trait, one we mark present or absent. Rather, it is a complex quantitative trait. Defining the interactions that occur among the genes that underlie the process of branching is essential to understanding its variability. When a gene that is critical to this developmental process mutates, the differentiating cells reprogram transcription via a cognate genetic circuit, ultimately altering the expression of the many genes beyond the mutated one. As such, the phenotype of a given genotype cannot be simply predicted by the sum of its component single-locus effects, but must take account of the almost certain epistasis in gene function . The most efficient way of doing this is by systems analysis, correlation modeling for hypothesis generation (e.g. [32, 73, 74]) and kinetic (mechanistic) modeling for hypothesis testing [46, 48, 49, 75, 76].
The molecular pathology of embryonic Tabby (EdaTa) SMGs is characterized by significantly (P < 0.01) downregulated Shh and Fgf8 gene expression, on the order of 50–75% less than WT. Typically, this would be expected to be associated with a severely hypoplastic to aplastic gland [38, 44], but this is not the case (; Figs. 1, 2). One obvious explanation for this is the many-fold upregulation of Egf/Tgfα/Egfr gene (Table 3) and protein (Fig. 4) expression. This dramatic robustness of the studied subcircuit (Fig. 3; see Additional file 1) is largely due to degeneracy, namely the ready availability of multiple parallel pathways at a key "choice point . Mechanistic modeling reveals that the primary linchpin of this degeneracy is not canonical NFκB as previously supposed ; rather, we provide extensive quantitative evidence that the Eda/Edar/NFκB cascade plays a minor role, apparently neither necessary nor sufficient for SMG development. We confirm this with an in vitro strategy that demonstrates that loss of NFκB function does not abrogate Eda-enhanced SMG branching and differentiation (Figs. 7, 8). These outcomes add support to the suggestion by Pispa et al  of an important NFκB-independent pathway downstream of Eda/Edar signaling in vivo.
The most parsimonious alternative explanation is a second Eda-activated transcription factor ("TFx") that regulates, alone or complexed with NFκB, the expression of all five target genes: Shh, Fgf8, Egf, Tgfα, Egfr (Fig. 10). Our initial in silico and in vivo investigations suggest that the identity of TFx is C/EBPα, a regulator of cell proliferation and differentiation . There are likely other, as yet unknown, TFs in the Eda/Edar pathway. Unmasking TFx with certainty, and delineating its participation in transcription factor complexes, could prove formidable ; delineating the stochastic variation in cognate gene expression more so .
The variable phenotypic expression of human hypohydrotic ectodermal dysplasia is considerable ([22, 80]; Melnick unpublished). This is particularly so for glandular phenotypes (salivary, sweat, sebaceous, lacrimal, mammary, and mucous). If branching epithelia use variations of a canonical developmental program, the experimental results presented here should be applicable to understanding the variable phenotypic expression of other branching organs affected by Eda (EDA) mutation.
Our prior studies of in vitro SMG branching morphogenesis suggested that Eda/Edar signaling largely regulates ontogeny through the canonical NFκB pathway . The present in vivo quantitative systems analyses indicate that this conclusion must be amended. The need to do so is inescapable because, for most NFκB-regulated genes, the observed time course of gene expression is nearly unchanged in Tabby mice as compared to wildtype mice (see Additional file 2), as is NFκB itself (Fig. 6). Importantly, a subset of genes is dramatically differentially expressed in the Tabby mouse (Edar, Fgf8, Shh, Egf, Tgfa, Egfr) (Figs. 5, 9), strongly suggesting the existence of an alternative Eda-mediated transcriptional pathway pivotal for SMG branching morphogenesis (Fig. 10). Experimental and in silico investigations have identified C/EBPα as a promising candidate (Fig. 11). Finally, it should be noted that upregulation of the Egf/Tgfα/Egfr pathway appears to mitigate the potentially severe abnormal phenotype predicted by the downregulation of Fgf8 and Shh. It has recently been suggested by Harris et al. that the buffering of the phenotypic outcome that is coincident with variant Eda signaling could be a common mechanism that permits viable and diverse phenotypes, including those we would consider normal . Our results support this proposition.
Wildtype (WT) mice were either B6CBACaF1-AW-J/A (AW-J) or B10.A/SgSn obtained from Jackson Laboratories (Bar Harbor, ME) and bred as previously described (Jaskoll and Melnick, 1999; Jaskoll et al., 2003); plug day = day 0 of gestation. Tabby breeding pairs [B6CBACA AW-J/A-EdaTa/O/J (Ta/0) and B6CBACA AW-J/A-EdaTa/Y (Ta/Y)] were obtained from Jackson Laboratories and kept by breeding Ta/0 females to Ta/Y males. All embryos from the cross were either Ta/0 or Ta/Ta females and Ta/Y males and displayed the Tabby phenotype. The AW-J wildtype SMGs were used for all comparisons with Tabby glands. E13-19 Tabby and WT pregnant females and newborn mice were sacrificed, SMGs were dissected in cold phosphate-buffered saline (PBS) and glands were collected for morphological, immunolocalization or quantitative RT-PCR. For histological analysis, SMGs were fixed for 4 hrs in Carnoy's fixative at 4°C or overnight in 10% neutral buffered formalin at room temperature, embedded in low melting point paraplast, serially-sectioned at 8 μm and stained with hematoxylin and eosin as previously described . For each embryonic day from gestation day 13 to 18, 5–15 WT and Tabby SMGs were analyzed.
For in vitro culture experiments, B10A/SnSg mice were mated and pregnant females were sacrificed on day 14 of gestation (E14) as previously described [14, 44] Embryonic SMGs were dissected under sterile conditions and in vitro experiments conducted as outlined below.
All animal studies were conducted with the approval of the appropriate committees regulating animal research. An Animal Review Board and a Vivaria Advisory Committee review all applications to ensure ethical and humane treatment.
E16-18 Tabby and WT (AW-J) glands were fixed for 4 hrs in Carnoy's fixative at 4°C or overnight in 10% neutral buffered formalin at room temperature, embedded in low melting point paraplast, serially sectioned at 8 μm and immunolocalization was conducted essentially as previously described [14, 38, 82] using the following affinity-purified antibodies: polyclonal rabbit anti Muc10 [82, 83], polyclonal rabbit anti C/EBPα (sc-61; Santa Cruz Biotechnology), polyclonal rabbit anti Egfr (sc-03; Santa Cruz Biotechnology), polyclonal goat anti Egf (sc-1343; Santa Cruz Biotechnology) and monoclonal mouse anti Tgfα (sc-36; Santa Cruz Biotechnology). In selected experiments, nuclei were stained with 4,6-diamidino-2-phenylindole (DAPI). For mucin protein localization, 3–4 WT and Tabby glands were analyzed for each embryonic day and 3–6 explants per treatment group were analyzed. For C/EBPα protein localization, 5 WT or Tabby glands were analyzed and 3–5 explants per treatment group were analyzed. For Egfr, Egf and Tgfα protein localization, 3–4 WT and Tabby glands were analyzed.
For analysis of gene expression, quantitative RT-PCR was conducted as previously described . WT and Tabby SMGs were pooled (WT-E15-17: 6–20 SMGs/sample; E18-NB: 2–4 SMGs/sample; Tabby-E16-17:15–20 SMGs/sample; E18-NB-4-6 SMGs/sample). For each embryonic day, in each of WT and Tabby, we performed quantitative RT-PCR on 9 or more independent samples. RNA was extracted and 1 μg RNA was reverse transcribed into first strand cDNA using ReactionReady™ First Strand cDNA Synthesis Kit: C-01 for reverse transcription (Superarray Biosciences, Frederick, MD). The primer sets used were prevalidated to give single amplicons and purchased from Superarray Biosciences: Eda (#PPM40603E); Edar (#PPM32386A); Edaradd (PPM32128E); Fgf8 (#PPM02962C); Fgf10 (PPM0345A); Fgfr2 (PPM03706E); PI3k (PPM03469A); Akt (PPM03377A); ERK1 (PPM03575A); ERK2 (PPM03571A); Stat3 (PPM04643E); Egf (PPM03703C); Tgfa (PPM03051A); Egfr (PPM03714E); Ikbkb (PPM03198A); Nfkb1 (#PPM02930A); Nfkb2 (PPM03204A); Rela (#PPM04224E); Relb (#PPM03202C); Tnf (#PPM03113E); Traf2 (PPM03083E); Traf6 (PPM03082A); Il6 (#PPM03015A); p53 (#PPM02931A); Shh (#PPM04516B); cyclin D1 (PPM02903A); Cdk1 (PPM02907A); Casp3 (PPM02922E); Myc (PPM02924A); Ikka (PPM03197A); Ikbkap (PPM37360A); Cebpa (PPM04674A). Primers were used at concentration of 0.4 microM. The cycling parameters were 95°C, 15 min; 40 cycles of (95°C, 15 sec; 55°C, 30–40 sec and 72°C, 30 sec). Specificity of the reactions was determined by subsequent melting curve analysis. RT-PCRs of RNA (not reverse transcribed) were used as negative controls. GAPDH was used to control for equal cDNA inputs and the levels of PCR product were expressed as a function of GAPDH. The relative fold changes of gene expression between the gene of interest and GAPDH, or between the WT and Tabby, were calculated by the 2-ΔΔCTCT method.
Kinetic modeling and hypothesis testing
Pathway and literature databases were mined for processes that interrelate a panel of pivotal genes in SMG development (Fig. 3; Table 1). We extracted relevant interactions and pathways from KEGG http://www.genome.jp, OMIM http://www.ncbi.nlm.nih.gov/sites/entrez?db=omim, GeneCards http://www.genecards.org, three models from the BioModels database http://www.ebi.ac.uk/biomodels/, iHOP http://www.ihop-net.org/UniPub/iHOP/, and the scientific literature accessed through the PubMedCentral database http://www.pubmedcentral.nih.gov/index.html. These data were manually curated and entered into the ProcessDB software http://www.integrativebioinformatics.com. In some cases pathways were available in SBML http://sbml.org/Main_Page format and these were programmatically imported from the BioModels.net database into the ProcessDB database for re-use. The resulting mechanistic system diagram contains 138 states (a state is a molecule or a complex in a physical place) in 5 cellular locations (cytoplasm, nucleus, secretory pathway, plasma membrane and extracellular fluid), and 217 processes (transport, chemical reaction, or binding). It is too large to be displayed legibly on a journal page, but is provided as a scalable PDF file (see Additional file 1).
Kinetic analysis and modeling
Experimental mRNA expression data on the genes listed in Table 1 were collected as described above and provided mRNA data on E15, E16, E17, E18, E19, and newborn (NB) for wild type (AW) embryos and on E16, E17, E18, E19, and NB for Tabby embryos whose glands are uniformly and significantly smaller. For each gene analyzed the data were normalized to the value obtained on E15 in the wild type animals. The defining characteristic of the Tabby mouse is an X-linked mutation that causes a major reduction in Eda mRNA. The underlying goal of the modeling work was to understand how this marked reduction in Eda propagates through the signal transduction and genetic networks to produce the observed temporal phenotypes for the measured panel of developmental genes.
A standard biophysical/bioengineering approach to kinetic analysis of a complex system containing many feedback and feedforward controls is to open the control loops by use of forcing functions [84, 85]. We took advantage of the rich data set of measured mRNA time courses by using them as forcing functions for protein synthesis in the large-scale mechanistic kinetic model. This approach accelerates model development by propagating data-determined (and presumably correct) time courses through the network even before a consistent model has been obtained. Because Tabby glands are too small at ED15 to be practical for RT-PCR-based measurements, some method of extrapolation back from ED16 to ED15 was required for the Tabby forcing functions. Two methods were possible: 1) assume the value at ED15 is equal to the value at ED16 or 2) assume the slope determined by values at ED16 and ED17 is the same as the slope between ED15 and ED16. Choice 1 was seen as more conservative and was implemented. Rate law parameters and initial abundances were chosen throughout the protein network to propagate the main features of each mRNA developmental transient as faithfully as possible. Default rate laws, based on mass action kinetics were created automatically by ProcessDB for each process. This automatically produces saturation behavior for processes limited by the abundance of one or more reactants. Enzyme catalyzed processes and processes activated or inhibited by other states in the model were supplied with rate laws based on classical rapid equilibrium enzyme kinetics .
Another useful modeling technique is, wherever feasible, to analyze multiple experiments simultaneously using the same mechanistic model. This could be implemented in the ProcessDB software, even though the present experiments are performed in two different mouse strains, by adopting the working hypothesis that WT and Tabby mice are identical except for the genetic defect in the Eda gene of the Tabby mutant. In other words, all differences between WT and Tabby developmental mRNA profiles are assumed to be secondary to the difference in Eda expression. In some cases fitting the WT and Tabby data sets required that states have different initial numerical values at ED15. This was allowed because the cytokine environment of the salivary gland in early embryonic development is undoubtedly different in the two strains. In all cases, however, parameter values and rate laws were taken as identical for the two strains. This is an extremely powerful modeling constraint. It does not mean, of course, that all unmeasured protein time courses are the same in the WT and Tabby models. Each is driven by the measured mRNA abundances so that known differences are propagated through common rate laws with the objective of testing the structure of the network – from cytokine production and secretion to receptor binding to activation of enzymes and transcription factors, to nuclear localization and activation or inhibition of cognate genes.
Importantly, the use of mRNA forcing functions in no way guarantees that the full model will simultaneously fit the mRNA data that serve as forcing functions. This is because the control of nuclear transcription factors is, in general, multifactorial and mechanistically distant from synthesis of the proteins involved. Consequently, it was possible to test the mechanistic model's ability to account for the experimental data by searching for transcription and mRNA degradation parameters able to fit the WT and Tabby data simultaneously. The full mechanistic model file including all differential equations, ancillary algebraic equations, rate laws and parameters is included in Additional files (see Additional file 1; Additional file 2; Additional file 3; Additional file 4; Additional file 5).
All differential equations and ancillary algebraic equations were formulated in ProcessDB and exported to the Berkeley Madonna http://www.berkeleymadonna.com solver for numerical integration and parameter optimization using standard methods. Assembled pathway diagrams were treated as hypotheses and tested against the WT and Tabby experimental data.
Transcription factor analysis
The mouse sequences of Fgf8, Shh, Eda, Tgfa, Egf, and Egfr were analyzed for the presence of putative TF DNA binding sites that are common to all 6 genes. We began with an unbiased search with AliBaba 2.1 http://www.generegulation.com/pub/programs.html, the most specific tool for predicting TF binding sites in an "anonymous" DNA sequence using the TRANSFEC database of TFs. The outcome was manually curated based upon mined pathway and literature databases (see Kinetic Modeling/Data mining above) and rank ordered. This analysis revealed SP1 and C/EBPα to be the most likely candidates. We then analyzed both candidates with greater stringency using MATCH 11.2 (proprietary; http://www.biobase-international.com) and P-MATCH 1.0 http://www.gene-regulation.com/pub/programs.html.
MATCH 11.2 is a weight matrix-based tool for searching putative TF binding sites in DNA sequences . MATCH 11.2 uses the matrix library collected in TRANSFAC 12.1 (proprietary; http://www.biobase-international.com). Multiple sets of optimized matrix cut-off values are built into the tool to provide a variety of search modes of different stringency. The matrix similarity is a score (0–1) that describes the quality of a match between matrix and an arbitrary part of the input sequence. Analogously, the core similarity score (0–1) denotes the quality of a match between the core sequence of a matrix (the five most conserved positions within a matrix) and a part of the input sequence. A match has to contain the core sequence of a matrix, i.e. the core sequence has to match with a score higher than or equal to the core similarity cutoff. In addition, only those matches which score higher than or equal to the matrix similarity threshold appear in the output. Cut-offs were chosen to minimize false positive and false negative outcomes. Of the two candidate TFs revealed by AliBaba 2.1, only C/EBPα survived the more rigorous analysis.
This outcome was confirmed using P-MATCH 1.0. P-MATCH combines weight-matrix and pattern matching analytical strategies, thus providing higher accuracy of recognition than either method alone . P-MATCH 1.0 uses the library of mononucleotide weight matrices from TRANSFAC 6.0 (public) along with the site assignments associated with the matrices. Comparisons with MATCH, show that P-MATCH generally provides enhanced recognition accuracy vis-à-vis lower false negative errors (i.e. high sensitivity).
EDA-A1 supplementation in vitro
To determine the optimal concentration of exogenous soluble human recombinant EDA-A1 that induces a significant increase in SMG branching and morphogenesis, we conducted a dose response study. This dose response study was conducted because the EDA-A1 recombinant set (EDA-A1 + enhancer molecule) (Alexis Biochemicals, Axxora, LLC, San Diego, CA) in the present set of experiments was different from the EDA-A1 peptide employed in our previous study . A stock solution of 1 Φg/ml EDA-A1 (10 Φl EDA-A1 + 5 Φl Enhancer + 985 Φl BGJb containing 1% BSA) was made following the manufacturer's protocol and then diluted to 100, 250 and 500 ng/ml. Paired E14 SMG primordia were cultured for 2 days (E14 + 2) and 5–7 days in the presence or absence of 100, 250 and 500 ng/ml EDA-A1. For E14 + 2 explants, Spooner branch ratios (epithelial bud number on day 2/bud number on day 0) were calculated for each explant, comparisons made between right and left glands (treated and control) from each embryo and mean Spooner ratios determined as previously described [14, 38, 44]. The data were arcsin transformed and compared by paired t-test for all embryos studied . Since branching morphogenesis is too complex to count in E14 explants cultured 5–7 days, the morphology of these explants were analyzed by routine hematoxylin and eosin histology; 4–6 explants per dose were analyzed for each day of culture. We determined the optimal dose to be 250 ng/ml EDA-A1 and this concentration was used in all subsequent experiments.
SN50 inhibition of canonical NF6B nuclear translocation
The cell permeable peptide SN50 (Biomol Research, Plymouth Meeting, PA) has been shown to inhibit translocation of the canonical NFκB1/RelA pathway into the nucleus [32, 56, 90]. To determine the effect of SN50 on E14 SMG morphogenesis, paired E14 SMG primordia were cultured for 2 days or 5–7 days in the presence or absence of 100 μg/ml SN50. This concentration was previously shown in our laboratory to be the optimal inhibitor of NFκB1/RelA translocation and interrupts SMG morphogenesis [32, 57]. The morphology of these explants was analyzed by routine hematoxylin and eosin histology; 4–6 explants per treatment group were analyzed.
In vitrorescue experiments
To determine if exogenous EDA-A1 rescues SN50-induced abnormal phenotypes, E14 SMGs were cultured for 7 days in the presence of 250 ng/ml EDA-A1, 100 μg/ml SN50 or 250 ng/ml EDA-A1 + 100 μg/ml SN50; control SMGs consisted of explants cultured in control medium. E14 + 7 explants were collected for histological analysis or quantitative RT-PCR as described above. The morphology of these explants was analyzed by routine hematoxylin and eosin histology; 20–25 explants per treatment group were analyzed. For RT-PCR, 9 independent samples per treatment of pooled glands (20–25 SMGs/sample/group) were analyzed.
Probabilistic neural network (PNN) analysis
We used PNN analyses to determine the contribution of each individual gene to the discrimination between experimental groups with 100% sensitivity and specificity. As such, PNN analyses identify the relative importance (0–1, with 0 being of no relative importance and 1 being relatively most important) of specific gene expression changes that discriminate between phenotypes. It is the contextual change in expression, not the direction of change that is important in defining the molecular phenotype. The foundational algorithm we used is based upon the work of Specht and colleagues [91–93]. The proprietary software designed by Ward Systems Group (Frederick, MD) formulates Specht's procedure around a genetic algorithm . A genetic algorithm is a computational method modeled on biologic evolutionary processes that can be used to find the optimum solution to a problem that may have many solutions . These algorithms have been found to be very powerful in solving optimization problems that appear to be difficult or unsolvable by traditional methods. They use a minimum of information about the problem and they only require a quantitative estimation of the quality of a possible solution. This makes genetic algorithms easy to use and applicable to most optimization problems.
Wessells NK: Tissue Interactions and Development. 1977, Menlo Park: Benjamin/Cummings
Davies JA: Do different branching epithelia use a conserved developmental mechanism?. Bioessays. 2002, 24: 937-48.
Ferguson BM, Brockdorff N, Formstone E, Ngyuen T, Kronmiller JE, Zonana J: Cloning of Tabby, the murine homolog of the human EDA gene: evidence for a membrane-associated protein with a short collagenous domain. Hum Mol Genet. 1997, 6: 1589-94.
Srivastava AK, Pispa J, Hartung AJ, Du Y, Ezer S, Jenks T, Shimada T, Pekkanen M, Mikkola ML, Ko MS, Thesleff I, Kere J, Schlessinger D: The Tabby phenotype is caused by mutation in a mouse homologue of the EDA gene that reveals novel mouse and human exons and encodes a protein (ectodysplasin-A) with collagenous domains. Proc Natl Acad Sci USA. 1997, 94: 13069-74.
Ezer S, Bayés M, Elomaa O, Schlessinger D, Kere J: Ectodysplasin is a collagenous trimeric type II membrane protein with a tumor necrosis factor-like domain and co-localizes with cytoskeletal structures at lateral and apical surfaces of cells. Hum Mol Genet. 1999, 8: 2079-86.
Ezer S, Schlessinger D, Srivastava A, Kere J: Anhidrotic ectodermal dysplasia (EDA) protein expressed in MCF-7 cells associates with cell membrane and induces rounding. Hum Mol Genet. 1997, 6: 1581-7.
Chen TC, Hinton DR, Sippy BD, Hofman FM: Soluble TNF-alpha receptors are constitutively shed and downregulate adhesion molecule expression in malignant gliomas. J Neuropathol Exp Neurol. 1997, 56: 541-50.
Elomaa O, Pulkkinen K, Hannelius U, Mikkola M, Saarialho-Kere U, Kere J: Ectodysplasin is released by proteolytic shedding and binds to the EDAR protein. Hum Mol Genet. 2001, 10: 953-62.
Courtney JM, Blackburn J, Sharpe PT: The Ectodysplasin and NFkappaB signalling pathways in odontogenesis. Arch Oral Biol. 2005, 50: 159-63.
Drew CF, Lin CM, Jiang TX, Blunt G, Mou C, Choung CM, Headon DJ: The Edar subfamily in feather placode formation. Dev Biol. 2007, 305: 232-245.
Mikkola ML: TNF superfamily in skin appendage development. Cytokine Growth Factor Rev. 2008, 19: 219-30.
Schmidt-Ullrich R, Tobin DJ, Lenhard D, Schneider P, Paus R, Scheidereit C: NF-kappaB transmits Eda A1/EdaR signalling to activate Shh and cyclin D1 expression, and controls post-initiation hair placode down growth. Development. 2006, 133: 1045-57.
Smahi A, Courtois G, Rabia SH, Döffinger R, Bodemer C, Munnich A, Casanova JL, Israël A: The NF-kappaB signalling pathway in human diseases: from incontinentia pigmenti to ectodermal dysplasias and immune-deficiency syndromes. Hum Mol Genet. 2002, 11: 2371-5.
Jaskoll T, Zhou YM, Trump G, Melnick M: Ectodysplasin receptor-mediated signaling is essential for embryonic submandibular salivary gland development. Anat Rec A Discov Mol Cell Evol Biol. 2003, 271: 322-31.
Blecher SR, Debertin M, Murphy JS: Pleiotropic effect of Tabby gene on epidermal growth factor-containing cells of mouse submandibular gland. Anat Rec. 1983, 207: 25-9.
Clarke A, Phillips DI, Brown R, Harper PS: Clinical aspects of X-linked hypohidrotic ectodermal dysplasia. Arch Dis Child. 1987, 62: 989-96.
Grüneberg H: Genes and genotypes affecting the teeth of the mouse. J Embryol Exp Morphol. 1965, 14: 137-59.
Grüneberg H: The glandular aspect of the tabby syndrome in the mouse. J Embryol Exp Morphol. 1971, 25: 1-19.
Nordgarden H, Jensen JL, Storhaug K: Oligodontia is associated with extra-oral ectodermal symptoms and low whole salivary flow rates. Oral Dis. 2001, 7: 226-32.
Nordgarden H, Johannessen S, Storhaug K, Jensen JL: Salivary gland involvement in hypohidrotic ectodermal dysplasia. Oral Dis. 1998, 4: 152-4.
Nordgarden H, Storhaug K, Lyngstadaas SP, Jensen JL: Salivary gland function in persons with ectodermal dysplasias. Eur J Oral Sci. 2003, 111: 371-6.
Smith DW: Recognizable Patterns of Human Malformation. 1982, Philadelphia: W.B. Saunders Company, 540-41.
Cui CY, Hashimoto T, Grivennikov SI, Piao Y, Nedospasov SA, Schlessinger D: Ectodysplasin regulates the lymphotoxin-beta pathway for hair differentiation. Proc Natl Acad Sci USA. 2006, 103: 9142-7.
Esibizione D, Cui CY, Schlessinger D: Candidate EDA targets revealed by expression profiling of primary keratinocytes from Tabby mutant mice. Gene. 2008, 427: 42-6.
Hong JW, Hendrix DA, L MS: Shadow enhancers as a source of evolutionary novelty. Science. 2008, 321: 1314-
Wray GA, Babbitt CC: Enhancing gene regulation. Science. 2008, 321: 1300-1301.
Michelson A, Kopan R: Differentiation and gene regulation: toward a holistic understanding of animal development: intercellular communication and transcriptional regulation are two sides of the same coin. Curr Opin Genet Dev. 2002, 12: 499-502.
Jaskoll T, Melnick M: Embryonic Salivary Gland Branching Morphogenesis. Branching Morphogenesis. Edited by: Davies J. 2005, Austin: Eurekah
Melnick M, Jaskoll T: Mouse submandibular gland morphogenesis: a paradigm for embryonic signal processing. Crit Rev Oral Biol Med. 2000, 11: 199-215.
Patel VN, Rebustini IT, Hoffman MP: Salivary gland branching morphogenesis. Differentiation. 2006, 74: 349-64.
Tucker AS: Salivary gland development. Semin Cell Dev Biol. 2007, 18: 237-44.
Melnick M, Chen H, Min Zhou Y, Jaskoll T: The functional genomic response of developing embryonic submandibular glands to NF-kappa B inhibition. BMC Dev Biol. 2001, 1: 15-
Benfey PN, Mitchell-Olds T: From genotype to phenotype: systems biology meets natural variation. Science. 2008, 320: 495-7.
Mitchell KJ: The genetics of brain wiring: from molecule to mind. PLoS Biol. 2007, 5: e113-
Gresik EW, Kashimata M, Kadoya Y, Mathews R, Minami N, Yamashina S: Expression of epidermal growth factor receptor in fetal mouse submandibular gland detected by a biotinyltyramide-based catalyzed signal amplification method. J Histochem Cytochem. 1997, 45: 1651-7.
Jaskoll T, Abichaker G, Witcher D, Sala FG, Bellusci S, Hajihosseini MK, Melnick M: FGF10/FGFR2b signaling plays essential roles during in vivo embryonic submandibular salivary gland morphogenesis. BMC Dev Biol. 2005, 5: 11-
Jaskoll T, Melnick M: Submandibular gland morphogenesis: stage-specific expression of TGF-alpha/EGF, IGF, TGF-beta, TNF, and IL-6 signal transduction in normal embryonic mice and the phenotypic effects of TGF-beta2, TGF-beta3, and EGF-R null mutations. Anat Rec. 1999, 256: 252-68.
Jaskoll T, Witcher D, Toreno L, Bringas P, Moon AM, Melnick M: FGF8 dose-dependent regulation of embryonic submandibular salivary gland morphogenesis. Dev Biol. 2004, 268: 457-69.
Jaskoll T, Zhou YM, Chai Y, Makarenkova HP, Collinson JM, West JD, Hajihosseini MK, Lee J, Melnick M: Embryonic submandibular gland morphogenesis: stage-specific protein localization of FGFs, BMPs, Pax6 and Pax9 in normal mice and abnormal SMG phenotypes in FgfR2-IIIc(+/Delta), BMP7(-/-) and Pax6(-/-) mice. Cells Tissues Organs. 2002, 170: 83-98.
Kashimata M, Gresik EW: Epidermal growth factor system is a physiological regulator of development of the mouse fetal submandibular gland and regulates expression of the alpha6-integrin subunit. Dev Dyn. 1997, 208: 149-61.
Melnick M, Chen H, Zhou Y, Jaskoll T: Embryonic mouse submandibular salivary gland morphogenesis and the TNF/TNF-R1 signal transduction pathway. Anat Rec. 2001, 262: 318-30.
Melnick M, Chen H, Zhou YM, Jaskoll T: Interleukin-6 signaling and embryonic mouse submandibular salivary gland morphogenesis. Cells Tissues Organs. 2001, 168: 233-45.
Neves SR, Iyengar R: Modeling of signaling networks. Bioessays. 2002, 24: 1110-7.
Jaskoll T, Leo T, Witcher D, Ormestad M, Astorga J, Bringas P, Carlsson P, Melnick M: Sonic hedgehog signaling plays an essential role during embryonic salivary gland epithelial branching morphogenesis. Dev Dyn. 2004, 229: 722-32.
Davidson EH, McClay DR, Hood L: Regulatory gene networks and the properties of the development process. Proc Natl Acad Sci. 2003, 100: 1475-1480.
Davidson EH, Rast JP, Oliveri P, Ransick A, Calestani C, Yuh CH, Minokawa T, Amore G, Hinman V, Arenas-Mena C, Otim O, Brown CT, Livi CB, Lee PY, Revilla R, Rust AG, Pan Z, Schilstra MJ, Clarke PJ, Arnone MI, Rowen L, Cameron RA, McClay DR, Hood L, Bolouri H: A genomic regulatory network for development. Science. 2002, 295: 1669-78.
Oliveri P, Tu Q, Davidson EH: Global regulatory logic for specification of an embryonic cell lineage. Proc Natl Acad Sci USA. 2008, 105: 5955-62.
Phair RD, Misteli T: Kinetic modelling approaches to in vivo imaging. Nat Rev Mol Cell Biol. 2001, 2: 898-907.
Smith J, Theodoris C, Davidson EH: A gene regulatory network subcircuit drives a dynamic pattern of gene expression. Science. 2007, 318: 794-7.
Viswanathan GA, Seto J, Patil S, Nudelman G, Sealfon SC: Getting started in biological pathway construction and analysis. PLoS Comput Biol. 2008, 4: e16-
Kitano H: Computational system biology. Nature. 2002, 420: 206-210.
Gardner TS, di Bernardo D, Lorenz D, Collins JJ: Inferring genetic networks and identifying compound mode of action via expression profiling. Science. 2003, 301: 102-5.
Ideker T, Galitski T, Hood L: A new approach to decoding life: systems biology. Annu Rev Genomics Hum Genet. 2001, 2: 343-72.
Lee PS, Shaw LB, Choe LH, Mehra A, Hatzimanikatis V, Lee KH: Insights into the relation between mRNA and protein expression patterns: II. Experimental observations in Escherichia coli. Biotechnol Bioeng. 2003, 84: 834-41.
Hirschberg K, Miller CM, Ellenberg J, Presley JF, Siggia ED, Phair RD, Lippincott-Schwartz J: Kinetic analysis of secretory protein traffic and characterization of golgi to plasma membrane transport intermediates in living cells. J Cell Biol. 1998, 143: 285-503.
Lin YZ, Yao SY, Veach RA, Torgerson TR, Hawiger J: Inhibition of nuclear translocation of transcription factor NF-kappa B by a synthetic peptide containing a cell membrane-permeable motif and nuclear localization sequence. J Biol Chem. 1995, 270: 14255-8.
Melnick M, Mocarski ES, Abichaker G, Huang J, Jaskoll T: Cytomegalovirus-induced embryopathology: mouse submandibular salivary gland epithelial-mesenchymal ontogeny as a model. BMC Dev Biol. 2006, 6: 42-
Walhout AJ: Unraveling transcription regulatory networks by protein-DNA and protein-protein interaction mapping. Genome Res. 2006, 16: 1445-1454.
Nerlov C: The C/EBP family of transcription factors: a paradigm for interaction between gene expression and proliferation control. Trends Cell Biol. 2007, 17: 318-24.
Stein B, Cogswell PC, Baldwin AS: Functional and physical associations between NF-kappa B and C/EBP family members: a Rel domain-bZIP interaction. Mol Cell Biol. 1993, 13: 3964-3974.
Shelest E, Kel AE, Göessling E, Wingender E: Prediction of potential C/EBP/NF-kappaB composite elements using matrix-based search methods. In Silico Biol. 2003, 3: 71-9.
Darwin C: The Variation of Animals and Plants Under Domestication. 1897, New York: D. Appleton Company
Lexner MO, Bardow A, Hertz JM, Almer L, Nauntofte B, Kreiborg S: Whole saliva in X-linked hypohidrotic ectodermal dysplasia. Int J Paediatr Dent. 2007, 17: 155-62.
Balleza E, Alvarez-Buylla ER, Chaos A, Kauffman S, Shmulevich I, Aldana M: Critical dynamics in genetic regulatory networks: examples from four kingdoms. PLoS ONE. 2008, 3: e2456-
Nykter M, Price ND, Aldana M, Ramsey SA, Kauffman SA, Hood LE, Yli-Harja O, Shmulevich I: Gene expression dynamics in the macrophage exhibit criticality. Proc Natl Acad Sci USA. 2008, 105: 1897-900.
Rämö P, Kesseli J, Yli-Harja O: Perturbation avalanches and criticality in gene regulatory networks. J Theor Biol. 2006, 242: 164-70.
Serra R, Villani M, Graudenzi A, Kauffman SA: Why a simple model of genetic regulatory networks describes the distribution of avalanches in gene expression data. J Theor Biol. 2007, 246: 449-60.
Shmulevich I, Kauffman SA, Aldana M: Eukaryotic cells are dynamically ordered or critical but not chaotic. Proc Natl Acad Sci USA. 2005, 102: 13439-44.
Kauffman SA: The origins of order. 1993, New York: Oxford University Press
Sha WC, Liou HC, Tuomanen EI, Baltimore D: Targeted disruption of the p50 subunit of NF-kappa B leads to multifocal defects in immune responses. Cell. 1995, 80: 321-330.
Schmidt-Ullrich R, Aebischer T, Hülsken J, Birchmeier W, Klemm U, Scheidereit C: Requirement of NF-kappaB/Rel for the development of hair follicles and other epidermal appendices. Development. 2001, 128: 3843-53.
Carlborg O, Haley CS: Epistasis: too often neglected in complex trait studies?. Nat Rev Genetics. 2004, 5: 618-25.
Bakal C, Aach J, Church G, Perrimon N: Quantitative morphological signatures define local signaling networks regulating cell morphology. Science. 2007, 316: 1753-6.
Janes KA, Albeck JG, Gaudet S, Sorger PK, Lauffenburger DA, Yaffe MB: A systems model of signaling identifies a molecular basis set for cytokine-induced apoptosis. Science. 2005, 310: 1646-53.
Bar-Yam Y, Harmon D, de Bivort B: Attractors and democratic dynamics. Science. 2009, 323: 1016-1017.
Kearns JD, Hoffmann A: Integrating computational and biochemical studies to explore mechanisms in NF-κB signaling. Journal of Biological Chemistry. 2009, 284: 5493-5443.
Pispa J, Pummila M, Barker PA, Thesleff I, Mikkola ML: Edar and Troy signalling pathways act redundantly to regulate initiation of hair follicle development. Hum Mol Genet. 2008, 17: 3380-91.
Lin LH, Lee HC, Li WH, Chen BS: A systematic approach to detecting transcription factors in response to environmental stresses. BMC Bioinformatics. 2007, 8: 473-
Ansel J, Bottin H, Rodriguez-Beltran C, Damon C, Nagarajan M, Fehrmann S, François J, Yvert G: Cell-to-cell stochastic variation in gene expression is a complex genetic trait. PLoS Genet. 2008, 4: e1000049-
Gorlin RJ, Cohen MM, Levin LS: Syndromes of the Head and Neck. 1990, New York: Oxford University Press
Harris MP, Rohner N, Schwarz H, Perathoner S, Konstantinidis P, Nusslein-Volhard C: Zebrafish eda and edar mutants reveal conserved and ancestral roles of ectodysplasin signaling in vertebrates. PLoS Genet. 2008, 4: e1000206-
Melnick M, Chen H, Zhou Y, Jaskoll T: An alternatively spliced Muc10 glycoprotein ligand for putative L-selectin binding during mouse embryonic submandibular gland morphogenesis. Arch Oral Biol. 2001, 46: 745-57.
Jaskoll T, Chen H, Denny PC, Denny PA, Melnick M: Mouse submandibular gland mucin: embryo-specific mRNA and protein species. Mech Dev. 1998, 74: 179-83.
Berman M: The formulation and testing of models. Ann N Y Acad Sci. 1963, 10: 182-94.
Jacquez JA: Compartmental analysis in biology and medicine. 1996, Ann Arbor, Michigan: BioMedWare Publishing
Segel IH: Enzyme kinetics: behavior and analysis of rapid equlibrium and steady state systems. 1975, New York: John Wiley and Sons Inc
Kel AE, Gössling E, Reuter I, Cheremushkin E, Kel-Margoulis OV, Wingender E: MATCH: A tool for searching transcription factor binding sites in DNA sequences. Nucleic Acids Res. 2003, 31: 3576-3579.
Chekmenev DS, Haid C, Kel AE: P-Match: transcription factor binding site search by combining patterns and weight matrices. Nucleic Acids Res. 2005, 33: W432-7.
Sokal R, Rohlf FJ: Biometry. 1981, New York: Freeman
Liu RY, Fan C, Olashaw NE, Wang X, Zuckerman KS: Tumor necrosis factor-alpha-induced proliferation of human Mo7e leukemic cells occurs via activation of nuclear factor kappaB transcription factor. J Biol Chem. 1999, 274: 13877-85.
Chen C, ed: Fuzzy logic and neural network handbook. 1996, New York: McGraw-Hill
Specht D: Probabilistic neural network for classification, mapping, or associative memory. Proceeding of the IEEE International Conference on Neural Networks. 1988, 1: 525-32.
Specht D, Shapiro P: Generalization accuracy of probabilistic neural networks captured with back-propagation networks. Proceeding of the IEEE International Joint Conference on Neural Networks. 1991, 1: 887-92.
Goldberg DE: Genetic Algorithms in Search, Optimization, and Machine Learning. 1989, Reading, Mass: Addison-Wesley
Holland JH: Adaptation in Natural and Artificial Systems. 1975, Ann Arbor: University of Michigan Press
We would like to thank Pablo Bringas for photographic assistance. This research was supported by NIH grant RO1 DE014535 (TJ/MM).
TJ and MM conceived and designed the study. MM was involved in and coordinated all experiments, as well as the modeling, and drafted the manuscript. TJ participated in RT-PCR analyses, analyzed the histological and localization data, conducted the in vitro experiments, and helped draft the manuscript. RP and SL participated in kinetic modeling and helped draft the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 1: . Diagram of the full mechanistic gene network model. Black arrows represent processes (chemical reactions, transport, or binding). Rectangles represent states. A state is a molecule or complex in a physiologic place. Places are represented by the background ivory or mauve bands of color and are labeled at the bottom of the diagram. Green and red dashed arrows represent, respectively, positive or negative regulation of processes by states. Processes with only starts or ends cross the boundary of the modeled system. This diagram and corresponding computational model were produced by ProcessDB software http://www.integrativebioinformatics.com. (PDF 50 KB)
Additional file 2: . Time course fits of the WT and Tabby experimental data by the full mechanistic model that are not included in Figures 4, 5, 8 of the main text. Note, Eda is not included among the genes because differences in Eda expression in WT and Tabby mice are not attributable to transcriptional control. The vertical axis is the relative abundance of mRNA, as presented in Tables 1 and 2. The lines labeled "WT data" are the quantitative RT-PCR derived mRNA data in wildtype SMGs; the lines labeled "WT model" are model simulated expected mRNA expression for wildtype SMGs. The lines labeled "TA data" are the quantitative RT-PCR derived mRNA data in Tabby SMGs; the lines labeled "TA model" are model simulated expected mRNA expression for Tabby SMGs mice. (PDF 560 KB)
Additional file 3: . Text file displaying the corresponding equations of the computational model shown in Additional file 1. (TXT 71 KB)
A Berkeley Madonna model file modified slightly from the model file generated automatically by ProcessDB
Additional file 4: . A Berkeley Madonna http://www.berkeleymadonna.com model file (MODE3680_20080326_final.mmd) containing all the equations and parameters for both WT and Tabby experiments, as well as the fits of the experimental data. This file can be run and examined using the free trial version of Berkeley Madonna. (MMD 417 KB)
Additional file 5: . Table showing Berkeley Madonna variable names corresponding to those of the model diagram shown in Additional file 1. (XLS 52 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Melnick, M., Phair, R.D., Lapidot, S.A. et al. Salivary gland branching morphogenesis: a quantitative systems analysis of the Eda/Edar/NFκB paradigm. BMC Dev Biol 9, 32 (2009). https://doi.org/10.1186/1471-213X-9-32
- Mucin Protein
- Canalicular Stage
- Fgf8 Gene Expression
- Quantitative System Analysis
- Pseudoglandular Stage