Evolutionary simulations reveal role for genomic recombination in the evolution of gene regulatory network complexity and robustness

Evolutionary simulations reveal role for genomic recombination in the evolution of gene regulatory network complexity and robustness | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Evolutionary simulations reveal role for genomic recombination in the evolution of gene regulatory network complexity and robustness View ORCID Profile Madison Chapel , View ORCID Profile Carl G de Boer doi: https://doi.org/10.1101/2025.08.28.672878 Madison Chapel 1 Bioinformatics Program, University of British Columbia , Vancouver, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Madison Chapel Carl G de Boer 1 Bioinformatics Program, University of British Columbia , Vancouver, Canada 2 School of Biomedical Engineering, University of British Columbia , Vancouver, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Carl G de Boer For correspondence: carl.deboer{at}ubc.ca Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract The gene regulatory networks (GRNs) of eukaryotes are dramatically more complex than the GRNs of prokaryotes, but we lack a complete picture of the selective pressures that have shaped this difference. Here, we use a biochemically informed model of gene regulation to simulate GRN evolution and explore the role that reproductive strategy plays in shaping regulatory complexity. We find that recombining and non-recombining populations converge to the same level of complexity, even in the absence of selection. However, recombination modifies the rate at which complexity emerges, accelerating convergence to the complexity plateau in changing environments while slowing the process in static environments. Our results suggest that, rather than being under direct selection, regulatory complexity may emerge as a byproduct of other evolutionary processes. These results highlight how reproductive strategy and environmental change interact to influence evolutionary trajectories. Introduction While many of the mechanisms underlying gene expression are similar in eukaryotes and prokaryotes, the topographies of their gene regulatory networks (GRNs) are fundamentally different 1 . Prokaryotes have TFs that are specific enough to bind only a single site in the genome 2 , and the vast majority of genes are regulated by few transcription factors (TFs), resulting in simple GRNs with few connections between TFs and genes. Meanwhile, the vast majority of eukaryotic TFs are much less specific than prokaryotic TFs. Regulation depends on combinatorial binding by multiple TFs to achieve specificity 2 , resulting in much more complex GRNs with many connections between TFs and genes 3 . There are a number of factors that may contribute to the differences in prokaryotic and eukaryotic GRN complexity. One possibility is that these differences are not adaptive, but arose through genetic drift associated with reduced effective population sizes in eukaryotic organisms, as has been proposed for other forms of genomic complexity 4 , 5 . Supporting this idea, simulations have shown that the gain and loss of TF binding sites is a slow process, particularly for longer motifs, unless selection is strong (as it is in large prokaryotic populations), which could prevent eukaryotes from evolving specific motifs (which are necessarily long) 6 . However, adaptive explanations have also been proposed. For example, simulations have revealed that weak cooperative interactions between gene products tend to evolve when organisms have to perform a larger number of biological tasks 7 . As even the simplest eukaryotes have more complex cellular architectures than prokaryotes 8 , complex regulatory networks may have evolved to enable the gene expression profiles necessary to complete these tasks. Finally, eukaryotes and prokaryotes differ in the physical organization of DNA within the cell. In eukaryotes, dense chromatin restricts the number of genomic regions accessible to TFs 9 , although this restriction is not sufficient to enable the same degree of specificity as in prokaryotes 10 , 11 . One possible explanation for differences in GRN complexity that has been underexplored is how differing reproductive strategies between eukaryotes and prokaryotes may influence the complexity of their regulatory networks. Prokaryotes typically reproduce asexually, resulting in clonal lineages, with mutation as the main driver of genetic variation in the population. While some additional variation can be obtained through horizontal gene transfer, most offspring are nearly identical to their parents. Eukaryotes favour sexual reproduction, with each offspring inheriting a mixture of genes from both parents through recombination. Under asexual reproduction, a tight coupling between TFs and their binding sites could evolve. In contrast, such specificity may be strongly disfavoured by a sexual reproductive strategy, as recombination will inevitably combine incompatible TF and binding site alleles. For example, a TF allele from one parent combines with a binding site allele from the other, resulting in dysregulation of the gene in the offspring. In a simple GRN, a differentially acting TF allele would change the expression of its target genes far more than in a complex GRN, where the TF is just one among many contributing to each gene’s expression ( Figure 1A ) . Download figure Open in new tab Figure 1. GRN model and simulation set-up. A) Recombination may introduce alternate TF alleles (i.e., alleles from the other parent). In a simple GRN, a differentially acting TF could alter target gene expression, while in a complex GRN, expression levels would be buffered by the activity of additional TFs. B) GRNs are initialized in one of three states: a high complexity state where each TF binds with minimal affinity (thin arrows), a high complexity state where each TF binds with moderate affinity (heavy arrows), and a low complexity state where a single TF binds with strong affinity while the remainder bind with weak affinity. C) During each generation of the simulation, individuals (squares) are selected with probability proportional to fitness (greyscale intensity). In populations with recombination (green), offspring are a hybrid of two parents, inheriting binding affinities (coloured arrows) from each. In populations without recombination (yellow), offspring are clones of parents. Offspring in both populations are subject to mutation (lightning bolt) before the process repeats. Previous work has explored the relationship between recombination and robustness, the ability to maintain a stable phenotype after perturbations such as mutation or recombination itself 12 . Across numerous studies, recombining populations achieve higher robustness than non-recombining ones 13 – 19 . The underlying mechanisms that drive this increased robustness are yet to be fully understood. Singhal et al. suggested two paths by which robustness may be achieved: increased physical linkage, so that recombination doesn’t disrupt regulatory interactions, or reduced epistasis, so that loci act more independently 13 . How recombination and regulatory complexity interact to influence robustness remains an outstanding question. We hypothesized that if complexity makes populations more robust to changes introduced during recombination, recombining populations should evolve higher complexity than non-recombining populations, and more complex GRNs should exhibit greater robustness than less complex ones. Here, we develop an in silico GRN model and perform evolutionary simulations to explore differences in regulatory complexity and robustness resulting from recombination. We demonstrate that recombining and non-recombining populations converge to similar complexity levels under a variety of conditions, including mutation rates, initial binding affinity strengths, and starting GRN topologies. However, recombining populations consistently demonstrate greater mutational robustness, indicating that robustness can evolve independently of GRN complexity. Finally, we reveal that, although all populations converge to the same end state of complexity, factors such as recombination and environmental changes influence the rate at which complexity emerges in populations. Results Using a biochemically informed GRN model to simulate evolution To explore how recombination influences GRN evolution, we developed a simulation framework ( Methods ). In this simulation, there are two gene types: target genes, the expression of which determines fitness, and TFs, which regulate the target genes and each other, but do not directly impact fitness. The 200 target genes and 20 TFs are randomly arranged on a single linear “chromosome”, and the binding affinities between them are captured in a 220 x 20 matrix. A gene’s expression is the sum of TF binding probabilities (determined by the combined effect of their affinity and expression), each weighted by their effect on expression (+1 for activators, −1 for repressors; n =10 of each). Each target gene has a randomly initialized ‘optimal’ expression; its fitness is normally distributed with respect to the log(expression), centred at this optimum. An individual’s fitness is the product of the fitness measures for all target genes, simulating a system where all target genes are equal and essential. In each generation, individuals are selected with probability proportional to their fitness and allowed to reproduce. Without recombination, offspring are a clonal replicate of a parent genome. With recombination, offspring are a hybrid of both parents, with a single randomly selected recombination site determining which genotypes are passed on. The TF affinities for each gene are inherited as a unit from a single parent, simulating inheritance of a single cis-regulatory region for each gene. For both with and without recombination, the TF affinities of offspring are subject to mutation. A new generation then begins, and the process repeats ( Figure 1C ). Each simulation was done in 10 replicates with a population size of 1,000 across one million generations. Regulatory complexity converges to similar plateaus across conditions Using our GRN simulation framework, we first asked how varying the mutation rate impacted GRN evolution. Higher mutation rates would lead to more genetic diversity between individuals, intensifying the selective pressure for greater recombinational robustness and potentially favouring the evolution of more complex GRNs. Higher mutation rates should also increase the mutational burden at each generation, resulting in different fitness trajectories. We tested both high (∼4 per individual per generation) and low (∼1 per individual per generation) mutation rates. Populations with higher mutation rates consistently plateaued at a lower fitness due to the greater mutational burden ( Figure 2 ), as expected 20 . Further, populations reproducing with recombination consistently reached higher fitness levels and evolved more rapidly than those reproducing without recombination, consistent with their increased ability to eliminate deleterious alleles 21 – 24 . Download figure Open in new tab Figure 2. Mutation rate and recombination alter mutation-selection balance. Fitness (y-axis) shown as mean ± SEM (n = 10 replicates) for 1 million generations (x-axis). Recombining (blue) and non-recombining (yellow) populations were evolved with (A) low or (B) high mutation rates. All populations were initialized with minimal binding affinity. We next tested how GRN initialization strategy influenced GRN evolution ( Methods ). Although previous work 14 initialized networks with random affinities, this approach may produce an initial population with substantial fitness variability. Consequently, early generations are dominated by the fittest individuals, resulting in a large ‘founder effect’ and increased variability in simulation outcomes. To avoid this, we initialized GRNs in three different states: two high complexity states, one with minimal initial binding affinities and one with moderate initial binding affinities for all TF-TF and TF-Target pairs, and a low complexity state where expression of each TF or Target is regulated by a single randomly-selected TF with maximal binding affinity while the remaining TFs have minimal affinity ( Figure 1B , Supplementary Figure S1 ). For sparse affinity initialization, the same initial GRN states were used to seed both recombining and non-recombining populations. Despite converging to similar complexity levels (more on this below), different GRN initialization strategies led to distinct evolutionary trajectories. Here, complexity for each gene is quantified as 1 – the Gini index of TF binding probabilities 3 , with binding probabilities determined by binding affinity and TF concentration 25 . GRN complexity is then calculated as the mean complexity across all genes ( Methods ). Populations initialized in the low complexity state, where each gene is strongly bound by only a single TF, initially maintained low complexity before gradually increasing as additional regulatory interactions accumulate ( Figure 3A ) . A similar trend was observed in populations initialized with minimal affinity GRNs ( Figure 3B ) . Although these began with maximal complexity due to equal input from all TFs, a low complexity state rapidly emerges as strong binding affinities become fixed in the population. From here, complexity increases similarly to the populations initialized in the low complexity state. Download figure Open in new tab Figure 3. Initial binding affinities alter how complexity emerges but similar plateaus are reached. Complexity (y-axis) shown as mean ± SEM (n = 10 replicates) for 1 million generations (x-axis). Recombining (blue) and non-recombining (yellow) populations were initialized with (A) sparse binding affinity, (B) minimal binding affinity, or C) moderate initial affinity. All populations were evolved with high mutation rates. In contrast, the populations initialized with moderate binding affinities followed a slightly different trajectory. Here, all activating and repressing TFs contribute moderately and equally to expression at the beginning of the simulation. Adaptive mutations are partially buffered by other regulating TFs, resulting in a gradual reduction in complexity as some edges are strengthened and others are pruned ( Figure 3C ). This buffering may reduce the individual impact of mutations, resulting in smaller fitness differences between individuals, and causing populations to reach the fitness plateau more slowly than populations initialized with minimal or sparse affinities ( Supplementary Figure S3 . Interestingly, while GRN complexity converged to a similar plateau regardless of GRN initialization, convergence occurred more rapidly in populations without recombination ( Figure 3 , Supplementary Figure S4 ) . Changing environments increase the rate at which complexity plateaus In our previous simulations, expression goals remained constant, allowing populations to optimize their regulatory networks under stable conditions. However, selective pressures are rarely static – environments change over time, populations migrate, and other organisms create competition. Any of these scenarios could change expression optima. We reasoned that a changing environment might result in increased regulatory complexity due to the presence of lingering regulatory interactions from past adaptations. The new expression optimum could be achieved by novel regulatory interactions that counteract existing ones. Additionally, complex regulatory networks may facilitate faster adaptation to shifting conditions through polygenic adaptation 26 , 27 . In a simple GRN, where a target gene is regulated by a single TF, matching a new optimum typically requires the fixation of a large-effect mutation, such as one altering the binding of a TF or one that alters expression of a regulating TF. In contrast, when a gene is influenced by multiple TFs, with each contributing relatively little, selection can act on many loci simultaneously and adaptation can proceed through small changes throughout the network that combine to have greater effects. Recombining populations are particularly well-suited to this strategy, as they can combine independently beneficial variants that may already exist as standing variation within the population 21 – 24 . Accordingly, we next tested whether environmental shifts could promote the evolution of regulatory complexity. We altered our simulation by introducing environmental shifts after populations have adapted to their initial environment (150,000 generations for high mutation rate populations, 300,000 generations for low mutation rate). From then on, the expression optima change for a subset of 20 randomly selected target genes every time the population approaches the fitness plateau ( Supplementary Figure S5 ). As expected, recombining populations adapted to new environments faster than non-recombining ones, requiring an average of 1,112 and 2,289 generations, respectively, to regain fitness after each shift ( Supplementary Figure S6 ). Despite repeated environmental changes, regulatory complexity ultimately plateaued at a similar level across all conditions, stabilizing at around ∼0.55 – the same value observed previously in the static environments. However, for populations that had not yet reached this level when environmental shifts were introduced, the changing environment accelerated the rate at which complexity increased. This was particularly pronounced in recombining populations, where a distinct inflection point can be seen when environmental changes are introduced ( Figure 4B , Supplementary Figure S7 , S8 ). Download figure Open in new tab Figure 4. Complexity emerges more rapidly during environmental changes, especially for recombining populations. Complexity (y-axis) shown as mean ± SEM (n = 10 replicates) up to 1 million generations (x-axis) for populations under static (green) or changing (purple) environmental conditions. Populations were evolved either (A) without recombination or (B) with recombination. All populations were initialized with minimal binding affinities and evolved with low mutation rates. Environmental changes begin at the generations marked by vertical dashed lines. Recombining populations achieve higher robustness We next investigated whether recombination increased robustness, and whether mutational robustness increased with increasing regulatory complexity. To test this, we perturbed individual GRNs by introducing mutations. ‘Robustness’ was measured by determining the mean Euclidean distance between the original expression profile and 100 mutated profiles ( Methods ). Despite converging at similar complexity levels, recombining populations consistently evolved greater mutational robustness than non-recombining ones ( Figure 5 , Supplementary Figure S9 ), as demonstrated in previous studies 14 – 19 . Furthermore, no correlation was found between complexity and robustness (Pearson’s r 2 <0.014 for all populations tested; Supplementary Figure S10 ). The trend of greater robustness in recombining populations was observed regardless of mutation rate or initialization conditions, although these variables did affect the robustness in different ways, especially at earlier generations ( Figure 5 ). Download figure Open in new tab Figure 5. Recombining populations consistently evolve higher robustness than non-recombining. Mutational robustness (y-axis) shown as mean ± SEM (n = 10 replicates) for 1 million generations (x-axis). Populations were evolved with (blue) or without (yellow) recombination. Panels show: (A, C, E) Low mutation rate, (B, D, E) high mutation rate, (A, B) sparse initial affinity, (C, D) minimal initial affinity, (E, F) moderate initial affinity. Download figure Open in new tab Figure 6. Complexity plateaus more rapidly in neutrally evolving populations. Complexity (y-axis) shown as mean ± SEM (n = 10 replicates) for 1 million generations (x-axis). Populations were evolved with (colors) or without (grey) selection. Panels show: (A, D) Sparse initial affinity, (B, E) minimal initial affinity, (C, F) moderate initial affinity, (A, B, C) no recombination, (D, E, F) recombination. All populations were evolved with high mutation rates. Similar regulatory complexity is produced by neutrally evolving populations To determine whether the consistently observed complexity plateau represents an adaptive optimum or a neutral byproduct, we evolved populations without selection. Instead of selecting individuals with probability proportional to their fitness, individuals were selected randomly at each generation. All other parameters (i.e., mutation rate, initial binding affinities, recombination) were matched to previous simulations. These simulations represent an evolutionary null, providing an expectation of regulatory complexity in the absence of selection. While the neutrally evolving populations converge to the same complexity as matched populations under selection, they did so more rapidly than with selection. In static environments, selection thus appears to delay this convergence, an effect that is more pronounced when recombination is present. Discussion In our simulations, recombining populations consistently evolved GRNs that were more robust to mutational perturbation. This finding aligns with previous work demonstrating smaller mutational effects in recombining finite populations 13 – 19 . Our simulations also captured the widely-demonstrated effect of recombining populations adapting to new environments faster than non-recombining ones 28 – 35 . Our simulations suggest that recombination and changing environmental pressures can influence the emergence of regulatory complexity. Reproductive mode alone did not drive the evolution of more complex GRNs, as recombining and non-recombining populations tended to plateau at similar levels of regulatory complexity. Even when selection was removed entirely, complexity still reached the same plateau. However, recombination did influence the rate at which complexity plateaued. In static environments, populations under selection plateaued more slowly than neutrally evolving populations, with this effect being most pronounced in recombining populations. When the environment was dynamic, this effect was reversed; recombination drove populations to the plateau more rapidly than in non-recombining populations. These results suggest that selection does not act directly on complexity. Instead, complexity arises as a byproduct of how recombination and selection shape the exploration of genotype space. In static environments, when populations are near the fitness optimum, most mutations will be neutral or deleterious. Recombination more effectively purges these deleterious variants from the population, maintaining a stable phenotype but also limiting the population’s ability to navigate the space of viable regulatory networks, slowing the accumulation of regulatory complexity. In contrast, when selection pressures change, recombination promotes rapid adaptation by combining beneficial variants from otherwise unfit backgrounds 36 , 37 . In our simulations, complexity emerged most rapidly for recombining populations in changing environmental conditions, suggesting that when selection pressures shift, recombination instead promotes exploration of genotype space. Thus, the changes in complexity may reflect changes in how recombination mediates adaptation rather than selection for complexity. The consistent complexity plateau across our simulations raises further questions: does this represent a biologically meaningful upper limit on regulatory complexity, or is it an artifact of the model’s constraints? Future work should aim to refine GRN models with additional regulatory mechanisms to more accurately reflect biological systems. This, along with testing a wider range of mutation rates, recombination frequencies, and other parameters, will provide further insights into how differences in regulatory complexity emerge and whether complexity itself provides an adaptive benefit. Consistent with previous simulations of gene regulatory evolution, there may be minimal selection on complexity, and instead the complexity of extant regulatory sequences may simply reflect sampling from the space of solutions that satisfies the expression constraints 3 , 38 . Instead, the different complexities of prokaryote and eukaryote GRNs could simply reflect the different gene regulatory hardware they utilize. While we simulated recombining and non-recombining populations with identical hardware, the reality is that prokaryotic TFs tend to be far more specific than eukaryotic TFs 2 . The question then becomes, why do prokaryotes have such specific TFs relative to eukaryotes? While our GRN model draws on established literature to simulate regulatory networks in a biologically informed manner 3 , 25 , several simplifications were made to streamline the design and improve the model’s interpretability. In our model, TFs functioned exclusively as activators or repressors. In reality, TFs can have both activating and repressing activity in different contexts, can be conditionally active only after receiving an appropriate signal 39 , and these activities can change over evolutionary time. Further, the number of TF and Target genes was fixed, preventing regulatory networks from evolving via gene duplication – a process partially responsible for the redundancy in eukaryotic binding motifs 40 . Mutations in TF coding sequences, which were not simulated here, could have even more widespread effects, altering the activity of that TF throughout the genome. To make our simulations computationally tractable, we tested only a single, constant, population size of 1000 individuals. Future simulations may explore whether drift in small populations or bottlenecks of reduced population size influence the emergence of regulatory complexity. Methods Modeling GRNs Binding affinity matrix Individual GRNs were modeled as a matrix of binding affinities bounded between ln( k ) = −5.0 and ln( k ) = 5.0, with 20 TFs and 200 target genes. Half of the TFs were designated as activators and half as repressors. TF and target genes were randomly distributed throughout a single linear chromosome. TFs regulate the expression of all target genes and each other, resulting in a 20 x 220 matrix whose elements ln( k i,j ) describe the binding affinity of the TF product of gene i to the promoter of gene j . Gene expression The binding probability R i,j of TF i to gene j was adapted from Granek and Clarke 25 : Where [TF i ] is the expression level of TF i and k i,j the binding affinity of TF i to gene j . The total regulatory activity R j of all TFs acting on gene j was calculated as the sum of the binding probabilities of all activators minus the sum of the binding probabilities for all repressors: Where Act is the set of all activating TFs while Rep is the set of repressing TFs. The regulatory activity was then passed into a sigmoidal function to determine the expression level E in log space: Expression was calculated iteratively until expression levels for all genes had stabilized (differences of <1 x 10 -4 between iterations). For the first iteration, log expression levels were randomly sampled uniformly between −1 and 1, with the same initialization applied to every individual in the population. Each of the 200 target genes was assigned an optimal expression level, . These values were drawn from a uniform distribution in log space: Fitness The fitness of an individual, F , was calculated as the product of the fitness of all target genes. The fitness of each target gene was defined by how closely its expression level, E j , matched the optimal expression level, , using a Gaussian function with a standard deviation of σ = 0.75. Expression levels were compared in log space: Complexity Complexity of a GRN was measured using the Gini coefficient 3 . While typically used in economics to quantify inequalities in wealth distribution, the Gini coefficient can also model other inequalities within a population, such as the degree of inequality in regulatory interactions of TFs in a regulatory network. A Gini coefficient of 0 represents a distribution where all members of the population have equal wealth, while a coefficient of 1 represents a distribution where all wealth is controlled by a single individual. Complexity of regulatory interactions is measured as 1 – Gini, and ranges from 0 (expression is controlled by a single TF) to 1 (all TFs contribute equally to expression). Complexity, C , for each gene was measured using the binding probabilities, R , sorted in ascending order, of n TFs acting on it: Complexity for each individual was measured as the mean complexity of all genes. Robustness Robustness was measured by introducing mutations to a GRN and measuring changes in expression. Mutation effect sizes were sampled from a Gaussian distribution with standard deviation of 2.0; the same range used for mutations during evolution simulations (see below). The change in expression was measured as the Euclidean distance between the vector of initial expression pattern and the expression pattern after mutation, compared in log space. Robustness, ρ, was measured every 50,000 generations across k = 100 mutations: Evolution simulations Evolution was simulated using a genetic algorithm (GA) implemented with the DEAP Python package 41 . For all simulations, the population size was 1000 and reproduction probability was 0.5. GRNs within a population were initialized in one of three states: a high complexity state with minimal initial binding affinities for all TF-TF and TF-Target pairs, a high complexity state with moderate initial binding affinities, and a low complexity state with sparse affinities where expression of each TF or Target is regulated by a single randomly-selected TF with maximal binding affinity while the remaining TFs have minimal affinity ( Supplementary Figure S1 ) . The moderate binding affinity ( k ) was defined as ln( k ) = 0.0. For the minimal affinity initialization, our goal was to evolve GRNs from a state without regulatory interactions. However, because affinities were represented in log space, it was necessary to determine a sufficiently low cut-off for a minimal binding affinity value. A value of 0 (no affinity) cannot be represented in log space, and choosing too low a value would result in GRNs evolving slowly or not at all (e.g., a mutation from ln( k ) = −10.0 to ln( k ) = −9.0 may not alter expression enough to confer a fitness advantage). The minimal initial binding affinity was set such that initial expression values varied minimally from the basal expression level of ln( E ) = 0.0. To achieve this, we tested a range of binding affinities and selected the largest value that yielded expression levels with a standard deviation ≤ 0.025 in log space. This value depended on the number of TFs in the network, and was determined to be ln( k ) = −(0.517 · ln ( n ) + 2.516), where n is the number of TFs (in our case of 20 TFs, ln( k ) = −4.06) ( Supplementary Figure S2 ) . For the high complexity initialization, a single maximal initial binding affinity was set as ln( k ) = (0.517 · ln ( n ) + 2.516), which, with 20 TFs, is ln( k ) = 4.06. For populations under selection, selection was performed using a roulette function, where the probability of selection is proportional to fitness. For neutrally evolving populations, selection was performed randomly where each individual had equal probability of being selected. For recombining populations, a custom one-point crossover function was used to independently generate two distinct recombinant offspring for each mating pair. Each offspring inherited all binding affinities prior to a randomly selected crossing-over point from one parent, and all binding affinities after the cross-over point from the other parent. The process was repeated with a new cross-over point to generate the second offspring. Mutation effect sizes were sampled from a Gaussian distribution with standard deviation of 2.0, allowing minimal or maximal expression to be achieved within two mutations, as previously observed for yeast promoters 3 . Mutation effects were applied to the binding affinities ( k i , j ) in logarithmic space. The mutation probability was set such that each new individual would receive an expected 4.4 ( p = 1.0 × 10 −3 , high mutation rate) or 1 ( p = 2.2727 × 10 −4 , low mutation rate) mutations. Each GA was run for 1 million generations. A total of 10 replicate populations were evolved for each set of conditions. Optimal expression goals were randomly initialized for each replicate and matched across populations with or without recombination. To simulate selection in changing environmental conditions, simulations were resumed after fitness levels had plateaued (150,000 generations for high mutation rate simulations, 300,000 generations for low mutation rate). The mean fitness of the population was determined as a percentage of the maximum possible fitness. 20 target genes were selected at random, and their optimal expression values shifted by values sampled from a Gaussian distribution with standard deviation of 0.25, simulating a change in the environment and selection for a new expression pattern. New optimal expression values were bounded between ln(−0.9) and ln(0.9) to prevent the emergence of difficult to attain expression goals (e.g., requiring maximum expression for all target genes). When the population’s mean fitness was within a percentage point of the original mean fitness, the environmental shift was repeated. Declaration of interests The authors declare no competing interests. Data availability Code used to run simulations and analyze results can be found at https://github.com/de-Boer-Lab/GRN-evolution/ . Supplementary Figures Download figure Open in new tab Supplementary Figure S1. Sample GRNs in various states. GRNs are represented as a matrix of log binding affinities (colour) for 20 TFs (y-axis) regulating all other TFs and target genes (x-axis). Shown here are GRN log binding affinities (color) initialized with minimal, moderate, and sparse binding affinities, and a GRN that has converged at the complexity plateau. When initialized with sparse binding affinity, a single TF per column binds with maximal affinity while the remainder bind with minimal affinity. Download figure Open in new tab Supplementary Figure S2. Determining minimal binding affinity initialization. (A) Relationship between sufficiently low binding affinity (ln( k )) values (y-axis) and number of TFs in a network (x-axis). A regression line was fit to cut-off values (points) for varying numbers of TFs. (B) Distributions of gene expression values (x-axis) for a network with 20 TFs initialized with varying binding affinities. (C) Relationship between binding affinity strength (x-axis) and standard deviation of expression (y-axis) for a network of 20 TFs. The threshold of σ = 0.025 is indicated with a dashed line, and the largest binding affinity value (ln( k ) = −4.10) not exceeding the threshold is indicated in orange. Download figure Open in new tab Supplementary Figure S3. Fitness trajectories for all simulations in static environments. Fitness (y-axis) shown as mean ± SEM (n = 10 replicates) for 1 million generations (x-axis). Populations were evolved with (blue) or without (yellow) recombination. Panels show: (A, C, E) Low mutation rate, (B, D, E) high mutation rate, (A, B) sparse initial affinity, (C, D) minimal initial affinity, (E, F) moderate initial affinity. Download figure Open in new tab Supplementary Figure S4. Complexity trajectories for all simulations in static environments. Complexity (y-axis) shown as mean ± SEM (n = 10 replicates) for 1 million generations (x-axis). Populations were evolved with (blue) or without (yellow) recombination. Panels show: (A, C, E) Low mutation rate, (B, D, E) high mutation rate, (A, B) sparse initial affinity, (C, D) minimal initial affinity, (E, F) moderate initial affinity. Download figure Open in new tab Supplementary Figure S5. Environmental changes cause fitness to drop and recover. Example of average fitness (y-axis) evolution for a single replicate population during environmental changes, beginning at generation (x-axis) 150,000. Population was initialized with moderate binding affinity and evolved with recombination and a high mutation rate. Download figure Open in new tab Supplementary Figure S6. Recombining populations recover more quickly after environmental changes. Mean generations between environment changes (y-axis) for recombining (green) and non-recombining (yellow) populations across different test conditions (x-axis). Download figure Open in new tab Supplementary Figure S7. Non-recombining populations: changing environments accelerate complexity convergence unless plateau is already reached. Complexity (y-axis) shown as mean ± SEM (n = 10 replicates) up to 1 million generations (x-axis) under static (green) or changing (purple) environments. Panels show: (A, C) Low mutation rate, (B, D) high mutation rate, (A, B) minimal initial affinity, (C, D) moderate initial affinity. Download figure Open in new tab Supplementary Figure S8. Recombining populations: changing environments accelerate complexity convergence unless plateau is already reached. Complexity (y-axis) shown as mean ± SEM (n = 10 replicates) up to 1 million generations (x-axis) under static (green) or changing (purple) environments. Panels show: (A, C) Low mutation rate, (B, D) high mutation rate, (A, B) minimal initial affinity, (C, D) moderate initial affinity. Download figure Open in new tab Supplementary Figure S9. Recombining populations tend to evolve higher robustness than non-recombining populations in changing environments. Mutational robustness (y-axis) shown as mean ± SEM (n = 10 replicates) for 1 million generations (x-axis). Populations were evolved with (blue) or without (yellow) recombination. Panels show: (A, C) Low mutation rate, (B, D) high mutation rate, (A, B) minimal initial affinity, (C, D) moderate initial affinity. Environmental changes begin at the generations marked by vertical dashed lines. Download figure Open in new tab Supplementary Figure S10. Complexity and robustness are uncorrelated. Pearson’s r 2 values (y-axis) between complexity and robustness across all test conditions (x-axis), calculated individually for each population of 1,000 individuals (dots). Correlation was tested every 50,000 generations for each of 10 replicate populations per condition (n = 200 tests). Acknowledgements We thank S. Otto and J. Dennis for helpful discussions. This research was supported by the Natural Sciences and Engineering Research Council of Canada (RGPIN-2020-05425), and the Canadian Institute for Health Research (PJT-180537). M.C. was supported by a UBC IGF. C.G.D. is a Michael Smith Health Research BC Scholar. Funder Information Declared Natural Sciences and Engineering Research Council, https://ror.org/01h531d29 , RGPIN-2020-05425 Canadian Institutes of Health Research, https://ror.org/01gavpb45 , PJT-180537 Footnotes https://github.com/de-Boer-Lab/GRN-evolution/ References 1. ↵ Struhl , K . Fundamentally Different Logic of Gene Regulation in Eukaryotes and Prokaryotes . Cell 98 , 1 – 4 ( 1999 ). OpenUrl CrossRef PubMed Web of Science 2. ↵ Wunderlich , Z. & Mirny , L. A . Different gene regulation strategies revealed by analysis of binding motifs . Trends Genet. TIG 25 , 434 – 440 ( 2009 ). OpenUrl PubMed 3. ↵ Vaishnav , E. D. et al. The evolution, evolvability, and engineering of gene regulatory DNA . Nature 603 , 455 – 463 ( 2022 ). OpenUrl CrossRef PubMed 4. ↵ Lynch , M . The frailty of adaptive hypotheses for the origins of organismal complexity . Proc. Natl. Acad. Sci. U. S. A . 104 , 8597 – 8604 ( 2007 ). OpenUrl Abstract / FREE Full Text 5. ↵ Lynch , M. & Conery , J. S . The Origins of Genome Complexity . Science 302 , 1401 – 1404 ( 2003 ). OpenUrl Abstract / FREE Full Text 6. ↵ Tuğrul , M. , Paixão , T. , Barton , N. H. & Tkačik , G . Dynamics of Transcription Factor Binding Site Evolution . PLOS Genet . 11 , e1005639 ( 2015 ). OpenUrl PubMed 7. ↵ Gao , A. et al. Evolution of weak cooperative interactions for biological specificity . Proc. Natl. Acad. Sci. U. S. A . 115 , E11053 – E11060 ( 2018 ). OpenUrl Abstract / FREE Full Text 8. ↵ Koonin , E. V . The origin and early evolution of eukaryotes in the light of phylogenomics . Genome Biol . 11 , 209 ( 2010 ). OpenUrl CrossRef PubMed 9. ↵ Guertin , M. J. & Lis , J. T . Mechanisms by which transcription factors gain access to target sequence elements in chromatin . Curr. Opin. Genet. Dev . 23 , 116 – 123 ( 2013 ). OpenUrl CrossRef PubMed 10. ↵ Kribelbauer , J. F. , Rastogi , C. , Bussemaker , H. J. & Mann , R. S . Low-Affinity Binding Sites and the Transcription Factor Specificity Paradox in Eukaryotes . Annu. Rev. Cell Dev. Biol . 35 , 357 – 379 ( 2019 ). OpenUrl CrossRef PubMed 11. ↵ Perkins , M. L. , Crocker , J. & Tkačik , G . Chromatin enables precise and scalable gene regulation with factors of limited specificity . Proc. Natl. Acad. Sci . 122 , e2411887121 ( 2025 ). OpenUrl CrossRef PubMed 12. ↵ de Visser , J. A. G. M. et al. Perspective: Evolution and detection of genetic robustness . Evol. Int. J. Org. Evol . 57 , 1959 – 1972 ( 2003 ). OpenUrl 13. ↵ Singhal , S. , Gomez , S. M. & Burch , C. L . Recombination drives the evolution of mutational robustness . Curr. Opin. Syst. Biol . 13 , 142 – 149 ( 2019 ). OpenUrl PubMed 14. ↵ Azevedo , R. B. R. , Lohaus , R. , Srinivasan , S. , Dang , K. K. & Burch , C. L . Sexual reproduction selects for robustness and negative epistasis in artificial gene networks . Nature 440 , 87 – 90 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 15. Misevic , D. , Ofria , C. & Lenski , R. E . Sexual reproduction reshapes the genetic architecture of digital organisms . Proc. R. Soc. B Biol. Sci . 273 , 457 – 464 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 16. Whitlock , A. O. B. , Peck , K. M. , Azevedo , R. B. R. & Burch , C. L . An Evolving Genetic Architecture Interacts with Hill–Robertson Interference to Determine the Benefit of Sex . Genetics 203 , 923 – 936 ( 2016 ). OpenUrl Abstract / FREE Full Text 17. Klug , A. & Krug , J . Conflicting effects of recombination on the evolvability and robustness in neutrally evolving populations . PLoS Comput. Biol . 18 , e1010710 ( 2022 ). OpenUrl PubMed 18. Szöllősi , G. J. & Derényi , I . The effect of recombination on the neutral evolution of genetic robustness . Math. Biosci . 214 , 58 – 62 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 19. ↵ Lohaus , R. , Burch , C. L. & Azevedo , R. B. R . Genetic Architecture and the Evolution of Sex . J. Hered . 101 , S142 – S157 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 20. ↵ Haldane , J. B. S . The rate of spontaneous mutation of a human gene . J. Genet . 31 , 317 – 326 ( 1935 ). OpenUrl CrossRef Web of Science 21. ↵ Felsenstein , J . The Evolutionary Advantage of Recombination . Genetics 78 , 737 – 756 ( 1974 ). OpenUrl Abstract / FREE Full Text 22. Fisher , R. A. The Genetical Theory Of Natural Selection . ( At The Clarendon Press , 1930 ). 23. Muller , H. J . THE RELATION OF RECOMBINATION TO MUTATIONAL ADVANCE . Mutat. Res . 106 , 2 – 9 ( 1964 ). OpenUrl CrossRef PubMed 24. ↵ Muller , H. J . Some Genetic Aspects of Sex . Am. Nat . 66 , 118 – 138 ( 1932 ). OpenUrl CrossRef Web of Science 25. ↵ Granek , J. A. & Clarke , N. D . Explicit equilibrium modeling of transcription-factor binding and gene regulation . Genome Biol . 6 , R87 ( 2005 ). OpenUrl CrossRef PubMed 26. ↵ Pritchard , J. K. , Pickrell , J. K. & Coop , G . The Genetics of Human Adaptation: Hard Sweeps, Soft Sweeps, and Polygenic Adaptation . Curr. Biol . 20 , R208 – R215 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 27. ↵ Milligan , W. R. , Hayward , L. K. & Sella , G. When should adaptation arise from a polygenic response versus few large effect changes ? 2025.05.15.654234 Preprint at doi: 10.1101/2025.05.15.654234 ( 2025 ). OpenUrl Abstract / FREE Full Text 28. ↵ Becks , L. & Agrawal , A. F . The Evolution of Sex Is Favoured During Adaptation to New Environments . PLoS Biol . 10 , e1001317 ( 2012 ). OpenUrl CrossRef PubMed 29. Rice , W. R. & Chippindale , A. K . Sexual Recombination and the Power of Natural Selection . Science 294 , 555 – 559 ( 2001 ). OpenUrl Abstract / FREE Full Text 30. Colegrave , N . Sex releases the speed limit on evolution . Nature 420 , 664 – 666 ( 2002 ). OpenUrl CrossRef PubMed Web of Science 31. Poon , A. & Chao , L . Drift Increases the Advantage of Sex in RNA Bacteriophage Φ6 . Genetics 166 , 19 – 24 ( 2004 ). OpenUrl Abstract / FREE Full Text 32. Goddard , M. R. , Godfray , H. C. J. & Burt , A . Sex increases the efficacy of natural selection in experimental yeast populations . Nature 434 , 636 – 640 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 33. Malmberg , R. L . THE EVOLUTION OF EPISTASIS AND THE ADVANTAGE OF RECOMBINATION IN POPULATIONS OF BACTERIOPHAGE T4 . Genetics 86 , 607 – 621 ( 1977 ). OpenUrl Abstract / FREE Full Text 34. Greig , D. , Borts , R. H. & Louis , E. J . The effect of sex on adaptation to high temperature in heterozygous and homozygous yeast . Proc. R. Soc. Lond. B Biol. Sci . 265 , 1017 – 1023 ( 1998 ). OpenUrl CrossRef PubMed 35. ↵ Luijckx , P. et al. Higher rates of sex evolve during adaptation to more complex environments . Proc. Natl. Acad. Sci . 114 , 534 – 539 ( 2017 ). OpenUrl Abstract / FREE Full Text 36. ↵ Pai , S. V. et al. Sex decreases the pleiotropic costs of local adaptation . 2025.03.28.646016 Preprint at doi: 10.1101/2025.03.28.646016 ( 2025 ). OpenUrl Abstract / FREE Full Text 37. ↵ McDonald , M. J. , Rice , D. P. & Desai , M. M . Sex speeds adaptation by altering the dynamics of molecular evolution . Nature 531 , 233 – 236 ( 2016 ). OpenUrl CrossRef PubMed 38. ↵ He , X. , Duque , T. S. P. C. & Sinha , S . Evolutionary origins of transcription factor binding site clusters . Mol. Biol. Evol . 29 , 1059 – 1070 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 39. ↵ Bondra , E. R. & Rine , J . Context-dependent function of the transcriptional regulator Rap1 in gene silencing and activation in Saccharomyces cerevisiae . Proc. Natl. Acad. Sci . 120 , e2304343120 ( 2023 ). OpenUrl CrossRef PubMed 40. ↵ Rosanova , A. , Colliva , A. , Osella , M. & Caselle , M . Modelling the evolution of transcription factor binding preferences in complex eukaryotes . Sci. Rep . 7 , 7596 ( 2017 ). OpenUrl CrossRef PubMed 41. ↵ Fortin , F.-A. , Rainville , F.-M. D. , Gardner , M.-A. , Parizeau , M. & Gagné , C . DEAP: Evolutionary Algorithms Made Easy . J. Mach. Learn. Res . 13 , 2171 – 2175 ( 2012 ). OpenUrl View the discussion thread. Back to top Previous Next Posted August 28, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Evolutionary simulations reveal role for genomic recombination in the evolution of gene regulatory network complexity and robustness Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Evolutionary simulations reveal role for genomic recombination in the evolution of gene regulatory network complexity and robustness Madison Chapel , Carl G de Boer bioRxiv 2025.08.28.672878; doi: https://doi.org/10.1101/2025.08.28.672878 Share This Article: Copy Citation Tools Evolutionary simulations reveal role for genomic recombination in the evolution of gene regulatory network complexity and robustness Madison Chapel , Carl G de Boer bioRxiv 2025.08.28.672878; doi: https://doi.org/10.1101/2025.08.28.672878 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41913) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13372) Ecology (19889) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40365) Molecular Biology (17163) Neuroscience (88540) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15136) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9818) Zoology (2269)