MATERIALS AND METHODS

Plant material: The RIL set derived from a cross between Cvi and Ler accessions was used for these studies (Alonso-Blancoet al. 1998b). Seeds of 162 RILs and the parents (CS22000), Col-1 (CS3176), and the Col-1 er-2 (CS3401) mutation were obtained from the Arabidopsis Biological Resource Center (ABRC) in Columbus, Ohio (http://www.arabidopsis.org) and used directly for hypocotyl measurements. Lan-1 (La ERECTA) was obtained from Carlos Alonso-Blanco. F₁ hybrids were made by reciprocal crosses for Ler × Cvi F₁ (Ler female) or Cvi × Ler F₁ (Cvi female).

Growth conditions: Seeds were sterilized in 1.5-ml microcentrifuge tubes for 10 min in 70% ethanol, 0.01% Triton X-100, followed by a 10-min wash with 95% ethanol, and then resuspended in 1 ml sterile water. After imbibition overnight at 4° in the dark, seeds were placed individually onto 0.7% phytagar plates containing ½ Murashige and Skoog salts using a Pipetman. Seedlings were spaced at a uniform density so that they did not shade each other. Plates were kept at 4° in the dark for another 3 days, followed by 4 hr of 120 μE m⁻² sec⁻¹ white light to induce germination. Further incubation was at 23°. Preliminary experiments in six conditions (all except dark) with Cvi, Ler, and most of the CvL RILs revealed substantial variation in hypocotyl lengths among lines and slight variation from week to week and from plate to plate (data not shown). For the results reported here, all light environments and all RILs were done in the same week to minimize week-to-week and week-to-plate variation. Furthermore the number of RILs per plate was increased to 12, providing better statistical control of plate-to-plate variation. Plates were rotated within each incubation chamber every 12 hr for the duration of the experiment to reduce variation among plates within each incubator. Ideally, the entire experiment would be replicated several times over to reduce the contribution of uncontrolled variation in the observed differences between RILs and light environments and to provide precise estimates of the magnitudes of the components of variation. Without such replication, differences due to uncontrolled variation between growth conditions in different incubators could be attributed to light environments and lead to spurious genotype-by-light environment associations. However, a single run of 162 RILs under seven light environments requires ~240 person-hours to perform. From the preliminary studies we believed that extraneous variation could be controlled sufficiently to obtain useful results from a single-week experiment alone. In all, 15–30 seedlings of each of 162 CvL RILs, Cvi, Ler, reciprocal F₁ hybrids, and photoreceptor mutants were arrayed in groups of 12 lines per plate across 15 plates. This was replicated for the seven environmental conditions.

Light/hormone conditions: Incubators used for all environments were Percival model E30B (Percival Scientific, Boone, IA). One incubator (Percival E30LED) was equipped with LED lights and used for the far-red environment. Neutral density screens were used to vary light fluence rate. Light measurements were made with a LI-1800 instrument (Li-Cor, Lincoln, Nebraska). We wanted to identify a fluence rate that would maximize the subtle variation in light sensitivity seen in natural populations. Pilot experiments showed that at high light fluence rates CvL RILs had relatively uniform, short hypocotyls and at low light fluence rates CvL lines were much longer but more variable. We chose intermediate light fluence rates, from a fluence response curve, for each light condition, where the broad-sense heritability was maximized for subsequent experiments. White light was provided by three 35-W cool white fluorescent bulbs and two 25-W incandescent bulbs. The photosynthetic active radiation (PAR, 400–700 nm) was 35 μE m⁻² sec⁻¹, the Pfr/P ratio was 0.72 (Kendrick and Kronenberg 1994, p. 268), and the R/FR ratio (655–665 nm)/(725–735 nm) was 1.3. The same light conditions were replicated in another incubator for the GA environment except that 30 μm GA₃ (Sigma, St. Louis) was added to the medium. Blue light (PAR = 4 μE m⁻² sec⁻¹) was provided by three 20-W cool-white fluorescent bulbs and a filter that blocked light above 550 nm. Red light (PAR = 35 μE m⁻² sec⁻¹) was provided by three 20-W Gro-Lux fluorescent bulbs (Osram Sylvania, Danvers, MA) and a red filter that blocked light below 600 nm. Far-red light (0.5 μE m⁻² sec⁻¹; 700–730 nm) was provided by LED lights. The same incubator was used for the dark and BRZ environments, and plates were wrapped in aluminum foil and received no further light after the 4-hr germination light pulse. A dose response curve, using different concentrations of the brassinosteroid biosynthetic inhibitor BRZ 91, identified 0.75 μm as the optimum concentration to maximize the heritability. A total of 0.75 μm BRZ (synthesized at RIKEN) was used unless otherwise indicated.

Hypocotyl length measurements: On day 2, poor germination was scored in white light as 1 or 0. Most lines had already germinated (and were scored as “0”), but 14 lines (CvL nos. 1, 3, 8, 15, 16, 22, 24, 27, 38, 39, 152, 185, 186, and 188) had not (and were scored as “1”). All lines except CvL 3 germinated by day 3 and were measured on day 7. The germination state seen in white light was representative of all conditions and likely reflected both environmental and genetic variation in the state of the seeds rather than light response. Therefore, germination was used as a covariate in subsequent analyses. Seedlings were transferred to acetate sheets containing moist tissue paper and scanned on a flat bed scanner. Hypocotyl lengths were measured in millimeters using National Institutes of Health (NIH) image version 1.62 (http://rsb.info.nih.gov/nih-image). The effect of the covariate germination on hypocotyl length and the average number of seedlings measured in each environment are shown in Table 1.

Statistical analysis: Data analysis was performed using the statistical package R (Ihaka and Gentleman 1996; http://www.R-project.org/). Hypocotyl length data are approximately normally distributed, so no transformation was needed. Hypocotyl length was fit using a statistical model that is described in detail below. Briefly, hypocotyl lengths were fit by a mixed-effects linear model with terms for germination, plate, and RIL. RIL and plate were modeled as random effects with RIL nested under plate; germination status was modeled as a fixed effect. Data for each light environment were fit separately. The variation due to plate, RIL, and residual variation is shown in Table 2. The relative variation explained by line means is an estimate of the broad-sense heritability. Best linear unbiased predictors (BLUPs) of RIL means under this model were used for QTL mapping; those lines showing poor germination had their means augmented by the coefficient for germination state. Note that these BLUPs of RIL means implicitly omit the plate effects, and this reduces the uncontrolled variation in the trait for QTL mapping. Calculations used the linear mixed effects (lme) function in the “nlme” package to R (Pinheiro and Bates 2000). The genetic correlation between environments (r_GE) was calculated using cov₁₂/(σ_L1σ_L2), where cov₁₂ is the covariance in line means, corrected for germination and plate effect, and σ_L1 and σ_L2 are square roots of among-line variances from the linear mixed-effects model (Robertson 1959). Here BLUPs of RIL means used a model with RIL effects as fixed effects; this leads to unbiased estimates of cov₁₂ when the reduction in degrees of freedom is accounted for The coefficient of genetic variation CV_G was calculated for each environment by dividing σ_L1 by the grand mean of line means and multiplying by 100. Confidence intervals for r_GE and CV_G were calculated, using a “leave one out” jackknife procedure on 161 lines (Efron and Tibshirani 1993).

The model for the dependence of hypocotyl length on germination status, RIL identity, and plate effects was yij=Xiβ+ZiΓRIL+WiΓplate+εij, (1) where y_ij are the measurements of hypocotyl length under a single light environment with i = 1, … , 161 indexing the 161 RILs and j = 1, … , n_i indexing the seedlings of the ith RIL. X_i = (1, germ_i), where germ_i = 1 if germination was poor and 0 otherwise. β = (β₀, β₁)′ is a column vector of two unknown coefficients; β₀ is the mean under good germination and β₁ is the increment due to poor germination. Z_i is a row vector of 161 elements all of which are zero except for the ith, which is 1. Γ_RIL ~ N(0, δRIL2I161×161 ) is a column vector of 161 elements. W_i is a row vector of 14 elements, all of which are zero except for the kth, which is 1 when RIL i was incubated on plate k such that ∑k(∑iZi′Wi)k=1161×14 (i.e., lines “nest” within plates). Γ_plate ~ N(0, δplate2I14×14 ) is a column vector of 14 elements and ε_{i j} ~ N(0, σ²) is the residual and independent of other model terms. The phenotypic means are taken as τ_i = β₀ + Z_iΓ_RIL and estimates are obtained by replacing the respective parameters with restricted maximum-likelihood (REML) estimates. These estimators are best linear unbiased, so phenotypic means based on τ_i are the BLUPs under this model.

QTL analysis: The CvL RILs had been previously genotyped (Alonso-Blancoet al. 1998b), using amplified fragment length polymorphisms (AFLP; Vos et al. 1995) and cleaved amplified polymorphic sequences (CAPS) markers (Konieczny and Ausubel 1993). Marker data and the genetic map were obtained from the web at the Nottingham Arabidopsis Stock Center (http://nasc.nott.ac.uk/). We used 163 of the 293 available markers that mapped to unique genetic loci and that had been genotyped on an average of 160 out of 162 RI lines. The BLUP data representing the line mean coefficients corrected for germination and plate effect were used as the phenotypic values for QTL mapping. The CIM (Zeng 1994) function of QTL Cartographer (http://statgen.ncsu.edu/qtlcart/cartographer.html) was used to map QTL. Background markers were chosen using the forward/backward stepwise multiple regression of SR map at a P value of 0.001. When SR map chose adjacent background markers in different environments, QTL models were tested where the same background marker was used in each environment to optimize the LOD score and minimize the LOD support interval. The numbers of background markers ranged from four to seven and are shown in Figure 3. Thresholds in each environment were set internally by running sets of 5000 permutations (Doerge and Churchill 1996). In each environment, a LOD score of 3.43–3.63 (depending on light environment) corresponded to an experiment-wise P value of 0.01 as determined by permutations. Instead of using seven different thresholds we used the largest one. A LOD of ~2.8 would correspond to P = 0.05. We used a window size of 1 cM because hypocotyl length has high heritability, many QTL were in tight linkage, and many lines and many markers were used in this experiment. QTL maps with larger window sizes (1–10 cM) gave broader QTL peaks; however, the two LOD support intervals were equivalent to the 1-cM window size map. Generally the width of QTL peak was defined by the flanking markers, at various window sizes.

Recombinant inbred line-by-environment testing: The RILs are nested within plates in each light environment, so there is no “RIL-by-light environment-by-replicate” term to use as the error term for the “RIL-by-light environment” interaction. To test this interaction, two approaches were used. One was the sequential F-test of the RIL-by-light environment interaction term in a model including terms for line, light environment, plate, germination, and the germination-by-light environment interaction. In the absence of spatial or other effects on plates that increase between-RIL variation on a plate without also increasing “within RIL” variation, this is a powerful and appropriate test. The other test refers Tukey's “1 d.f. for interaction” statistic to its distribution under permutation of RIL interaction terms within plates; this yields a correct P value even in the presence of uncontrolled variation within RIL on a plate, but generally has limited power. For the first approach, Equation 1 and definitions above are extended as yijk=X∼ikβ∼+ZiΓ∼RIL+Z∼ikΓRIL.light+W∼ikΓ∼plate+εijk, (2) where k indexes the light environment, X∼ik=Xi⊗Lk , Z∼ik=Zi⊗Lk , and W∼ik=Wi⊗Lk with L_k being a 1 × 7 vector with a 1 in the kth element and zeroes elsewhere. ε_ijk ~ N(0, σk2 ) are independent residuals. Other terms on the right-hand side are suitably sized vectors of coefficients following the normalizations ∑i(Γ∼RIL)i=0 , ∑i(Γ∼RIL⋅light)7(i−1)+k=0 , ∑k(Γ∼RIL⋅light)7(i−1)+k=0 , summing i over 1, … , 161 and k over 1, … , 7. The mean square for the RIL-by-light environment interaction, Z∼ikΓ∼RIL⋅light , has 876 d.f. after accounting for the other terms and the F-statistic takes the residual mean square to be the error. The Tukey 1 d.f. for interaction statistic (Scheffé 1959, section 4.8) decomposes the sum of the squared RIL-by-light environment terms into two parts, SSG=(∑i=17∑j=1161α^i(Γ∼̂RIL)j(Γ∼̂RIL⋅light)7(j−1)+i)2∑i=17α^i2∑j=1161(Γ∼̂RIL)jSSres=∑i=17∑j=1161(Γ∼̂RIL⋅light)7(j−1)+i2−SSG, where α^i=β∼̂2i−1−17∑i=17β∼̂2i−1. (3) The statistic is F = SS_G/(SS_res/d.f.), taking d.f. as 875. A permutation test can be constructed by permuting (Γ∼̂RIL⋅light)7(j−1)+i with respect to the index j and calculating F under each permutation. Possible plate effects should be preserved in the reference distribution, so permutations of j must respect the assignment of RILs to plates. This type of permutation honors the normalizations above.

Multienvironment QTL mapping:The multitrait CIM (mCIM) program JZmapqtl in QTL Cartographer was used (Jiang and Zeng 1995). mCIM mapping calculates a joint likelihood to detect QTL in multiple environments and a genotype-by-environment likelihood to determine if QTL are specific to certain environments. Light and hormone interactions were tested separately by including four light environments (white, blue, red, and far-red) in one analysis, white light and GA in a second analysis, and dark and BRZ in a third. A common set of background markers was used for each analysis (Figure 3, Table 4), to avoid problems of overparameterization. When SRmap chose adjacent background markers closer than 3 cM apart in different environments, a common marker was chosen that was selected in the majority of environments. An experiment-wise P = 0.01 threshold for both the joint likelihood (main effects) and G × E likelihood was determined separately for each analysis by performing 5000 permutations (Doerge and Churchill 1996; G × E LOD = 5.7 for four light environments and 3.6 for each of the hormone comparisons). These routines were provided by Chris Basten and are available upon request. The likelihood of the G × E test at each QTL was compared to the threshold to determine if that QTL showed a significant G × E interaction. This occurs when the estimated QTL effect is different from the joint effect in at least one of the tested environments.

Tests of epistasis: We tested interactions between QTL and then performed an exhaustive search for pairwise marker interactions using BQTL. Each environment was analyzed separately and included main effect QTL as background markers (Figure 3, Table 4) and the covariate germination. A total of 43,956 pairwise tests were done between 296 loci. These included 163 actual markers and 133 pseudomarkers, at marker intervals <2 cM, creating an ~2-cM walking speed. The test statistic is the LOD score difference between a model with only additive effects and one that included an epistatic term. Thresholds for statistical tests used a sequential permutation procedure (Nettleton and Doerge 2000) to ensure that enough permutations were performed to assert that each test attained (or failed at) P = 0.05. This is discussed in detail below. A total of 5000 permutations were done in the white light environment. Each permutation tested 43,956 pairs of markers. A LOD score of ~4.6 corresponded to an experiment-wise threshold of P = 0.05 and was similar across light environments. The effect of the epistatic interaction is shown as 4i (Table 4), which represents the difference between the homotypic and heterotypic means (Juengeret al. 2000). The interacting loci on chromosome 5 were also detected as the pair with the largest test statistic from the preliminary experiment in white light.

The maximum-likelihood method provides good power for detecting epistasis (Kao 2000) but requires more time for computation than linear methods for QTL mapping (Knappet al. 1990; Haley and Knott 1992). To obtain correct P values with a model with covariates and additional loci, a permutation procedure is needed—further increasing the computational burden. Hybrid procedures for QTL mapping that are linear with respect to some, but not all, loci (Jansen 1993; Zeng 1994) are widely used and provide some of the benefits of a full maximum-likelihood approach at a reduced computational cost. Such a hybrid approach and an associated permutation test for scanning for epistasis were implemented as follows. The log-likelihood used for scanning for epistasis is L(α,β,γ1,γ2,γ12,θ,σ2;y,x,l1,l2)=∑log∑j=12∑k=12ϕ(yi;μijk,σ2)πZ∣M(Zl1l2=(zj,zk)∣M=mi), (4) where μijk=α+xiβ+zjγ1+zkγ2+zjzkγ12+EZ∣M(zl3,…,zln∣mi)θ, (5) yi=τ^i , α is the intercept, x_i is zero if the ith line had good germination and one otherwise with β as the coefficient for germination, j and k index the parental lines of the alleles at loci l₁ and l₂ being tested for epistasis, z_j = 2(j − 1.5), z_k = 2(k − 1.5), γ₁ and γ₂ are the main effects at those alleles and γ₁₂ is the epistatic effect, E_Z|M(z_l3, … , z_ln|m_i) is the expectation of vector of z's corresponding to n − 2 other loci given the marker information for subject i and φ is a vector of coefficients, π_Z|M(Z = (z_j, z_k)|M = m_i) is the probability that the two loci are in states j and k given the marker information for subject i, σ² is the residual variance, and φ(y; μ, σ²) is the normal density function. The maximum-likelihood solution of (4) and (5) with respect to all of the coefficients is carried out. In addition, the maximum-likelihood solution under γ₁₂ = 0 (no epistasis) is found. The log-likelihood ratio statistic X²(l₁, l₂) = 2(sup_ν L(ν; y, x, l₁, l₂) − sup_ν0 L(ν₀; y, x, l₁, l₂)) is formed for all pairwise combinations of loci, l₁ = 1, … , 295, l₂ = l₁ + 1, … , 296, taking ν as the vector of free parameters and ν₀ as that vector with γ₁₂ fixed at zero. Statistical significance is ascertained via permutation testing using the “residual empirical threshold” method (Doerge and Churchill 1996). Predicted values and residuals are formed using a model in which γ₁ =γ₂ =γ₁₂ = 0; i.e., only the germination effect and the effects of the n − 2 loci used in all models are included. A new vector of trait values is formed by adding the fitted values to a permutation of the residuals from that model. The log-likelihood ratio statistic is found as above for every combination of loci, and the maximum of these is found for each permutation. Attained P values are found as the fraction of permuted maxima that equals or exceeds X²(l₁, l₂). This produces genome-wide P values that are nominally correct under the null hypothesis of no epistasis anywhere on the genome. However, the randomness in the procedure is considered objectionable especially when claiming to have attained a fixed significance level. This can be overcome by following the recommendations of Nettleton and Doerge (2000), requiring that the 95% confidence interval for P values exclude 0.05 and 0.01. Calculations were performed by the function “bqtl” available in the R package, BQTL (http://hacuna.ucsd.edu/bqtl and http://cran.r-project.org).

QTL effect estimation: QTL effects were estimated by applying the method of maximum likelihood to the QTL model (Kao 2000). This model included main effect QTL identified by CIM, significant epistatic loci, and the covariate germination. QTL effects are estimated using a likelihood analogous to (4) and (5). Since the number of loci is manageably small and there are only a few models to be fit, no use of linearized terms E_Z|M is needed, and full maximum-likelihood fitting is used. The modifications required are to replace μ_ijk by μij1,…,jk=α+xiβ+∑k=1Kzjkγk+∑l=1Lzjnzjslγrlγsl , where j₁, … , j_K epistatic index K loci included in the QTL model, terms are included with r_l and s_l indexing the main effects upon which they depend. Obvious modifications are made to the summation and to π_Z|M, the joint allele state probabilities, in (4). These calculations were also performed by the function bqtl available in the R package. The additive effect is shown as 2a, the difference between homozygous classes. The percentage of change caused by a single QTL is the effect in millimeters (2a) divided by the average RIL hypocotyl length for that environment (Table 4) multiplied by 100. The percentage of variance explained for each QTL was determined by squaring the coefficient (a) and by dividing the residual variance in a null model without genetic loci (σrN2 ). Total variance explained was determined as 1−(σrG2)∕(σrN2) , where (σrG2 ) is the residual variance in the model with all genetic terms.

Near isogenic lines: The LIGHT1 near isogenic line (NIL) was derived from line N42 created to map EDI (Alonso-Blancoet al. 1998a). N42 contains only 35 cM of Cvi from the top of chromosome 1 in a Ler background determined by selection against other markers throughout the rest of the genome (gift from Carlos Alonso-Blanco and Maarten Koornneef). N42 was crossed to Ler, and F₂ plants that had the Ler allele at EDI and were heterozygous at marker g2395 were selected. The LIGHT1 NIL is an F₃ line (F₃-77), derived by selfing, that is heterozygous for the AFLP marker GD143L-Col and the CAPS marker m235 at 22 and 34 cM on chromosome 1, respectively. After hypocotyl lengths in white light were measured, individual seedlings were genotyped at g2395 as a marker for the LIGHT1 QTL. The LIGHT2 and HYP2 NILs were made by crossing the RIL CvL 125 to Ler. The CvL125 × Ler F₂ cross segregates Cvi DNA from 34 to 63 cM on chromosome 2 containing the PHYB and ERECTA loci, as well as from 84 to 107 cM on chromosome 5. For LIGHT2, 100 F₂ plants were measured in white light and genotyped at PHYB and GPA1 as markers for the LIGHT2 QTL. Interval mapping was done between these two markers. For HYP2, F₂ plants were measured in the far-red and BRZ environments and genotyped at BAS1 (Neffet al. 1999). The additive and dominance effects of each marker were assessed using linear regression.

Genotyping:Genotyping was done using CAPS markers. GPA1, g2395, and m235 information was from TAIR (http://www.arabidopsis.org). PHYB oligonucleotide primers were 5′ CTGCTGACGAGAACACG 3′ and 5′ GAAAGTTGGCTTAAATGG 3′; Ler has a PstI restriction site absent from Cvi. BAS1 oligonucleotide primers were 5′ ATATAATAGGCGTTCATCTAATG 3′ and 5′ CTCGGAGTTCGTACATG 3′; Cvi has an AccI restriction site absent in Ler. The BAS1 marker is 170 kb from the ERECTA gene, on the same bacterial artificial chromosome (BAC) T9J22.

Data and statistical routines are available on our web page (http://naturalvariation.org).