Genetic architecture of plant stress resistance: multi‐trait genome‐wide association mapping

Summary Plants are exposed to combinations of various biotic and abiotic stresses, but stress responses are usually investigated for single stresses only. Here, we investigated the genetic architecture underlying plant responses to 11 single stresses and several of their combinations by phenotyping 350 Arabidopsis thaliana accessions. A set of 214 000 single nucleotide polymorphisms (SNPs) was screened for marker‐trait associations in genome‐wide association (GWA) analyses using tailored multi‐trait mixed models. Stress responses that share phytohormonal signaling pathways also share genetic architecture underlying these responses. After removing the effects of general robustness, for the 30 most significant SNPs, average quantitative trait locus (QTL) effect sizes were larger for dual stresses than for single stresses. Plants appear to deploy broad‐spectrum defensive mechanisms influencing multiple traits in response to combined stresses. Association analyses identified QTLs with contrasting and with similar responses to biotic vs abiotic stresses, and below‐ground vs above‐ground stresses. Our approach allowed for an unprecedented comprehensive genetic analysis of how plants deal with a wide spectrum of stress conditions.


Fig. S1
Narrow sense heritability for Arabidopsis thaliana resistance to abiotic and biotic stresses.            Genes in linkage with SNPs with -log 10 (P) score above 4 (20 kb half-window size) in the contrastspecific GWA mapping of parasitic plants and aphids on the one hand vs fungus, caterpillar, thrips and drought on the other hand 2 Table S5 Candidate genes in linkage with SNPs with -log 10 (P) score above 4 (20 kb half-window size) that have common effects on plant response to parasitic plants and aphids on the one hand vs fungus, caterpillar, thrips and drought on the other hand Table S6 Candidate genes in linkage with SNPs with -log 10 (P) score above 4 (20 kb half-window size) that have common effects on biotic and abiotic stress responses Methods S1 Salt.
Methods S11 Screening of T-DNA lines.
Methods S12 Simulations to compare power for full MTMM, contrast MTMM and univariate analysis. 3

Phenotyping the HapMap population of 350 Arabidopsis thaliana accessions
Methods S1 Salt.

Traits
Salt_1: root response to mild salt stress (75 mM NaCl), in terms of a combination of main root vector length (MRVL), number of lateral roots per main root (noLR) and straightness (main root length divided by MRVL) Salt_2: root response to severe salt stress (125 mM NaCl), in terms of MRVL Salt_3: root response to severe salt stress (125 mM NaCl), in terms of noLR Salt_4: root response to severe salt stress (125 mM NaCl), in terms of straightness.  Silique zero belongs to the flower that opened first on the day of the treatment. 2 Stress did not disappear when watering with PEG-containing nutrient solution stopped, because PEG is not evaporating.

