Introduction

Environmental canalization and phenotypic plasticity

An organism may show a range of phenotypes in different environments. A phenotypically plastic genotype shows a broad range of phenotypes (Schmalhausen, 1949). Homeostasis, or canalization, refers to the opposite response, when the individual is insensitive to changes in the environment (Waddington, 1942). Thus, environmental canalization and phenotypic plasticity can be viewed as the two possible extremes of response to changes in environment (see Flatt, 2005). Because different genotypes of the same species may display a wide range of variation in the level of their plasticity response (Breitkreutz et al., 2003; Ungerer et al., 2003; Ros et al., 2004; Turner, 2004), there must be a genetic basis of response to environmental change. Alleles at loci that affect variation in a phenotype across environments should therefore be considered as determinants of plasticity and canalization.

Three mechanisms were originally proposed to explain the genetic basis of phenotypic plasticity (see for synthesis Scheiner and Lyman, 1989). In the overdominance model, the degree of stability (or homeostasis) of a genotype is proportional to its degree of heterozygosity for relevant genes: an individual homozygous at these loci shows higher plasticity (Gillespie and Turelli, 1989). In the pleiotropic model, it is the differential expression across environments of alleles at a locus that contributes to phenotypic plasticity (Via and Lande, 1985). In the epistatic model (Scheiner and Lyman, 1989), plasticity genes interact epistatically with, and thus regulate, the ‘constitutive’ loci that determine the mean value of the trait.

Is QTL detection useful for understanding the genetic basis of plasticity?

One approach to testing the genetic models of phenotypic plasticity is to examine the genetic basis of quantitative trait variation, with a special focus on the genetic mechanisms shaping phenotypic variation across environments. QTL analysis allows dissection of quantitative variation by detecting the underlying genetic determinants (Stuber et al., 1992). These determinants are characterized in terms of their effects, modes of action and interactions with other genetic loci (QTL epistasis) and the environment (QTL × E) (Jansen et al., 1995; Korol et al., 1998). Two basic patterns of QTL × E interaction are recognized, they are: (i) inversion in the ranking of allelic effects across environments (the crossover effect, CO) and (ii) variation in the magnitude of allelic effects without inversion of ranks (the scale effect, SC) (Yadav et al., 2002).

QTL showing QTL × E interaction should be implicated in the control of phenotypic plasticity because they are related to variable phenotypic expression across environments, as predicted in the pleiotropic model. Provided with a suitable measure of phenotypic plasticity as a quantitative trait, QTL detection tools can therefore be used to map plasticity genes. Measures have been proposed to ‘quantify’ the intensity of phenotypic plasticity or to ‘qualify’ the direction of phenotypic plasticity according to the shapes of reaction norms. In the early literature, however, measures of canalization and plasticity had been attributed to two alternative categories, static and dynamic stability, as defined by Lin et al. (1986) and Becker and Léon (1988). Static stability is analogous to canalization: a stable genotype shows a constant phenotype across environments. Dynamic stability refers to a genotype that responds the same as the average genotype of the population. Becker and Léon (1988) define three types of stability parameters. The first category relates to static stability. It includes the phenotypic variance across environments (S2), that is the variance of the genotype's phenotypic trait across environments. Therefore, the most stable genotypes, based on static stability, are those that show canalization—a constant phenotype across environments. The dynamic concept of stability is based on the regression slope, according to Finlay and Wilkinson (1963). This involves estimating the regression of genotypic values across environments on phenotypic mean values within environments and taking the slope of the regression as a stability parameter (Si). The highest stability is manifested as Si=1. Other measures are used to characterize dynamic stability (Becker and Léon, 1988). These include Wricke's (1962) ecovalence that is calculated for each genotype i as

where Rij is the observed quantitative trait response, mi and mj correspond to the mean values for genotype i and environment j, respectively, and m is the general mean. The greatest dynamic stability is achieved when W2=0. Another proposed measure relates to additive main effects and multiplicative interaction (AMMI) model decomposition (Zobel et al., 1988). This decomposition is based on principal component analysis of the G × E matrix resulting from analysis of variance with declared genotype and environment main effects. The various scores of each genotype for each principal component axis can be taken as dynamic stability parameters. Eberhart and Russell (1966) proposed the estimated mean square of genotype deviations from Finlay–Wilkinson's regression (sd2) as another stability parameter. This measure was later proposed as a third type of stability parameter by Lin et al. (1986), as it allows a choice of which type of stability (static or dynamic) is the most suitable to the data set.

