Heterosis Is a Systemic Property Emerging From Non-linear Genotype-Phenotype Relationships: Evidence From in Vitro Genetics and Computer Simulations

Heterosis, the superiority of hybrids over their parents for quantitative traits, represents a crucial issue in plant and animal breeding as well as evolutionary biology. Heterosis has given rise to countless genetic, genomic and molecular studies, but has rarely been investigated from the point of view of systems biology. We hypothesized that heterosis is an emergent property of living systems resulting from frequent concave relationships between genotypic variables and phenotypes, or between different phenotypic levels. We chose the enzyme-flux relationship as a model of the concave genotype-phenotype (GP) relationship, and showed that heterosis can be easily created in the laboratory. First, we reconstituted in vitro the upper part of glycolysis. We simulated genetic variability of enzyme activity by varying enzyme concentrations in test tubes. Mixing the content of “parental” tubes resulted in “hybrids,” whose fluxes were compared to the parental fluxes. Frequent heterotic fluxes were observed, under conditions that were determined analytically and confirmed by computer simulation. Second, to test this model in a more realistic situation, we modeled the glycolysis/fermentation network in yeast by considering one input flux, glucose, and two output fluxes, glycerol and acetaldehyde. We simulated genetic variability by randomly drawing parental enzyme concentrations under various conditions, and computed the parental and hybrid fluxes using a system of differential equations. Again we found that a majority of hybrids exhibited positive heterosis for metabolic fluxes. Cases of negative heterosis were due to local convexity between certain enzyme concentrations and fluxes. In both approaches, heterosis was maximized when the parents were phenotypically close and when the distributions of parental enzyme concentrations were contrasted and constrained. These conclusions are not restricted to metabolic systems: they only depend on the concavity of the GP relationship, which is commonly observed at various levels of the phenotypic hierarchy, and could account for the pervasiveness of heterosis.

Heterosis, the superiority of hybrids over their parents for quantitative traits, represents a crucial issue in plant and animal breeding as well as evolutionary biology. Heterosis has given rise to countless genetic, genomic and molecular studies, but has rarely been investigated from the point of view of systems biology. We hypothesized that heterosis is an emergent property of living systems resulting from frequent concave relationships between genotypic variables and phenotypes, or between different phenotypic levels. We chose the enzyme-flux relationship as a model of the concave genotype-phenotype (GP) relationship, and showed that heterosis can be easily created in the laboratory. First, we reconstituted in vitro the upper part of glycolysis. We simulated genetic variability of enzyme activity by varying enzyme concentrations in test tubes. Mixing the content of "parental" tubes resulted in "hybrids," whose fluxes were compared to the parental fluxes. Frequent heterotic fluxes were observed, under conditions that were determined analytically and confirmed by computer simulation. Second, to test this model in a more realistic situation, we modeled the glycolysis/fermentation network in yeast by considering one input flux, glucose, and two output fluxes, glycerol and acetaldehyde. We simulated genetic variability by randomly drawing parental enzyme concentrations under various conditions, and computed the parental and hybrid fluxes using a system of differential equations. Again we found that a majority of hybrids exhibited positive heterosis for metabolic fluxes. Cases of negative heterosis were due to local convexity between certain enzyme concentrations and fluxes. In both approaches, heterosis was maximized when the parents were phenotypically close and when the distributions of parental enzyme concentrations were contrasted and constrained. These conclusions are not restricted to metabolic systems: they only depend on the concavity of the GP relationship, which is commonly observed at various levels of the phenotypic hierarchy, and could account for the pervasiveness of heterosis.

