Gene network simulations provide testable predictions for the molecular domestication syndrome

The domestication of plant species lead to repeatable morphological evolution, often referred to as the phenotypic domestication syndrome. Domestication is also associated with important genomic changes, such as the loss of genetic diversity compared to adequately large wild populations, and modifications of gene expression patterns. Here, we explored theoretically the effect of a domestication-like scenario on the evolution of gene regulatory networks. We ran population genetics simulations in which individuals were featured by their genotype (an interaction matrix encoding a gene regulatory network) and their gene expressions, representing the phenotypic level. Our domestication scenario included a population bottleneck and a selection switch mimicking human-mediated directional and canalizing selection, i.e., change in the optimal gene expression level and selection towards more stable expression across environments. We showed that domestication profoundly alters genetic architectures. Based on four examples of plant domestication scenarios, our simulations predict (i) a drop in neutral allelic diversity, (ii) a change in gene expression variance that depends upon the domestication scenario, (iii) transient maladaptive plasticity, (iv) a deep rewiring of the gene regulatory networks, with a trend towards gain of regulatory interactions, and (v) a global increase in the genetic correlations among gene expressions, with a loss of modularity in the resulting coexpression patterns and in the underlying networks. We provide empirically testable predictions on the differences of genetic architectures between wild and domesticated forms. The characterization of such systematic evolutionary changes in the genetic architecture of traits contributes to define a molecular domestication syndrome.