Weber and Scheiner (1992) were the first to report mapping of plasticity genes. They were able to detect chromosome regions having quantitative effects on plasticity for thorax size in Drosophila melanogaster. Flies were raised in two environments and the difference in thorax size across environments was taken as a measure of plasticity. Wu (1998) mapped stability QTL and asked the question ‘do QTL displaying QTL × E interaction show patterns of crossover or scale effects?’ The implicit assumption was that all QTL showing QTL × E interaction would also affect phenotypic plasticity. Stratton (1998) proposed a different measure of plasticity—the linear and quadratic components of the reaction norm. Kleibenstein et al. (2002) used data from a population of recombinant inbred lines of Arabidopsis thaliana phenotyped in two environments, followed the example of Weber and Scheiner (1992), and used the difference between phenotypic trait scores in two environments as a direct measure of plasticity. Ungerer et al. (2003) extended this approach to data from more than two environments. Their measure of phenotypic plasticity was based on quantifying the difference in a phenotypic trait for each genotype across each pair of environments. Then, for each genotype, the trait differences between each pair of environments were summed and divided by the average of the trait across genotypes and environments. Kraakman et al. (2004) proposed a descriptor of phenotypic plasticity based on the slope of Finlay and Wilkinson's (1963) regression (Si for genotype i) and the deviation from that regression (si2) as defined by Eberhart and Russell (1966). Recently, Emebiri and Moody (2006) proposed using the ecovalence of Wricke (1962) and the AMMI first components genotypic scores, Si, and si2, as defined previously.

Multiple descriptors of phenotypic plasticity are, therefore, considered in the literature. The static stability parameters tend to display higher repeatability and heritability than the dynamic measures of stability (Léon and Becker, 1988; Lin and Binns, 1991). Surprisingly, none of the reports involving QTL mapping of plasticity used static stability parameters as measures of plasticity. Moreover, no estimates of heritability of plasticity measures were provided.

It is important, from both theoretical and practical reasons, to understand the genetic basis of phenotypic plasticity. The relationship between QTL × E interaction and plasticity genes is not clear. In this study, we address the following questions regarding the genetic basis of plasticity:

  1. 1)

    What is the relationship between loci showing QTL × E interaction and those affecting plasticity?

  2. 2)

    Do various measures of plasticity identify the same QTL?

  3. 3)

    How do the various measures of plasticity relate to the different concepts of plasticity and stability?

To address these questions, we used two data sets based on different barley mapping populations phenotyped in multiple environments. For each data set, two measures of plasticity were computed. This allowed us, through QTL detection, to test the specificities of the two measures of plasticity according to different patterns of QTL × E interaction. We then followed analysis of different sets of simulated data with loci showing different patterns of QTL × E interactions.

Materials and methods

Data sets

