Transcriptional interpretation of Shh morphogen signaling: computational modeling validates empirically established models

In an article published in Development titled ‘A theoretical framework for the regulation of Shh morphogen-controlled gene expression’, the Briscoe laboratory presented a mathematical model of the transcriptional interpretation of graded Shh signaling and bifunctional Gli transcription factors (TFs) in neural patterning (Cohen et al., 2014). This model makes the following major predictions: (1) gene expression boundaries will be shifted as a function of the affinity properties of Gli-binding sites (GBSs) in cis-regulatory modules (CRMs); (2) target gene expression is determined by the combinatorial input of transcriptional effectors which are integrated in individual target gene CRMs. In addition to Gli TFs, such input includes (3) uniformly expressed TFs and (4) morphogen-regulated target genes that are dynamically regulated downstream of the morphogen and which comprise repressor or activator functions.

While overall we do not dispute the theoretical model itself, we are troubled that the study largely overlooks previous experimental work on Shh-regulated CRMs by the Ericson (Oosterveen et al., 2012, 2013) and McMahon laboratories (Peterson et al., 2012), which collectively arrived at conclusions we believe to be conceptually indistinguishable from those of Cohen et al. In these studies, functional analyses of endogenous CRMs firmly establish that transcriptional interpretation of graded Shh signaling and Gli TFs is critically reliant on the cooperative activity of uniformly expressed SoxB1 TFs. Oosterveen et al. (2012) find that SoxB1 and Gli TFs operate synergistically at CRMs, rendering Gli-mediated gene activation a largely concentration-independent event. In a subsequent paper (Oosterveen et al., 2013), it was shown that the activity of SoxB1 can be extended and applied to morphogen signals other than Shh, as well as that morphogen-regulated transcriptional networks underlying neural patterning are functionally recapitulated in limb bud tissue in response to forced SoxB1 expression and morphogen signaling. Moreover, the Oosterveen studies provide direct evidence that morphogen-regulated activator and repressor inputs directly influence the output of Shh-regulated CRMs, and this is also supported by studies from the Matise laboratory (Lei et al., 2006; Wang et al., 2011). Furthermore, a relationship between GBS affinity properties and regional expression of Shh target genes was reported by the Oosterveen et al. (2012) and Peterson et al. (2012) studies.

Cohen et al. (2014) cite the Oosterveen and Peterson studies, but none of the conceptual conclusions presented in those papers, as briefly outlined above, are properly introduced, acknowledged or discussed at any point in their text. Instead, when the concepts of GBS affinity, multiple transcriptional inputs or their integration into Shh-regulated CRMs are being introduced, non-vertebrate model systems are typically discussed and the aforementioned studies either disregarded or misrepresented. Given that the study by Cohen et al. (2014) explicitly models the transcriptional interpretation of Shh signaling in the vertebrate neural tube, we believe that the failure to accredit previously established concepts in the system actually being modeled and only referring to related concepts defined in unrelated morphogen systems is completely inappropriate. This is particularly important, as the conclusions drawn by Cohen et al. (2014) show a high degree of overlap with previously established models of Shh interpretation.

With respect to the role of GBS affinity properties, the authors predict a theoretical ‘neutral point’ around the boundary of Nkx2.2 and Olig2 expression, which results in opposite effects on the range of expression for target genes expressed below or above this point in response to alterations of GBS affinity. This is strikingly reminiscent of the ‘mechanistic differences between local and long-range interpretation of Shh’ outlined by Oosterveen et al. (2012). In that study, by interfering in an unbiased manner with the ability of both GliA and GliR to bind their DNA-binding sites, it was empirically established that genes expressed with a dorsal limit above or below the Nkx2.2-Olig2 boundary interpret Gli input differently. Consequently, this position was defined as discriminating between local and long-range interpretation of Shh signaling (see figs 4 and 7 in Oosterveen et al., 2012). Based on these data, GBS-swapping experiments and other functional analyses of CRMs, it was proposed that genes induced at long range require GliR, together with CRM-specific repressive input, to prevent default (permissive) activation by GliA and SoxB1 proteins at ectopic positions (see figs 3-5,7 in Oosterveen et al., 2012). Accordingly, the Oosterveen study established that genes activated at long range (i.e. above the neutral point) become derepressed when GBS affinity (or GliR) is lowered, in agreement with the Cohen computational model.