Growing conditions
Four types of stress were studied in four different experiments. Seeds were sown in Petri dishes on wet filter paper. After 4 d of cold treatment, they were placed in the light at room temperature for 1.5 d to germinate.
Germinated seeds were placed on rockwool blocks saturated with nutrient solution (Hyponex,

Experimental design
In the salt experiment three blocks received control treatment and three blocks received treatment conditions.
For the drought experiments the Phenopsis phenotyping platform was used, preventing position related differences in plant growth within the climate chamber (Granier et al., 2006;Bac-Molenaar et al., 2015). Plants were grown under similar conditions during the first 3 wk. Drought stress was imposed by withholding water for 7 d while the rest of plants were watered every 2 d with 1 l of water per tray. Botrytis inoculation was carried out 24 h before Pieris inoculation. Plants were 4 wk old when they were exposed to stress by P. rapae as single or combined stress. Plants were inoculated with two newly hatched L1 and the larvae were allowed to 9 feed for 5 d until harvesting. Phenotypic measurements were taken 24 h and 5 d after inoculation with P. rapae as a single and combined stress. After 24 h, the number of leaves damaged and number of feeding sites was counted in plants exposed to P. rapae, drought and P. rapae, and Botrytis and P. rapae. After 5 d, fresh weight was measured for the five treatments.

Experimental design
Plants were screened in rounds of 37 accessions. Three control accessions were present in all rounds (Col-0, Tsu-0, Fei-0). In each round, six replicates were included per accession-treatment combination. Treatments were screened simultaneously. Plants were randomly allocated in trays (28 accessions per tray). Plant positions within a tray were recorded (X pos and Y pos ). The chamber where the experiment were conducted consist of six racks, and each rack contained four shelves. Positions of trays within shelves within racks were also recorded.

Genotypic means
For each of the three traits (shoot fresh weight, number of leaves damaged and number of feeding sites) we fitted the following mixed model: where TRT is the treatment factor (Control and four treatment levels), GEN is genotype (accession) and GEN:TRT is the genotype by environment interaction. The terms GEN, TRT and GEN:TRT were fitted as a fixed effect and all others as random. For each of the three traits, significance of each model term was assessed and only significant terms were retained for the estimations of genotypic means. Genotypic means (BLUEs) were calculated for shoot fresh weight, number of leaves damaged and number of feeding sites, for each accessiontreatment combination.

Definition of target traits
Target traits were defined based on the genotypic means for shoot fresh weight, number of leaves damaged

Growing conditions
Arabidopsis seeds were put on filter paper in the dark at 4 o C for 2 d. Then, Arabidopsis seeds were sown on river sand (with a thin layer of soil on the top of river sand). Arabidopsis plants were grown for 2 wk on river sand at 21°C, 60% RH, 100 µmol m -2 s -1 light intensity, 12 h : 12 h, light : dark photoperiod. After 2 wk, Arabidopsis seedlings were surface-sterilized with 70% ethanol for 5 s and washed with sterile demi-water. The rhizotron system was prepared by cutting a hole at the side of 14.5 cm diameter round Petri dish, putting successively a piece of round rock wool slice (14.5 cm diameter, 1.5 cm in thickness) at the bottom of Petri dish, a piece of 12 cm diameter glass-fibre filter discs and a piece of 14.5 cm diameter nylon mesh on top. The rhizotron system was supplied with sterile ½ Hoagland liquid medium. Sterile seedlings were then moved to prepared rhizotron system by fitting the plant in the hole of the Petri dish. Leaves and shoots of the seedlings were kept outside of Petri dishes. The roots were carefully separated and organized on the top of nylon mesh by forceps. Arabidopsis seedling were grown in rhizotron system at 21°C, 60% RH, 100 µmol m -2 s -1 light intensity, 12 h : 12 h, light : dark photoperiod for another 2 wk.
Sterile Phellipanche ramosa seeds were spread on 5 cm diameter glass-fiber filter discs (Whatman GF/A paper) which were wetted with 0.8 ml sterilized demi-water and placed in 9 cm diameter Petri dishes. The Petri dishes were sealed with parafilm and then kept in dark in a growth chamber at 20°C for a 12 d precondition period. Preconditioned seeds on a glass-fiber filter disc were dried and treated with 0.8 ml strigolactone analog GR24 at the concentration of 3.3 x 10 -3 µM for 1 d under dark at 25°C. GR24 treatment triggered the initial germination of P. ramosa. After 1 d, GR24 was immediately washed off the P. ramosa seeds by sterile demi-water.
Pre-germinated P. ramosa seeds were spread along 4-wk-old Arabidopsis seedlings in the rhizotron system with painting brushes. The rhizotron Petri dish were sealed with tape and covered by aluminium foil.
Plant were grown at the same condition for the following 4 wk. Pictures of P. ramosa-infested roots in the rhizotron system were taken 4 wk after infection with Canon camera EOS 60D DSLR (with EF-S 18-135 mm IS Lens).

Experimental design
The 359 accessions were screened in two rounds (the first 200 accessions, the second with 160 accessions, two accessions were used for control in both rounds). Rhizotron Petri dishes were randomly arranged in trays.
Positions of trays and Petri dishes were rearranged randomly every 3 d. Pictures of rhizotrons were taken after 4 wk. Image analysis was done with the ImageJ software (Schneider et al., 2012). The number of attachment organs was counted. The total number of pre-germinated P.ramosa seeds was recorded as a co-variable.

Genotypic means and definition of target traits
Values for diameter of attachment organs, number of attachment organs and number of pre-germinated seeds were log-transformed for normality and were averaged over technical replicates where present. Since there was significant correlation between the two variables of interest and the number of pre-germinated P.ramosa seeds, we used the residuals from the regression on the number of pre-germinated seeds for further analysis.

Trait
Nematode: Plant response to nematodes (Melodoigyne incognita, 180 stage-2 juveniles), in terms of number of M. incognita egg masses per plant.

Growing conditions
Seeds were vapor-sterilized for 5 h and transferred to a six-well plate with MS20 (5% gelrite). After four nights in the dark at 4°C the plates were transferred to 24°C in 12 h light. At the age of 1 wk the seedlings were transferred individually to a well of a six-well plate. Melodoigyne incognita infection was induced with 180 juveniles stage-2 added to 2-wk-old seedlings. Six-well plates with nematodes and seedlings were incubated in the dark at 24°C for 6 wk. Plants were grown for 2 wk 24°C, 12 h : 12 h, light : dark, then 6 wk 24°C, dark. Egg masses were quantified manually.

Experimental design
Plants were screened in rounds of 20 accessions. Each round included a six-well plate with one Col-0 plant as reference.

Genotypic means
'Nematode' was defined as the number of egg masses, after arcsine-square root transformation. Genotypic means were calculated using the following mixed model: where GEN is genotype (accession) and RND is a random effect for round.

Traits
Whitefly_1: Plant response to whitefly (Aleyrodes proletella, five females), in terms of whitefly survival Whitefly_2 Plant response to whitefly (Aleyrodes proletella, five females), in terms of number of eggs

Growing conditions
Plants were grown for 5 wk at 20 o C, 70% RH, 100 µmol m -2 s -1 light intensity and 10 h : 14 h, light : dark photoperiod. One leaf of each accession was infested with five female whiteflies (placed in clip cages) that were allowed to feed and oviposit. Seven days after infestation, the number of living and dead females was counted as well as the number of eggs. From this, we calculated the survival (number of living flies divided by the total number of flies) and oviposition rate (eggs laid per female d -1 ).

Experimental design
Accessions were screened in three blocks of 120 accessions with five reference accessions (Col-0, Ler-1, WS-0, Cvi-0, Kin-0) in each block. The whole experiment was repeated five times.

Genotypic means
Genotypic means were calculated with Genstat 15th edition (Payne, 2009), using the following mixed model: where GEN is genotype (accession), REP denotes complete replicates (experiments with three blocks) and

REP:BLOCK is a random effect for incomplete blocks within replicates.
14 Methods S7 Aphids.

Aphid_1
Plant response to aphids (Myzus persicae), in terms of behavior at t1 Aphid_2 Plant response to aphids (M. persicae), in terms of behavior at t2 Aphid_3 Plant response to aphids (M. persicae), in terms of aphid reproduction Methods wingless aphid was released on the leaf disc and cling film was used to cover the arena. Twenty Arenas were recorded simultaneously with a mounted camera. Aphids were observed for 85 min on two time points: (1) immediately after introducing the aphids into the arenas; and (2) 4.5 h after the start of the first observation.
Motion analysis was performed with EthoVision XT ® 8.5 software (Noldus Information Technology bv, Wageningen, The Netherlands). Aphids are phloem-feeding insects and probe with their piercing mouthparts between plant cells to feed from the plant sap. Start time and duration of probes were registered with automated video-tracking (Kloth et al., 2015). Aphid survival was checked 24 h after recording. Subject detection was checked on four time points within each movie. Samples with no survival, low subject detection or with less than five replicates were excluded from analysis. Probes were categorized into short (<3 min) probes and intermediate (<15 min) probes, both associated with penetration of the plant epidermis or 15 mesophyll, and long (≥15 min) probes, putatively associated with phloem uptake. Response variables expressed in seconds were arcsin or logit transformed to approach a normal distribution.
Aphid reproduction was measured in a whole-plant assay. Each 2-to-3-wk-old plant was inoculated with one 0-to-24-h-old nymph. Two weeks after infestation, aphid population size was measured per plant.
Plants were placed in a Petri dish in trays with a water barrier to prevent aphids to move between plants. Each tray contained 20 plants, none of the aphids developed wings.

Experimental design
Automated video tracking of aphid behaviour was performed in an incomplete block design with each complete replicate consisting of 18 incomplete blocks of 20 accessions. One replicate of the complete Hapmap collection was acquired in 6 d, 60 plants were screened each day across three batches. An alpha design was generated with Gendex (http://designcomputing.net/gendex/) to assign accessions to blocks. For each accession five to six replicates were acquired.
Phenotyping of aphid reproduction was performed in an incomplete block design with seven incomplete blocks. Blocks were defined according to the position in the climate cell. Each replicate consisted of three to four blocks and plant genotypes were randomized across blocks between replicates. For each genotype two to three replicates were acquired.

Genotypic means and definition of target traits
Genotypic means were calculated using the following linear mixed model:

Experimental design
Plants were screened in five rounds (complete replicates) of 360 accessions, using an incomplete block (alpha) design. Within each round plants were randomly allocated to 18 blocks of 20 accessions, the blocks representing plants being screened in one recording. One sampling day consisted of five blocks (100 accessions), with the exception of the last day (three blocks, 60 accessions).

Genotypic means
Genotypic means were calculated using the following linear mixed model: light intensity, and 10 h : 14 h SD photoperiod. Projected leaf area, rosette feret and rosette perimeter were measured using ImageJ software in the P. rapae and PEG8000 combined treatment group (Schneider et al., 2012). Data were recorded at three time points: T1, before applying P. rapae; T2, before applying PEG8000 treatment; T3, after 7 d PEG8000 treatment. Rosette fresh weight was measured from both control and combinatorial stress treatment groups at T3. We used 332 Arabidopsis accessions, grown on rock wool blocks in a climate controlled growth chamber. Plants were first treated with P. rapae L1 larvae for 24 h, and then irrigated with nutrient solution that containing PEG8000 for 7 d (P. rapae and PEG8000 combined stress treatment). In additional, plants were grown without any stress treatment (control).

Experimental design
All traits were measured in a randomized complete block design with two complete blocks (replicates) under treatment conditions and 2 complete blocks under control conditions.

Genotypic means
The square root transformation was first applied to the area traits. For each treatment, genotypic means were calculated, using a linear model with a fixed effect for block.

Definition of target traits
We regressed the genotypic means of fresh weight at T3 after P. rapae and PEG8000 treatment on the means of fresh weight at T3 under control conditions; these residuals represent the effect of P. rapae and PEG8000 treatment on fresh weight (Caterpillar_&_osmotic_2). Similarly, we regressed each of the rosette area related traits (projected leaf area, rosette feret, and rosette perimeter) observed at T2 (only P. rapae treatment) on the corresponding means of these traits measured at T1 (Caterpillar_1). The resulting residuals represent the effect of P. rapae treatment on rosette area related traits. Finally we performed the regression of rosette area related traits at T3 on the values measured at T1 as well as T2, whose residuals represent the combined effect of P. rapae and PEG8000 treatment on rosette area related traits (Caterpillar_&_osmotic_1). In both the Caterpillar_1 and Caterpillar_&_osmotic_1 group, the three traits were replaced by the first principal component.

19
Methods S10 Fungus -combinatory stress. Lesions that did not exceed the size of the droplet, (5 µl) were scored as zero, whereas a spreading lesion was scored as a one.

Experimental design
Plants were screened in rounds of 35 accessions. Col-0 was present in all rounds as a control. The three treatments were screened simultaneously.

Genotypic means and definition of target traits
An arcsine transformation was applied to the proportion of leaves with spreading lesions, that is, for each observed count k = 0,1,2,3,4,5,6, the transformed phenotype was defined as arcsin√((k)/6). Before transformation, counts equal to zero or 6 were replaced by respectively 1/4 and 5.75=6-1/4. The transformed phenotypic observations were corrected for round effects by subtracting from each observation the mean of the round it was contained in, and genotypic means were calculated based on the round corrected phenotypes.
Differential sensitivity of each double stress was calculated as the residuals obtained from the linear regression of the double stress on the single stress phenotype.
Methods S11 Screening of T-DNA lines.
T-DNA lines were ordered and screened for homozygosity, using primers described in Methods Table M6.
Seeds from homozygous mutants were harvested and grown and screened individually by consortium partners (Methods Table M7). The top 300 coexpressed genes of RMG1 were retrieved from Atted-II version 8.0 (Obayashi et al., 2014). GO enrichment analysis was performed with the application BiNGO in Cytoscape (Cline et al., 2007;Maere et al., 2005). Methods S12 Simulations to compare power for full MTMM, contrast MTMM and univariate analysis.

Methods
To compare the performance of the full and contrast MTMM, we repeatedly simulated 30 traits for all of the n = 350 Hapmap accessions, for two scenarios. In scenario A, SNP-effects had the same sign within two predefined groups of 15 traits, whereas in scenario B, each SNP-effect was given a randomly chosen sign. For both scenarios, each of the 1000 simulations was performed as follows: • The 30 x 30 matrix containing the genetic covariances was simulated using a first order factor analytic model: • Given and , we simulated the 350 x 30 matrices and , containing respectively the genetic and residual effects. followed a matrix normal distribution with row-and column covariance K and , where K is the 350 x 350 genetic relatedness matrix. Equivalently, the length 10500 vector ( ) obtained by stacking all columns of followed a multivariate normal distribution with covariance ⨂ . Similarly, ( ) followed a zero mean normal distribution with covariance ⨂ .
• We randomly selected one of the 214051 available SNPs, under the restriction that its minor allele frequency was at least 0.4. This restriction is not essential, and only serves to avoid extra variation in the simulation results due to varying allele frequencies (which affect power).
• We defined a 350 x 30 matrix S by multiplying the vector of SNP-scores (x) with trait specific SNPeffects: the jth column of S was defined as x β j . The magnitude of the effects β j was chosen such that the explained variance for the jth trait was 1% of the polygenic variance of that trait (i.e. the jth diagonal element of ). In scenario A, the effects were made negative for the first 15 traits and positive for traits 16,...,30, whereas in scenario B each SNP-effect was given a randomly chosen sign.
• The matrix of phenotypes was defined as Y = G + E + S, which was used for all subsequent analyses.
• We fitted a first-order analytic model without marker effects, using the 37 x 37 compressed kinship matrix used for the full MTMM in the main text (for computational reasons the latter was kept constant throughout all simulations).
• As in the MTMM analyses in the main text, we tested the significance of marker effects conditional on the variance components estimated in the previous step. Using the Wald test, we tested the

22
In both scenarios, the power of the full MTMM, contrast MTMM and univariate analysis was estimated by the proportion of simulations where the -log 10 (P) value was above a certain threshold. Narrow sense heritability estimated using the 'heritability' R package.    Table 2b) in plants exposed to biotic or abiotic stress factors, relative to control conditions. (a) Shoot tissues and (b) root tissues. Black bars represent abiotic stresses, red bars represent biotic stresses and purple bars represent phytohormonal treatments. Expression data from Arabidopsis eFP browser (http://bbc.botany.utoronto.ca). ) for the common response are clustered according to trait-specific effects estimated from the full MTMM. If there was another SNP in LD that had a higher effect size, this SNP was used as representative for the LD block. Negative effect sizes (blue) were cases where the rare allele was associated with a detrimental effect on the plants, positive effect sizes (yellow) were cases where the rare allele was associated with increased resistance to the stress. The rare alleles of the top 9 SNPs are associated with enhanced resistance to abiotic and biotic stresses; the bottom 11 SNPs are associated with reduced resistance to abiotic and biotic stresses. Stresses are clustered according to SNP effect size, using Ward's minimum variance method. If SNPs were located within a 20 kb half-window of each other, only the SNP with the highest absolute cumulative effect size was included. The key shows the frequency distribution of SNPs across effect sizes.

46
Fig. S8 Expression data of three candidate genes (resulting from MTMM, see Table S6) in plants exposed to biotic or abiotic stress factors, relative to control conditions. (a) Shoot tissues and (b) root tissues. Black bars represent abiotic stresses, red bars represent biotic stresses and purple bars represent phytohormonal treatments. Expression data from Arabidopsis eFP browser (http://bbc.botany.utoronto.ca).

Fig. S9
Genetic associations specific for plant responses to either below-or aboveground stress. Genetic associations were estimated with a contrast analysis using MTMM. Significant SNPs (P ≤ 10 -4 ) for the belowground-aboveground contrast are clustered according to trait-specific effects estimated from the full MTMM. If there was another SNP in LD that had a higher effect size, this SNP was used as representative for the LD block. Negative effect sizes (blue) were cases where the rare allele was associated with a detrimental effect on the plants, positive effect sizes (yellow) were cases where the rare allele was associated with increased resistance to the stress. The rare alleles of the top 12 SNPs are associated with enhanced resistance to aboveground stresses and reduced resistance to belowground stresses; the bottom 8 SNPs show the inverse.
Stresses are clustered according to SNP effect size, using Ward's minimum variance method. If SNPs were located within a 20 kb half-window of each other, only the SNP with the highest absolute cumulative effect size was included. The key shows the frequency distribution of SNPs across effect sizes. ) for the common response are clustered according to trait-specific effects estimated from the full MTMM. If there was another SNP in LD that had a higher effect size, this SNP was used as representative for the LD block. Negative effect sizes (blue) were cases where the rare allele was associated with a detrimental effect on the plants, positive effect sizes (yellow) were cases where the rare allele was associated with increased resistance to the stress.
The rare alleles of the top 5 SNPs are associated with enhanced resistance to above-and belowground stresses; the bottom 11 SNPs are associated with reduced resistance to above-and belowground stresses. Stresses are clustered according to SNP effect size, using Ward's minimum variance method. If SNPs were located within a 20 kb half-window of each other, only the SNP with the highest absolute cumulative effect size was included.
The key shows the frequency distribution of SNPs across effect sizes. univariate analysis (red) as a function of P-value thresholds, in case of contrasting SNP-effects (Scenario A) and SNP-effects with random sign (Scenario B). Power was estimated based on 1000 simulations, which were performed as described in Methods S12. shown. Markers within 20 kb half window size of each other were considered to be in LD.