(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



System drift in the evolution of plant meristem development [1]

['Pjotr L. Van Der Jagt', 'Sainsbury Laboratory', 'University Of Cambridge', 'Cambridge', 'United Kingdom', 'Department Of Genetics', 'Steven Oud', 'Renske M. A. Vroomans']

Date: 2026-04 

Developmental system drift (DSD) is a process where a phenotypic trait is conserved over evolutionary time, while the genetic basis for the trait changes. DSD has been identified in models with simpler genotype-phenotype maps (GPMs), such as RNA folding, however the extent of DSD in more complex GPMs, such as developmental pattern formation, is debated. To investigate the occurrence of DSD in complex developmental GPMs, we constructed a multi-scale computational model of the evolution of gene regulatory networks (GRNs) governing plant meristem (stem cell niche) development. We found that, during adaptation, some regulatory interactions became essential for the correct expression of stem cell niche genes. These regulatory interactions were subsequently conserved for thousands of generations. Nevertheless, we observed that these deeply conserved regulatory interactions could be lost over an extended period of stabilising evolution. These losses were compensated by changes elsewhere in the GRN, which then became conserved as well. This gain and loss of regulatory interactions resulted in a continual cis-regulatory rewiring in which accumulated changes caused changes in the expression of several genes. Using two publicly available datasets we found frequent changes in conserved non-coding sequences across six evolutionarily divergent plant species, and showed that these changes do not correlate with changes in gene expression patterns, demonstrating the occurrence of DSD. These findings align with the results from our computational model, showing that DSD is pervasive in the evolution of complex developmental systems.

A key open question in evolution of development (evo-devo) is the evolvability of complex phenotypes. Developmental system drift (DSD) contributes to evolvability by exploring different genotypes with similar phenotypic outcome, but with mutational neighbourhoods that have different, potentially adaptive, phenotypes. We investigated the potential for DSD in plant development using a computational model of developmental evolution. We found that the regulatory interactions between genes changed extensively, resulting in the continual rewiring of the gene regulatory network underpinning development. Even regulatory interactions that were essential for correct development were replaced over long evolutionary time scales. Using plant genome and gene expression data from two publicly available datasets, we found high turnover of conserved non-coding sequences, which often contain regulatory sequences, occurring at both short and long time scales. This did not correlate consistently with gene expression changes in plant tissue, supporting the prevalence of DSD as predicted by our model.

Data Availability: The code and scripts for running and analysing the evolutionary simulations, and the bioinformatic analysis can be found at: https://gitlab.developers.cam.ac.uk/slcu/teamrv/publications/vanderjagt_2025 . The publicly available datasets we used for our bioinformatic analysis are the Schuster dataset: https://github.com/schustischuster/evoGE/tree/master , and the Conservatory Project dataset: https://conservatorycns.com/dist/pages/conservatory/analysis.php .

To investigate the potential for DSD in plant SAMs, we developed a computational model of gene regulatory network evolution, extending an evo-devo approach previously applied to animal development [ 36 – 38 ]. In our evolutionary simulations, a small number of regulatory interactions become highly conserved as GRNs evolve to generate a functional tissue pattern. Surprisingly, we found that even these deeply conserved, essential regulatory interactions can diverge over longer evolutionary timescales, resulting in concomitant shifts in gene expression patterns and DSD. To validate these theoretical findings we performed a bio-informatics analysis on two publicly available datasets: one on conserved non-coding sequences (CNSs) in plants [ 39 ], and another on organ-specific RNA expression across plant species at different evolutionary distances [ 40 ]. By combining the CNS data with the cross-species RNA expression, we found that entirely different sets of CNSs can still have similar gene expression patterns, providing empirical evidence for a many-to-one GP mapping underlying plant gene expression and regulatory rewiring. Altogether these findings highlight the prevalence of DSD and its role in shaping developmental evolution in plants.

DSD in plants has received little attention. One multicellular structure which displays evidence of DSD in plants is the shoot apical meristem (SAM), which are multicellular structures containing a stem cell niche that generate all above-ground plant tissues [ 29 ]. Differences in both morphology and gene expression across various vascular plant lineages suggest that SAMs originated independently multiple times before the emergence of the angiosperm clade [ 30 – 33 ]. Within angiosperms however, SAMs are generally regarded as homologous structures due to their structural similarity and strong overlap in associated genes. Some key regulators of SAM stem cells, such as the CLAVATA3/Embryo Surrounding Region-Related (CLE) peptide family, are deeply conserved among land plants [ 34 ]. Nevertheless, the precise CLE peptides governing SAM function differ between angiosperm species, and their expression patterns also vary [ 35 ], suggesting that DSD plays a role in SAM evolution.

However, it remains poorly understood how complex phenotypes are distributed throughout genotype space within highly complex GPMs: do neutral paths continue to percolate through genotype space, or do complex phenotypes occur in isolated genotype islands, limiting the extent of DSD? Theoretical and computational models of high-dimensional GPMs, such as RNA folding [ 13 , 19 – 21 ], protein folding [ 22 , 23 ], gene regulatory networks (GRNs) [ 24 , 25 ] have revealed extensive neutral paths that allow for significant genotype variation while maintaining the same phenotype. For other GPMs, such as those resulting from multicellular development, it has been suggested that complex phenotypes are sparsely distributed in genotype space, and have low potential for DSD because the number of neutral mutations anti-correlates with phenotypic complexity [ 26 , 27 ]. On the other hand, theoretical and experimental studies in nematodes and fruit flies have shown that DSD is present in a phenotypically complex context [ 3 , 11 , 28 ]. It therefore remains debated how much DSD actually occurs in species undergoing multicellular development.

Phenotypic traits are often conserved between related species, even when the developmental process that generates them diverged significantly [ 1 – 4 ]. This phenomenon is called developmental system drift (DSD) or phenogenetic drift [ 5 – 7 ]. It can result from compensatory mutations after adaptive change in a pleiotropic gene [ 8 , 9 ], or from neutral mutations that change the genotype but not the phenotype [ 10 , 11 ]. The potential for DSD in the latter case depends on the number of genotypes resulting in the same phenotype and how they are mutationally connected – aka how long is the neutral path in genotype-phenotype space [ 12 – 14 ]. DSD can drive speciation [ 15 ], may accelerate adaptation [ 16 , 17 ], and is a possible evolutionary mechanism giving rise to the developmental hourglass [ 11 , 18 ].

Results

Model overview We developed a computational model of gene expression evolution in shoot apical meristems (SAMs). We modelled a population of SAMs, each undergoing a developmental process encoded by a heritable genome. A genome consists of a string of genetic elements representing genes (each of which encodes a transcription factor (TF)), and transcription factor binding sites (TFBSs) which regulate the expression of the downstream gene (Fig 1A) [41]. The genome therefore encodes a gene regulatory network (GRN), which governs gene expression in the cells of a two-dimensional tissue representing a longitudinal cut through the SAM (Fig 1B–1D). A subset of TFs exhibit specific properties, such as the ability to diffuse, form dimers, or mediate cell-cell communication with directly neighbouring cells. Gene expression and protein production are subject to a small amount of molecular noise, modelled with stochastic differential equations. Development therefore consists of spatiotemporal changes in protein distribution within the tissue due to gene expression, diffusion and noise (Fig 1D). PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 1. Overview of the computational evo-devo model. (A) A schematic representation of an in silico genome including genes and Transcription Factor Binding Sites (TFBSs). (B) Such a genome can be translated into a gene regulatory network (GRN) which is used to simulate the development. (C) The different functional areas in the simulated tissues. The L1 layer is shown as red-lined cells, the CZ is shown in the tip of the tissue, the OC is in the center, and their regions overlap in a single row of cells. A cell being in the CZ/OC area is determined by their centroid being in the respective bounding box (indicated with dotted lines). (D) An example of the expression pattern of the CZ gene at the end of development in a fit individual. (E) Selection is based on the expression of fitness genes, shown here are examples of a high fitness, medium fitness and low fitness CZ gene expression pattern (left-to-right). (F) The different types of structural genome mutations possible during the simulations. In addition to these mutations, mutations are possible in the binding constants of TFBSs, and the maximum transcription rates of genes. https://doi.org/10.1371/journal.pgen.1012089.g001 The development of each individual begins with one TF uniformly distributed throughout the tissue, while a second diffusible TF is constitutively expressed in the epidermal layer (L1) (Fig 1C). From this initial condition, individuals have a fixed amount of time to express other genes based on the interactions encoded by the genome. At the end of this period, each individual is assigned a fitness score based on the protein concentration of two target genes in specific regions of the SAM: one in the central zone (CZ), and one in the organizing center (OC) (Fig 1C). This fitness score determines the probability that the individual will produce offspring in the next generation (Fig 1E). During reproduction, the parent’s genome is inherited by the offspring with random mutations (Fig 1F). This cycle of development and reproduction with mutation is repeated for 50 000 generations.

Conservation of regulatory interactions during evolution of developmental programs We ran 20 simulations, each with a constant population size of 1 000 individuals. Out of the 20 populations, 15 evolved to correctly express both the OC and CZ genes, resulting in a fitness (f) score of at least 75 out of 100 (Fig 2A, examples of high (≥75) and low (<75) fitness patterns in B). To investigate how the developmental programs evolved that generated the pattern, we tracked how long each regulatory interaction between TFs was conserved in the ancestral lineage of each high-fitness population. Most interactions persisted for only a few hundred to a few thousand generations (Fig 2C), in the same pattern as a control simulation (indicated with C) run with random reproduction and no selection for a pattern. However, a small but significant subset of interactions, which we call “conserved interactions,” remained consistently present for more than 5 000 generations, which did not occur in the control simulation (Fig 2C and 2D). These conserved interactions emerged as individuals in a population achieved higher fitness (Figs 2E and S1, control in S2). This suggests that deeply conserved interactions are a consequence of selection. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 2. Fitness increase is correlated with GRN conservation. (A) Fitness of all individuals in the 20 populations at generation 50 000 (outliers not shown). Dotted line denotes the cutoff fitness of 75. (B) Protein pattern of the two fitness genes for the fittest individuals of simulations 10 and 11, with a fitness of 51.9 and 89.4, respectively. (C) The number of generations each regulatory interaction was conserved in the ancestor trace of a population. The box plot indicated with C is a control simulation without selection, resulting in the absence of highly conserved interactions. Dotted line indicates the ’conservation cutoff’ determined by the maximum conservation time of interactions within the control simulation (7 100 generations), which we rounded down to 5 000 for the rest of this work. (D) The sorted conservation times of all regulatory interactions within the ancestry trace of simulation 2. In red are shown all interactions with a conservation time greater than the conservation cutoff of >5 000 generations. (E) Number of conserved (>5 000 generations) interactions and fitness of ancestor trace for the first 50 000 generations of simulation 2. (F) Fitness distribution of 10 000 randomly generated offspring of the fittest individual from simulation 2 at generation 50 000. Orange: offspring inheriting genomes without any mutations n = 7 503 IQR = 5.69; Blue: offspring with mutation(s) n = 2 497 IQR = 13.24. https://doi.org/10.1371/journal.pgen.1012089.g002 To assess the potential for neutral evolution and DSD after the target expression pattern evolved, we created 10 000 offspring of the highest-fitness individual of population 2 at generation 50 000 (Fig 2F). As mutations are probabilistic, the majority of offspring inherit the parental genome without any mutations. Any variation in fitness of these non-mutated offspring results from variation in cell connectivity in the SAM and noisy gene expression. We found that fitness variation in these non-mutated individuals follows a tight distribution around a high mean fitness, showing that the developmental mechanism is robust to these sources of noise (Fig 2F, orange). In offspring which inherited a genome with mutations, the fitness distribution instead followed a U shape (Fig 2F, blue). The majority of mutations was (near) neutral, indicating a high degree of mutational robustness and redundancy in the GRN, while a smaller set was very deleterious, resulting in a near-zero fitness. Mutations resulting in an ‘in between’ fitness were more rare, consistent with previous findings on fitness landscapes [42]. Overall this shows that the evolved genotypes are both developmentally and mutationally robust which indicates the existence of neutral areas within the GPM around this complex expression phenotype.

Developmental system drift in evolved gene regulatory networks DSD occurs when the phenotype remains conserved along the ancestral lineage while the underlying regulatory architecture changes. Our simulations matched the target phenotype by 30 000 generations, after which they entered a fitness plateau where evolution was mostly driven by stabilising selection (S1 Fig). To study DSD at this fitness plateau, we selected the eight populations which reached the highest maximum fitness, created five clones of each, and evolved these for an additional 50 000 generations (Fig 3A). The median fitness increase remained indistinguishable from background variability during this period (S3 Fig), indicating predominantly neutral evolution (we excluded the clones of 2/8 populations due to more significant fitness increase, S3 Fig, S1 Appendix). PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 3. Neutral evolutionary change in regulatory mechanisms. (A) At generation 50 000, populations that reached a fitness plateau are cloned into 5 separate populations that continue to evolve independently but share a common ancestor before generation 50 000. (B) Fitness distributions of mutated offspring (mutational robustness) of related ancestors from different generations. (C) Fitness distributions of clones (developmental robustness) of related ancestors from different generations. (D) Divergence of interactions between lineages of each clonal population. Divergence is shown for lineages of the full GRNs with that of the common ancestor of all lineages (Full-CA); only the conserved set of interactions with the common ancestor (Core-CA, > 5000 generations); and the conserved set of interactions between the populations (Core). Shown are the medians and IQRs. Divergence between networks is calculated as described in Methods. (E) Functional networks and expression patterns of the five fittest individuals at generation 100 000 from simulations 5 0 to 5 4 , and their common ancestor (CA) at generation 49 800. The top row shows the functional networks, that descend from the CA in different cloned populations. The expression patterns of genes 1,2,8 and 11 are shown to illustrate different levels of phenotypic divergence. For instance, expression of gene 2 diverged significantly between some lineages, whereas expression of gene 1 is very conserved. Different TF types indicated by symbols. https://doi.org/10.1371/journal.pgen.1012089.g003 Next, we investigated whether developmental or mutational robustness increased over this time period, which could explain small fitness gains of the population over longer periods of time. During the fitness plateau, mutational robustness and developmental robustness both fluctuate between generations, without either of their distributions becoming consistently higher or lower compared to ancestral states, indicating drift rather than adaptation of robustness (Fig 3B and 3C; S2 Appendix). We quantified genetic divergence of individuals from their ancestor before the cloning at generation 50 000, using a divergence score based on the adjacency matrices of their GRNs. We found that the full GRNs diverged rapidly from the common ancestor (CA) (Figs 3D and S4, orange line), which was likely due to changes in the redundant or non-functional parts in the GRN, which can be mutated without any phenotypic consequence (see Fig 2F). Therefore, we also measured divergence by including only the conserved interactions in the adjacency matrix (Figs 3D and S4, blue line). Strikingly, these conserved interactions gradually diverged from the CA, indicating turnover despite their initial conservation. Since conserved interactions emerge during adaptation (Fig 2E) and are therefore likely functionally important, their turnover may indicate that the regulatory dynamics that generate the target pattern have changed as well. To assess whether conserved interactions follow similar evolutionary trajectories across independent lineages, we performed a pairwise comparison of evolved GRNs between cloned populations (Figs 3D and S4, green line). We found that conserved regulatory interactions diverged between populations at a similar rate as each lineage diverged from the CA. This indicates that overall, the different lineages follow different paths of divergence from the CA. To investigate whether this divergence can be explained by indirect selection on developmental robustness, we ran simulations without gene expression noise. In these simulations, GRNs still diverge after a fitness plateau has been reached (S5 Fig). Although this does not exclude developmental robustness playing a role in the GRN divergence of the noisy simulations, it does show it is not necessary for divergence. Mutational robustness exhibits a similar drift-like pattern of alternating increases and decreases to that observed in the noisy simulations. To observe the divergence in GRNs more closely, we pruned the full GRNs to remove non-functional and redundant genes and interactions, and compared these pruned GRNs between individuals from different clonal populations, as well as with their CA (Fig 3E). As expected, some regulatory interactions in these networks were highly conserved, e.g., the regulation of gene 3, which is consistently activated by gene 11 and inhibited by gene 1 across all GRNs. In contrast, other interactions were rewired, e.g., the activation of gene 8, which is driven by gene 0 in networks CA, 5 1 , 5 3 but by gene 1 in 5 0 and 5 4 , and gene 5 in 5 2 . Since these GRNs were pruned to eliminate redundancy and all descend from the same CA, differences in their regulatory interactions reflect functional divergence through neutral evolution, indicating the occurrence of DSD. To understand the effects of this regulatory turnover on development, we examined the pattern of the proteins that were not subject to selection for a specific pattern. Indeed, we found that the patterns of some free proteins remained conserved, but the expression of other genes underwent significant change (Figs 3E, S6 and S7). For example, the expression of protein 1 was nearly identical among different populations, while proteins 2, 8, and 11 displayed varying degrees of divergence in pattern compared to the CA. Interestingly, protein 2 had varying patterns, but was absent in the pruned networks of CA, 5 1 , 5 3 . In line with our findings in Fig 3D, this shows that different lineages explored different parts of the neutral evolutionary space, where in some cases TFs were recruited to regulate the expression of the genes under selection, whereas in others the TF remained nonfunctional and accumulated neutral changes. As only 2 out of the 14 genes are under selection for a target pattern, a large number of genes is “free” to evolve, which might contribute to the necessary redundancy for rewiring. Nevertheless, we still observe network divergence in simulations with only 6 instead of 12 “free genes,” showing that redundancies can be created even in more constrained GRNs (S8 Fig). Taken together, we found that DSD can drive functional divergence in the underlying GRN resulting in novel spatial expression dynamics of the genes not directly under selection.

Network redundancy creates space for rewiring To understand how conserved and functional interactions can diverge without disrupting fitness, we traced conserved regulatory interactions over evolutionary time in two of the cloned populations. Consistent with our earlier observations, we found that over time, interactions were lost and new conserved interactions arose (Fig 4A and 4B). To investigate the conditions under which an interaction can be rewired, we examined a single interaction ( ) which is conserved in the lineage of one population but lost in the lineage of a related population, as indicated by the in Fig 4A. We measured in each generation how much its removal affected fitness: an importance value of 1 indicates a total loss of fitness after removal, while a value of 0 indicates no change in fitness. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 4. Rewiring of conserved interactions through loss of function. (A) Presence of conserved interactions in two related lineages, every bar representing a specific conserved interaction. The green bars indicated with a show the interaction 0 activates 8. The lineages are 5 1 and 5 2 , respectively. (B) Importance of conserved interactions that upregulate gene 8 expression over evolutionary time of the respective simulations shown in A. The legend shows the TF regulating 8. Indicated under the plot is the persistence time of each interaction. (C) Importance of conserved interactions upregulating gene 12 (fitness gene). (D) Expression of subset of genes (1,4,7,8) upregulating expression of gene 12 at generations 20 000, 40 000, 60 000, 80 000 and 100 000. Pattern is transparent if there is no interaction between the respective gene and gene 12 at that time point. https://doi.org/10.1371/journal.pgen.1012089.g004 In both lineages, the importance of the interaction ( ) remained high for the first 60 000 generations. However, in lineage 5 2 , a new interaction ( ) emerged that activated gene 8, which coincides with a drop in the importance of . This reduction in importance provided an opportunity for the deletion of without a significant loss of fitness. This example shows the general mechanism by which functional redundancy enables the turnover of regulatory interactions, causing neutral GRN divergence over long evolutionary timescales. We even observed this process among interactions regulating the fitness genes (Fig 4C), which could be rewired multiple times in relatively quick succession. The TF taking over regulation did not necessarily have the same expression pattern as the original: e.g., gene 8, the new regulator of gene 12 at generation 60 000, has a different pattern from gene 4. However, they were both expressed around the center, which is where gene 4’s regulation of gene 12 was active. This shows how the evolution of GRNs exploits overlaps in expression patterns to generate redundancies, allowing the rewiring of regulatory interactions. Finally, we tracked each rewiring event from 4 simulations to investigate more closely how these redundancies emerge in the first place. A general intuition is that gene duplications give rise to redundant but functional copies which can diverge to perform a novel function [43]. We therefore compared the copy number of genes where a conserved interaction was rewired, to the copy number of genes where a non-conserved interaction was rewired; while we would expect that rewiring of conserved interactions is more likely for duplicated genes, we did not find such a bias (S9 Fig).

[END]
---
[1] Url: https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1012089

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/