By contrast, for local genes, Gli interference experiments resulted in a notable ventral retraction of expression of Nkx2.2, demonstrating that locally restricted Shh target genes critically require instructive input by GliA to overcome default repression (see fig. 4 in Oosterveen et al., 2012). It was further shown that the induction of local genes (Nkx2.2) is less dependent on CRM contextual co-activator input (compared with long-range genes; see fig. 7 in Oosterveen et al., 2012) and that high-affinity GBSs examined in isolation were able to largely recapitulate the expression of endogenous local CRMs in vivo (see fig. 3E in Oosterveen et al., 2012). Based on these and other data, it was proposed that local gene activation is centered around individual high-affinity GBSs associated with these genes, and that the role of GliR is to restrict gene activation by GliA at these sites (see fig. 2B and fig. 3 in Oosterveen et al., 2012). These empirically based conclusions are in line with the effect of reducing Gli affinity for genes below the neutral point, as subsequently outlined by Cohen et al. (2014). Thus, the differential responses of local and long-range genes to Gli proteins and alterations of Gli affinity defined by Oosterveen et al. correspond precisely to the predicted output of genes expressed below and above the neutral point of Cohen et al. Importantly, Cohen and colleagues define the neutral point as being ‘determined only by the concentrations of GliA and GliR and the strength of their cooperative binding with polymerase (Eqn 6) and is independent of the basal level of gene expression’. It is notable that the term ‘basal level of expression’ corresponds to the ‘CRM-specific contextual [non-Gli] input’ described and functionally examined in Oosterveen et al. (2012). In light of this, it seems obvious that the Oosterveen model outlining mechanistic differences regarding short- and long-range interpretation of Shh signaling bears strong conceptual resemblance to most aspects of Cohen's theoretical model, albeit employing different terminology.

Oosterveen et al. (2012) reported a broadly inverse correlation between the range of gene expression and GBS affinity (see fig. 3 in Oosterveen et al., 2012). That said, Oosterveen et al. (2012) did not identify the positive correlation between GBS affinity and regional expression pattern for genes expressed below the neutral point, but this was clearly outlined and empirically validated by a number of criteria, including GBS affinity-swapping experiments in CRMs for Nkx2.2 and Foxa2, in Peterson et al. (2012) (see fig. 5D-G and Discussion of that study). However, when Cohen and colleagues address these studies in their introduction, they incorrectly state ‘…analysis of GBSs within enhancers of Shh target genes failed to find a positive correlation between binding site affinity and range of gene induction (Oosterveen et al., 2012; Peterson et al., 2012)’. This could easily be interpreted to mean that no correlation has been identified, let alone the positive correlation of the Peterson study. Inconsistent with this introductory statement, however, they later affirm that the Oosterveen and Peterson studies describe correlations between GBS affinity and expression.

It is notable that, a year and a half before the present study, Cohen and Briscoe extensively reviewed the findings of the Oosterveen and Peterson studies (Cohen et al., 2013). In this review, they presented a schematic summary of the Gli-affinity data in the Oosterveen and Peterson studies that is remarkably reminiscent of Cohen et al.'s present theoretical model, including the positioning of a presumptive ‘neutral point’ close to the Nkx2.2/Olig2 boundary – apparently without any need for computational modeling (see fig. 1 in Cohen et al., 2013, reproduced as Fig. 1 here). We have not found any reference that the input parameters for the model in Cohen et al. (2014) are based on empirical observations, and these instead appear to have been arbitrarily selected. Thus, given that their modeling so comprehensively matches the existing empirical models that they reviewed in Cohen et al. (2013) (see Fig. 1), we find it astonishing that they consider it irrelevant to discuss or even mention these similarities in their more recent paper (Cohen et al., 2014).