INTRODUCTION
Domestication is a process of rapid evolution over successive generations of anthropogenic selection, leading to adaptation to habitats created by humans and acquisition of profitable traits for them. Such innovations originate from genetic and plastic variation sustaining phenotypic shifts in domesticates compared to their wild counterparts (Fuller et al. 2010). In plants, traits targeted by those shifts alter architecture (more compact morphology), life-history (loss or partial loss of seed dispersal and seed dormancy, increased synchronicity of germination and ripening), as well as production-and usage-related traits (taste, increase of harvestable organs). They are often associated with convergent phenotypic changes across species (Larson et al. 2014), and collectively referred to as the phenotypic domestication syndrome.
The discovery of the genetic bases underlying variation of domesticated traits has been the focus of ample empirical work. Dozens of domestication genes have been discovered, most of which are transcription factors (Martínez-Ainsworth and Tenaillon 2016;Fernie and Yan 2019) embedded into complex gene regulatory networks (GRNs). Perhaps the most emblematic example is provided by the Tb1 gene, which together with other genes controls maize branching architecture via hormone and sugar signaling (Doebley et al. 1997;Whipple et al. 2011;Dong et al. 2017Dong et al. , 2019. It is responsible for the strong apical dominance phenotype, i.e. repression of axillary bud outgrowth (Clark et al. 2006). Interestingly, in contrast to the maize allele, the Tb1 allele from its wild ancestor (teosinte) confers a responsiveness to light when introgressed into a maize background (Lukens and Doebley 1999). It therefore appears that domestication has triggered the selection of a constitutive shade avoidance phenotype in maize (Studer et al. 2017), that has translated into a loss of phenotypic 3 55 60 65 70 plasticity. Along this line, recent results indicate that reduced Genotype-by-Environment (GxE) interactions may be a general consequence for traits targeted by human selection. For example, genomic regions displaying footprints of selection explain less variability for yield GxE than "neutral" regions (Gage et al. 2017). Decreased phenotypic plasticity during domestication likely results both from the stability of human-made compared with wild habitats and from selection for stable crop performance across environments; it has yet to be characterized in other crops and for a broader range of traits.
In addition to the phenotypic domestication syndrome, genome-wide sequencing data have revealed the outlines of a molecular domestication syndrome. This molecular syndrome includes a loss of genetic diversity through linked selection and constriction of population size due to sampling effects . Severity of those genetic bottlenecks as estimated by nucleotide diversity loss, ranges from 17% to 49% in annuals while often no loss is observed in perennial fruit crops (reviewed in Gaut et al. 2015). The combined effect of bottlenecks, increased inbreeding (Glémin and Bataillon 2009) and linked selection in domesticates translates into shrink in effective population size, which in turn reduces the efficacy of selection against deleterious mutations (Moyers et al. 2018). Although increased recombination rate in domesticates compared with their wild relatives may partially compensate this effect (Ross-Ibarra 2004), fixation of deleterious mutations in domesticates and a resulting genetic load is often observed as exemplified in African rice (Nabholz et al. in rice, cotton (Liu et al. 2019), beans (Bellucci et al. 2014); a significant gain as in tomato (Sauvage et al. 2017); or no substantial change as in soybean (Liu et al. 2019), olives (Gros-Balthazard et al. 2019), and maize (Swanson-Wagner et al. 2012). In the latter, however, reduced variation in expression was observed at domestication candidate genes, indicating that selection primarily acts on cis-acting regulatory variants ): most evolutionary-relevant mutations affecting the evolution of gene expressions are located in (or in the close vicinity) of the domestication genes. This result was further confirmed in F1 hybrids from maize / teosinte crosses where large differences in expression were primarily caused by cisdivergence, and correlated with genes targeted by selection during domestication (Lemmon et al. 2014).
Beyond quantitative measures of gene expression, domestication is also associated with gene network rewiring. A pioneer work in maize indeed indicates that 6% of all genes display altered co-expression profiles among which, genes targeted by selection during domestication and/or breeding are over-represented (Swanson-Wagner et al. 2012). Interestingly, networks encompassing domestication targets display greater connectivity in wild than in domesticated forms as if selection had triggered connection loss to/from these genes. In beans, coexpression networks at the genome level revealed a global excess of strong correlations in domesticates compared with wild, the latter being sparser with more isolated nodes and smallest connected components than the former (Bellucci et al. 2014). In contrast to maize, little qualitative difference was reported as for networks surrounding selected and neutral contigs.
While population genetic tools have been broadly used to estimate domestication bottlenecks and associated genetic load in plants (Eyre-Walker et al. 1998 interpreted as an allele, i.e., the set of regulatory sites in the promoter of the gene.
The model considered discrete regulatory time steps, and the expression of the n genes, stored in a vector P , changes during the development of an individual as ., x n ) applies a sigmoid scaling function f ( x ) to all elements to ensure that gene expression ranges between 0 (no expression) and 1 (full expression). We used an asymmetric scaling function as in Rünneburger and Le Rouzic (2016); Odorico et al. (2018): . This function is defined such that a=0.2 stands for the constitutive expression (in absence of regulation, all genes are expressed to 20% of their maximal expression).
The kinetics of the gene network was simulated for 16 time-steps in each individual, starting from P 0 =( a ,... , a ) . The simulation program reports, for each gene i , the mean p i and the variance V i of its expression level over the four last time steps. A non-null variance characterizes networks that have not reached equilibrium at 16-4=12 time-steps, either because of slow network dynamics or because the network is unstable (cyclic pattern). In addition to this traditional framework, we considered that one of the network genes was a "sensor" gene influenced by the environment. This makes it possible for the network to react to an environmental signal, and evolve expression plasticity. In practice, the environmental signal at generation g was drawn in a uniform distribution e g ∼ U ( 0,1) and the value of the 7 155 160 165 170 sensor gene was e g at each time step (the sensor gene had no regulator and was not influenced by the internal state of the network).

Population model
The Individual fitness w was calculated as the product of two components, w =w U ×w S . The first term w U corresponds to the penalty for networks that have not reached being the strength of selection on gene expression variance (i.e., selection against expression instability). The second term w S corresponds to a Gaussian stabilizing selection component, which depends on the distance between the expression phenotype and a selection target θ : , s i standing for the strength of stabilizing selection on gene i. As detailed below, some genes were not selected (in which case s i =0 ), some genes were selected for a stable optimum θ i ("stable" genes), while a last set of genes were selected for optima that changed at every generation g ("plastic genes"), half of them being selected for θ i g =e g , and the other half for θ i g =1−e g .
Mutations occurred during gametogenesis with a rate m , expressed as the mutation probability per haploid genome. A mutation consists in replacing a random element plastic genes ( Figure S2). At the onset of domestication, we modified the selection regime to mimic increased environmental stability and, in turn, decreased plasticity (12 unselected, 10 stable and 2 plastic genes, Figure S2). The mutation rate was set to m=10 − 3 / gamete/ generation, which, given the estimated mutational target of 24 genes of 1kb (average estimated length of enhancers from Oka et al. 2017;Ricci et al. 2019) roughly corresponded to a per-base mutation rate of 3×10 − 8 par generation, close to the maize estimate (Clark et al. 2004).
In addition to the maize default domestication scenario, we considered three additional domestication scenarios (African rice, pearl millet, and tomato In addition to the four default domestication scenarios (maize, African rice, pearl millet, tomato) described above, we explored control simulations to disentangle the contribution of the bottleneck and the selection switch in emerging patterns, based on the maize scenario ('Default'): a scenario with no bottleneck, and a scenario with no selection switch. Given the importance of these three scenarios (Default, no bottleneck, no selection switch) in the analysis of the results, simulations were run with a longer burn-in ( T a =24,000 ) to ensure that the network was close to the mutation-selection-drift equilibrium at the onset of domestication. We further assessed the sensitivity of our results for the maize default domestication scenario to independent changes in parameters values by (1) increasing the number of genes of the GRN, from 24 to 48, and doubling the number of selected genes and the mutation rate per genome accordingly; (2) setting the mutation rate to 0 at the time of domestication to evaluate selection response from standing variation only; (3) modulating selection intensity both through a decrease in selected genes count (by two-fold), and through a modification of the fitness function to simulate stronger selection; (4) dissociating selection switch from a loss of plasticity, either by maintaining the selection for plasticity over genes during domestication, or by keeping the same number of plastic genes before and after domestication ; (5) testing the effect of a harsher (10 times less individuals) bottleneck. We also assessed the sensitivity of the model to arbitrary parameters influencing the gene network dynamics, such as the number of time steps, or the selection on network instability. All scenarios were replicated 1000 times; unless specified otherwise, the reported variables were averaged over all individuals from the population; figures report the grand mean over the replicates; colored areas stand for the 10% -90% quantiles over the simulation replicates.

Model output and descriptive statistics
For each simulation run, summary statistics were computed every 100 generations.
The output includes the population mean and variance of (i) the absolute fitness w , (ii) gene expressions p i , (iii) gene regulations W i j for all pairs of genes. In addition, the environmental index e g and all selection optima θ g were recorded.
Effective population sizes were estimated as and E w being the population variance and the population mean of the absolute fitness, respectively). When computed over a time interval (e.g. over the duration of the bottleneck), the harmonic mean effective size N e =σ / ∑ t =1 σ 1 / N et was reported.
A proxy for neutral molecular variance around gene i was obtained by reporting the average population variance of the W i j for a subset of genes j which expression was very low (p j <0.1 over the whole simulation), as regulatory sites sensitive to nonexpressed transcription factors are expected to evolve neutrally.
Environmental reaction norms (gene expression plasticity) were estimated for each gene i by regressing the average expression p i over the environmental index e g , taken over a sliding window of 10 consecutive measurements (1000 generations).
The effect of gene regulations W i j being quantitative (and thus, never exactly 0), the presence/absence of a connection in the network was determined by the following procedure: the expression phenotypes P and P i j 0 were calculated both 13 300 305 310 315 from the full W matrix, and from each of the n ² possible W i j 0 matrices in which W i j was replaced by 0. The regulation W i j was considered as a meaningful connection when the Euclidean distance d ( P , P i j 0 ) exceeded an arbitrary threshold of 0.1. Using other thresholds shifted the number of connections upward or downward, but did not affect the results qualitatively.
Genetic correlation matrices were estimated directly from the population covariances in gene expressions (hereafter called G matrices, although they reflect here all genetic components and not only additive (co)variances). The evolution of G matrices was tracked by computing the distance between consecutive matrices G g and G g +500 in the simulation output. In practice, genetic covariances were turned into genetic correlation matrices, and then into genetic distance matrices d =√ 2(1−r) .
The difference between both genetic distance matrices was calculated from their element-wise correlation, as in a Mantel test (function mantel.rtest in the R package ade4, Dray and Dufour 2007). Network topological features, including the number of clusters used as an index of modularity, were measured with the package igraph (Csardi and Nepusz 2006).

Implementation and Data availability
The simulation model was implemented in C++ and compiled with gcc v-7.5.0. The simulation software is available at https://github.com/lerouzic/simevolv. Simulation runs were automated via bash scripts, and simulation results were analyzed with R version 4.0 (R Core Team 2020). All scripts (simulation launcher, data analysis, and figure generation) are available at https://github.com/lerouzic/domestication. Demographic parameters were inspired by four documented domestication histories, two outcrossers (maize and pearl millet) and two selfers (African rice and tomato).
We modeled the selection switch associated with domestication both as a change in the gene expression optima and a partial loss of plastic responses. The default maize domestication scenario was compared with simulations without bottleneck (albeit a selection switch), and simulations without selection switch (albeit a bottleneck), and we also explored independent variation of parameters values to explore the sensitivity of our results. The whole simulation approach is summarized in Figure S3.

Adaptation during habitat shift
The We simulated the loss of plasticity during domestication as a change in selection regime for four plastic genes (out of 6) towards stable selection or neutrality ( Figure S2). The speed at which the gene network evolved increased by a factor of roughly 13 when the selection regime shifted ( Figure 3A). The selection switch translated into an abrupt change in reaction norm for genes that became selected for a flat reaction norm (plastic → stable in Figure 3B). We indeed observed a rapid loss of plasticity, showing that it was an evolvable feature that responded to selection.
More surprisingly, however, the loss of plasticity also affected genes that (i) were no longer under direct selection (Plastic → Neutral in Figure 3B), and (ii) were supposed to remain plastic (Plastic → Plastic in Figure 3B), albeit to a lower extent. This shortterm maladaptive evolution highlighted the genetic constraints during the rewiring of the network caused by the selection switch. Immediately after it, plastic genes were still tightly connected to genes that were selected to evolve a flat reaction norm, and the first stage of this evolutionary change involved a maladaptive trade-off. showing that it resulted from underlying constraints of the network, where selectiontriggered changes in reaction norms at some genes affect the evolution of plastic genes.

Molecular variation is affected by both demography and selection
Regulation strength was modeled as a quantitative variable directly affected by mutation ( Figure 1). For any given gene in the network, we measured its neutral molecular variance among individuals of the population as the average variance of the regulation strength at regulatory sites that had no influence on gene expression.
Hence, molecular variance is an analogous measure of neutral nucleotide genetic diversity of the genes of the network.
Based on empirical evidence, the first signal that we expected was a loss of neutral genetic diversity. The variance indeed dropped at the beginning of the domestication ( Figure 4A). Such variance drop was driven by genetic drift, that increased during the bottleneck. The maximum observed drop in genetic diversity was ~30% loss during the bottleneck for the default scenario. Recovery was slow and still ongoing at the end of the simulations.
In addition to change in molecular variance, we investigated the evolution of phenotypic (expression) variance during the domestication. Our results showed that in contrast to the neutral genetic variance, phenotypic variance may increase during domestication ( Figure 4B). Expression variance bursts, absent from the simulations 17 395 400 405 410 without selection switch, can be associated with ongoing adaptation: they corresponded to the segregation of selected variants that brought the phenotype closer to the new optimum. Domestication was thus associated with an increase in the gene expression variance, as a result of the balance between the selection switch (which increased temporarily the variance, Figure 4B) and the bottleneck (which slightly reduced the variance, Figure 4B). In case of a stronger bottleneck, however, the expression diversity was reduced at the selection switch showing that the net effect on phenotypic diversity strongly depends on the details of the domestication scenario ( Figure S4D).

Domestication is associated with the rewiring of gene networks
Genetic ( Figures 6A and S7A). This rewiring was solely due to the selection switch, as there was no effect of the bottleneck alone on the network evolution. The rewiring was associated with a systematic excess of gained connections over lost connections, i.e., domestication caused an increase in the total number of connections ( Figure   6A). As a consequence of the gain of new connections, the number of clusters decreased (some connections appeared between previously independent modules, Figure 6B). New connections appeared to be distributed evenly across the network ( Figure S7B and S8). in all scenarios, and raises again after the bottleneck. As predicted from the maize scenario, the direction of the evolution of the variance in gene expression was sensitive to the demographic scenario, and depends on a complex balance between drift, selection, and selfing rate. The evolution of plasticity was very similar in maize, African rice, and pearl millet. Interestingly, however, the extremely small ancestral population size of tomato hampered the evolution of plasticity. In all cases, the number of network connections evolved similarly as in maizenetwork rewiring at the onset of domestication, with an excess of connection gains vs. losses. Patterns of pleiotropy in the network (average expression correlation) were consistent with an overall strong short-and mild long-term increase, except for the recent marked expansion of population size in tomato that translated into a decrease in pleiotropy below the initial level ( Figure 7).

DISCUSSION
Domestication is a complex process, involving deep modifications of the demographic, environmental, and selective context in which populations evolve.
Here, we explored the consequences of domestication-like changes on the evolution of gene regulatory networks underlying domestication traits, combining a population bottleneck, directional selection and phenotypic canalization, simulated as the evolution of selection pressure towards decreased plasticity i.e., environmental stability of phenotypes.

Adaptive dynamics under domestication
We observed that the bottleneck had a substantial effect on genetic diversity, including (i) a substantial loss of neutral genetic (molecular) diversity ( Figure 4A) Less expected perhaps was the fact that even such a mild bottleneck penalized substantially the response to anthropic selection ( Figure 2). This may be due to less frequent occurrence of adaptive mutations during the bottleneck or a diminished efficiency of selection, or a combination of both. In the simulations, the bottleneck was associated with a burst of segregating adaptive alleles ( Figure 4B), which suggests a two-stage domestication scenario: (i) during the bottleneck, the adaptive alleles that segregate (either from the standing genetic variation and/or from new mutations) increased the population fitness, but tend to have suboptimal effects (e.g. negative side effects on well-adapted genes are illustrated by plastic genes whose reaction norm diminishes while they are continuously selected to be plastic, Figure   3B); (ii) after the end of the bottleneck, a new set of adaptive alleles can invade the population (because more mutations are available and selection is more efficient), fine-tuning genetic effects, e.g. on reaction norms ( Figure 3B). Hence, we expect mutations segregating during the first stage and surviving to drift to display greater effects than those segregating during the second stage. In line with this prediction, early work on maize domestication has identified several quantitative trait loci (QTLs) 21 495 500 505 510 with large effects, some of which were fine-mapped down to individual genes such as Tb1 (Doebley et al. 1997;Studer et al. 2011) and Tga1 (Wang et al. 2005).
Examples of early mutations with large effects have also been recovered in tomato (Frary et al. 2000), in wheat (Simons et al. 2006), in rice (Konishi et al. 2006;Li et al. 2006), in barley (Komatsuda et al. 2007)  Most phenotypic changes associated with domestication are controlled by mutations in transcription factors, and therefore involve a re-orchestration of gene networks (Martínez-Ainsworth and Tenaillon 2016) as described in cotton (Rapp et al. 2010), maize ), bean (Bellucci et al. 2014) and tomato (Sauvage et al. 2017).
Consistently, in our simulations, the gene network was deeply rewired, as the rate of gain/loss connections increased by more than one order of magnitude ( Figure 6A).
This effect was solely due to the shift in the selection regime. Before domestication, the population was well-adapted to an arbitrary wild type fitness landscape, involving genes which expression was constant and genes which expression was selected to track the environment. The structure of the underlying network evolved so that expressions of genes of the same type were genetically correlated, suggesting direct or indirect regulatory connections. When the fitness landscape changed, some genes that were previously correlated were forced to become independent. The results suggest that this was easier to achieve by adding connections rather than removing them, illustrating evolution by genetic tinkering instead of re-engineering.
Interestingly, there was no apparent cost to this additional complexity, as the fitness after domestication reached similar levels as before domestication.

