Modeling the zebrafish segmentation clock's gene regulatory network constrained by expression data suggests evolutionary transitions between oscillating and nonoscillating transcription

Abstract

During segmentation of vertebrate embryos, somites form in accordance with a periodic pattern established by the segmentation clock. In the zebrafish (Danio rerio), the segmentation clock includes six hairy/enhancer of split-related (her/hes) genes, five of which oscillate due to negative autofeedback. The nonoscillating gene hes6 forms the hub of a network of 10 Her/Hes protein dimers, which includes 7 DNA-binding dimers and 4 weak or non-DNA-binding dimers. The balance of dimer species is critical for segmentation clock function, and loss-of-function studies suggest that the her genes have both unique and redundant functions within the clock. However, the precise regulatory interactions underlying the negative feedback loop are unknown. Here, we combine quantitative experimental data, in silico modeling, and a global optimization algorithm to identify a gene regulatory network (GRN) designed to fit measured transcriptional responses to gene knockdown. Surprisingly, we find that hes6, the clock gene that does not oscillate, responds to negative feedback. Consistent with prior in silico analyses, we find that variation in transcription, translation, and degradation rates can mediate the gain and loss of oscillatory behavior for genes regulated by negative feedback. Extending our study, we found that transcription of the nonoscillating Fgf pathway gene sef responds to her/hes perturbation similarly to oscillating her genes. These observations suggest a more extensive underlying regulatory similarity between the zebrafish segmentation clock and the mouse and chick segmentation clocks, which exhibit oscillations of her/hes genes as well as numerous other Notch, Fgf, and Wnt pathway genes.

Keywords: gene regulatory network; her/hes; segmentation clock; simulated annealing; zebrafish.

Figure 1

Design of computational model for her/hes negative feedback network. (A) Schematic showing autoregulation of a generic her gene by Her dimers, including key biochemical parameters. The her genes are targets of Notch signaling, and the Notch intracellular domain (NICD) activates transcription along with the cofactor Su(H). Processes of transcription, translation $\left(\beta \right),$ dimerization $\left(k^{\mathrm{-}},k^{\mathrm{+}}\right),$ and degradation $\left({\delta }^m,\delta ,{\delta }^D\right)$ are modeled explicitly. Maximum transcription rate is given by ${\beta }^m$ and repression by Her/Hes dimers is a function of the dimer concentration, $f\left(D\right).$ Simulations consider reactions within a single cell over time. (B) Example of simulation results showing oscillations of multiple mRNA species over time. For the first half of the simulation the posterior PSM network is used, and the anterior PSM network is used during the second half of the simulation. (C) Schematic showing known Her/Hes dimers and their respective expression domains. In the anterior PSM, her1, her7, and her11 are transcribed, and two DNA-binding dimers can form: Her1/Her1 and Her11/Her11. In the posterior PSM, her1, her7, her12, hes6, and her15 are expressed; six DNA-binding dimers and three non-DNA-binding dimers can form. B and C use the same gene/protein color coding.

Figure 2

Derepression of hes6 transcription in response to her/hes knockdown. (A) Expression of her1 (green) and hes6 (red) in a wild-type 8- to 10-somite stage tailbud. Stripes of her1 are indicative of oscillating expression, while hes6 expression is graded. (B) Response of hes6 transcription rate to all single and pairwise her/hes gene knockdowns. Knockdowns are shown along the x-axis, and the y-axis indicates average fold change of hes6 transcription rate relative to wild type with standard error. Highlighted in pink are the seven treatments resulting in greater than twofold increase in hes6 transcription: her15 MO, her1+hes6 MO, her7+hes6 MO, her11+hes6 MO, her12+hes6 MO, her12+her15 MO, and hes6+her15 MO.

Figure 3

Transcription rate data and validation of simulated annealing. (A) Normalized qPCR data for all six her/hes genes. Each column shows the normalized log-fold change of transcription relative to wild type across all single and pairwise knockdowns. Values are color coded according to the scale on the right. Red boxes indicate greater than average response and blue boxes indicate less than average response, where the average response refers to all her/hes genes within a given knockdown treatment. Note the abundance of red squares for her12 and blue squares for hes6. (B) Validation of simulated annealing, using synthetic networks. Results are shown for 50 trials each of five different in silico networks. Each network was of equal complexity to that of the in vivo network, and noise was added to in silico “qPCR” data. Network predictions are color coded for accuracy and sorted along the x-axis by Jensen–Shannon divergence and along the y-axis by energy score. Here, accuracy corresponds to the proportion of network connections that match the synthetic input GRN. Note the cluster of red and orange (high accuracy) networks near the origin. For the five different networks, the best predictions ranged from 64% to 97% accuracy.

Figure 4

Performance of the predicted GRN for her/hes negative feedback. (A) BioTapestry (Longabaugh et al. 2005) diagram of the predicted negative feedback network showing regulatory connections between all DNA-binding dimers and their targets in either the posterior or the anterior PSM. Thin lines indicate “weak” repression, and thick lines indicate “strong” repression. For clarity, activation of genes and non-DNA-binding dimers are not shown. (B) Comparison of prediction and experiment for three novel triple knockdowns. Pink circles indicate prediction made using the GRN model, and green bars are standard deviations for three experimental replicates. (C) Results of simulation under wild-type conditions showing oscillations of mRNA within a single cell over time. Each gene is color coded as indicated in the key. Note the oscillations of hes6 (pink). (D) The predicted GRN also recapitulates synergy seen between her1+her7 and her1+hes6. mRNA oscillations cease or are highly damped under both knockdown conditions in the posterior and anterior PSM. The only remaining cyclic gene is her11 in the anterior PSM, although these oscillations are still damped compared to the wild-type simulation. mRNA traces are color coded as in C. Some gene traces may not be visible since they overlap completely with others.

Figure 5

Oscillations of hes6 can be minimized by fitting biochemical parameters. (A) Simulation results using one sample parameter set that results in minimized hes6 oscillations. Traces for each mRNA are color coded as indicated in the key on the right. (B) Bar graph showing the coefficient of variation for eight biochemical parameters describing hes6 and repressor strength across all parameter sets (yellow) or only those with minimized hes6 oscillations (pink). Parameters are along the x-axis. Parameters with average values significantly different between the two groups are marked with an asterisk below the x-axis (see also Table S2). (C) Bar graph showing the linear pairwise correlation of qPCR response across all 21 knockdowns between hes6 (pink), sef (green), and axin2 (blue) and the her/hes genes. Significant correlations are marked with an asterisk above the bar. Correlation was calculated using Pearson’s linear correlation on the log-fold change qPCR response. All qPCR values are provided in File S1.