INTRODUCTION
The relationship between genetic polymorphism and phenotypic variation is a major issue in many branches of biology. It is usually referred to as the genotype-phenotype (GP) relationship, GP correlation or GP map (Lewontin, 1974). Strictly speaking, the genotype corresponds to the DNA level, whereas the phenotype can correspond to any trait observed or measured at any organizational level from transcription rate to fitness. In fact, following S. Wright (Wright, 1934), the GP relationship can also refer to the relationship between genetically variable phenotypic traits (e.g., enzyme activities) and more integrated phenotypic traits, i.e., traits measured at higher levels of the biological hierarchy (e.g., fluxes through networks).
Almost inevitably, the relationship between parameters at different organizational levels is non-linear, because the cellular processes that shape biological functions and structures are highly intertwined. The kinetics of biochemical and molecular reactions, transport and signaling mechanisms, positive and negative feedbacks, growth kinetics, etc., are intrinsically nonlinear. Thus, from allosteric regulation to fitness landscapes, from genetic regulatory circuits to developmental processes, there are myriads of examples showing complex relationships between genotype and phenotypes, or between phenotypic levels (Kauffman, 1993;Heinrich and Schuster, 1996). These GP relationships display a large diversity of shapes (Rendel, 1962;Pryciak, 2008;Felix and Barkoulas, 2015), but two prevail to a large extent: (i) The sigmoidal relationship ( Figure 1A). Cooperativity and synergy underlie such behavior (Carey, 1998;Kazemian et al., 2013), which is classically modeled using Hill functions. For instance in Drosophila, transcription of the hunchback gene in response to the gradient of Bicoid transcription factor concentration follows a typical S-shaped curve (Perry et al., 2012). Examples at other phenotypic levels exist, such as the relationship between the expression level of various genes and growth rate, a proxy for fitness, in Escherichia coli (Perfeito et al., 2011;Rodrigues et al., 2016) and in Saccharomyces cerevisiae (Rest et al., 2013). (ii) The concave relationship (Figure 1B), the archetype of which is the enzyme-flux relationship. The longstanding Metabolic Control Theory (MCT; reviewed in Fell, 2007) is probably the most elaborated corpus to describe this particular, but essential, case of genotypephenotype relationship. Both theoretical developments and experimental observations show that increasing from zero the value of an enzyme parameter usually results first in an increase of the flux value, then in attenuation of the response and finally in saturation. The phenotypic response resembles a rectangular hyperbola, and at the asymptotic upper limit there is robustness of the phenotype toward genotypic variation.
Such concave relationships are common at other organizational levels. For instance, it is observed during transcriptional regulation when transcription factors do not bind to DNA cooperatively (Giorgetti et al., 2010). In yeast, the rate of protein synthesis in response to variation in abundance of some 20 different translation factors is concave in most cases (Firczuk et al., 2013), and can be modeled using a deterministic model (Poker et al., 2014;Zarai et al., 2014). Other examples show the effects of variable gene/protein expression on the fitness of an organism. Izard et al. (2015), by engineering the control of RNA polymerase expression in E. coli, consistently observed a saturation curve for growth rate in diverse media. In mitochondrial diseases, the relationship between the proportion of mtDNA carrying pathogenic mutations and various cellular responses-mitochondrial protein synthesis rate, respiratory chain complex activity, mitochondrial respiration or ATP synthesis and growth rate-is concave, with a large plateau. Thus, an effect will only be detectable above a threshold value of the proportion of deleterious mtDNA (reviewed in Rossignol et al., 2003). In Parkinson's disease, the symptoms do not begin to appear until about 80% of neurons in the substantia nigra, which produce dopamine, have died. This is consistent with the hyperbola-shaped relationship between the proportion of alive neurons and the concentration of extra-cellular dopamine (Nijhout et al., 2015). In a wide range of species, allometry also provides relevant examples of non-linear relationships, such as between organismal mass and whole-organism metabolic rate, but also life span, tree trunk diameter, blood circulation time, etc. (West et al., 1999). Allometric equations are power laws whose scaling exponents are lower than 1, resulting in concave relationships.
In the above-mentioned studies, only one genotypic (or phenotypic) parameter is assumed to be variable. In fact, most integrated phenotypic traits are polygenic, i.e., they are potentially affected by a very large number of parameters, giving rise to complex and rugged phenotypic landscapes. However, from the previous considerations, one can expect global concavity of these multi-dimensional surface responses. Indeed, within the MCT framework, Dykhuizen et al. (1987) showed that the fitness response of E. coli to variation of both β-galactoside permease and β-galactosidase activity is a twodimensional hyperbolic-like surface. Similarly, Nijhout et al. (2003) analyzed quantitatively the pairwise effects of various components of the MAPK signaling cascade on the MAPK-PP output and observed largely concave response surfaces with saturation. MacLean (2010) also found a concave relationship between growth rate and the joint variation of transcription and translation mediated by streptomycin and rifampicin, respectively, in Pseudomonas aeruginosa. Another interesting example is the allometric relationship described in a population of recombinant inbred lines of Arabidopsis thaliana between, on the one hand, photosynthetic rate and, on the other hand, dry mass per area and age at flowering (Vasseur et al., 2012).
Whether the phenotypic response is sigmoidal or concave, there is often a horizontal asymptote which accounts for the strong asymmetry of the phenotypic response, depending on whether the value of the genotypic variable is increased or decreased. Thus, in metabolic engineering it is usually much easier to decrease than to increase flux and metabolite concentrations (Moreno-Sanchez et al., 2008). At the genomewide level, this asymmetry was recognized early on through the FIGURE 1 | The two most common types of GP relationships with possible inheritance. (A) The S-shaped relationship. There is convexity below the abscissa of the inflection point (red point) and concavity above. The three phenotypes P1, F1, and P2 are associated with the three genotypes A 1 A 1 , A 1 A 2 , and A 2 A 2 , respectively. MP: mid-parental value. In the convex region of the curve (mauve background), the low allele is dominant over the high allele (negative dominance), while in the concave region of the curve (yellow background) the high allele is dominant over the low allele (positive dominance). (B) The hyperbolic concave relationship. Whatever the genotypic values, there is positive dominance. In all cases there is additive inheritance of the genotypic parameters.
classical work of Lindsley et al. (1972) on the effects of segmental aneuploidy in Drosophila. Using Y-autosome translocations with defined breakpoints, the authors produced a large number of small duplications and deficiencies spanning chromosomes 2 and 3. Examining the phenotypic effects of having the same chromosomal region present in one, two or three copies, the authors found that monosomy had a much greater effect on viability than trisomy, revealing a strong asymmetry of the response to gene dosage. "The pervasive robustness in biological systems" (Felix and Barkoulas, 2015), which has been observed at all phenotypic levels (transcript, protein and metabolic abundance, metabolic flux, chemotaxis, cell cycle period, cell fate patterning, cell survival, etc.), most likely stems from this saturation effect. In some cases the mechanisms of robustness have been unraveled (Barkai and Leibler, 1997;Eldar et al., 2004;Wagner, 2013), and a widespread link is thought to exist between robustness and network complexity (Martin and Wagner, 2008;Whitacre, 2012), even though robustness could be partly adaptive (Ho and Zhang, 2016).
In addition to robustness, the "diminishing return" curves could account for various fundamental and apparently unrelated observations. The best-known example is the dominance of "high" over "low" alleles for metabolic genes, which was first pointed out by Wright (1934), then formalized in detail by Kacser and Burns (1981) in a seminal paper. Due to concavity, the phenotypic value of the heterozygote exceeds the homozygote mean value, and alleles with large deleterious homozygous effects are more recessive than alleles of weaker effect (Figure 1). The metabolic concave relationship also has major implications regarding the evolution of selective neutrality (Hartl et al., 1985), the variability of enzyme activity in populations under mutationselection balance (Clark, 1991), the relationship between flux and fitness (Dykhuizen et al., 1987), the epistasis between deleterious mutations (Szathmary, 1993), the response to directional selection (Keightley, 1996), the distribution of QTL effects (Bost et al., 1999, the dynamics of retention/loss of metabolic genes after whole genome duplications (Gout et al., 2009) and the evolution of genes in branched pathways (Rausher, 2013).
The model of positive dominance of Wright, Kacser and Burns is not restricted to metabolic genes, since it is valid whenever there is concavity of the GP relationship. In the case of a sigmoidal response, dominance of deleterious mutations may occur, depending on the allelic values relative to the inflection point of the curve (Gilchrist and Nijhout, 2001;Bost and Veitia, 2014) (Figure 1A). It would be quite a laborious task to try to estimate the respective proportions of positive and negative dominance in the literature. Nevertheless, the fact remains that there is a marked bias toward the dominance of "high" alleles. As early as 1928, Fisher (1928) noted that, in Drosophila, 208 out of 221 autosomal or sex-linked deleterious mutations were recessive. Similar proportions have been observed in all species studied to date, including in artificial diploids of the normally haploid Chlamydomonas (Orr, 1991). More recently, genomewide estimates from collections of deletion strains in yeast have been published (Steinmetz et al., 2002a;Deutschbauer et al., 2005). Whereas, 891 genes (20% of the genome) contributed to slow growth as homozygotes, only 184 (3% of the genome) were haploinsufficient (Deutschbauer et al., 2005). The mean values of the dominance coefficients were found to be h = 0.25 or less, which means that null alleles were on average recessive (Agrawal and Whitlock, 2011;Manna et al., 2011). Finally, predominant recessivity of deleterious mutations is fully consistent with common observations in population genetics, such as inbreeding depression (Charlesworth and Willis, 2009).
In this paper, we show that generalizing this dominance model to the multilocus case directly leads to a simple model of heterosis, i.e., a model that accounts for the common superiority of hybrids over their parents for quantitative traits (Shull, 1914;Jones, 1918;Gowen, 1952). Described more than 250 years ago (cited in Roberts, 1965) and studied in detail by Darwin (1876), hybrid vigor is of great importance for food production because it affects complex traits such as growth rate, biomass, fertility, disease resistance, drought tolerance, etc. In crops, it is currently used to derive hybrid varieties with exceptional performances. For instance in maize, the hybrid grain yield can exceed by more than 100% the mean parental yield (Sprague, 1983;Duvick, 2001). In rice, heterosis is about 20% for grain yield and about 10% for plant height (Xiao et al., 1995). In livestock and farms animals, such as cattle, swine, sheep and poultry, heterosis is also present and commonly exploited (Dickerson, 1973;Yang et al., 1999). This phenomenon is not restricted to species of agronomic interest. In the model plant A. thaliana, there is heterosis for biomass, rosette diameter and flowering date (Barth et al., 2003;Rohde et al., 2004;Ni et al., 2009). In animals, it has been described for diverse traits in drosophila (Houle, 1989), frog (Hotz et al., 1999), Pacific oyster (Launey and Hedgecock, 2001), owned dogs (O'Neill et al., 2013), etc. Heterosis has also been found in microorganisms, for instance for growth rate in Neurospora (Emerson et al., 1952) and S. cerevisiae (Steinmetz et al., 2002b). In the latter case it may be spectacular, with growth rate one order of magnitude higher in the hybrid than in the best parent.
Theoretical models of heterosis have been proposed based on the functioning of gene circuits (Omholt et al., 2000) or metabolic systems Fiévet et al., 2010). In the latter case the authors unraveled the links between heterosis, dominance and epistasis within the MCT framework. In support of this model, we show here that heterosis can easily be created in the laboratory. We relied on two complementary approaches: (i) Test-tube genetic experiments. We crossed in vitro "parents" that differed in the concentration of four successive glycolysis enzymes, and obtained "hybrids" displaying heterosis for flux through the pathway; (ii) In silico genetics. We crossed parents that differed in the concentration of 11 enzymes from the carbon metabolism network in S. cerevisiae, and we computed the input flux and two output fluxes in parents and hybrids using a system of differential equations. Heterotic fluxes were frequently observed. The conditions for having heterosis were analyzed using computer simulations and mathematically using the theory of concave functions. Small phenotypic difference between parents, together with contrasted and constrained parental enzyme concentrations, proved to be reliable predictors of heterosis. These conclusions are qualitatively valid beyond metabolic networks because they only depend on the shape of the GP relationship. Thus, the so-called "mysterious" phenomenon of heterosis appears to be a systemic property emerging from the non-linearity that exists at various levels of the biological hierarchy.
A common index of heterosis, the potence ratio (Mather, 1949): directly gives the type of inheritance: We frequently used in this study the index of best-parent heterosis: -If H WP ≥ 0, there can be -MPH, ADD, +MPH or BPH. The mean of the negative H WP values was noted H − WP .
Finally, we used in some instances the index of mid-parent heterosis: Note that in this context, various authors wrongly refer to heterosis as overdominance, complete dominance, partial dominance, and use various circumlocutions for negative heterosis. Strictly speaking, the words "dominance" and overdominance only apply to monogenic traits.
The Concave (Resp. Convex) Genotype-Phenotype Relationship Results in Positive (Resp. Negative) Heterosis , defined as the function obtained by fixing x j = x • j , ∀j = i, and letting x i be variable. Assume that ∀i,z =f (x i ) is strictly concave. Following the standard concavity argument, we can write: where 0 < t < 1. The additivity of genotypic variables corresponds to t = 1/2 for all i, so we have: Hence which means that +MPH or BPH necessarily holds, and that: Note that if the parents have the same phenotypic value (z 1 = z 2 ), it comes from Equation (1) that there is necessarily BPH. If ∀i,z =f (x i ) were strictly convex, Equation (1) would become: and there would be necessarily -MPH or WPH.

Application of These Theoretical Models to the Case of the Enzyme-Flux Relationship
Since Wright (1934), the enzyme-flux relationship is considered the archetype of the concave genotype-phenotype (GP) relationship. Using a simple system of differential equations, he relied on the concavity of the GP curve to propose his famous model of dominance. Kacser and Burns (1981) refined his approach using a formalism in which the relationship between enzyme concentrations and the flux through a metabolic pathway is described by a heuristic equation which is at the heart of the metabolic control theory (Kacser and Burns, 1981;Heinrich and Schuster, 1996;Fell, 2007). Similarly, but more generally, Fiévet et al. (2006Fiévet et al. ( , 2010 expressed the relationship between the overall flux J (M −1 .s −1 ) through a metabolic network and n enzyme concentrations {E j } j = 1,...,n in the following way: where E j is the concentration of enzyme j, X is a positive constant, A j and d j are systemic parameters (unit M −1 .s −1 ) that account for the kinetic behavior of enzyme E j in the network, and E tot is the total enzyme amount: Because in the standard case A j and d j are positive, the flux is an ascending function of E j , for all j. We will consider that A j and d j are not genetically variable, so that variation in enzyme activity will only be due to the variability of the enzyme concentration E j . Nevertheless formal developments with variable A j 's would also be possible since A j and E j play a similar role in the relationship.

Flux heterosis
The variables and parameters of Equation (2) are all positive, and the function J = f (E 1 , ..., E n ) has a negative partial second derivative with respect to E j , for all j. Therefore, the function is strictly concave, and the previous developments apply. In particular, consider two individuals P1 and P2 that differ in their distribution of enzyme concentrations, respectively {E j1 } j = 1,..., n and {E j2 } j = 1,..., n . Assuming additive inheritance of enzyme concentrations, the distribution of enzyme concentrations of the P1*P2 hybrid is: and its flux is: Due to concavity, we have for all j: Frontiers in Genetics | www.frontiersin.org Thus, either +MPH or BPH are expected in this GP relationship, and if the parents have the same flux value (J 1 = J 2 ) there is necessarily BPH.

Conditions for best-parent heterosis
In the particular case where there is only one variable enzyme, BPH cannot occur, as can be shown from Equation (2). Otherwise, with J 2 > J 1 , the condition for having BPH is J 1 * 2 > J 2 , i.e., After rearrangements, we get: where .
Note that ∀j, a j and b j are positive.
If the total enzyme concentration does not vary between individuals due to global cellular constraints of space and energy, we have E tot1 = E tot2 , and the condition for BPH is simplified: In this case, enzyme concentrations do not vary independently, and this covariation will increase heterosis because the concavity of the flux-enzyme relationship is higher. This can be shown by comparing the flux control coefficient (C J k ) (Kacser and Burns, 1981) to the combined response coefficient (R J k ) (Lion et al., 2004). Both coefficients measure the sensitivity of the flux to variations in enzyme concentration k, but the former applies when concentrations vary independently whereas the latter applies when there is a constraint on the total enzyme concentration. Following (Lion et al., 2004), we have: is the redistribution coefficient, i.e., the ratio of an infinitesimal change in concentration E j to an infinitesimal change in concentration E k . As ∀j, C J j > 0 and ∀(k, j), α kj < 0, we have R J k < C J k , which means that constraint on E tot increases concavity.

Optimal distribution of enzyme concentrations
The optimal distribution of enzyme concentrations is the one that maximizes the flux through the network for a given E tot value. The relevant unit of measure here is g.L −1 and not M.L −1 , because it is the total weight that matters in terms of energy and crowding in the cell. So we re-wrote Equation (2): where A ′ j = A j /M j , with M j the molecular weight of enzyme j, and where g E j is the concentration of enzyme j in g.L −1 .
The constraint is on We derived the optimal distribution of enzyme concentrations, { g E * i } i = 1,..., n , using the Lagrange multiplier method (see Fiévet et al., 2006 for details): In Vitro and in Silico Genetics
In order to mimic genetic variability in enzyme activity, we varied the concentrations of PGI, PFK, FBA and TPI in test tubes to create "parents, " each parent being defined by a particular vector of concentrations. HK concentration "Hybrids" were produced by mixing 1:1 the content of the parental tubes, which corresponds to additivity of enzyme concentrations. In order to have hybrids derived from a large range of parental values, we defined four classes of parental flux values (< 5 µM.s −1 , 5-8 µM.s −1 , 8-11 µM.s −1 and > 11 µM.s −1 ), and we performed "crosses" within and between classes. We created 61 hybrids using 36 different parents (Table  S1). For each of the 97 genotypes (61 hybrids + 36 parents), three flux measurements were performed (Fiévet et al., 2006), with the exception of two parents and four hybrids for which there were only two replicates (Table S2).
Two statistical tests were performed for every cross: (i) a Student's t-test to compare hybrid and mean parental fluxes; (ii) a one-way ANOVA with the genotype as factor (three levels: parent 1, parent 2 and hybrid) followed by a Tukey's test to classify observed inheritance as either BPH, +MPH, ADD, -MPH or WPH. For both, p-values were adjusted for multiple tests using the conservative Holm's method.

Computer Simulations of the in Vitro System
We simulated the in vitro system using the parameter values published by Fiévet et al. (2006), who estimated XA j and Xd j by hyperbolic fitting of the titration curves obtained by varying the concentration of each enzyme in turn. The values were XA PGI = 499.4 s −1 , XA PFK = 115.5 s −1 , XA FBA = 22.5 s −1 , XA TPI = 22940 s −1 , Xd PGI = 0, Xd PFK = 0, Xd FBA = 0 and Xd TPI = 21.9 s −1 . Thus, for any set of E j values, i.e., for any virtual genotype, the flux value could be computed from the equation: We analyzed the effects of three factors on heterosis: (i) mean parental enzyme concentrations: either equidistributed or centered on their optima; (ii) constraint on total enzyme concentrations: either free E tot or fixed E tot ; (iii) changes in enzyme concentrations, by varying their coefficients of variation. For equidistributed enzyme concentrations, the means of the gamma distributions were µ = 25.475 mg.L −1 (i.e., E ϕ tot /4) whatever the enzyme. For the optimum-centered distributions, the means were µ = 15.768 mg.L −1 , µ = 25.162 mg.L −1 , µ = 59.497 mg.L −1 and µ = 1.473 mg.L −1 for PGI, PFK, FBA and TPI, respectively, computed from Equation (7). In all cases we varied the coefficients of variation from c v = 0.1 to c v = 1.2, within the range of observed c v 's for these enzymes (Albertin et al., 2013). The shape and scale parameters of the gamma distributions were respectively k = 1/c 2 v and θ = µc 2 v . To apply a strict constraint on the total amount of enzymes allocated to the system, we assumed that for each parent i the total enzyme amount, j g E ji = g E toti , was 101.9 mg.L −1 , namely the physiological total enzyme concentration. After drawing enzyme concentrations under the different conditions described above (equidistributed/centered on the optima for different c v values), the concentration g E ki of enzyme k in individual i was computed as: where x ki was the result of the draw. Thus, the total enzyme concentration was 101.9 mg.L −1 for every individual. It is worth noting that the constraint on E tot changes the c v 's in an inverse relation to enzyme concentrations: the c v 's of the most (resp. the less) abundant enzymes will decrease (resp. increase) under constraint (see Appendix S1).  For every simulation condition, we drew two series of 10,000 parents which were "crossed" to get 10,000 hybrids. Parent and hybrid fluxes were computed according to Equation (6), assuming additivity of enzyme concentrations.

Computer Simulations of the Glycolysis and Fermentation Network in Yeast
We modeled the glycolysis and fermentation network in yeast using a simplified topology in which there was one input flux, the rate of glucose consumption, and two output fluxes, the rate of glycerol production and the rate of transformation of pyruvate into acetaldehyde (Figure 3). The system of ordinary differential equations we used was derived from Conant & Wolfe's model (Conant and Wolfe, 2007) obtained from the Biomodels database (http://www.ebi.ac.uk/biomodels-main/). This model is a slightly different version of a previously published model (Teusink et al., 2000). We modified this model in order to increase the range of enzyme concentrations that can lead to a steady state. First, we excluded the reactions producing trehalose and glycogen because, due to their constant rates, they caused incompatibility when the entry of glucose decreased. Second we suppressed mitochondrial shuttling, which is a composite reaction, because it is not associated with any enzyme. Thus, our model included 17 reactions and 22 metabolites, including 6 cofactors (ATP, ADP, AMP, NAD, NADH, and Fru-2,6-BP) (Figure 3). For all simulation conditions (see below), we solved the system of ordinary differential equations numerically using solvers (deSolve package) of the Runge-Kutta family (Hindmarsh, 1983) implemented within the R software tool (R Core Team, 2013). We considered that the steady state was reached when metabolite concentrations varied less than 10 −4 between two successive time steps.
Conant & Wolfe's model uses the maximal enzymatic rate V max , i.e., there is no explicit expression of enzyme concentration and catalytic constant k cat . To simulate variation in enzyme concentrations we needed to know the k cat values and molecular weights of the enzymes. We found most k cat values in the Brenda database (http://www.brendaenzymes.info/). When they were not known in yeast (in two cases), we used data from other organisms. Thus, we got k cat values for 11 enzymes (Table S3): HK, hexokinase (S. cerevisiae); PGI, phosphoglucose isomerase (S. cerevisiae); PFK, phosphofructokinase (Schizosaccharomyces pombe); FBA, fructose-1,6-bisphosphate aldolase (S. cerevisiae); PGK, 3phosphoglycerate kinase (S. cerevisiae); PGM, phosphoglycerate mutase (S. cerevisiae); ENO, enolase (S. cerevisiae); PYK, pyruvate kinase (S. cerevisiae); PDC, pyruvate decarboxylase (S. cerevisiae); ADH, alcohol dehydrogenase (Flavobacterium frigidimaris) and G3PDH, glycerol 3-phosphate dehydrogenase (Oryctolagus cuniculus). For these enzymes, we determined the reference concentrations E ref of the model from the equation: where M r is the molecular weight of the enzyme (Table S3). As expected, the flux response to variation in concentration of each enzyme, with the other 10 concentrations being maintained at their reference values, displayed saturation curves in most cases. However, for certain enzymes the reference concentration was far or very far along the plateau, so that varying the concentration was almost without effect on the fluxes. Therefore, we tried to arbitrarily decrease the reference values of these enzymes, but as a consequence the rate of successful simulations decreased dramatically. In the end, we solved this problem by decreasing the 11 reference concentrations by the same factor, which was empirically fixed at 5. Thus, we simulated the genetic variability of parental enzyme concentrations by drawing values from gamma distributions with parameters k = 1/c 2 v and θ = E ref c 2 v /5. Variances were set to have seven c v 's ranging from 0.1 to 0.7. With higher c v 's the steady state was frequently not reached, and c v 's smaller than 0.1 produced results similar to those obtained with c v = 0.1. For each c v value, we drew 10,000 pairs of parents. The enzyme concentrations of their hybrids were computed assuming additive inheritance. We considered two situations, free E tot and fixed E tot . For free E tot , we used the concentration values drawn as described above. For fixed E tot , the concentrations were scaled to make their sum equal to the sum of the reference enzyme concentrations (31833.97 mg), as x ki j x ji . We checked that the input flux (glucose) was equal to the sum of the output fluxes (glycerol and acetaldehyde), indicating that there was no "leak" in the system. For each parent-hybrid triplet we computed the fluxes of glucose, glycerol and acetaldehyde at the steady state. There were two possible sources of bias on the a posteriori c v 's. First, if one member of the triplet did not reach a steady state, the triplet was excluded. Thus, to get 10,000 simulations, ≈ 10, 060 to ≈ 36, 200 simulations were necessary, depending on the c v and the presence of constraint on E tot . Second, even if a steady state was reached, the glucose flux could be null, and we obviously eliminated the corresponding triplets. An analysis of the resulting biases is presented in Appendix S2. Overall, the differences between a priori and a posteriori c v values were insignificant or moderate when initial c v 's were small to medium, and were significant only for certain enzymes with high initial c v 's. In any case, these differences maintained the wide range of c v 's, making it possible to study the effects of variability in enzyme concentrations on heterosis.

Flux Differences and Enzymatic Distances Between Parents
The relationship between parental flux values and heterosis was analyzed using the absolute value of the flux difference between parents i and i ′ : The relationship between parental enzyme concentrations and heterosis was analyzed using the weighted normalized Euclidean distance between parents i and i ′ : where n is the number of enzymes. With normalization, we have 0 ≤ D enz ii ′ ≤ 1.
To examine the relationships between the heterosis indices and the distances, we used -when applicable -the Pearson correlation coefficient (r) and/or the coefficient of determination (R 2 ). Otherwise we used dedicated graphical representation methods (see Results section).

RESULTS
In order to analyze the way the genotype-phenotype (GP) relationship shapes heterosis, we used theoretical developments, in vitro genetics and computer simulations of the glycolytic/fermentation flux.

A Geometrical Approach: The Concavity of the GP Relationship Results in Heterosis
As an archetype of the concave GP relationship, we first chose a hyperbolic function relating the flux J i of individual i to enzyme concentrations in a n-enzyme metabolic system: where X is a positive constant, A j and d j are systemic positive parameters accounting for the kinetic behavior of enzyme E j in the network, and E ji and E tot i are respectively the concentration of enzyme E j and the total enzyme concentration in individual i (Fiévet et al., 2006(Fiévet et al., , 2010Fell, 2007). Because we assumed in this study that there was additivity of enzyme concentrations, the flux of the hybrid of parents P1 and P2, with fluxes J 1 and J 2 respectively, was written: whereĒ j = (E j1 + E j2 )/2 andĒ tot = (E tot1 + E tot2 )/2. Since the relationship between enzyme concentrations and flux is a rectangular multivariate hyperbolic function, the concavity of the surface results in positive heterosis for the flux (J 1 * 2 > J 1 + J 2 2 ), i.e., either positive Mid-Parent Heterosis (+MPH) or Best-Parent Heterosis (BPH) (see Figure S1 for definitions of the types of heterosis). This has been shown analytically (see Theoretical background), and is illustrated geometrically on the two-dimensional hyperbolic flux response surface obtained with two variable enzymes ( Figure 4A). Interestingly, it is possible to show that there is BPH, i.e., J 1 * 2 > J 2 (for J 2 > J 1 ), whenever: where a j and b j are positive terms (see Theoretical background). This equation predicts that heterosis will be observed only if a sufficient number of enzymes have a lower concentration in the parent displaying the highest flux than in the other parent, making (E j1 − E j2 ) positive, and this effect is strengthened if E tot is smaller in the parent that has the highest flux. These conditions are more likely to arise when parental fluxes are similar ( Figure 4A). In other words, a condition for observing BPH is that parents display contrasted distributions of enzyme concentrations, i.e., with complementary "high" and "low" alleles.
In the case of constrained resources, i.e., if the total enzyme amount allocated to the system is limited, which is biologically realistic, the flux-enzyme relationship is no longer a multivariate hyperbolic function but exhibits a hump-shaped response surface because enzyme concentrations are negatively correlated ( Figure 4B). Increasing the concentration of certain enzymes leads to decreasing the concentration of other enzymes (Heinrich and Schuster, 1998;Lion et al., 2004;Malmstroem et al., 2009). The constraint on E tot increases the concavity of the surface, which reinforces the incidence of heterosis (see Theoretical background). If E tot is the same in the two parents, the condition for BPH is simply written as: Note that when E tot is constant, BPH cannot be observed if a parent has optimal enzyme concentrations, because no offspring can exceed the maximal flux value.

In Vitro and in Silico Heterosis in a Four-Enzyme Pathway
In order to assess these predictions, we performed in vitro genetic experiments. We reconstructed the first part of glycolysis (Figure 2), and created 61 "hybrid" tubes by mixing 1:1 the content of "parental" tubes that differed in their distribution of the concentrations of four enzymes, PGI, PFK, FBA and TPI, with total enzyme concentration being constant (see Materials and Methods).
Inheritance was clearly biased toward high hybrid values (Figure 5). Among the 61 hybrids, 27 (≈ 44.3 %) exhibited a significantly higher flux value than the mid-parental value: seven (≈ 11.5%) displayed BPH and 20 (≈ 32.8%) displayed +MPH (p-values < 0.05). The maximum BPH value, quantified with the index H BP = J 1 * 2 −max(J 1 ,J 2 ) max(J 1 ,J 2 ) , was H BP = 0.37, i.e., the hybrid displayed a flux that was 37% higher than the best parental flux. By contrast, only eight hybrids had a flux value that was significantly lower than their mid-parental value (negative Mid-Parent Heterosis, or -MPH), with no case found of significant Worst-Parent Heterosis (WPH) (Figure 5 and Table S1).
Examining enzyme concentrations in the parents confirmed that high values of BPH were likely to occur when parental distributions were contrasted and when neither parent had enzyme concentrations close to the optimal distribution (Figures 5D,E). However, as shown in Figure 6A, this relationship is not straightforward, with large distances being necessary but not sufficient to get high H BP values, which resulted in a non-significant correlation between H BP and D enz , the Euclidean distance between parental enzyme concentrations.
To test the prediction that a cross between parents with similar fluxes is more likely to give heterosis than a cross between parents with contrasted fluxes, we computed the correlation coefficient between H BP and D flux , the difference between parental fluxes (| J 1 − J 2 |), and found a highly negative significant value (r = − 0.33, p < 0.01) ( Figure 6B).
Interestingly, even though H BP and D enz were not correlated, the heterosis value H reg , computed from the equation of the multiple linear regression performed with both D enz and D flux as predictor variables, was better correlated with H BP than D flux alone: r = 0.46 (p < 1.8.10 −4 ) ( Figure 6C).
Finally, in order to know if the instances of heterosis we observed in vitro were consistent with our heuristic model of flux-enzyme relationship, we computed from Equations (9) and (10) the theoretical fluxes of the 61 hybrids and their parents, using the values of enzyme concentrations of the parental tubes and the values of parameters XA j and Xd j estimated in an independent experiment (see Materials and Methods). Then, for each virtual cross, we calculated the heterosis that would be expected from the model, H BPmod . As shown in Figure 6D, the correlation between H BPmod and H BP is very highly significant (r = 0.77, p < 2.5.10 −13 ). Red and blue curves show two types of crosses between parents P1 and P2, or P3 and P4, respectively (from Fiévet et al., 2010). The positions of hybrids F1 (light blue and orange points) were determined assuming additivity of enzyme concentrations. (A) In the "red" cross, where parent P2 has a flux close to the maximum (high concentration of both enzymes), there is positive mid-parent heterosis (+MPH) for the flux. In the "blue" cross, where parents have low flux values due to low concentrations of enzyme 1 (parent P3) or enzyme 2 (parent P4), the hybrid displays best-parent heterosis (BPH). (B) If there is a constraint on the total enzyme amount, the "red" cross also displays BPH because the concavity of the surface is increased. (9) and (10) resulted in quite reliable predictions of flux heterosis, we used them to simulate a large series of crosses. The parental enzyme concentrations were drawn to be centered on their optimum or equidistributed, and their coefficients of variation (c v ) varied from 0.1 to 1.2. In every case we considered two situations, free and fixed E tot . For every simulation condition, 10,000 crosses were performed and the parental and hybrid values were computed.

As Equations
The percentage of BPH and the mean of positive H BP values (H + BP ) were compared over the range of c v 's. Figure 7 shows that: (i) in all conditions, the values of these two variables increased, usually to a large extent, when enzyme variability increased; (ii) by far the highest % BPH (about 50%) was observed when enzyme concentrations were constrained (fixed E tot ) and centered on their optimum ( Figure 7A); (iii) H + BP was highly dependent on the c v values (varying from ≈ 0 to ≈ 1.4), but was barely affected by the other simulation conditions ( Figure 7B).
The typical influence of D flux on H BP is illustrated Figure 8 for c v = 0.6. BPH was observed when the parental fluxes were close to each other, as the hybrids displaying BPH were distributed around the diagonal of the parental flux space. This pattern is consistent with the negative relationship between D flux and H BP (Figure 9A). The highest D flux never resulted in BPH, while small D flux could result in high H BP values. The weakest correlation between D flux and H BP was observed when both enzyme concentrations were centered on their optimum and E tot was fixed (r = −0.35 instead of r > −0.66 in other conditions).
The graphs showing the relationship between D enz and H BP revealed a typical triangular relationship, regardless of the simulation conditions: BPH was never observed with the smallest distances, and large distances were necessary but not sufficient to get high H BP values ( Figure 9B).
To assess the joint influence of D flux and D enz on heterosis, we divided the crosses into 16 classes according to a 4×4 grid of D flux and D enz values, and computed for each class the % BPH and H + BP .
In all conditions, the highest values of these two variables were observed with both small D flux and large D enz (Figures 9C,D).
As expected, constraint on E tot resulted in more heterosis, with almost 100% BPH when: (i) enzyme concentrations were centered on their optimum; (ii) D enz values were maximal; (iii) D flux values were minimal ( Figure 9C). In these conditions, H + BP increased to 0.37, i.e., hybrids had on average a flux that was 37% higher than the best parental flux ( Figure 9D).

In Silico Heterosis in the Yeast Glycolysis/Fermentation Network
To simulate heterosis in a larger system, we modeled the glycolysis and fermentation network in yeast using a simplified topology in which there was one input flux, the rate of glucose consumption, and two output fluxes, the production rates of glycerol and acetaldehyde (Figure 3). We used a system of differential equations derived from Conant and Wolfe (2007) and varied the concentration of 11 enzymes. We first examined how varying the concentration of one enzyme at a time, the other concentrations being maintained at their reference values (see Material and methods), affected flux variation. When E tot was free, the enzyme-flux relationships were concave with a horizontal asymptote in most cases, although a slight convexity was visible in two instances ( Figure S2): (i) sigmoidal curves were observed upon variation of HK, an allosteric enzyme that is by far the most abundant in the system (Table S3); (ii) due to the balance between the acetaldehyde and glycerol branches, increasing the concentration of an enzyme in one branch convexly decreases the flux of the other branch (except for enzyme at low concentrations). The decrease of the glucose flux upon G3PDH increase is likely to be due to the related reduction of the acetaldehyde flux, which dampens ATP production and in turn reduces the glucose flux.
When E tot was fixed, these observations remain qualitatively valid, however we no longer found horizontal asymptotes:   due to this constraint, flux values decreased when enzyme concentrations were high ( Figure S3). This effect was marked for the most abundant enzymes (particularly HK and PYK) but was imperceptible for less abundant ones (e.g., PFK). Note that in all situations the three fluxes are null when enzyme concentrations are null, even for enzymes downstream of the bifurcation. Cofactor availability probably explains this observation.
We then simulated 10,000 crosses between parents that differed in their concentrations of the 11 variable enzymes. Parental concentrations were drawn from gamma distributions, the means of which were the reference values. Variances were chosen within the range of seven c v 's, from 0.1 to 0.7, with free or fixed E tot . Hybrid fluxes and inheritance were computed for each cross (see Materials and Methods). Unlike the previous four-enzyme model, where the enzyme-flux relationship was concave by design and could only produce BPH or positive MPH, we observed the four possible types of inheritance, BPH, +MPH, -MPH and WPH. For the glucose and acetaldehyde fluxes, inheritance was very biased toward positive heterosis whatever the c v , since the sum of the percentages of +MPH and BPH varied from 69.4 to 96.9%, and there was no or very few cases of WPH. The highest percentage of WPH, observed for acetaldehyde when c v = 0.7 with free E tot , was only 0.09%, with a very small mean of the negative H WP values (H − WP = −0.024) (Tables 1, 2 and Figure 10). The one-step glycerol branch, which is convexly related to the six enzymes of the acetaldehyde branch, had higher percentages of negative heterosis (from ≈ 65% when c v =0.1 to ≈ 43% when c v = 0.7 with fixed E tot ), but WPH remained quite low: 1.23 and 3.02% for free and fixed E tot , respectively, with H − WP ≈ −0.04 in both cases (c v = 0.7). For all fluxes, constraint on E tot decreased the percentage of -MPH in favor of BPH and/or +MPH, depending on the c v value (Tables 1, 2 and Figure 10).
We hypothesized that negative heterosis can occur only when the enzyme-flux relationship is convex, as exemplified for negative dominance (Figure 1A). To verify this, we ran simulations where HK concentration was maintained at its reference value and/or the glycerol branch was deleted. When both sources of convexity were removed, no case of negative heterosis was observed (Table S4 and Figure S4).
As previously, both % BPH and H + BP increased as the c v increased (Tables 1-3, Figures 11A,B). This effect was larger for glucose and acetaldehyde than for glycerol. For instance, for glucose the % BPH varied from 3.76 to 35.12% with free E tot and from 11.96 to 42.46% with fixed E tot , while for glycerol it varied from 2.92 to 18.84% with free E tot and from 3.71 to 24.10% with fixed E tot . For the three fluxes, constraint on E tot increased the % BPH, as expected, but decreased H + BP . The latter observation is explained by the fact that this constraint decreases the variance of the concentrations, which in turn decreases the flux variance and hence both H BP variance and H + BP (see Appendix S1 and Appendix S2).
Whatever the c v , BPH was usually observed when parental fluxes were similar, with heterotic hybrids roughly distributed around the diagonal of the parental flux space, and this was valid for the three fluxes ( Figure S5). The correlations between H BP and D flux were consistent with these observations. They were negative and highly significant, although slightly smaller when E tot was fixed (r = −0.60, −0.69 and −0.61 for glucose, glycerol and acetaldehyde, resp.) than when E tot was free (r = −0.55, −0.59 and −0.54, resp.) (Figures 12A,C).
The triangular relationship between H BP and D enz was slightly less apparent with free E tot than in the four-enzyme model, but was obvious with fixed E tot . Again high H BP values were never observed with small distances, and the highest H BP values could only be observed with medium or large distances (Figures 12B,D).
As previously, we assessed the joint influence of D flux and D enz on % BPH and H + BP using 4 × 4 grids (Figure 13, for c v = 0.4). Results were very similar: both variables were positively related to D enz and negatively related to D flux , whatever the flux. High D enz associated with small D flux could result in about 50% BPH and H + BP ≈ 0.25 when E tot was free and almost 100% BPH and H + BP ≈ 0.40 when E tot was fixed. The positive effect of the constraint on H + BP is not in contradiction with the above-mentioned argument: H + BP was computed here for each square of the grid, and its  Figure S6 shows that contrasted distributions of enzyme concentrations can lead to BPH whereas close distributions favor additivity. Indeed, the cross that is the most heterotic for both glucose and acetaldehyde (H BP = 4.23 and H BP = 5.55, resp.) had divergent enzyme concentrations (D enz = 0.68) ( Figures  S6A,B); this was also observed for the most heterotic cross for glycerol (D enz = 0.66) (Figures S6C,D). These enzymatic distances are close to the highest distance observed (D enz = 0.71), which also corresponded to large heterosis for glucose and acetaldhyde (H BP = 0.69 and H BP = 1.09, resp.) (Figures S6E,F). By contrast, the cross with the minimum enzymatic distance (D enz = 0.16) exhibited additivity or near additivity for the three fluxes (H MP = 0, H MP = 0.02 and H MP = −0.01 for glucose, glycerol and acetaldehyde, resp.) (Figures S6G,H).
The geometric models of Figure 14 and Figure S7 show a relationship between D enz , D flux and H BP that is fully consistent with all our observations. For a given value of D enz , H BP can be positive or negative depending on the position of the parents in the enzyme concentration space, with quite complex relationships between these two variables. BPH is observed when parents have complementary enzyme concentrations, with the highest H BP values observed when D flux is small. Of note, D flux and H MP are roughly positively related, because the reference for this index is the mid-parent and not the best parent (Figures 14E,F). Regarding the constraint on E tot , comparison of Figure 14A and Figure S8 shows how greater surface curvature increases the incidence of heterosis.

The So-Called Mystery of Heterosis
Nearly all papers on heterosis mention in their introduction that is it an elusive, enigmatic or even magical phenomenon, whose molecular (or physiological, or biochemical) bases are not known or not well understood, and which would be the "lasting mystery  in biology." Yet, this "fascination with the idea that standard genetic models are not sufficient to explain heterosis" (Kaeppler, 2011) is quite surprising given the literature on this phenomenon. First, regardless of any genetic or molecular basis, it was recognized early on (Richey, 1942;Grafius, 1961) that multiplicative traits, i.e., traits determined as the product of two or more components, display commonly predictable levels of heterosis. For instance, heterosis for plant height in bean is quite well explained by multiplying internode length by internode number (Coyne, 1965); similarly, yield heterosis in barley is predicted by the product of ears/plant, kernels/ear and average kernel weight (Grafius, 1959), etc. It is worth noting that this cause of heterosis does not even require non-additive inheritance of the components of the trait (reviewed in Schnell and Cockerham, 1992).
Second, heterosis has been explained at the genetic and/or molecular level for various traits in various species. Paradoxically, the best understood cases of heterosis are certainly the least common, since they involve only one locus, which makes them easier to study. Here, heterosis is due to overdominance, i.e., the inherent superiority of the heterozygote over both homozygotes (East, 1936;Hull, 1946). Relative to the innumerable examples of complete or partial dominance, well-attested overdominance remains quite rare, with probably less than 20 cases described over more than one century of genetics. The bases of overdominance are diverse: (i) Pleiotropy. The individual traits controlled by the gene harbor only dominance, but the direction of dominance is reversed between the traits, which results in increased viability or fitness of heterozygotes compared to the parental homozygotes (Allison, 1954;Vrebalov et al., 2002;Lippman and Zamir, 2007); (ii) Dosage effects. When there is an optimal expression level of the gene controlling the trait, the hybrid will be heterotic whenever its expression level will be closer to the optimum than that of its parents (Hall and Wills, 1987;Flintham et al., 1997;Krieger et al., 2010); (iii) Formation of favorable oligomers in the hybrids (Tahiri-Alaoui et al., 2006;Singh et al., 2013).
Understanding heterosis when several-possibly a wealth ofpolymorphic genes control the trait is challenging. However, this has been done when a limited number of genes are involved. As early as 1910, an example of heterosis accounted for by complementary dominance at two loci L (controlling internode length) and T (controlling stem thickness) was found in garden pea (Keeble and Pellew, 1910), in line with Davenport's model (Davenport, 1908), which stated that heterosis was due to recessive unfavorable alleles being masked at loci in repulsion. Since then, similar examples have been described (e.g., Warner et al., 1969). It is worth noting that many alleged overdominant QTLs (Quantitative Trait Loci) could actually correspond to pseudo-overdominance, where two (or more) linked dominant QTLs are in repulsion (Eshed and Zamir, 1995;Graham et al., 1997;Lariepe et al., 2012;Martí-Raga et al., 2017). Epistasis also accounts for heterosis through favorable intergenic interactions created in the hybrids (Powers, 1944). The textbook case (Sinnott et al., 1950), which involves both dominance and epistasis, was recognized very early on (Bateson and Gregory, 1905) and was illustrated in various plants: crossing two individuals with colorless flowers may produce hybrids with colored flowers because the parents each have a mutation inactivating a particular step of the pigment biosynthesis pathway, which disrupts the flux, and the flux is restored in the hybrid because both steps are active (e.g., Dooner et al., 1991 for anthocyanins in maize). Other cases of heterosis that have been elucidated at the molecular level also combine more than one genetic effect. For instance in yeast, heterosis for growth at high temperature involved both dominance and epistasis at three genes MKT1, END3, and RHO2 (Steinmetz et al., 2002b;Sinha et al., 2006), and heterosis for fermentation performance proved to be partially explained by dominance and pseudo-overdominance at four genes VMA13, PDR1, PMA1, and MSB2 (Martí-Raga et al., 2017). Similar complex effects were found in the circadian clock genes CCA1 and LHY, which play a central role in heterosis for chlorophyll and starch content in A. thaliana (Ni et al., 2009). In tomato, epistasis and overdominance at the SFT gene account for yield heterosis (Krieger et al., 2010).
Third, given the pervasiveness of non-additive genetic effects (Veitia, 2006;Phillips, 2008), it is the lack of heterosis that would be a mystery. Whenever there is dominance, and/or overdominance and/or epistasis at loci governing a quantitative trait, heterosis can arise. Any one of these three genetic effects is necessary, but also sufficient, to generate heterosis. Besides, their respective importance in producing heterosis is variable, and depends on the mating system, the traits under study and the genetic material used, as shown from QTL/association mapping analyses and other approaches. In maize, dominance seems to be the main factor producing heterosis for grain yield and its components, with a modest role of overdominance and possibly epistasis (Garcia et al., 2008;Guo et al., 2014;Mezmouk and Ross-Ibarra, 2014). In rice, dominance was first invoked as the main factor for heterosis of grain yield and its various components (Xiao et al., 1995), but subsequently overdominance and additive x additive or dominance x dominance epistasis have been repeatedly found for various yield-related traits (Yu et al., 1997;Li et al., 2001;Luo et al., 2001;Hua et al., 2003;Mei et al., 2005;Garcia et al., 2008;Huang et al., 2016). In Arabidopsis the main factor was claimed to be epistasis for seven growthrelated traits (Melchinger et al., 2007), however, for biomass the three genetic effects seemed to be involved (Meyer et al., 2010), and a recent study showed the role of dominance and overdominance for flowering date, rosette diameter and rosette biomass (Seymour et al., 2016). In upland cotton, the three effects contribute to heterosis for yield and yield components (Shang et al., 2016). In yeast, the three effects also play a role in heterosis for growth rate measured in five different culture media (Shapira et al., 2014), while overdominance and synergistic epistasis seem mainly involved in growth at high temperature (Shapira and David, 2016). It is worth noting that in most cases it is not possible to distinguish true overdominance from pseudo-overdominance, even if indirect arguments suggest true overdominance (Semel et al., 2006).
Despite the multiple cases where the bases of heterosis have been clarified, some authors have been searching for a "unifying principle" (Birchler et al., 2003), such as genomewide changes in DNA methylation (Tsaftaris and Polidoros, 1999;Shen et al., 2012), small RNA expression and epigenetic regulation Groszmann et al., 2011;Chen, 2013), reduced metabolic cost of protein recycling in hybrids (Goff, 2011), gene dosage effects in macro-molecular complexes (Veitia and Vaiman, 2011) , enhanced metabolic efficiency due to weak co-aggregation of allozymes (Ginn, 2017), mitochondrial complementation (McDaniel andSarkissian, 1966;Srivastava, 1981) and phytohormonal expression (Rood et al., 1988). These molecular processes may indeed distinguish heterozygotes from homozygotes, but if they were the general hidden causes of heterosis, correlations between levels of heterosis for different traits would be observed, and this has never been reported so far (Flint- Garcia et al., 2009;Kaeppler, 2012). In fact, searching for the mechanistic bases of heterosis of a certain trait in a certain cross is equivalent to searching for the bases of its underlying genetic effects, and dominance, epistasis and overdominance have never been considered "mysterious" phenomena.

Curvature of the Genotype-Phenotype Relationship and Heterosis
The "unifying principle" proposed in this paper is not mechanistic. It is a general framework built on the common observation of the concavity of the GP relationship. The two genetic consequences of this non-linearity -dominance and epistasis -have been known for a long time in the context of the metabolic control theory (Wright, 1934;Kacser and Burns, 1981). More recently, modeling gene and metabolic networks have led to the idea that heterosis could also be a systemic property emerging from non-linear processes in the cell (Omholt et al., 2000;Fiévet et al., 2010). We provide strong arguments in support of this proposition by taking advantage of theoretical developments, and using both experimental approaches and computer simulations.
Experimentally, we performed in vitro 61 crosses between "parents" that differ in the concentrations of four enzymes from the upstream part of the glycolysis pathway, and observed inheritance that was biased toward positive heterosis despite additivity of enzyme concentrations in the hybrid tubes: about 44% of hybrids displayed positive MPH or BPH, whereas only 13% displayed negative MPH and none WPH. In addition, the H BP index could increase up to ≈ 0.37, which means that with enzyme concentrations half-way between the parents, the hybrid had a flux that was 37% higher than the best parental flux. From a geometrical standpoint, these results indicate that the GP hypersurface in this in vitro system is globally concave, accounting for the predominance of positive heterosis, the cases of negative heterosis revealing local unevenness of the surface. We did not explore further the possible causes -metabolic and/or technical -of these irregularities.
Computer simulations based on the previous system and using a hyperbolic GP relationship allowed us to perform a large number of crosses. Enzyme concentrations were drawn, varying means and coefficients of variation, and total enzyme concentrations were either free or fixed. In all conditions only +MPH and BPH were observed, which was expected due to the concavity of the GP relationship. The results of the simulations together with our theoretical developments allowed us to identify the factors favoring heterosis: -A contrast between parental enzyme concentrations. A high positive H BP could only be observed when D enz was large. Accordingly high c v 's resulted in the highest % BPH and H + BP . However, high D enz is necessary but not sufficient to get BPH: the position of the parents in concentration space is also a key factor (Figure 14 and Figure S7). If a parent contains most of the high alleles and the other most of the low alleles, BPH is unlikely because the hybrid can hardly exceed the best parent, while if the parents are complementary for the high and low alleles, the hybrid is expected to display BPH. -Phenotypic proximity between parents. In the case of strict concavity, BPH is very likely when the parents have close phenotypic values, and the % BPH decreases as the parental fluxes diverge. If a parental value is close to the maximum, MPH alone will be observed (Figure 14). -Constraint on total enzyme concentration allocated to the system. All other things being equal, the % BPH is higher when E tot is fixed, due to increased concavity of the flux response surface ( Figure 4B).
Combining fixed E tot , a large enzymatic distance and a small phenotypic difference makes BPH almost inevitable, particularly when mean enzyme concentrations are centered on their optimal values. This is consistent in geometric terms: when these conditions are met, the hybrids are closer to the "top of the hill" than their parents ( Figure S8). We performed similar simulations from a branched network based on glycolysis and fermentation in yeast, modeled with a system of differential equations. Three fluxes were considered: glucose input and glycerol and acetaldehyde outputs. Again BPH could only arise when the parents were phenotypically close to each other, but due to some convex enzyme-flux relations in this system, negative heterosis was also observed. For glucose and acetaldehyde fluxes, positive MPH and BPH prevailed to a large extent, since the percentage of positive heterosis varied from 69.5 to 96.9%, and for all fluxes and conditions the % BPH increased with the c v , up to ≈ 44% when E tot was fixed. For glycerol, positive heterosis was predominant only for medium and high c v 's. In all cases WPH remained very rare. Thus, the three fluxes, even though they are pleiotropically related, do not have the same type of inheritance. However, the factors favoring BPH are the same as before for the three fluxes, namely constraint on E tot and a large enzymatic distance together with a small phenotypic divergence. Therefore, these factors do not depend on a particular flux-enzyme function, but only on the shape of the GP relationship. In addition we see that heterosis is a systemic property which enables a better exploitation of cell resources, since with a distribution of enzyme concentrations closer to the optimum distribution than both parents, a hybrid can "do more with less, " i.e., exceeds the best parent with a lower total enzyme concentration. In epistemological terms, it is an emergent property, in the sense that the individual properties of the molecular components (transcription/translation factors, mRNAs, proteins/enzymes, metabolites, etc.) cannot alone account for the phenomenon.
In order to avoid confounding effects, in particular in the glycolysis/fermentation system, we did not consider in this study possible non-additive inheritance of enzyme concentrations. Non-additivity of genetic variables has logical consequences: if it is positive, there is a higher occurrence of positive heterosis, while negative non-additivity can result in negative heterosis even in the case of a strictly concave GP relationship (Fiévet et al., 2010;Blein-Nicolas et al., 2015). In actual fact, transcripts and proteins display mostly additive inheritance (Stupar et al., 2008), and in the case of non-additivity it is usually biased toward positive values (Swanson-Wagner et al., 2006;Hedgecock et al., 2007;Shen et al., 2012;Chen, 2013;Blein-Nicolas et al., 2015;Xing et al., 2016;Guo et al., 2017), which will reinforce positive heterosis. For certain gene/protein categories that are found to be non-additive (histone modifications, small RNAs, protein and carbon metabolism, ribosome proteins, photosynthesis, stress response and energy production), a causal link with heterosis has been suggested or evidenced (Chen, 2013;Xing et al., 2016). Incidentally, positive non-additivity at the molecular level also reflects concave GP relationships: transcript, protein and metabolite abundances are polygenic traits (Damerval et al., 1994;Brem et al., 2002;Kliebenstein, 2009) that behave like any quantitative trait and can display non-additive inheritance.
How can our results be applied to whole-cell metabolism? Crowding (Ellis, 2001) and energy cost (Dekel and Alon, 2005;Vilaprinyo et al., 2010) create constraints on the total content of cell components. For instance in yeast, Albertin et al. (2013) showed that the enzymatic pool allocated to the glycolysis, glycerol, acetate, and ethanol pathways was invariant regardless of the culture medium and strain considered. When some proteins are overexpressed, this constancy can sometimes lead to burden effects, which causes a decline in flux (Snoep et al., 1995), as we found under constrained E tot . However, due to the huge number of molecular species in a cell, the increase of particular components does not usually impair the whole system. Similarly, situations of convexity within cell networks are expected to be buffered at the scale of the whole flux of matter and energy through the system. Thus, despite global constraints and local convexity, plateaued GP relationships, either S-shaped or more often concave, are commonly observed at integrated levels of cell organization, up to fitness, and most likely account for the pervasiveness of positive heterosis. Therefore, heterosis could almost be viewed as the measurement of the curvature of the GP relationship.
Our results may have practical consequences for heterosis prediction, a major issue in plant and animal breeding. A classical predictor of heterosis is the genetic distance between the parents, but the quality of the prediction depends markedly on the traits, on the genetic history of the material and on the range of genetic distances considered (Moll et al., 1965;Melchinger et al., 1992;Charcosset and Essioux, 1994;Flint-Garcia et al., 2009;Seymour et al., 2016). Our results show that better predictions of heterosis could be achieved by combining parental information at two (or more) phenotypic levels, namely the level of the trait of interest and the level (or levels) of underlying components. Omics methods allow thousands of transcripts, proteins and metabolites to be identified and quantified. Thus, we now have access to a plethora of putative predictors (Riedelsheimer et al., 2012) to be incorporated in non-linear statistical models. The high-dimensionality is challenging, but relevant variables could be selected by using both biological knowledge of the trait under study (e.g., Miller et al., 2015) and recently developed regularization methods (Dalalyan and Tsybakov, 2012;Heinzl and Tutz, 2014).