Model approximations
Gene network models based on Wagner (1994)  For the sake of realism, and to connect the model results to quantitative genetics theory, we proposed several changes to the original framework from Wagner (1996).
We adopted the setting used in e.g. Siegal and Bergman (2002), in which gene expression was considered as a quantitative character, with a continuous scaling function between 0 (no expression) and 1 (maximal expression), instead of the traditional on/off binary setting (Wagner 1996;Ciliberti et al. 2007). We used an asymmetric sigmoid scaling (as in Rünneburger and Le Rouzic 2016) to ensure that a non-regulated gene has a low constitutive expression (here, 20% of the maximal expression). Our model allows for the possibility to evolve a plastic response. We starting state of the network, mimicking developmental plasticity (Masel 2004), or trans-generational plasticity (Odorico et al. 2018).
Computational constraints limited the population size to a maximum of N=20,000.
Estimates of the effective population sizes of both maize and teosinte vary roughly between 10 5 and 10 6 depending on data/methods/models (Eyre-Walker 1998; Tenaillon et al. 2004;Beissinger et al. 2016;Wang et al. 2017), suggesting that genetic drift before and after the bottleneck was substantially larger in the simulations than expected in a realistic domestication scenario. Domestication scenarios were also greatly simplified, with a single bottleneck. Refinements of this initial setting could include multiple expansion waves of semi-domesticated forms, as well as rapid population growth and gene flow with wild relatives post-domestication (Beissinger et al. 2016;Kistler et al. 2018). For simplicity, we parameterized the bottlenecks by setting the census size ( N ) to the effective population sizes ( N e ) documented in the literature. Because of selection, N e < N , which made bottleneck slightly stronger than expected. Yet, the difference was modest (< 10%, Figure 5), and was unlikely to affect the results. Larger population size would raise the neutral diversity, but is unlikely to impact general outcomes. Due to computational constraints, we also had to limit the number of generations prior to domestication ( T a ) for some simulations; as a consequence, "wild" populations were not necessarily at mutation-selection-drift equilibrium. However, the effect remains limited compared to the strong effects due to domestication (e.g., Figure 5A, T a =24,000 , vs. Figure 5C, T a =12,000 prior to domestication).
Likewise, network size also had to be limited to n=24 genes, as the complexity of the gene network algorithm increases with the square of the number of genes. Defining In the default scenario, most of the response to selection was due to new mutations, as the standing genetic variation alone could not explain more than about 20% of the (log) fitness recovery ( Figure S4). The contribution of standing genetic variation to the response to selection is a complex function of the mutational variance, the strength of stabilizing selection before domestication, and the strength of directional selection during domestication (Stetter et al. 2018 the phenotypic diversity of the wild ancestor. We also considered that the expression level of only half of the network genes was under direct selection pressurethis would happen if half of the TFs were regulating directly key enzymes or growth factors. Simulating twice less selected genes did not affect the qualitative outcomes of the model ( Figure S5).

The molecular syndrome of domestication
Simulations confirm that the domestication process is expected to be associated with several characteristic signatures (S) at the molecular level; S1: a decrease of allelic diversity, S2: a change in gene expression variance, S3: the rewiring of the gene regulatory networks, and S4: less modularity of coexpression patterns.
The loss of genetic diversity (S1) was both due to the bottleneck (genetic drift removed rare alleles from the population) and to the selection shift (selective sweeps decreased the genetic diversity at linked loci), it is thus expected to be a general signature of domestication ( Figure 4A). Empirically, a loss of genetic diversity is indeed always associated with domestication, although its amplitude may vary (reviewed in Gaut et al. 2015).
The direction and magnitude of the evolution of gene expression variance (signature S2) depends on the balance between selection and drift; bottlenecks tend to reduce diversity, while a shift in the selection regime tends to increase it transiently (segregation of adaptive variants). Given our simulation parameters, inspired from the maize domestication scenario featuring a mild bottleneck, expression variance increased ( Figure 4B). This was not necessarily the case with all parameter combinations, as a stronger bottleneck as in African rice led to a decrease in both 26 625 630 635 640 molecular and expression variance (Figure 7). The strength and the pattern of selection also affect the speed and the nature (soft vs. hard) of the selective sweeps, which may differ across species. As a consequence, domestication is not expected to be associated with a systematic evolution of gene expression variance: it may increase when the bottleneck is moderate, as in maize, or decrease in species where the bottleneck was drastic and/or associated with an autogamous mating system, such as rice, cotton (Liu et al. 2019), and beans (Bellucci et al. 2014).
Genetic networks were rewired (signature S3) and evolved towards less modularity ( Figure 5B), as a consequence of swapping the selection pattern among genes (shift in the optimal expression for stable genes, and loss of plasticity for others). The network was less plastic after domestication, which was a consequence of a modelling choice (domestication was associated with a decrease in the number of genes expected to respond to the environmental cue). New connections occurred among previously isolated modules, but former connections were not all eliminated.
As a result, the rewiring of regulatory connections lead to a moderate increase in gene coexpressions (signature S4), associated with a loss of structure in the coexpression network (uncorrelated genes became correlated, and strongly correlated genes became more independent). This illustrates a realistic evolutionary scenario towards non-adaptive complexity, where the final network structure is not the more efficient one, but rather results from the accumulation of successive beneficial mutations in an existing, constrained genetic background. Empirically, we therefore predict that connections involving genes targeted by domestication should pleiotropy, and less independent modules. Interestingly, the general increase in genetic correlations was associated with a trend towards homogenization, i.e. strong correlations tended to weaken whereas uncorrelated genes became slightly correlated ( Figure 5B). Empirical comparisons at 18 domestication-related traits between two independent populations of offspring generated by the intermating of multiple parents from a teosinte population and from a maize landrace, revealed several interesting features in line with our observations: only a subset of genetic correlations (33 out of 153) were conserved between teosinte and maize, teosinte correlations were more structured among trait groups (Yang et al. 2019). Investigating carefully the transcriptome evolution for several pairs of domesticated / ancestral populations will be necessary to assess the predictive power of our theoretical model.
The genetic diversity available in modern cultivated species is often considered as a limitation to further response to artificial selection. Controlling recombination has been proposed as crucial for plant breeders to engineer novel allele combinations and reintroduce diversity from wild crop relatives (reviewed in Taagen et al. 2020).
Yet, if the domestication syndrome was also associated with changes in the pleiotropy of the genetic architecture, genetic progress might also be limited by undesirable genetic correlations among traits of interest (Yang et al. 2019).
Understanding how genetic constraints evolved under anthropic selection and whether it is possible to avoid or revert them requires a better understanding of the complex non-linear mapping between domestication genes and phenotypes.
thank Clémentine Vitte for useful literature suggestions on enhancers. Simulations were performed on the core cluster of the Institut Français de Bioinformatique (https://www.france-bioinformatique.fr/ifb-core/).

FUNDING INFORMATION
This work received financial support from two grants overseen by     Average genetic correlation for each pair of genes, at generation -9000 (just before the onset of domestication) below the diagonal, and at generation 0 (last generation of the simulations), above the diagonal. Gene selection status is indicated (n: nonselected, s: stable, p:plastic); capital letters indicate genes whose selection status changed during domestication. Red lines delimitates gene categories, before and after domestication. The two groups of plastic genes correspond to genes selected to correlate positively and negatively with the environmental index.