The data sets used in this study were obtained from Graingenes (http://wheat.pw.usda.gov/GG2/maps.shtml) and consist of genotypic and multi-environment phenotypic data for two doubled haploid mapping populations developed by the US Barley Genome Mapping Project. Each population was derived from a cross between two inbred parental lines in the following: The Harrington × TR306 (HT) population (n=145 doubled haploids) was the subject of prior QTL mapping reports (Hayes et al., 1993; Tinker and Mather, 1994; Tinker et al., 1996; Mather et al., 1997; Spaner et al., 1998). The Steptoe × Morex (SM) population (n=150 doubled haploids), first reported by Kleinhofs et al. (1993), was also the subject of prior QTL reports (Hayes et al., 1993; Hayes and Iyamabo, 1994).

Genetic maps

The marker data were used for linkage map construction as described by Mester et al. (2003) using MultiPoint software (http://www.multiqtl.com). Owing to the relatively small population sizes (150 per population), markers in high-density regions were removed to establish higher certainty of map order. The final maps for each population have an average marker interval of 10 cM. Marker orders in the newly generated maps are consistent with those available at GrainGenes. However, some genomic regions in the HT population are sparsely populated, leading to gaps of 30 cM or more. As a result, the linkage groups corresponding to chromosomes 2H and 5H were each split into two parts.

Quantitative traits

For SM, the traits selected for QTL analyses were grain yield and grain protein content (GPC). These traits were measured in seven environments throughout the United States. Grain yield and thousand kernel weight (TKW) data from 22 environments across the United States were used for HT. The SM and HT data trace to different locations and years.

QTL × E interactions

The quantitative trait data from different environments were analyzed simultaneously within the multi-environment framework of MultiQTL (http://www.multiqtl.com) using the fixed effect model proposed by Jansen et al. (1995) as follows:

Consider a simplified situation when the trait of interest (x) depends on a single QTL with two alleles q1/q2. For an arbitrary genotype j of the doubled haploids population, the trait measurement in the ith environment can be presented as

where μi is the mean trait value for ith environment (i=1, 2,…, I), gk is either +1 (k=1 corresponds to q1q1 genotypes) or –1 (k=2 corresponds to q2q2 genotypes), ai is the effect of allele substitution at putative QTL on trait in the ith environment and eij is a random normally distributed variable with zero mean and variance σi2. Assume that the putative QT locus q1/q2 resides in an interval flanked by two marker loci, M1/m1 and M2/m2, with recombination rates r1 and r2 in intervals M1/m1q1/q2 and q1/q2M2/m2. If we assume no interference (recombination in one segment does not affect the chances of recombination at an adjacent segment), then r=r1+r2–2r1r2, where r is the recombination rate between M1/m1 and M2/m2. On the basis of marker scores and measurements of the quantitative trait of interest for individuals from the mapping population, we should test whether or not variation of x indeed depends on the interval M1/m1M2/m2, and, if so, identify the corresponding locus q1/q2. The expected distribution of the trait in each of the four marker groups, U m 1 m 2 (x)= U 1 (x), U M 1 m 2 (x)= U 2 (x), U m 1 M 2 (x)= U 3 (x) and U M 1 M 2 (x)= U 4 (x), can be written as

with

The specification of the densities f q 1 q 1 (x) and f q 2 q 2 (x) depends on the assumptions made about the genetic control of the trait (in our case Equation (1)) and the residual variation (usually, normal distribution is used, although other possibilities also exist (Hackett and Weller, 1995; Kruglyak and Lander, 1995)).

Assuming that locus q1/q2 belongs to interval M1/m1M2/m2, the log likelihood for a sample of measurements xijs in marker groups with sizes Js (s=1,…, 4) can be written as

where vector ϑ=ϑ1={r1,μ,a,σ2} of unknown parameters specifies the putative QTL position in the considered interval, population means (μ1, μ2,…, μn), QTL effects (a1, a2,…, aI) and residual variances {σ12, σ22,…, σI2} in each environment. Clearly, vector ϑ1 defines the hypothesis H1e that there is a fluctuating QTL effect across environments (QTL × E interaction) associated with the considered marker interval. The assumption of a QTL with no such interaction can be parameterized by a vector ϑ=ϑ10={r1,μ,a*,σ2}, where a=a* fits the condition a1=a2=…=aI. For each pair of flanking markers, maximization of the likelihood function is conducted using sliding scanning across ‘trial’ positions rtrial of the putative QTL between the flanks, resulting in estimates of parameters μi, ai and σi2 for each rtrial. The likelihood for the H10 hypothesis is calculated similarly, but the QTL effect is considered constant across environments, ai=constant. To reduce the residual variation that is not accounted for by the mapping model, the QTL analysis was conducted using the multiple interval mapping approach (Kao et al., 1999).

The existence of a QTL was confirmed by permutation tests (Doerge and Chirchill, 1996). To detect a QTL showing QTL × E interaction, the hypothesis of constant substitution effect of the QTL across environments was compared with the more general hypothesis of varying effect across environments. For this comparison, interval QTL analysis allowing for QTL × E interaction (H1e model) was performed, followed by fitting a model assuming no QTL × E interaction (H10 model). We then simulated 1000 data sets using the parameter values estimated from the H10 model. For each of these sets, a QTL analysis was performed with and without an assumption of QTL × E. The difference in LOD score between the two models was used as a criterion for the significance of QTL × E interaction; that is we counted the proportion of cases where the difference in LOD score between the two models using simulated data was higher than the difference obtained using real data. This proportion gives the significance of QTL × E interaction. To take into account multiple tests, the obtained significance values were corrected based on false discovery rate control (Benjamini and Hochberg, 1995).

The plasticity measures

To quantify phenotypic plasticity, two measures were computed that relate (respectively) to the concepts of static and dynamic plasticity: (a) phenotypic variance across environments (VAR) and (b) the slope of regression based on Finlay and Wilkinson's regression (SLOPE). For each quantitative trait (for example, yield), we calculated SLOPE and VAR for each doubled haploid genotype in each mapping population. We then analyzed jointly the two VAR or two SLOPE measures for each set of traits (for example, yield VAR and VAR–GPC in the SM population) as implemented in MultiQTL (two-trait model).

Repeatability of phenotypic plasticity measures

To evaluate the repeatability of the plasticity measures, re-sampling of multi-environment data was performed on the HT data set. The HT data were chosen because it included the most environments. A total of 40 sub-samples, each including 4 of 22 environments, were generated at random using jackknifing. The sub-samples were generated separately for yield and TKW. For each re-sampling run, VAR and SLOPE were computed. To assess the repeatability of each of these scores, a fixed effects analysis of variance was performed. The repeatability was then computed as the ratio of the mean square of genotypic effects to the mean square for total variation.

Simulations

To understand the specificity of each plasticity measure relative to plasticity QTL, we performed simulations. Each simulation was based on the same genetic map, generated with MultiQTL, corresponding to a population of 300 doubled haploids. Only one QTL was simulated, with specific patterns of QTL × E interaction in each simulation (Figures 1 and 2). The simulations reflected the type of QTL × E interaction: crossover (CO) or scale (SC) interactions versus symmetry or dissymmetry of QTL effects across environments. For each simulated data set of two environments, we computed VAR and SLOPE. Using MultiQTL, we then conducted separate QTL analysis for the two plasticity measures.

Figure 1
figure 1

Different cases of QTL × E interactions and corresponding phenotypic variance across environments. QT refers to the quantitative trait. SC and CO represent scale and crossover patterns of QTL × E interaction effects, respectively. Sym and disym correspond to symmetrical and dissymmetrical allelic effects at the given QTL position. All1 and All2 represent the reaction norms of individuals carrying one or the alternative allele at the given QTL. V1 and V2 are the phenotypic variances across environments associated with each of the alternative alleles.

Figure 2
figure 2

Different cases of QTL × E interactions and corresponding slope of Finlay–Wilkinson's regression. QT refers to the quantitative trait. SC and CO represent scale and crossover patterns of QTL × E interaction effects, respectively. Sym and disym correspond to symmetrical and dissymmetrical allelic effects at the given QTL position. All1 and All2 represent the reaction norms of individuals carrying one or the alternative allele at the given QTL. m1 and m2 are environmental phenotypic mean values. b11 and b12 are phenotypic values associated with allele 1 in environments 1 and 2, respectively. SLOPE1 is the slope of Finlay–Wilkinson's regression for individual carrying allele 1.

Results

Real data analysis

Reaction norms profile

As an example of the shapes of the reaction norms for the quantitative traits, TKW data from the HT population are shown in Figure 3.

Figure 3
figure 3

Reaction norms of the 145 doubled haploids of Harrington × TR306 population for quantitative trait thousand kernel weight (TKW).

Trait correlations

In SM population, yield and GPC were negatively correlated (Pearson's correlation) in six of seven environments. The average correlation was −0.25; individual environment values ranged from −0.44 to 0.04. Yield and TKW, in the HT data set, were positively correlated in all but two environments, with an average of 0.33 and a range of −0.06 to 0.74. In light of our findings that trait QTL and plasticity QTL tend to coincide (see below), it was of primary interest to test whether traits and plasticity measures also tend to correlate. We therefore calculated the correlation between the plasticity scores and the traits scores in each environment by regressing the 150 (SM) or 145 (HT) plasticity values for the different genotypes with the individual trait scores. Overall, the plasticity traits (VAR and SLOPE) were not highly correlated with their corresponding individual quantitative environments in HT (Table 1). However, for SM, plasticity measures for yield were highly correlated with one environment of yield data (Table 2, Environment MonI91, R=0.83 for yield VAR and R=0.89 for yield SLOPE). The doubled haploid lines with higher yields in this environment also had higher plasticity measures for this trait. Overall, the correlations were higher for SM than for HT and higher for SLOPE than for VAR.

Table 1 Correlations between plasticity measures and individual quantitative traits typed in the several environments for HT dataset
Table 2 Correlations between plasticity measures and individual quantitative traits typed in the several environments for SM dataset

QTL main effects and QTL × E interactions

Twelve and 11 QTL were detected for TKW and yield, respectively, in the HT population (Figure 4, Table 3). Considering the linkage map positions of these 23 QTL, 14 distinct QTL peaks were defined; and of these, 12 QTL were coincident for the two traits. Of the 14 QTL, all but 4 showed significant QTL × E interaction and 9 of these showed QTL × E interactions for both traits.

Figure 4
figure 4

QTL analysis for Harrington × TR306 data set for yield, TKW and plasticity trait (phenotypic variance across environments). TKW indicates thousand kernel weight; dashed, QTL showing significant QTL × E interaction; SLOPE, slope of Finlay–Wilkinson's regression; LG, linkage group.

Table 3 QTL detected in HT dataset

A total of 16 QTL were detected using the SM data set: 8 for each trait, yield and GPC (Figure 5, Table 4). There were nine clearly defined QTL positions for at least one of the two traits; at eight of these positions, there were coincident QTL effects for the two traits. Two of the GPC and four of the yield QTL did not show QTL × E interaction.

Figure 5
figure 5

QTL analysis for Steptoe × Morex data set for yield, GPC and plasticity trait. GPC, grain protein content; dashed, QTL showing significant QTL × E interaction; VAR, phenotypic variance across environments; SLOPE, slope of Finlay–Wilkinson's regression.

Table 4 List of QTL detected in SM dataset

Plasticity and environmental canalization

In the HT data set, there were two QTLs for SLOPE and all were coincident with yield and/or TKW QTL. No significant QTL were found for VAR (Figure 4, Table 5). In all cases, the same parent was the source of the higher plasticity alleles for both traits. In the SM population, seven QTL were detected, which had significant effects on VAR of the both traits (yield and GPC) analyzed simultaneously (Figure 5). All seven QTL were colocalized with main effect QTL for one or both traits. The four QTL detected for SLOPE were also coincident with main trait QTL and with QTL for the VAR measure of plasticity. Figure 6 illustrates the difference in the slope of reaction norms for alternative alleles of the plasticity QTL located on chromosome 3H. At four of the six plasticity QTL, the same parent contributed to higher plasticity (as measured by VAR and SLOPE) for yield and the lower plasticity allele for GPC.

Table 5 List of plasticity QTL detected in both datasets
Figure 6
figure 6

Example of a yield QTL in SM data set affecting the slope of Finlay–Wilkinson's regression. The two trend lines show the mean difference in the slopes of reaction norms of individuals carrying alternative alleles at the plasticity QTL located on chromosome 3H. SM, Steptoe × Morex.

Repeatability of the plasticity measures obtained by resampling

The repeatabilities for the plasticity of TKW in the HT data set were 0.07 for VAR and 0.15 for SLOPE. For yield, the repeatabilities were lower, that is 0.05 and 0.06 for VAR and SLOPE, respectively.

Simulations

For the four main simulated situations (Figures 1 and 2), the results are specific to the type of QTL × E interaction pattern: for crossover QTL × E interaction and symmetrical QTL effects, no significant VAR QTL were detected. On the contrary, as one could expect from Figure 2, there were significant SLOPE QTL (LOD=16.9). For the crossover and dissymmetry cases, significant QTL were detected for both VAR and SLOPE with LOD scores of 3.4 and 35.3, respectively. For the scale QTL × E interaction situation, the results were consistent with the previous cases. VAR only detected QTL with dissymmetrical effects (LOD=2.9), whereas SLOPE detected QTL at LOD scores of 2.98 for symmetrical effect and LOD=15.2 for the dissymmetrical case. Therefore, VAR and SLOPE reveal different types of plasticity QTL according to their QTL × E interaction patterns. Overall, SLOPE revealed plasticity QTL in more situations than VAR.

Discussion

With a large number of environments, mixed model formulation, with mean QTL effects per interval taken as fixed effects and QTL × E interactions as random effects, is widely considered to be the most suitable approach for QTL analysis (Piepho, 2000). However, we considered environments to be a fixed factor, which is more appropriate when environments have been purposefully selected and when the analysis focuses on the detection of QTL × E interaction. Fixed model formulation ignores genetic correlations between environments: even if the major genetic effects are accounted for in the fixed effect model, some genetic correlation across environments will remain (Korol et al., 1998). Such correlations cause ‘residual’ QTL and QTL × E effects that have not been accounted for by the fixed part of the mapping model to reduce the power of the QTL detection tests. This is one of the main reasons for considering environment-specific effects as random effects (Piepho, 2000).

In a more suitable framework for QTL × E analysis, that is, mixed model formulation, QTL × E interaction is considered in terms of random effects. Instead of estimating and testing variation of the QTL additive effects ai across environments, it is the difference in variance σ2(ai) associated with each QTL allele across environments that is estimated and tested (Wang et al., 1999). With such an approach, the predicted environment-specific effects ai should be estimated using the best linear unbiased prediction (BLUP) method. However, the corresponding tests for detecting QTL × E interaction effects show low detection power (Wang et al., 1999; Xing et al., 2002). Thus, if the focus of the study is to detect QTL × E interaction, then despite the drawbacks noted in recent studies, fixed model (with fixed QTL effects and fixed QTL × E interactions) can be applied for detecting QTL × E interaction with higher power than mixed models (Jansen et al., 1995; Tinker and Mather, 1995; Romagosa et al., 1996; Korol et al., 1998; Wang et al., 1999). One of the possibilities for improving the suitability of fixed models for multiple environment situations is to reduce the residual genetic variation that is not accounted for by the fixed part of the mapping model. This can be achieved by using multiple interval mapping analysis (Kao et al., 1999). Our previous analysis (Korol et al., 1998) based on the application of the Eberhart and Russell (1966) approach for modeling QTL-by-environment interaction shows that the combination of multiple environment and multiple QTL analyses is indeed very promising for addressing this issue (see also Piepho, 2000). This combination is available in MultiQTL and we used it for the analyses in this study.

Comparison of QTL across quantitative traits and mapping populations

There were coincident yield and TKW QTL at 12 of 14 locations in HT and 8 of 9 QTL for yield and GPC were coincident in SM. This coincidence of QTL supports the need for measuring as many potentially related phenotypes as possible in QTL mapping experiments. The relationship of yield and the yield component TKW is obvious because grain yield is the result of developmental processes affecting kernel number from sowing to flowering and kernel weight after flowering (Slafer, 2003). All QTL affecting TKW should therefore affect yield. Yield and GPC are usually negatively correlated, as GPC is directly related to nitrogen sink/source relationship. This relationship, in turn, is related to the number of grains per spike formed and the deposition of carbohydrate during grain filling (Martre et al., 2003).

Seventy nine percent (11/14) of QTL significantly affecting quantitative traits in the HT population, and 55% (5/9) in the SM population, showed QTL × E interaction. This is not surprising, given the large number of environments sampled and the large differences between environments in terms of available moisture, growing conditions, and so on. It is quite improbable that a QTL affecting yield will have an expression independent of the environment. Indeed, the elaboration of yield is dependent on the number of spikelets per spike and on post-fertilization grain filling processes. These mechanisms are known to be dependent on environmental factors (Mozafar and Oertli, 1990; Loss and Siddique, 1994; Voltas et al., 1999; Savin et al., 1996; Savin et al., 1997; Schelling et al., 2003; Samarah, 2005). Therefore, the expression of QTL involved in the determination of yield will be highly dependent on environmental factors; and hence, the additive effect at the QTL will vary across environments. Ten of the coincident grain yield and TKW main effects QTL in the HT data set also showed significant QTL × E interaction, as did nine of the coincident grain yield/GPC QTL in SM.

The relationship of QTL × E interaction and phenotypic plasticity

In our analysis, QTL affecting environmental plasticity were coincident with trait QTL. This finding lends support to the arguments of Via (1993) and de Jong (1995), who reasoned that phenotypic plasticity may evolve as a consequence of selection for plasticity in underlying traits. Via et al. (1995) proposed two genetic models for phenotypic plasticity, and we have adopted their simplified terminology in this study. In the allelic sensitivity model (which is equivalent to the aforementioned pleiotropic model), constitutive genes are directly sensitive to the environment, and it is their differential expression across environments that is responsible for plasticity. The gene regulation model (equivalent to the epistatic model) hypothesizes that regulatory genes, sensitive to the environment, mediate the expression of constitutive genes that actually determine the trait. Our results support the allelic sensitivity rather than the gene regulation model (no cases favoring the gene regulation model were observed in our analysis). Therefore, it seems that a plasticity QTL should simultaneously affect the trait itself. However, in most prior reports, the authors conclude that some QTL affect only plasticity. For example, Wu (1998) reported that 13 out of 15 QTL for plasticity showed clear indications of gene regulation. Ungerer et al. (2003) found that 55–83% of the plasticity QTL colocalized with QTL having an effect on the respective quantitative traits. In the study by Kraakman et al. (2004), three of five QTL for plasticity co-located with QTL having an effect on the quantitative trait itself; thus, two of the QTL for plasticity were good candidates for gene regulation. Overall, reports in the literature indicate a greater prevalence of gene regulation than we observed. The main limitation of the QTL detection to identify QTL affecting only plasticity resides in the uncertainty of QTL location. In the literature, the use of colocalization for differentiating between gene regulation vs allelic sensitivity is considered a weak line of evidence. Therefore, conclusions regarding the relative importance of the two phenomena can, at this point, only be of a very limited value.

The main difference between prior reports and our report is the number of environments sampled. In our case, we have many more environments, particularly with HT. Wu (1998) sampled only 2 different environments and all 17 of the QTL detected showed gene regulation. Our findings are more similar to those of Kleibenstein et al. (2002), where only two environments were considered and yet all of the plasticity genes showed allelic sensitivity. Kraakman et al. (2004) also found a low rate of genes showing gene regulation: 2 QTL out of 18 showed gene regulation in an experiment with 15 environments. Therefore, on the basis of the literature and our results, we can hypothesize that gene regulation is more detectable as the number of environments decreases. Moreover, we estimated commonly used phenotypic plasticity parameters for each individual and then mapped QTLs for these parameters. In doing so, we introduced two different sources of noise. First, the estimation of phenotypic plasticity parameters has a sampling error (estimation of yield in each environment is, for instance, dependent on the harvesting conditions). Second, QTL mapping with these estimated parameters has sampling errors (position and QTL effects). However, by increasing the number of environments, we also increase the precision of the estimates of the slope of Finlay–Wilkinson's regression, which would be closer to the true value that one could obtain by sampling all ‘possible’ types of environments.

Alternative plasticity measures: results of simulations

Our results indicate that SLOPE is more repeatable than VAR. This may be due to the SLOPE measure being less sensitive to changes in environmental composition and its computational basis (regression). VAR and SLOPE are able to discriminate between different types of QTL affecting plasticity as is evidenced by the fact that in the simulations, the two plasticity measures reveal different plasticity QTL. The distinction between the two measures corresponds to the kind of simulated QTL × E interaction. Each of the measures relates to one of the two concepts of plasticity: one static (VAR) and the other dynamic (SLOPE). A ‘VAR’ plasticity QTL affects, in a quantitative manner, the range of variation of the reaction norms and relates to the definition of static plasticity. Plasticity QTL detected by SLOPE affects phenotypic plasticity ‘qualitatively,’ as the two QTL alleles show alternative shapes of reaction norms (see Figure 6 for illustration). As this measure of plasticity depends on the environmental mean of the genotype set, it relates to the dynamic concept of plasticity. Our analysis suggests that there may not be a universal measure of plasticity that will identify all plasticity QTL. Rather, there are different types of plasticity QTL related to the different concepts of plasticity and the two measures will show different degrees of repeatability, with SLOPE being more ‘heritable.’

Kumar et al. (1998) reached the same conclusion, but these authors, like others (Annicchiarico, 1997; Ortiz et al., 2001), reported higher repeatability values for the two measures of plasticity than we did. This may be due to higher heterogeneity among the environments in our data sets. The number of plasticity QTL varied with data set. More plasticity QTL were detected in SM than in HT. This may be due to the smaller number of environments in the former (7 vs 22). It is likely that more plasticity QTL will be detected when a smaller number of more homogeneous environments are sampled.

The importance of homogeneity vs heterogeneity of environments

As described in the previous section, the repeatability values of the plasticity measures are very low. The dependence of the results on the set of environments is crucial, which is intuitively obvious. Indeed, the regulatory pathways affecting the traits chosen for study may vary from one environment to another due to unique features of each environment. Therefore, with an increasing number of environments, one will eventually reach a situation where it is impossible that a single ‘consistent’ regulatory pathway determines the pattern of plasticity across the entire set of environments. On the contrary, when considering smaller and/or more homogeneous samples of environments, it is more likely that some common regulatory pathways may be involved in the control of phenotypic plasticity. Our results suggest that the regulatory genes controlling plasticity are different from one subgroup of environments to another. Thus, for a better understanding of the genetic control of plasticity, one should take into account the homogeneity of the environments that are sampled.

Trait correlations and epistasis

In the HT data set, VAR and SLOPE were not highly correlated with the corresponding individual quantitative traits (for example, yield VAR) in the different environments. However, in SM, plasticity measures for yield were correlated with yield per se in one environment (MonI91). This suggests that this environment is different from the other environments in that it has a deep influence on the plasticity traits. The trends for correlations between plasticity traits and the quantitative traits themselves were similar for yield and GPC in the SM data set and for yield and KTW in the HT data set. This suggests that QTL detection results should have been similar across quantitative traits for each data set. In SM, however, the alleles at four of the seven QTL giving higher plasticity for yield gave lower plasticity for GPC. In HT, on the other hand, alleles giving higher plasticity for yield also gave higher plasticity for TKW. We believe that this difference is attributable to the physiological bases of these traits and that it reflects the phenotypic correlations—positive for yield and TKW and negative for yield and GPC.

It is therefore essential to take into account the correlation between related quantitative traits to understand the evolution of plasticity, as suggested by Schlichting (1986). In the epistatic model (Scheiner and Lyman, 1989), plasticity genes interact epistatically with, and thus regulate, the ‘constitutive’ loci that determine the mean value of the trait. In a limited number of models—Drosophila melanogaster (Scharloo, 1991), tomato (Eshed and Zamir, 1996) and RNA viruses (Burch and Chao, 2004)—there is evidence that canalization and phenotypic plasticity may be explained by epistatic mechanisms. The role of epistasis in the elaboration of phenotypic plasticity should be a subject of future studies.

Evolutionary perspectives

Phenotypic plasticity has long been recognized by evolutionary biologists and ecologists as an important phenomenon. When the environments with which a genotype is confronted are heterogeneous, it is not likely that having the same phenotype across the spectrum of environments will lead to the highest fitness. Therefore, the capacity to modulate the phenotype in response to the environment should be evolutionarily advantageous. One way an individual can adapt to environmental change is through a phenotypically plastic response (Scheiner and Lyman, 1989). Clearly, numerous examples could be cited where phenotypic stability may be of high value.

By extending the definition of plasticity to the level of gene expression, it is now possible to test for variability of gene expression across environments. In principle, plasticity/canalization at the gene expression level could be considered a molecular manifestation of plasticity/canalization at the level of quantitative traits. The doubled haploid barley mapping populations and reference data used in this study, together with tools for transcriptome profiling (Druka et al., 2006), would provide new insights into the genetic control of genetic and environmental plasticity.