Pleiotropy accelerates tooth phenotypic and genomic evolution - An in silico study under the lens of development

Pleiotropy accelerates tooth phenotypic and genomic evolution — An in silico study under the lens of development | 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 Pleiotropy accelerates tooth phenotypic and genomic evolution — An in silico study under the lens of development View ORCID Profile Pascal Felix Hagolani , View ORCID Profile Marie Sémon , View ORCID Profile Guillaume Beslon , View ORCID Profile Sophie Pantalacci doi: https://doi.org/10.1101/2025.04.11.648404 Pascal Felix Hagolani a Laboratoire de Biologie et Modelisation de la Cellule, Ecole Normale Superieure de Lyon, CNRS, UMR 5239, Inserm, U1293, Universite Claude Bernard Lyon 1 , 46 allee d’Italie F-69364 Lyon, France c Centre for Applied Mathematics and Bioinformatics, Department of Mathematics and Natural Sciences , Gulf University for Science and Technology, Kuwait Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Pascal Felix Hagolani For correspondence: pfhagolani{at}gmail.com Marie Sémon a Laboratoire de Biologie et Modelisation de la Cellule, Ecole Normale Superieure de Lyon, CNRS, UMR 5239, Inserm, U1293, Universite Claude Bernard Lyon 1 , 46 allee d’Italie F-69364 Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Marie Sémon Guillaume Beslon b INRIA, Lyon, France CITI, INSA Lyon , Villeurbanne, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Guillaume Beslon Sophie Pantalacci a Laboratoire de Biologie et Modelisation de la Cellule, Ecole Normale Superieure de Lyon, CNRS, UMR 5239, Inserm, U1293, Universite Claude Bernard Lyon 1 , 46 allee d’Italie F-69364 Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sophie Pantalacci Abstract Full Text Info/History Metrics Preview PDF Abstract Pleiotropy, which can occur when a gene affects multiple traits, is a central property of living organisms, influencing their response to mutations and their evolutionary trajectories. Despite many studies and discussions, it remains very difficult to reconcile molecular, developmental and quantitative evolutionary genetics viewpoints on pleiotropy and appreciate how much it puts constraints on genetic evolution, phenotypic evolution and adaptation. Here, we revisit this question by simulating evolution in silico. Our model captures multiple levels of integration observed in complex organisms: from genes to development to phenotype to fitness, additionally allowing to remove pleiotropy to directly test its effect. We focus on the pleiotropic interactions between two organs, specifically teeth. Fitness is determined from the functional interaction between these two teeth, which are produced by an in silico model of tooth morphogenesis. The developmental parameters are produced by a genome consisting of pleiotropic genes which may be influenced by one tooth-specific transcriptional regulator per tooth. Our simulations with and without pleiotropy confirm several acknowledged consequences of pleiotropy. It reduced genetic and phenotypic exploration, and facilitated the rapid accumulation of adaptive mutations in the modular cis-regulatory regions of pleiotropic genes, which are specific to each tooth. Unexpectedly, pleiotropy promoted fitness improvements, morphological complexity, and the accumulation of genetic divergence. Mutations in pleiotropic genes contributed significantly to adaptation, and removing pleiotropy did not increase the proportion of adaptive mutations. Thus pleiotropy does not act as a conservative force, but a channeling force promoting genetic divergence. Significance Statement The notion of “Gene pleiotropy”, just as the notion of “gene”, crosses different fields of biology, but appreciations vary on how frequent it is and how much it favors conservation or divergence at the molecular, developmental or phenotypic scale. One reason is the discrepancies of viewpoints between these three scales, another the difficulty to manipulate pleiotropy experimentally. Here, we performed in silico evolution experiments with a computational toy model focusing on the pleiotropy between two organs while bridging the molecular, developmental and phenotypic scales. Moreover, we compared evolutionary scenarios where the two organs shared or not pleiotropic genes. Our results reconcile observations and show that sharing pleiotropic genes accelerates genetic divergence, even though non-pleiotropic genes diverge faster than pleiotropic genes. Introduction Pleiotropy is the phenomenon whereby a genetic unit (a gene, a mutation) affects multiple traits. It is a central property of the genetic organization of organisms, influencing their response to mutations and their evolutionary trajectories. Fisher’s geometric model, the basis for many studies related to pleiotropy ( 1 – 7 ) predicts that pleiotropy poses such a great disadvantage that adaptation is severely reduced ( 3 ). Since then, many more models have been proposed which moderate this view. They showed that pleiotropy does not scale with complexity as in this original model ( 8 ), and that intermediate levels of pleiotropy could be advantageous ( 9 – 15 ). Despite these efforts, pleiotropy remains a puzzling phenomenon in complex organisms, whose homeostasis and development integrate multiple traits in multiple organs ( 16 ). Pleiotropy has been defined and studied from multiple viewpoints, e.g. evolutionary and functional genomics, quantitative genetics or evolutionary developmental genetics, which are often difficult to reconcile ( 17 – 21 ). Intriguingly, pleiotropy is often considered pervasive in molecular genetics and evolutionary genomics ( 16 ), whereas quantification of pleiotropy at the phenotypic level by quantitative genetics suggest it is limited ( 22 ). These different fields also highlight different consequences of pleiotropy on evolution. In the fields focused on molecular evolution (e.g. genomics, evolutionary or developmental genetics), it is widely assumed that pleiotropy favors genetic conservation ( 23 – 25 ) and modularity enables adaptation ( 26 ). Both pleiotropy and modularity are observed in the functional organization of the genome. Most gene products are used repeatedly in different organs or time points during lifetime, hence most coding sequences are considered pleiotropic, whereas their regulatory sequences are considered more modular, decoupling the functional regulation of gene products. The conservation of coding sequences is positively correlated with the level of gene pleiotropy, which is measured, for example, as a function of the different organs in which a gene is expressed ( 27 ). In contrast, CIS regulatory sequences, which regulate pleiotropic genes, evolve more rapidly and are considered to be the main drivers of phenotypic evolution of complex traits such as morphology ( 26 , 28 , 29 ). However, this view of pleiotropy as a strong constraint is challenged by fields related to quantitative genetics, where pleiotropy, quantified from genetic crosses, is limited. Here, models and experimental findings emphasize that correlated changes between traits can be advantageous and that, in this case, pleiotropic mutations could facilitate adaptation ( 9 , 10 , 12 – 15 ). It has even been argued that pleiotropy could promote genetic divergence by favoring the accumulation of mutations that compensate for each other ( 30 , 31 ). The reason for all these discrepancies is twofold. First, mutations in pleiotropic genes can have limited or even no pleiotropic effects. This may depend on epistatic interactions with other loci, and more generally on the structure of the genotype-phenotype (GP) map ( 32 ). Second, it is not possible to remove pleiotropy experimentally to test its influence. Here we built an in silico evolutionary framework to circumvent these limitations and study how pleiotropy influences adaptation, phenotypic and genotypic evolution. We have developed a complete Genotype-Phenotype-Fitness (GPF) map inspired from real organisms. We use simplified organisms with only two organs — one upper and one lower teeth — that interact to determine fitness. These teeth acquire their shape during development through the same pleiotropic genes (gray genes in Fig. 1A ), but modular regulation is possible thanks to tooth-specific transcriptional regulators (green and orange regulators in Fig. 1A ). After tooth development with an in silico model of tooth morphogenesis, fitness is determined by the way in which tooth shapes fit together, measured by the contact surface between them (occlusion). As detailed below, with this setting, we can focus on the pleiotropy linking the phenotypes of the two teeth, but selection still operates at the integrated level corresponding to the whole system composed of the two interacting teeth. Moreover, mutations in the pleiotropic genes could have either a pleiotropic or a non-pleiotropic effect depending on their interaction with the GP map, including its genetic background context. Download figure Open in new tab Fig. 1. Pleiotropic framework and tooth model. (A) The model includes a simplified genome containing developmental genes (gray rectangles), each with its own constitutive regulation (gray circle). The activity of the developmental genes is further controlled by two Hox-like genes, one for the upper tooth (in green) and one for the lower tooth (in orange). The Hox activity combines Cis and Trans effects. Hence, the genome contains pleiotropic regions (developmental genes, in which a mutation could affect both teeth) and modular regions (Hox-like genes and their target regions, in which a mutation would affect only one tooth) (B) In the with-pleiotropy setting, an individual is defined with 42 parameters (21 for each tooth, grouped in vectors A upper and A lower Sfor the upper and lower teeth respectively). Each parameter is characterized by its constitutive level (grouped in vector in gray, identical for the two teeth) added to the combination of the Cis and Trans effects of the Hox-like genes. Cis effects are local to the developmental genes (hence corresponding to the vectors and for the upper and lower teeth respectively) while Trans effects correspond to the activity of the Hox-like genes (hence corresponding to two real values and respectively). These parameters are combined to compute A upper and A lower . (C) In the without-pleiotropy setting, the constitutive levels are different for the upper and lower teeth ( and respectively). Hence, the values of A upper and A lower are computed by using totally independent (modular) parameters. (D) We use an evolutionary model with a generational scheme. We run the simulations for 200 generations. At each generation the genome is decoded to compute A upper and A lower and the tooth developmental model is used to compute the shape of the upper and lower teeth. Once the teeth have develop, we measure the occlusion between both teeth and select the best 50% of the occluding individuals to contribute to the next generation. During replication all genomic elements (constitutive genes, Hox-like genes and regulating elements) mutate independently (see methods for more details on the mutations rates, developmental model and measurement of occlusion). This framework has several advantages. First, despite its simplicity, it reflects the property of complex organisms, where the genetic independence of body parts is ensured by transcriptional regulators, such as Hox and Hox-related genes. These genes are specifically expressed in certain regions of the body and introduce part-specific regulation of a shared developmental toolkit. We have previously demonstrated that during mouse molar morphogenesis, Pou3f3 is expressed specifically in upper molars, while Nkx2-3 is expressed specifically in lower molars. ( 33 ). Second, our framework reproduces the complexity and degeneracy of the GPF maps, an important ingredient in biological systems. Indeed, teeth shape in this work is the product of a complex GP map, since we used a model of tooth morphogenesis that includes not only gene-gene interactions, but also biomechanical interactions that influence each other during tooth development ( 34 ). This means that this model takes into account many properties of development that can influence evolution, such as non-linearities in the GP map or developmental bias, as demonstrated previously ( 35 – 37 ). Furthermore, in our framework, fitness is determined by the interaction of two tooth shapes. Hence it is an emergent property of the phenotype-fitness (PF) map where nonlinear effects and degeneracy can occur. Third, our framework uses an explicit functional partitioning of the genome into pleiotropic genes versus organ-specific regulators of tooth development and their target region in pleiotropic genes. This way, we can follow independently the effect of mutations of these two functional categories. However, and importantly, some mutations in the pleiotropic genes could be non-pleiotropic, what we call “silent pleiotropy”. They can either have no effect on the phenotype or affect one of the two molars only. For other mutations, pleiotropy can be “effective”: they will have an effect on the phenotype of both molars. Their phenotypic effect on tooth shape and possible advantageous/deleterious/neutral effect on fitness is an emergent property of the GPF map. This allows us to 1) compare the rate of evolution of pleiotropic versus organ-specific genes, but also 2) assess the effect of gene pleiotropy on molecular and phenotypic evolution. Specifically, we can compare evolutionary outcomes of this model where teeth share pleiotropic genes possibly regulated by tooth-specific gene, with a control model where each of the two teeth owns its own set of developmental genes (corresponding to a fully modular setting, with no pleiotropy between the upper and lower teeth). We found that, although gene pleiotropy actually reduces the phenotypic and genotypic variation observed during evolution, it does not hinder evolution, since the collaboration of pleiotropic and modular parts of the genome can enhance the fitness space exploration capabilities through a mechanism akin to the path of least resistance ( 38 , 39 ), resulting in higher accumulation of adaptive mutations. Results Modeling framework to evaluate how the pleiotropy linking two teeth impacts the evolutionary dynamics We wished to develop two teeth from the same genome, but taking into account that genome activity during development differs in the two teeth, resulting in different lower and upper molar phenotypes. The tooth developmental model that is at the core of our GP map includes 21 parameters that interact to form a tooth. Therefore, here we consider that each individual owns a genome with 23 genes ( Fig 1A ): 21 pleiotropic genes, one for each developmental model parameter (colored in gray in Fig. 1A ) and 2 genes, coding Hox-like transcription factors, one for each tooth (in green and orange in Fig. 1A ). The 21 pleiotropic genes are under the influence of both pleiotropic cis-regulatory regions (in gray) and Hox-targeted cis-regulatory regions (in green and orange in Fig. 1A ). During development, the Hox genes are specifically expressed in one tooth or the other, therefore the activity of the 21 genes differs in the two teeth. Different part of the genome contribute to these activities: the pleiotropic regions of the 21 genes, the modular cis-regulatory sequences targeted by the Hox gene, and the Hox gene itself. These different regions of the genome mutate with different rates, proportional to the fraction they would represent in nature: large for the pleiotropic part, smaller for the Hox-targeted cis-regulatory regions, very small for the Hox gene. In order to reflect this partitioning and function of the genome, we model the activity of the 21 genes in each tooth as two vectors A upper and A lower for upper and lower teeth, respectively. Each of the 21 developmental genes has a specific constitutive activity (all grouped in the vector ) and is additively regulated by a product of a Hox gene activity in trans ( or ) with a gene-specific cis-regulatory vector ( or ). As shown in equation 1 and Fig. 1B , each parameter in this framework is influenced by regulatory inputs from pleiotropic regions as well as modular factors ( and ). Each of these factors will mutate independently of the others, allowing us to evaluate the contributions of pleiotropic and modular factors to adaptation. This setup is referred to as the “with-Pleiotropy” setting. To better understand the effect of pleiotropy, we compared the “with-Pleiotropy” setting to a fully modular model, called the “without-Pleiotropy” setting ( Fig. 1C ). In this setting, the vector is replaced by two tooth-specific vectors, and (in green and orange in Fig. 1C ), such that, the model inputs for the two teeth become: In this case, all factors are modular, and comparing the two settings allows us to test the impact of pleiotropy on evolutionary dynamics. A upper and A lower are used to parameterize the tooth developmental model ( 34 ) to calculate the shape of the upper and lower teeth, see Methods: developmental model. We then compute occlusion as the optimal contact surface between the upper and lower teeth. To do this, we transform the 3D tooth morphology into a 2D buco-lingual representation (see methods Occlusion for details), accurately capturing tooth shape as they primarily develop along a single line of cusps. We used the 2D outlines to identify the position with the best possible contact surface without overlap, where a larger surface indicates better occlusion, increasing fitness and selection likelihood. The model runs in an evolutionary loop ( Fig. 1D ) with 200 individuals per generation, selecting the top 50% to replicate. Each replication mutates genes and regulatory sequences. These different regions of the genome mutate with different rates, proportional to the fraction they represent: large for the pleiotropic part A constitutive , smaller for the tooth-specific cis-regulatory part hoxcis , very small for the tooth-specific protein hoxtrans (see Methods: Mutations). Evolutionary dynamics will drive teeth toward improved occlusion. As the two settings differ in the number of A constitutive genes (21 for the with-Pleiotropy setting vs. 42 for the without-Pleiotropy), their mutational landscapes vary. To address this, we adjusted the mutation rate in the without-Pleiotropy setting to match the average mutation count in the with-Pleiotropy setting. To ensure results are not solely due to a low mutation rate, we added a third setting, “without-Pleiotropy-Double-Mutation”, where mutation rates remain unchanged. For each setting — with-Pleiotropy, without-Pleiotropy, and without-Pleiotropy-Double-Mutation — we ran simulations lasting 200 generations, all starting with single-cusp teeth and biologically informed initial regulatory values for 10 independent starting conditions and a total of 100 simulations per setting (see Methods: Making initial conditions). At each generation, individuals were selected for better occlusion, and we recorded fitness and other key phenotypic traits, such as number of cusps in each tooth. By comparing the final occlusion and the fitness trajectories in the three settings, we analyze the effects of pleiotropy and modularity on evolution. Since the with-Pleiotropy setting contains both types of parameters, we also examined the specific contributions of pleiotropic and modular factors to evolutionary dynamics, which revealed several advantages of pleiotropy. Populations simulated with-pleiotropy reach better occlusion values and achieve more changes in morphology During the simulations, occlusion generally increases, with notable changes in tooth morphology and cusp count over time ( Fig. 2A-B ). Download figure Open in new tab Fig. 2. Populations simulated with-pleiotropy reach better occlusion values and achieve more changes in morphology. (A) Example simulation in the with-Pleiotropy setting. Top: 2D profiles of the upper and lower teeth in the first (left) and last (right) generations, positioned for optimal occlusion. Note the increase in cusp count and improvement in occlusion. Bottom: Fitness (occlusion) over time, with the top 50% of individuals (light blue) contributing to the next generation, while the rest (black) do not. The inset shows cusp count changes over generations. (B) Same as (A) but for a simulation in the without-Pleiotropy setting. Note the distinct morphologies of the upper and lower teeth in the final generation (top panel). (C) Fitness of the ten best individuals in the last generation for all simulations across settings. The with-Pleiotropy setting achieves significantly better occlusion on average than the without-Pleiotropy settings ( P < 2.2 e − 16). Doubling the mutation rate in the non-pleiotropic setting also improves average fitness ( P < 2.2 e − 16). (D) Morphological distance between the first and last generations for the ten best occluding individuals in all simulations across settings, shown separately for upper and lower teeth. Distances are measured using the Euclidean Morphological Distance (see Methods: EMD). The with-Pleiotropy setting shows significantly greater morphological change than both the without-Pleiotropy and without-Pleiotropy-Double-Mutation settings ( P < 2.2 e − 16). (E) Number of cusps per tooth in the last generation. In the with-Pleiotropy setting, the prevalence of teeth with three cusps is significantly higher than that of single-cusp teeth and greater than in either of the without-Pleiotropy settings. On average, the with-Pleiotropy populations achieve significantly higher occlusion than those without pleiotropy, outperforming even the without-Pleiotropy-Double-Mutation setting ( Fig. 2C ). Morphology changes more in the with-Pleiotropy setting, frequently resulting in multi-cusped teeth ( Fig. 2D-E ). In contrast, the without-Pleiotropy settings often retain single-cusp teeth, similar to the initial generation ( Fig. 3A , example 1). Download figure Open in new tab Fig. 3. Interaction between cusp gain dynamics and pleiotropy during evolution results in better than average fitness. (A) Left panel: Each point corresponds to a single simulation. The x-axis corresponds to the average fitness of the top 10 occluding individuals in the last generation. The y-axis shows the correlation of cusp change (gain or loss) between the upper and lower teeth along evolution (measured as the Spearman correlation between the average upper and lower cusp count of the top 50% occluding individuals for each generation). The with-Pleiotropy setting shows a higher correlation ( r = 0.6359), than the without-Pleiotropy ( r = 0.3488) and the without-Pleiotropy-Double-Mutation ( r = 0.1475) settings. Right panel: four examples of occlusions at the final generation (numbers correspond to labels on the left panel). (B) For the three settings, we show how the number of cusps correlate between the upper and lower teeth (for all teeth along all generations). Colors show the amount of teeth. To understand why these differences arise, we studied the strategies leading to better occlusion across settings. Two main strategies emerged: increasing cusp number to expand contact surface ( Fig. 3A , example 2) and increasing tooth size ( Fig. 3A , example 4). A combination of these strategies can optimize occlusion, as seen in ( Fig. 3A , example 3). The latter strategy, found in the without-Pleiotropy-Double-Mutation setting, achieved the best occlusion observed but appeared in only one simulation. Because pleiotropy links the evolution of parameters, morphological evolution may be more correlated between upper and lower teeth in the with-Pleiotropy setting. Taking teeth from all generations together, we found that cusp number is indeed better correlated in pleiotropic conditions than without it ( Fig. 3B-D ). However, this correlation in cusp number between upper and lower teeth is far from perfect. Thus parallel cusp change in the two teeth occurs, but not systematically, and the success of the with-Pleiotropy setting is not easily explained. We measured the correlation between cusp increases in upper and lower teeth throughout evolution and related it to the final fitness at the end of each simulation ( Fig. 3A ). A clear cluster shows that low correlation aligns with poor occlusion, indicating cases where tooth morphology remains largely unchanged. We observe a positive correlation between coordinated cusp changes and final occlusion in the with-Pleiotropy setting ( r = 0.6359). While a high correlation in cusp gain is often advantageous for occlusion in simulations with pleiotropy, it does not guarantee high occlusion. Some tooth pairs with a strong correlation in cusp number changes between upper and lower teeth can still exhibit poor occlusion. Moreover, other occlusion-enhancing strategies, such as those shown in Fig. 3A , do not necessarily rely on high cusp-gain correlation. Notice that the best occlusion observed ( Fig. 3A , example 3) appeared in the without-Pleiotropy-Double-Mutation setting and not in the with-pleiotropy setting. Fitness steadily increases in the with-Pleiotropy setting Previous results show that, in our model, pleiotropy does not impede adaptation and that the morphologies of upper and lower teeth are more correlated in the with-Pleiotropy setting than in the two settings evolving without pleiotropy. However, they also show that this correlation does not fully explain the advantage of pleiotropy. To understand this effect, we turn our investigations to analyzing how the different pleiotropic settings traverse the fitness landscape in terms of their rate of fitness gain during evolution. We examined the cumulative acceleration of fitness improvement over time. Acceleration of fitness allows us to explore whether the variations of fitness are more or less gradual or saltatory in a quantitative way, which in turn allows to reveal differences in how these settings evolve. For each evolutionary simulation, we first constructed a smooth curve of fitness values, derived from the average fitness of individuals contributing to the subsequent generation (black curves in Fig. 4 insets). We then computed a total acceleration value per simulation over its whole trajectory, and displayed it against the final fitness level achieved ( Fig. 4 ). The three settings clearly segregate between them. For the same values of final fitness, the with-Pleiotropy setting exhibits much lower total acceleration, showing that fitness improves much more steadily with pleiotropy than without. Download figure Open in new tab Fig. 4. The exploration of the fitness landscape is smoother in the with-Pleiotropy setting. In order to study how pleiotropy affects the exploration of the fitness landscape we measure, for each simulation of the three settings, the total amount of acceleration in fitness change during evolution and plot it against the final fitness achieved by the simulation. Insets show two examples of fitness trajectories. On the left, a trajectory in the without-Pleiotropy setting. On the right, a trajectory in the with-Pleiotropy setting. In both insets the black lines show the average fitness of the top 50% of occluding individuals at each generation while the gray curves show the cumulative absolute value of acceleration, computed from the fitness trajectories. For similar final fitness values, simulations in the with-Pleiotropy setting show significantly less total accelerations ( P < 2.2 e − 16). Parameters in the with-pleiotropy setting evolve more while exploring a smaller fraction of the parameter space So far we compared the different settings on the basis of their fitness trajectories. To better understand the difference between the with-Pleiotropy and without-Pleiotropy settings, we analyzed the genomic evolution by examining changes over time in the different functional parts of the genome (i.e. the pleiotropic part of the 21 pleiotropic genes coded in A constitutive and the modular part of the 21 pleiotropic genes coded in hox cis , which includes both and ). First, we measured the change from the initial to the final generation (Δ begin − end ) for each gene in each setting. For the 21 pleiotropic genes, it implies assessing the change in the A constitutive and hox cis separately. Since different fitness gains could be associated with different levels of gene adaptation, we used a linear model to adjust for fitness differences when comparing ranges in the with-Pleiotropy and without-Pleiotropy settings (see methods). Fig. 5A provides an example of the A constitutive part of the gene that regulates epithelial proliferation. Each point represents a simulation, showing Δ begin − end for this gene on the y -axis (normalized between 0-1) and final fitness on the x -axis. For each setting, the linear model is shown with trend lines and 95% confidence intervals. Simulations with pleiotropy (in red) consistently show larger ranges than those without (in blue and purple), regardless of fitness. By comparing the intercepts across the three settings, we can estimate the average parameter change independently of fitness. Download figure Open in new tab Fig. 5. In the with-pleiotropy setting, parameters evolve more but explore a smaller fraction of the parameter space. For each simulation, we calculated two ranges of variation for each gene considering its two functional parts separately ( A constitutive and hox cis ): the variation over time (difference between parameter values in the last and first generations Δ begin − end ) and the overall parameter exploration (difference between the highest and lowest values observed Δ min − max ). (A) Linear model for Epithelial Proliferation gene, showing their Δ begin − end for the A constitutive part. Each point represents an evolutionary simulation with the different colors specifying to which setting it corresponds. (B) We then compared the parameter ranges of the two functional parts of pleiotropic genes ( A constitutive and hox cis ) between the different settings. For this we used the intercepts calculated with the linear models to estimate the differences in parameter ranges between With- and Without-pleiotropy settings independently of fitness (see SFig. 3 for the results using the without pleiotropy double mutation setting). Points above zero (white circles) show a larger range in the with-Pleiotropy setting, while points below zero (black) indicate a larger range in the without-Pleiotropy setting; gray points represent no significant difference between settings. Results for genes are not shown in the plot, but for both ranges, only the genes vary significantly more ( P < 1.4248 e − 06 and P = 0.0184 for Δ begin − end and Δ min − max respectively) in the with-pleiotropy setting. Fig. 5B summarizes the results for the 21 genes. To see the results for each specific parameter of the model, see SFig. 1 ). Positive values (white circles) indicate larger ranges for the with-Pleiotropy setting, while negative values (black points) indicate larger ranges in the without-Pleiotropy setting. Gray points signify no significant difference. Fig. 5B shows that Δ begin − end is generally larger in the with-Pleiotropy setting, suggesting greater accumulation of genetic variation during evolution. To see if this greater accumulation of mutations also concerns adaptive mutations, we performed an additional analysis. Adaptive mutations for each gene were estimated by correlating how changes in gene values correlated to improvements in fitness (see Methods: Adaptive change during evolution). The with-Pleiotropy setting shows significantly more adaptive changes (see Fig. 2 ). This accumulation of adaptative mutations in the with-Pleiotropy setting indicates a higher retention of beneficial mutations over time, while the without-Pleiotropy setting shows a consistently lower capacity to find or maintain such mutations. Finally, we calculated the range between each gene’s minimum and maximum values (Δ min − max ) to assess overall parameter space exploration within the whole population along the 200 generations. Fig. 5B shows that without-Pleiotropy simulations generally explore a larger range of genetic values, while with-Pleiotropy simulations exhibit lower Δ min − max values, except for three parameters. Together, this suggests that pleiotropy may limit overall parameter exploration, while allowing larger cumulative variations, in coherence with the steady fitness dynamics that we observed. Single mutations in the modular parameters in the pleiotropic setting result more often in adaptive changes It is widely speculated that morphological adaptation mainly relies on non pleiotropic mutations that arise in cis regulatory regions of the genome ( 40 ). To investigate this, we examined what proportion of the total mutations increased fitness and were consequently selected to propagate to the next generation. We focused on mutations affecting one parameter at a time to clearly identify which one improved fitness. Therefore we calculated the proportion of single mutations that improved fitness. In Fig. 6 , we see that both the hox cis and A constitutive part of pleiotropic genes contribute to adaptations. However, in the with-Pleiotropy setting (depicted in red), mutations in the modular part ( hox cis ) contributed 24% significantly more ( P < 2.2 e − 16) to fitness improvements than mutations in A constitutive part. On the opposite, in the without-Pleiotropy setting (illustrated in blue), no significant differences ( P = 0.8960) were observed between the contributions of hox cis and A constitutive parts (an expected result as, in the without-Pleiotropy setting, both types of genes play the same role). The results including the Without pleiotropy double mutation setting can be found in SFig. 4 . Download figure Open in new tab Fig. 6. Mutations in the hox cis part of pleiotropic genes result more often in adaptive changes, but substantial adaptive changes also occur in the A constitutive part. We calculated the fraction of mutations that resulted in a fitness improvement and contributed to the next generation. For simplicity, we only consider mutation events that changed one parameter at a time and combined the results for upper and lower teeth. In the y-axis we calculate, for each simulation, the total amount of mutations and the proportions of those mutations that produce an improvement in fitness. In the with-pleiotropy setting, the proportion of mutations in the hox cis part of pleiotropic genes that result in an increase in fitness is significantly higher (stars; P < 2.2 e − 16) than this proportion in the A constitutive part. In the without-pleiotropy setting, proportions in hox cis and A constitutive parts are not significantly different (n.s.; P = 0.8526). This was expected, since both parts have similar effects in this setting. In the with-pleiotropy setting, favorable mutations of the hox trans genes are significantly more frequent (stars; P < 2.2 e − 16) than in the without-pleiotropy setting. An unexpected result was to find that the fraction of adaptive mutations in A constitutive was equally high ( P = 0.8531) in the with-Pleiotropy setting as in the without-Pleiotropy setting, indicating that pleiotropy did not limit the adaptive potential of this functional part of the genome. Discussion A model with “Silent pleiotropy” and “effective pleiotropy” as an emergent property of the GPF map Our model represents a new way to tackle the influence of pleiotropy on evolution with simulations. In previous efforts, NK models have been extensively used, to account for multiple traits, pleiotropy and epistasis ( 41 – 44 ). In these models, pleiotropy and epistasis are embedded together, making it easy to compare results with estimates obtained from quantitative genetics, albeit not allowing the specific study of the effect of pleiotropy ( 45 ). However, more recent models have made it possible to implement pleiotropy and epistasis as independent features, enabling to track their evolution separately. These studies have revealed that pleiotropy and epistasis evolve together ( 43 ). However, pleiotropy and epistasis remain one specific parameter of the model, while in biological systems, pleiotropy and epistasis are the product of developmental mechanisms (GP map) and their evolution ( 18 , 46 , 47 ). Therefore it continues to be challenging to connect these modeling results with experimental findings or genomic evolution. Complex phenotypes such as morphology are produced by development with complex GP maps. Different organizational levels (molecules, cells, tissues) influence each other through developmental time, which means that pleiotropy occurring at the molecular level may not necessarily translate into pleiotropic effects on the phenotype. Rather, properties of the GP map determine how, in a given genetic background (epistasis), the effect of mutations propagate and affect one or more traits. The PF map also matters for predicting the effect of pleiotropy, as integrated functions and not phenotypes are relevant for selection ( 48 ). Our modeling strategy is based on the molecular understanding of pleiotropy and modularity in developmental systems, and a development-explicit GPF map (see Table 1 ). This effort was meant to bridge different views on pleiotropy. In our particular setting, we introduce pleiotropy at the molecular level in an explicit manner, by splitting mutations in pleiotropic genes according to their molecular effect: pleiotropic or organ-specific, if they change how the Hox-like regulator regulate the gene. However, pleiotropic mutations at the molecular level may ultimately have a pleiotropic effect or not at the phenotypic level (effective or silent pleiotropy), and they may have a neutral, positive or negative effect on fitness ( Table 1 ). All of this is not predefined by us, but depends on organ-specific epistatic interactions determined by the Hox and the properties of the GPF map (which exhibit some level of degeneracy at both G-P and P-F levels). Importantly, it also implies that, as the Hox and pleiotropic genes evolve with each generation, effective pleiotropy and silent pleiotropy are both free to evolve as well. View this table: View inline View popup Download powerpoint Table 1. Partitioning of mutations and their effect in the “with pleiotropy” setting Generalization of the findings obtained by focusing on pleiotropy between two serial organs Both the number and the serial nature of organs could in principle limit the generalization of our findings. We examine each of them below. First, there are only two organs in our model, with a single organ-specific transcriptional regulator. Fisher’s model predicts that the cost of pleiotropy increases in an organism with more than two morphological traits. But evolutionary simulations have shown that modularity evolves to scale with the number of traits ( 48 ). Accordingly, in nature, we observe that organ development is regulated by combinations of transcriptional regulators rather than single regulators. Thus, in our model both organ number and modularity have been simplified in parallel. Second, our model examines serial organs, i.e. organs of similar type, and in this context, we found that coordinated evolution can be particularly advantageous. This advantage appears in our pleiotropy setting, though we should note that coordinated evolution itself is neither necessary nor sufficient for achieving high fitness. This advantage for coordinated evolution is likely reduced when considering unrelated organs. However, modularity also scales with organ relatedness in nature. This is visible in the sharing of regulatory regions, which is almost total between serial organs, but only partial between unrelated organs ( 49 , 50 ). Modularity, measured as a diversity of transcription factor and combinatorial regulatory logic, also appears to scale with organ and cell type diversity in metazoans ( 51 – 53 ). Thus, in our model, high organ relatedness also scales with the simple modularity brought by only one regulator per organ. Given these observations, we believe our findings can be extended to understand evolution across entire organisms, taking this scaling effect into account. Morphological adaptation in modular versus pleiotropic regions For complex traits such as morphology, it has long been postulated that selection should more frequently act in cis-regulatory sequences than in coding sequences ( 26 ). This prediction relied on two notions. First, morphological changes are dependent on changes in gene expression throughout time and space during development. Second, cis-regulatory sequences being supposedly highly modular, changes in these sequences have significantly less pleiotropic effects than those in coding sequences. Consequently, it was assumed that evolution mostly acts through non pleiotropic mutations in highly modular cis-regulatory regions. Two decades of genetics of morphological evolution now highlight a slightly different picture. The predominant role of cisregulatory mutations in morphological evolution is now demonstrated in insects, thanks to a compilation of two decades of experimental findings ( 54 ). But case studies typically focus on one trait and do not usually ask whether the cis-regulatory mutation also generated variation in other parts of the organism, with rare exceptions ( 14 ). In fact, whether morphological evolution occurs predominantly with pleiotropic or non pleiotropic mutations in these cisregulatory regions is, more than ever, an open question for two reasons. First, there has been a paradigm shift in molecular developmental biology, showing that cis-regulatory sequences are often pleiotropic ( 55 ). Second, an increasing number of case studies show that mutations with pleiotropic effects are involved in morphological adaptation both at a micro and macro-evolutionary scale ( 9 , 14 , 15 , 31 , 56 ). Our study brings quantification to this framework, with the limitations of relying on simulations in a simple organism with two serial organs (discussed above). Changes downstream of the Hox factor (in cis) strongly dominate over changes in the Hox factor itself (in trans). This is consistent with the idea that the latter are often counter-selected because they have a large effect size, but can occasionally be advantageous. We observe that adaptive mutations occur more frequently in the non-pleiotropic regions of the genome ( hox cis + hox trans ) than in pleiotropic regions. Yet mutations in the pleiotropic regions remain substantial during adaptation and take a considerable part in genomic evolution. Our study therefore adds to the growing realization that pleiotropic mutations play an important role in morphological adaptation ( 9 , 11 , 14 , 15 , 57 – 61 ). Pleiotropy, genetic and developmental divergence The traditional view in molecular evolution and evo-devo has regarded pleiotropy as a conservative force that constrains evolutionary change and genomic evolution. However, our understanding has evolved significantly with the introduction of the “selection, pleiotropy, compensation” (SPC) model by Pavlicev and Wagner ( 30 ). This model proposes that when selection drives adaptive changes in one organ, it triggers subsequent selection for compensatory mutations that minimize pleiotropic effects on other organs. Rather than promoting conservation, pleiotropy would facilitate developmental system drift (DSD) and the evolution of the genetic basis for development. Sequentiality is however a problem in this model: pleiotropic mutations with beneficial effect in a single organ could be eliminated before mutations compensating for their deleterious effect in other organs could arise. Several lines of evidence in the literature support the idea that pleiotropy promotes divergence. In our previous research on rodent molar adaptation, we found that developmental changes driving shape evolution in one tooth (upper molar) were also present in the development of the other tooth (lower molar), despite its phenotype remaining largely conserved ( 62 ). This demonstrated how pleiotropic changes could be effectively compensated, allowing developmental drift while maintaining phenotypic stability. Additional support comes from modeling studies of trait evolution under directional and stabilizing selection, which showed rapid developmental system drift when pleiotropic genetic loci responded to selection ( 63 ). Furthermore, recent analyses of gene expression evolution in mammalian development have revealed that highly pleiotropic genes show similar levels of expression change as less pleiotropic genes ( 64 ). Evidence for pleiotropy-driven compensatory evolution between organs was recently found at the level of promoters and enhancers of highly pleiotropic genes. Their sequence evolves faster than that of more organ-specific regulatory regions, even though their activity is maintained throughout evolution through compensatory evolution of transcription factor binding sites ( 30 ). Our model shows that mutations in pleiotropic regions have comparable chances of being adaptive in both pleiotropic and non-pleiotropic settings. This suggests either that the burden of pleiotropy is minimal (very few of the mutations with a pleiotropic effect at the molecular level are pleiotropic at the phenotype level, “silent pleiotropy”), or alternatively, that it is efficiently balanced by compensatory mutations (effective pleiotropy being converted to silent pleiotropy by epistatic compensatory mutations), as in the SPC model. A third, more nuanced mechanism is however possible. As the relationship between molecular pleiotropy and its phenotypic effects is continuously shaped by selection, new mutations with pleiotropic effects at the molecular level could arise in favorable epistatic contexts silencing pleiotropy at the phenotypic level. In this scenario, purifying selection alone would be sufficient to maintain the appropriate epistatic context. This mechanism would explain how higher genetic divergence could accumulate, potentially enhancing long-term evolvability. Future studies using a framework that allows individual mutations to be traced will be crucial to test this hypothesis and better understand how pleiotropy acts as a dynamic property that promotes genetic divergence. Materials and Methods Developmental model The developmental model simulates developmental dynamics seen in teeth by integrating gene dynamics and tissue mechanics ( 34 ). The model starts with a flat layer of epithelium and mesenchyme under it, resembling the initial stage of tooth development. Cells in this model can produce signaling molecules which can diffuse over the cells in the model. These signaling molecules can regulate some cell behaviors such as cell division, differentiation or cell adhesion for example, additionally they can also regulate the expression of other gene products. The model includes three types of signaling molecules. An activator, which promotes its own synthesis and the differentiation of epithelial cells into enamel knots, which are signaling centers key in the development of teeth ( 65 ). These enamel knots do not divide and produce the two other molecules present in the model. An inhibitor, which inhibits the activator and precludes the formation of enamel knots close to each other, depending on how the activator and inhibitor diffuse and decay, the distance of enamel knots can be affected. Finally, the secondary signal, which is promoted by the activator and promotes differentiation, inhibits epithelial growth and activates mesenchymal growth. The model also includes growth biases in the antero-posterior and bucco-lingual axis, which can greatly affect the tooth morphology. In total the model has 21 parameters that account for how the different signaling molecules interact, diffuse, are promoted and decay, and also different biomechanical aspects of the model, such as the pressure that the stellate reticulum is affecting on the epithelium or the adhesion between cells. Mutation rates Our framework to model pleiotropy includes three types of parameters: hox trans , hox cis and A constitutive . To account for their biological representation and to be able to consider differences in how they evolve, we assigned different mutation rates for each of these categories. A constitutive represent a bigger amount of genes and gene network, they summarize a broad range developmental dynamics in a single parameter, therefore we assign them the highest mutational rate. The next category are the hox cis , they represent enhancers and cis regulatory regions for specific genes, in general they represent smaller regions in the genome, so with this naive approach, we simply assume that they have less chances to mutate and assign them a smaller mutation rate. Finally, hoxtrans represents single gene products and therefore, following the same logic as before, we assigned to it the lowest mutational rate. The mutation rates for each of the categories is: For A constitutive µp = 0.025; hox cis µm = 0.0025; hoxtrans µt = 0.00125. These mutation rates indicate the chance of each parameter to mutate, this means that multiple mutations can occur in an individual. These mutation rates are the ones used for the with-pleiotropy setting. The without-pleiotropy setting has double the amount of A constitutive parameters, since instead they have and , therefore we adjusted the mutation rate to be µp = 0.0125. The without-pleiotropy double mutation setting has the same mutation rates as the with-pleiotropy setting. Mutations magnitude Mutation affects each parameter by s i ( N (0, 1 )), where s i is a specific value for each parameter i and N (0, 1 ) is a random number in the normal distribution. To establish the magnitude of each mutation s i we followed two steps. First we ensure that the mutational steps were big enough to be able to increase ten times its initial value after 200 generations. Ten times the original value is enough to cover the parameter space within certain limits that would result in anomalous developmental behavior of the model, which leads to extremely non biological morphologies or causes the model to crash. We also ensured that at least with the initial parameter values, the average mutational step would have a small effect on the morphology. Due to the complex interactions intrinsic to the developmental model, this can only be assessed for a reduced number of parameter combinations, which is why we focused on the starting parameter values of the simulations. This was done by measuring the euclidean morphological distance (EMD) between the original tooth and the mutated tooth. An EMD = 0 means that no change between the morphologies occurred. If EMD > 0 we consider the mutational step to be significant enough, otherwise we increased s i . If the average mutation changes the morphology too much ( EMD > 0.2) the mutational step is reduced. This way we ensure that at least in the initial stages of evolution, mutations can produce enough changes in the tooth morphologies in order for evolution to be able to proceed, but not such big changes that every mutation could only be fixed if it produces a hopeful monster. Making initial populations We created ten different initial individuals in order to start at different locations in the parameter space. Each of these individuals was created by finding random parameters that would result in teeth with a single cusp; this way, each individual has an upper and a lower tooth with different parameters values, although both are single cusped. Initially these teeth had only pleiotropic parameters. To introduce the modular parameters, we subtracted from the pleiotropic parameters so that upper and lower teeth had the exact same values, the subtracted values were assigned to the modular parameters. When calculating these differences, in order to introduce some asymmetry in how the upper and lower teeth could evolve, we assigned some parameters to be specifically regulated by the upper tooth and some to the lower. In the upper teeth, modular parameters initially regulated the parameters related to epithelial dynamics (epithelial growth, adhesion and repulsion of epithelial cells). In the lower teeth, modular parameters initially regulated the parameters related to mesenchymal dynamics (mesenchymal growth, buoyancy and border growth). This is a simplified representation of what we find in mouse molar development, the modular regulatory protein specific to the upper tooth (Pou3f3) is expressed in the mesenchymal part of the tooth, while the one specific to the lower tooth (Nkx2-3) is expressed in the epithelial part of the tooth ( 33 ). The rest of the parameters were randomly assigned to either the upper or the lower modular parameters. This way we created 10 different individuals. For each of them we performed 200 generations of stabilizing evolution, selecting to maintain the shape of the initially single cusped morphology for the upper and the lower teeth. After the conservative selection, each of the 10 populations underwent directional selection to improve their occlusion. These was repeated 10 times for each initial population for each of the 3 different settings. Therefore the total number of evolutionary simulations is: 10 × 10 × 3 = 300 Euclidean minimal distance (EMD) This measure allows to compare morphologies made of different numbers of cells and without having to arbitrarily pre-select landmarks or special morphological features ( 35 ). This is highly convenient for our study because teeth can have different numbers of cells. EMD is defined as the mean distance from one cell in a morphology to the closest cell in another morphology. This is calculated as: Where n 1 and n 2 are the number of cells in tooth 1 and 2 respectively, d k ,min( k ,2) is the distance between cell k in tooth 1 and its closest cell in tooth 2. d j ,min( j ,1) is the distance between cell j in tooth 2 and its closest node in tooth 1. Occlusion To establish the fitness of an individual, we calculate the surface of contact between the upper and lower teeth of an individual. This surface of contact is what we define as the occlusion value; the higher the occlusion, the more fit we consider the individual to be. Although occlusion in real organisms is a more complex subject than we present here, this measurement serves as a proxy to calculate how two teeth evolve to improve occlusion. The teeth that we model here can accurately be reduced to two dimensions. In order to reduce the tooth morphology to a two dimensional profile, we first find the extreme cells in the anteroposterior dimension (x-axis). Then we find the cusps of the tooth, this is simply done by finding the enamel knot cells present in the tooth. In most cases, each cusp only has one enamel knot cell, but occasionally they can have more, to avoid counting more cusps than are actually present, the cartesian coordinates of enamel knots that are close to each other ( d < 1) are averaged together and counted as a single enamel knot. Next, we establish the tooth shape. We start with the cusp tips and the extreme antero-posterior points. In the following step, we locate the first cusp along the x-axis. From this cusp, we find the highest neighboring cell with a lower x-axis value and repeat until no such cells remain or until we reach the lowest x-axis point. This defines one side of the cusp profile, descending gradually by selecting only cells with decreasing x-axis values to avoid backtracking. We repeat the process for the other side, selecting cells with increasing x-axis values. By saving the coordinates of each selected cell and dropping the y-axis coordinate, we create a 2D profile of the cusp. With the profiles established, we calculate the contact surface between the teeth to determine the occlusion value. We begin by identifying the highest point C (the peak cusp) and the point with the highest x -axis value E in the lower tooth profile. These points define the search area As , a rectangle where: - The bottom edge is set by Ey , the y -coordinate of E . - The top edge is set by Cy , the y -coordinate of C . - The right edge is set by Ex , and the left edge by Cx . Within As , we translate the upper tooth profile to find the maximum contact surface, avoiding overlaps by using John M. Zerzan’s algorithm to check for intersections between the tooth profiles (each as polygons). If overlap is detected, that translation is excluded. For each non-excluded translation of the upper profile, we calculate the contact surface between the upper and lower tooth profiles. Points within 0.5 model units are considered in contact. We then measure the distance between consecutive contact points on each profile, sum these distances for both the upper and lower profiles, and add them together. This total gives the occlusion value, representing the contact surface between the teeth. Adaptive change during evolution To measure adaptive change during evolution, we analyze correlations between parameter changes and fitness, in addition to using the range method. Fitness changes may result from multiple parameters shifting, or from parameters hitchhiking with others that improve fitness. Multiple parameters might also need to change together to enhance fitness. By measuring correlations across simulations, we observe how often each parameter links to fitness improvements, helping to reduce the effects of hitchhiking and random dynamics. We measure correlations in rolling windows of 50 generations, moving from generation 1-50, then 2-51, and so on until generation 200. To evaluate how adaptive each parameter is, we count windows where the Pearson correlation between parameter and fitness exceeds r = 0.7. This gives a count of fitness-related correlations for each parameter across simulations. To assess differences in parameter evolution between with-Pleiotropy and without-Pleiotropy simulations, we perform a GLM analysis with a binomial family, categorizing windows as either above or below r = 0.7. We control for fitness differences by including fitness as an independent variable and compare intercept values across settings to identify significant differences between pleiotropic and non-pleiotropic conditions. Supplementary Information Download figure Open in new tab Supplementary figure 1. For each gene, we used the intercepts calculated with the linear models to estimate the differences in parameter ranges between With- and Without-pleiotropy settings independently of fitness, see section “Parameters in the with-pleiotropy setting evolve more while exploring a smaller fraction of the parameter space” for more details. Here we see the difference for each specific parameter. As the difference is calculated as Intercept with Pleiotropy - Intercept without pleiotropy, positive values show that the range of that parameter is larger in the with Pleiotropy setting. If there are no significant differences we show the bar in gray. To distinguish the upper from the lower teeth we use green for the upper and orange for the lower. Download figure Open in new tab Supplementary figure 2. Using the same strategy shown in Fig. 5a and explained in more detail in Methods: Adaptive change during evolution, we calculate how adaptive the different parameters are in our simulations. As the difference is calculated as Intercept with Pleiotropy - Intercept without pleiotropy, positive values show that the range of that parameter is larger in the with Pleiotropy setting. This is true for most parameters in both A constitutive and hox cis . Download figure Open in new tab Supplementary figure 3. For each gene, we used the intercepts calculated with the linear models to estimate the differences in parameter ranges between With- and Without-pleiotropy double mutation settings independently of fitness. Points above zero (white circles) show a larger range in the with-Pleiotropy setting, while points below zero (black) indicate a larger range in the without-Pleiotropy double mutation setting; gray points represent no significant difference between settings. Download figure Open in new tab Supplementary figure 4. Mutations in the hox cis parameters in the with-pleiotropy setting result more often in adaptive changes. We calculated the fraction of mutations that resulted in a fitness improvement and contributed to the next generation. For simplicity, we only consider mutation events that changed one parameter at a time and combined the results for upper and lower teeth. In the y-axis we calculate, for each simulation, the total amount of mutations and the proportions of those mutations that produce an improvement in fitness. ACKNOWLEDGMENTS This work was supported by a grant from IDEX Lyon ANR-16-IDEX-0005, and an European Council Research grant (ERC 2022 COG PLEIOTROPY 101088398). Salaries were supported by the Centre National de la Recherche Scientifique, the Ecole Normale Supérieure de Lyon and the INSA de Lyon. HPC resources where provided by the CBPsmn (PSMN, Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon. We also thank Lisandro Milocco, Benjamin Audit, Sabrina Renaud, Teemu Häkkinen and the Beagle’s team for their considerable advice and valuable comments. Funding European Research Council https://ror.org/0472cxd90 ERC 2022 COG PLEIOTROPY 101088398 Footnotes ↵ 2 marie.semon{at}ens-lyon.fr ↵ 3 sophie.pantalacci{at}ens-lyon.fr SP, PF, MS and GB designed the experiments; PF wrote all codes and performed experiments; SP, PF, MS and GB drafted and revised the manuscript The authors declare no competing interest. References 1. ↵ RA Fisher , The genetical theory of natural selection: a complete variorum edition . ( Oxford University Press ), ( 1999 ). 2. HA Orr , The population genetics of adaptation: the distribution of factors fixed during adaptive evolution . Evolution 52 , 935 – 949 ( 1998 ). OpenUrl CrossRef PubMed Web of Science 3. ↵ HA Orr , Adaptation and the cost of complexity . Evolution 54 , 13 – 20 ( 2000 ). OpenUrl CrossRef PubMed Web of Science 4. JJ Welch , D Waxman , Modularity and the cost of complexity . Evolution 57 , 1723 – 1734 ( 2003 ). OpenUrl CrossRef PubMed Web of Science 5. NH Barton , Pleiotropic models of quantitative variation . Genetics 124 , 773 – 782 ( 1990 ). OpenUrl Abstract / FREE Full Text 6. PD Keightley , WG Hill , Quantitative genetic variability maintained by mutation-stabilizing selection balance in finite populations . Genet. Res . 52 , 33 – 43 ( 1988 ). OpenUrl CrossRef PubMed Web of Science 7. ↵ SP Otto , Two steps forward, one step back: the pleiotropic effects of favoured alleles . Proc. Royal Soc. London. Ser. B: Biol. Sci . 271 , 705 – 714 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 8. ↵ J Lourenço , N Galtier , S Glémin , Complexity, pleiotropy, and the fitness effect of mutations . Evolution 65 , 1559 – 1571 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 9. ↵ DJ Rennison , CL Peichel , Pleiotropy facilitates parallel adaptation in sticklebacks . Mol. Ecol . 31 , 1476 – 1486 ( 2022 ). OpenUrl CrossRef PubMed 10. ↵ RA Martin , AT Tate , Pleiotropy promotes the evolution of inducible immune responses in a model of host-pathogen coevolution . PLoS computational biology 19 , e1010445 ( 2023 ). OpenUrl CrossRef 11. ↵ C Knibbe , O Mazet , F Chaudier , JM Fayard , G Beslon , Evolutionary coupling between the deleteriousness of gene mutations and the amount of non-coding sequences . J. Theor. Biol . 244 , 621 – 630 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 12. ↵ JM Collet , K McGuigan , SL Allen , SF Chenoweth , MW Blows , Mutational pleiotropy and the strength of stabilizing selection within and between functional modules of gene expression . Genetics 208 , 1601 – 1616 ( 2018 ). OpenUrl Abstract / FREE Full Text 13. M Pavličev , JM Cheverud , Constraints evolve: context dependency of gene effects allows evolution of pleiotropy . Annu. Rev. Ecol. Evol. Syst . 46 , 413 – 434 ( 2015 ). OpenUrl CrossRef 14. ↵ O Nagy , et al. , Correlated evolution of two copulatory organs via a single cis-regulatory nucleotide change . Curr. Biol . 28 , 3450 – 3457 ( 2018 ). OpenUrl CrossRef PubMed 15. ↵ SL Archambeault , L. Bärtschi , AD Merminod , CL Peichel , Adaptation via pleiotropy and linkage: association mapping reveals a complex genetic architecture within the stickleback eda locus . Evol. letters 4 , 282 – 301 ( 2020 ). OpenUrl CrossRef 16. ↵ TF Mackay , RR Anholt , Pleiotropy, epistasis and the genetic architecture of quantitative traits . Nat. Rev. Genet . 25 , 639 – 657 ( 2024 ). OpenUrl PubMed 17. ↵ GP Wagner , M Pavlicev , JM Cheverud , The road to modularity . Nat. Rev. Genet . 8 , 921 – 931 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 18. ↵ GP Wagner , J Zhang , The pleiotropic structure of the genotype–phenotype map: the evolvability of complex organisms . Nat. Rev. Genet . 12 , 204 – 213 ( 2011 ). OpenUrl CrossRef PubMed 19. AB Paaby , MV Rockman , The many faces of pleiotropy . Trends genetics 29 , 66 – 73 ( 2013 ). OpenUrl CrossRef PubMed 20. DL Stern , V Orgogozo , The loci of evolution: how predictable is genetic evolution? Evolution 62 , 2155 – 2177 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 21. ↵ J Zhang , Patterns and evolutionary consequences of pleiotropy . Annu. review ecology, evolution, systematics 54 , 1 – 19 ( 2023 ). OpenUrl CrossRef 22. ↵ A Porto , R Schmelter , JL VandeBerg , G Marroig , JM Cheverud , Evolution of the genotype-to-phenotype map and the cost of pleiotropy in mammals . Genetics 204 , 1601 – 1612 ( 2016 ). OpenUrl Abstract / FREE Full Text 23. ↵ HB Fraser , AE Hirsh , LM Steinmetz , C Scharfe , MW Feldman , Evolutionary rate in the protein interaction network . Science 296 , 750 – 752 ( 2002 ). OpenUrl Abstract / FREE Full Text 24. MW Hahn , AD Kern , Comparative genomics of centrality and essentiality in three eukaryotic protein-interaction networks . Mol. biology evolution 22 , 803 – 806 ( 2005 ). OpenUrl CrossRef 25. ↵ RR Masalia , AJ Bewick , JM Burke , Connectivity in gene coexpression networks negatively correlates with rates of molecular evolution in flowering plants . PLoS One 12 , e0182289 ( 2017 ). OpenUrl CrossRef PubMed 26. ↵ SB Carroll , Evo-devo and an expanding evolutionary synthesis: a genetic theory of morphological evolution . Cell 134 , 25 – 36 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 27. ↵ X He , J Zhang , Toward a molecular understanding of pleiotropy . Genetics 173 , 1885 – 1891 ( 2006 ). OpenUrl Abstract / FREE Full Text 28. ↵ K Mattioli , et al. , Cis and trans effects differentially contribute to the evolution of promoters and enhancers . Genome biology 21 , 1 – 22 ( 2020 ). OpenUrl CrossRef PubMed 29. ↵ SA Signor , SV Nuzhdin , The evolution of gene expression in cis and trans . Trends Genet . 34 , 532 – 544 ( 2018 ). OpenUrl CrossRef PubMed 30. ↵ M Pavlicev , GP Wagner , A model of developmental evolution: selection, pleiotropy and compensation . Trends Ecol. & Evol . 27 , 316 – 322 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 31. ↵ M Sémon , et al. , Phenotypic innovation in one tooth induced concerted developmental evolution in another . bioRxiv pp. 2020 – 04 ( 2023 ). 32. ↵ M Pavlicev , EA Norgard , GL Fawcett , JM Cheverud , Evolution of pleiotropy: epistatic interaction pattern supports a mechanistic model underlying variation in genotype–phenotype map . J. Exp. Zool. Part B: Mol. Dev. Evol . 316 , 371 – 385 ( 2011 ). OpenUrl 33. ↵ S Pantalacci , et al. , Transcriptomic signatures shaped by cell proportions shed light on comparative developmental biology . Genome biology 18 , 1 – 25 ( 2017 ). OpenUrl CrossRef PubMed 34. ↵ I Salazar-Ciudad , J Jernvall , A computational model of teeth and the developmental origins of morphological variation . Nature 464 , 583 – 586 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 35. ↵ I Salazar-Ciudad , M Marín-Riera , Adaptive dynamics under development-based genotype–phenotype maps . Nature 497 , 361 – 364 ( 2013 ). OpenUrl CrossRef PubMed 36. L Milocco , I Salazar-Ciudad , Is evolution predictable? quantitative genetics under complex genotype-phenotype maps . Evolution 74 , 230 – 244 ( 2020 ). OpenUrl CrossRef PubMed 37. ↵ R Zimm , F Berio , M Debiais-Thibaud , N Goudemand , A shark-inspired general model of tooth morphogenesis unveils developmental asymmetries in phenotype transitions . Proc. Natl. Acad. Sci . 120 , e2216959120 ( 2023 ). OpenUrl CrossRef PubMed 38. ↵ D Schluter , Adaptive radiation along genetic lines of least resistance . Evolution 50 , 1766 – 1774 ( 1996 ). OpenUrl CrossRef PubMed Web of Science 39. ↵ PS Ungar , L Hlusko , The evolutionary path of least resistance . Science 353 , 29 – 30 ( 2016 ). OpenUrl Abstract / FREE Full Text 40. ↵ SB Carroll , Evolution at two levels: on genes and form . PLoS biology 3 , e245 ( 2005 ). OpenUrl CrossRef PubMed 41. ↵ P Mehra , A Hintze , The role of pleiotropy and epistasis on evolvability and robustness in a two-peak fitness landscape . Biology 13 , 1003 ( 2024 ). OpenUrl CrossRef PubMed 42. B Østman , Effects of epistasis and pleiotropy on fitness landscapes . Evol. Biol. Exobiol. Evol. Mech . pp. 37 – 46 ( 2013 ). 43. ↵ B Østman , A Hintze , C Adami , Impact of epistasis and pleiotropy on evolutionary adaptation . Proc. Royal Soc. B: Biol. Sci . 279 , 247 – 256 ( 2012 ). OpenUrl 44. ↵ LM Chevin , G Martin , T Lenormand , Fisher’s model and the genomics of adaptation: restricted pleiotropy, heterogenous mutation, and parallel evolution . Evolution 64 , 3213 – 3231 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 45. ↵ AJ Petit , J Guez , A Le Rouzic , Correlated stabilizing selection shapes the topology of gene regulatory networks . Genetics 224 , iyad065 ( 2023 ). OpenUrl CrossRef PubMed 46. ↵ MA Félix , A Wagner , Robustness and evolution: concepts, insights and challenges from a developmental model system . Heredity 100 , 132 – 140 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 47. ↵ JM Cheverud , Phenotypic, genetic, and environmental morphological integration in the cranium . Evolution pp. 499 – 516 ( 1982 ). 48. ↵ H Le Nagard , L Chao , O Tenaillon , The emergence of complexity and restricted pleiotropy in adapting networks . BMC evolutionary biology 11 , 1 – 15 ( 2011 ). OpenUrl CrossRef PubMed 49. ↵ Z Kliesmete , et al. , Evidence for compensatory evolution within pleiotropic regulatory elements . Genome Res . 34 , 1528 – 1539 ( 2024 ). OpenUrl Abstract / FREE Full Text 50. ↵ DJ McKay , JD Lieb , A common set of dna regulatory elements shapes drosophila appendages . Dev. cell 27 , 306 – 318 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 51. ↵ BM Degnan , M Vervoort , C Larroux , GS Richards , Early evolution of metazoan transcription factors . Curr. opinion genetics & development 19 , 591 – 599 ( 2009 ). OpenUrl CrossRef 52. A de Mendoza , et al. , Transcription factor evolution in eukaryotes and the assembly of the regulatory toolkit in multicellular lineages . Proc. Natl. Acad. Sci . 110 , E4858 – E4866 ( 2013 ). OpenUrl Abstract / FREE Full Text 53. ↵ T Ravasi , et al. , An atlas of combinatorial transcriptional regulation in mouse and man . Cell 140 , 744 – 752 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 54. ↵ V Courtier-Orgogozo , The loci of insect phenotypic evolution . Curr. Opin. Insect Sci . 60 , 101134 ( 2023 ). OpenUrl CrossRef PubMed 55. ↵ EPB Noon , et al. , Comprehensive analysis of a cis-regulatory region reveals pleiotropy in enhancer function . Cell reports 22 , 3021 – 3031 ( 2018 ). OpenUrl CrossRef PubMed 56. ↵ G Rice , M Rebeiz , Evolution: how many phenotypes do regulatory mutations affect? Curr. Biol . 29 , R21 – R23 ( 2019 ). OpenUrl CrossRef PubMed 57. ↵ AK Greenwood , MG Mills , AR Wark , SL Archambeault , CL Peichel , Evolution of schooling behavior in threespine sticklebacks is shaped by the eda gene . Genetics 203 , 677 – 681 ( 2016 ). OpenUrl Abstract / FREE Full Text 58. T Hämälä , et al. , Gene expression modularity reveals footprints of polygenic adaptation in theobroma cacao . Mol. biology evolution 37 , 110 – 123 ( 2020 ). OpenUrl CrossRef 59. JJ Lewis , et al. , Parallel evolution of ancient, pleiotropic enhancers underlies butterfly wing pattern mimicry . Proc. Natl. Acad. Sci . 116 , 24174 – 24183 ( 2019 ). OpenUrl Abstract / FREE Full Text 60. MG Mills , AK Greenwood , CL Peichel , Pleiotropic effects of a single gene on skeletal development and sensory system patterning in sticklebacks . EvoDevo 5 , 1 – 10 ( 2014 ). OpenUrl CrossRef PubMed 61. ↵ SD Smith , Pleiotropy and the evolution of floral integration . New Phytol . 209 , 80 – 85 ( 2016 ). OpenUrl CrossRef PubMed 62. ↵ M Sémon , et al. , Comparative transcriptomics in serial organs uncovers early and pan-organ developmental changes associated with organ-specific morphological adaptation . Nat. Commun . 16 , 768 ( 2025 ). OpenUrl CrossRef PubMed 63. ↵ NA Johnson , AH Porter , Evolution of branched regulatory genetic pathways: directional selection on pleiotropic loci accelerates developmental system drift . Genetica 129 , 57 – 70 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 64. ↵ M Cardoso-Moreira , et al. , Gene expression across mammalian organ development . Nature 571 , 505 – 509 ( 2019 ). OpenUrl CrossRef PubMed 65. ↵ J Jernvall , I Thesleff , Reiterative signaling and patterning during mammalian tooth morphogenesis . Mech. development 92 , 19 – 29 ( 2000 ). OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted April 13, 2025. Download PDF 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 Pleiotropy accelerates tooth phenotypic and genomic evolution — An in silico study under the lens of development 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 Pleiotropy accelerates tooth phenotypic and genomic evolution — An in silico study under the lens of development Pascal Felix Hagolani , Marie Sémon , Guillaume Beslon , Sophie Pantalacci bioRxiv 2025.04.11.648404; doi: https://doi.org/10.1101/2025.04.11.648404 Share This Article: Copy Citation Tools Pleiotropy accelerates tooth phenotypic and genomic evolution — An in silico study under the lens of development Pascal Felix Hagolani , Marie Sémon , Guillaume Beslon , Sophie Pantalacci bioRxiv 2025.04.11.648404; doi: https://doi.org/10.1101/2025.04.11.648404 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 Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7636) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41967) Biophysics (21460) Cancer Biology (18600) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13384) Ecology (19909) Epidemiology (2067) Evolutionary Biology (24326) Genetics (15613) Genomics (22512) Immunology (17740) Microbiology (40423) Molecular Biology (17193) Neuroscience (88645) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4825) Physiology (7647) Plant Biology (15159) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)