When introducing the concept of multiple transcriptional inputs in CRMs, the authors neglect the rather extensive body of data provided by the Oosterveen and Peterson studies, and instead cite different non-vertebrate morphogen systems, only vaguely stating later in their introduction: ‘The pan-neuronal transcriptional activator Sox2 provides neural specificity to these Shh target genes (Bailey et al., 2006; Graham et al., 2003; Oosterveen et al., 2013)’. By handling previously published work in this manner, the reader is left with the impression that very little is known about the transcriptional interpretation of Shh signaling at the CRM level. Subsequently, without having clearly introduced previously established Shh interpretation models, the authors propose at the end of the introduction, citing only their own work: ‘An alternative is that the dynamics of the transcriptional network, which is composed of Gli proteins, uniformly expressed TFs and TFs downstream of Shh signalling, explains the spatial pattern of gene expression in the neural tube (Balaskas et al., 2012)’. It is noteworthy that the Balaskas study did not include any CRM analyses nor examine uniformly expressed TFs.

Throughout the paper, we noted a tendency for the authors to cite their own work while overlooking contributions by others when outlining network dynamics or well-established repressive interactions amongst morphogen-regulated TFs. One such example relates to the phenomenon of hysteresis, or cellular memory of Shh morphogen exposure, on which both we and the Briscoe laboratory have presented related but distinct models (Lek et al., 2010 and Balaskas et al., 2012, respectively). In Cohen et al. (as well as other publications from the Briscoe laboratory), only the Balaskas study is cited, despite the fact that both describe mechanisms of adaptation in which a cell, once exposed to Shh, will not respond in the same manner to subsequent exposure. Whereas to our knowledge no studies have been published that would disqualify the earlier model by Lek and colleagues, both mechanisms are compatible with recent quantitative analyses by Junker et al. (2014).

In their Discussion, Cohen et al. conclude, ‘Three distinct classes of inputs can be defined: the MR-TF, the activity of which is determined by the distribution of the morphogen in the tissue; uniformly expressed TFs that are active throughout the tissue; and morphogen-controlled target genes that are dynamically regulated downstream of the morphogen. Each of these inputs can comprise multiple individual TFs with either inhibitor or activator function’. This bears strong resemblance to conclusions made in Oosterveen et al. (2013): ‘SoxB1 and Gli proteins therefore appear to define the central node of a neural-specific GRN required to translate graded Shh signaling into regional gene expression patterns…Although SoxB1 and Gli proteins are sufficient to trigger activation of this network…many genes cooperatively activated by SoxB1 and Gli proteins…encode transcriptional activators or repressors that themselves are integral components of the network…Such proteins are likely to act in a more CRM context-dependent manner to influence the regional expression pattern of Shh-regulated genes within the neural tube’. Considering the similarity between these conclusions, we find it remarkable that the Oosterveen and Peterson studies are not cited on a single occasion in their Discussion.

The conceptual conclusions of the Oosterveen studies regarding Shh-regulated transcriptional inputs are partly based on analyses of CRMs of genes that are not directly modelled in Cohen et al. However, we would argue that this does not justify ignoring them, given that Cohen et al. draw conceptually equivalent conclusions and likewise stress the general applicability of their model. In our opinion, the failure to properly accredit work by others risks leaving an impression that the model described in the Cohen study is entirely novel, despite the fact that all of its major conclusions primarily confirm published models of Shh interpretation. By publishing this Correspondence (and the associated response) we hope to protect the integrity of the scientific record by bringing these issues to light. Nevertheless, while we have strived to be factual in our criticism, we obviously cannot be considered to be unbiased and would therefore recommend that engaged readers establish their own opinion on this subject.