Fig. 1

Experimental overview. Conceptual depiction of the three GRN-based models used in this work: A A pure GRN model where the (two-trait) phenotype is the steady-state concentration of two arbitrary genes. B GRN + Multilinear model, where each phenotypic trait is calculated as the weighted sum of all the elements within the steady-state GRN. C Lattice model, where the phenotype is conceptualised as the steady-state expression pattern of one of the constituent genes (Gen 5 in this example) along a one-dimensional row of cells that can communicate between them through cell–cell signalling. In all of these models, phenotypic variation is created by perturbing one or several elements in the core GRN: Perturbations can be introduced in the strength of gene–gene interaction (i.e. as genetic mutations, D); in some environmental cue that may regulate some environmentally-sensitive gene (E); or in the (maternally inherited) initial concentrations of each GRN element (F). Perturbations on each of these three different sources of phenotypic variation (one element of the GRN perturbed at a time) will produce a collection of two-trait phenotypes (i.e., hind- and fore-limb lengths). If these phenotypes are plotted in a two-trait (T₁-T₂) morphospace, they can reveal the structure of the parameter-to-phenotype maps (D–F, right panels). The linear slopes of these maps can be used as a coarse description of these maps, allowing for map-to-map comparisons of random (Fig. 2) and evolved GRNs (Figs. 3, 4, 5)