Introduction

Under natural conditions, plants encounter many different stresses individually, sequentially or simultaneously. To withstand these stresses, plants make use of constitutive and induced defenses to ensure survival. Induced plant responses are often regulated by phytohormones, such as jasmonic acid (JA), salicylic acid (SA), ethylene (ET) and abscisic acid (ABA) that accumulate in the plant upon infection by a pathogen or infestation by an insect herbivore (Pieterse et al. 2012). Depending on the type of stress, plants activate different hormone-regulated responses, in which JA and ABA signaling is typically engaged upon herbivory or wounding, JA and ET signaling upon infection by necrotrophic pathogens, and SA signaling after infection by biotrophic pathogens. Cross-talk between these hormone signaling pathways is thought to play an important role in the ability of plants to swiftly adjust to a variety of biotic and abiotic environmental conditions, and has, therefore, been the focus of many studies on plant stress adaptation (Yamaguchi-Shinozaki and Shinozaki 2006; Pieterse et al. 2012; Vos et al. 2013; Caarls et al. 2015).

Plant responses to single-stress factors are subject to complex regulatory mechanisms. Environmental conditions and additional stress factors can greatly influence the outcome of a plant’s response to a single stress (Suzuki et al. 2014). In nature, plants co-evolved with a large variety of biotic and abiotic stresses. Hence, genetic variation amongst naturally occurring accessions can provide insight into naturally evolved plant adaptive responses. To study the genetics of these responses, genome-wide association (GWA) mapping has been used in many studies, revealing genes with important functions in diverse processes of plant growth and survival, whether or not under stress conditions (Atwell et al. 2010; Baxter et al. 2010; Thoen et al. 2017; Proietti et al. 2018).

One of the most notorious necrotrophic plant pathogens is Botrytis cinerea, which is responsible for yield losses in a wide range of crop plants (Dean et al. 2012). Primary necrotic lesion formation by B. cinerea is affected by the ability of the fungus to overcome plant constitutive defenses. For the outgrowth and success of B. cinerea infection at the stage of secondary (or spreading) lesion formation, inducible plant defenses need to be overcome (Van Kan 2006). At the primary stage of infection, B. cinerea produces a large array of enzymes and toxins that target plant cell walls and facilitate pathogen invasion (Blanco-Ulate et al. 2014). Breakdown of pectin in the plant cell wall results in the release of plant self-molecules (e.g. oligogalacturonides) that act as damage-associated molecular patterns (DAMPs). Together with pathogen-associated molecular patterns (PAMPs), such as chitin from the fungal cell wall, these elicitors contribute to the onset of the plant’s primary defense response. Fungal chitin is recognized via the LYSIN MOTIF (LysM) RECEPTOR KINASE 5 (LYK5) and the receptor-like kinase (RLK) CHITIN ELICITOR RECEPTOR KINASE 1 (CERK1) complex. Association of this complex with the kinases BOTRYTIS-INDUCIBLE KINASE 1 (BIK1) and PBS1-LIKE 27 (PBL27) leads to activation of downstream mitogen-activated kinases (MAPKs), ultimately resulting in the activation of plant defense responses (Liu et al. 2012; Cao et al. 2014; Le et al. 2014; Shinya et al. 2014; Yamada et al. 2016). B. cinerea-elicited defenses consist of the production of reactive oxygen species (ROS), nitric oxide (NO), callose, glucanases, chitinases, and anti-microbial phytoalexins, which collectively help the plant to limit pathogen ingress (Glazebrook 2005; Van Kan 2006; Laluk and Mengiste 2010). Important regulators of B. cinerea-induced defenses are the phytohormones JA and ET, which accumulate upon infection and in turn activate defense-related genes such as PLANT DEFENSIN1.2 (PDF1.2) (Thomma et al. 1998). Plant defenses against B. cinerea also involve the plant hormones ABA and SA (Audenaert et al. 2002; Ferrari et al. 2003; Liu et al. 2015; Vos et al. 2015). Their signaling pathways cross-communicate with those of JA and ET, therewith modulating plant resistance against B. cinerea infection (Thomma et al. 1998; El Oirdi et al. 2011; Vos et al. 2015). Hence, interactions between the different hormone-controlled signaling pathways may be decisive in the outcome of the B. cinerea resistance response during multi-stress conditions.

Attempts to unravel the genetics and transcriptional dynamics behind plant responses to necrotrophs, such as B. cinerea, revealed many genes associated with plant resistance to B. cinerea (Denby et al. 2004; AbuQamar et al. 2006; Rowe and Kliebenstein 2008; Davis et al. 2009; Anuradha et al. 2011; Windram et al. 2012; Hickman et al. 2017). They involve components of structural barriers, general stress signals such as ROS, different hormonal biosynthesis and signaling pathways (predominantly JA and ET), glucosinolates and phytoalexins (e.g. camalexin), the plant’s nutrient status, and the circadian rhythm (Audenaert et al. 2002; Bessire et al. 2007; Ferrari et al. 2007; Manabe et al. 2011; Ingle et al. 2015). However, still little is known about their role in B. cinerea resistance in multi-stress environments, such as frequently occurs in natural and agricultural environments.

In a previous study, we investigated the effect of prior herbivory by Pieris rapae and drought stress on the dynamics of the transcriptional response of Arabidopsis thaliana to B. cinerea infection (Coolen et al. 2016). We showed that when a plant is encountering different stresses in sequence, it is well capable of swiftly adapting its transcriptional response to the last stress encountered. The transcriptional response to herbivory or drought stress was largely overruled by the plant’s response to B. cinerea infection, although transcriptional signatures of the first stress were still detectable in the second stress transcriptional profile (Coolen et al. 2016). In the present study, we aimed to gain insight into the genetic basis of the plant’s ability to adapt to these sequential stresses in naturally occurring A. thaliana accessions. To this end, we studied the genetic basis of the phenotypic differences that we observed in B. cinerea disease severity among 346 globally collected A. thaliana accessions during single- and sequential double-stress conditions. The preceding biotic stress factor studied here is herbivory by P. rapae caterpillars, one of the most destructive plant pests on cruciferous plants. Plant defenses against this attacker are regulated by the phytohormones JA and ABA (Vos et al. 2013). The preceding abiotic stress factor that was used in this study is drought stress, one of the most-encountered abiotic stresses in agriculture. The phytohormone ABA plays an important role in the adaptation of plants to this stress (Yamaguchi-Shinozaki and Shinozaki 2006). The combination of these prior stresses with B. cinerea infection is particularly interesting because the plant’s adaptive response to these stresses shares components of the JA and ABA signaling pathways (Pieterse et al. 2012; Caarls et al. 2015). Hence, this study could potentially reveal novel players in cross-communication between the involved stress pathways. To study the genetics behind plant adaptive responses to B. cinerea infection under single- and multi-stress conditions, we performed GWA mapping on the 346 phenotyped A. thaliana accessions. These accessions are thought to have evolved under different selection pressures depending on their habitat, potentially leading to genetic variation that can be mined for naturally occurring plant adaptive responses to multiple stresses. To aid in the selection of candidate genes, we additionally performed quantitative trait loci (QTL) mapping and fine mapping using available full genome sequences of diverse A. thaliana accessions. Our results reveal several candidate genes with potential roles in the plant’s adaptive response during sequential double stresses.

Materials and methods

Plant material

In this study, 346 A. thaliana (L.) Heynh. accessions (Suppl. Tables S1–S3) of the haplotype map (HapMap) population (http://bergelson.uchicago.edu/wp-content/uploads/2015/04/Justins-360-lines.xls) were included, which are genotyped for ~ 250,000 bi-allelic single-nucleotide polymorphisms (SNPs) (Baxter et al. 2010; Platt et al. 2010; Chao et al. 2012). After quality control and imputation, this SNP set was reduced to a set of 214,051 SNPs (Thoen et al. 2017). Furthermore, a set of 431 multiparent advanced generation inter-cross (MAGIC) recombinant inbred lines (RILs) (Suppl. Table S6) was included that are genotyped for 1260 SNP markers (Kover et al. 2009; Kover and Mott 2012).

Plant growth conditions

Seeds of the A. thaliana accessions were sown in cultivation containers filled with autoclaved river sand. Sand was supplied with half-strength Hoagland solution containing sequestrene as described (Van Wees et al. 2013). To attain a high relative humidity (RH) for germination, cultivation containers were enclosed in a tray with water and covered with a transparent lid. Seeds were stratified for 2 days at 4 °C in the dark to ensure a homogeneous germination, after which the tray was moved to a growth chamber with an 8-h day/16-h night rhythm, a temperature of 21 °C, and a light intensity of 100 µmol m−2 s−1. After 8 days, the lids of the trays were slightly opened and gradually removed over a 2-day period to adjust to the 70% RH present in the growth chamber. Ten-day-old seedlings were transplanted into individual pots containing an autoclaved mixture of river sand and potting soil (1:1, v:v). Pots were supplied with water from the bottom up three times per week. At an age of 3 weeks, the plants were supplied once with half-strength Hoagland solution.

Rearing of P. rapae caterpillars

Pieris rapae caterpillars were reared on cabbage plants (Brassica oleracea convar. capitata var. alba) under greenhouse conditions (24 °C, with natural daylight). Butterflies were supplied with flowering plants such as Lantana camara for their (nectar) food. When flowers were scarce, additional food (solution of 20% honey and 10% sucrose) was offered to the butterflies. Inbreeding of the population was minimalized by adding wild butterflies and caterpillars from the Dutch Flevopolder to the existing population. First-instar (L1) larvae were placed on A. thaliana leaves using a fine paint brush as described by Van Wees et al. (2013).

Cultivation of B. cinerea

Botrytis cinerea strain B05.10 (Staats and Van Kan 2012) was grown on half-strength Potato Dextrose Agar (PDA; Difco BD Diagnostics, Franklin Lakes, NJ, USA) plates containing penicillin (100 µg ml−1) and streptomycin (200 µg ml−1) for 2 weeks at room temperature. Spores were collected, filtered through glass wool, and re-suspended in half-strength Potato Dextrose Broth (PDB; Difco BD Diagnostics) to a final density of 105 spores ml−1. After a 3-h incubation period, the spores were used for inoculation by applying 5-µl droplets on A. thaliana leaves as described (Van Wees et al. 2013).

Stress treatments and experimental design

Single- and sequential double-stress treatments were applied according to the schedule shown in Fig. 1. Pre-treatment with drought was started when plants were 4 weeks old by withholding water for 7 days, after which plants were re-watered and allowed to recover for 1 day before plants were inoculated with B. cinerea. P. rapae pre-treatment was started 1 day prior to B. cinerea inoculation by placing a single P. rapae L1 caterpillar on each plant and allowing it to feed for 1 day, after which it was removed. All plants were simultaneously inoculated with B. cinerea and were kept under 100% RH for the remaining time period until disease development was registered after 3 days of B. cinerea infection. As disease measurement, the occurrence of spreading lesions was recorded. Lesions that did not exceed the size of the inoculation droplet of 5 µl were scored as non-spreading. In total, six leaves per plant were inoculated and the three treatments were screened simultaneously. For our GWA experiments, each treatment was performed on six plants (biological replicates) of each accession. Furthermore, plants were screened in 11 experiments of approximately 35 randomly chosen accessions and Col-0 was present in all experiments as a reference. The B. cinerea disease severity in Col-0 in all 11 experimental rounds is given in Suppl. Fig. S1. For the QTL analysis experiments, 431 RILs of the MAGIC collection were screened individually in two experiments.

Fig. 1
figure 1

Experimental design of the timing of the single- and sequential double-stress treatments. Single-stress B. cinerea-treated plants (Bc) were inoculated on six leaves with B. cinerea spores when 4 weeks old. In the P. rapae/B. cinerea sequential double stress (Pr + Bc), a single P. rapae larva was allowed to feed for 24 h on a single leaf per plant (day 7) prior to inoculation of six leaves with B. cinerea (day 8). In the drought/B. cinerea sequential double stress, plants were withheld from watering for 7 days (day 0–7), after which they were re-watered and allowed to recover for 1 days (day 7) prior to inoculation of six leaves with B. cinerea (day 8). The percentage of leaves with spreading lesions was scored 3 days after inoculation (day 11). The treatments were synchronized so that all treatments were inoculated with B. cinerea simultaneously

Full size image

GWA mapping

GWA mapping was performed on genotypic means that were calculated from the raw disease severity data (Suppl. Table S3), as described for this dataset in a meta-analysis by Thoen et al. (2017). The original disease severity scores per plant of each of the 346 accessions (n = 6 plants per accession) were arcsine-square root transformed, where 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 0 or 6 were replaced by 1/4 and 5.75 = 6 − 1/4, respectively. Subsequently, the data were corrected for the variation within the 11 experiments (with ~ 35 randomly selected accessions each). This was done by subtracting the mean of each experiment from the scored values in that experiment. The corrected data were used for calculating genotypic means using the best linear unbiased estimator (BLUE). Residuals were calculated by regressing the disease severity scores of the prior stress/B. cinerea treatments on the predicted disease severity scores of the single stress (B. cinerea), which were obtained from a linear regression (Suppl. Table S3). The predicted prior stress/B. cinerea scores minus the predicted observations from the B. cinerea single-stress data represent the residual data.

GWA mapping was performed by employing factored spectrally transformed linear mixed models (FaST-LMM) software in R (Lippert et al. 2011; Thoen et al. 2017). A minor allele frequency (MAF) of > 0.05 was used together with an arbitrary threshold with a logarithm of odds (LOD; − log10(p)) score of 4.0 to determine SNP-trait associations of interest. Linkage disequilibrium (LD) was taken into account by including all SNPs within a 50-kb window surrounding the SNP of interest. Narrow sense heritability (h2) was estimated using the ‘heritability’ R package (Kruijer et al. 2015).

QTL mapping

QTL mapping was performed on transformed raw disease severity data of 431 tested MAGIC line RILs. The original mean disease severity scores per plant were divided by 100 and transformed using an arcsine-square root transformation where 0.00 was replaced by 0.04 (Suppl. Table S7). Residuals of the sequential double-stress data regressed on the B. cinerea single stress were calculated to obtain the relative effect of the sequential double stress on the single stress. QTL mapping was performed using publically available R-scripts for performing QTL mapping with the MAGIC line RILs (http://mtweb.cs.ucl.ac.uk/mus/www/magic/). Significant QTLs exceeded the genome-wide permutation threshold with p value < 0.1 (Kover et al. 2009; Kover and Mott 2012).

Fine mapping

Fine mapping was performed using full genome sequences of 164 available accessions from the 1001 genomes project (http://signal.salk.edu/atg1001/3.0/gebrowser.php, Suppl. Table S9). Genome sequences surrounding the SNP of interest with a 50-kb window and gene sequences (for mapping SNPs leading to amino acid changes within specific genes) were downloaded and aligned using Jalview (http://www.jalview.org/) (Waterhouse et al. 2009). Fine mapping was performed using GWA mapping phenotypic input data (Suppl. Table S3). Furthermore, a MAF of > 0.05 and a Kruskal–Wallis test was used for obtaining significant, false discovery rate (FDR)-corrected, SNP-trait associations using R and the ‘p. adjust’ function with the Benjamini and Hochberg (BH) method (Benjamini and Hochberg 1995).

Results

Experimental approach to study sequential stress responses

To study the genetic basis of the effect of prior abiotic and biotic stresses on the plant response to B. cinerea infection, 346 naturally occurring accessions of the A. thaliana HapMap collection (Baxter et al. 2010; Platt et al. 2010) were screened for susceptibility to infection by B. cinerea, either as a single stress, or when preceded by P. rapae herbivory or drought stress. These preceding stresses were chosen because the plant response to these stresses shares important signaling components with the defense response to B. cinerea, theoretically resulting in antagonistic or synergistic effects on the level of B. cinerea resistance during sequential stress conditions. For herbivory by P. rapae, a short 24-h period of exposure to herbivory was chosen (Fig. 1), because plant defense responses to this attacker are triggered very rapidly (Coolen et al. 2016). First effects of herbivore feeding on defense-related gene expression are typically detectable within 3 h after introduction of the caterpillar onto the leaf. At 12–24 h after infestation, the herbivory- and JA-responsive marker gene LIPOXYGENASE2 (LOX2) peaks in expression in this experimental setup (Coolen et al. 2016).

For drought stress, a moderate drought treatment was chosen in which plants were withheld from watering for 7 days after which they were re-watered and allowed to recover for 24 h (Fig. 1). The 7-day drought treatment resulted in visible growth retardation, dark coloration of the leaves, and first signs of wilting. By comparing leaf areas from normally watered and drought-stressed plants, it is clear that the 346 tested A. thaliana accessions respond differently to water withhold, leading to phenotypic variations that range from highly drought sensitive to highly drought tolerant (Suppl. Fig. S2 and Suppl. Table S1). Induced expression of the drought-responsive gene RAD18 (AT5g66400) is typically observed in this setup from day 5 of the water withhold (Coolen et al. 2016).

Plant responses to B. cinerea are detectable at the gene expression level within 6 h after application of the inoculum (Coolen et al. 2016). The infection itself becomes clearly visible at 48 h after inoculation with the start of watery lesions, rapidly followed by expansion of the lesions, and clearly visible cell death 3 days after inoculation (Fig. 2). As a readout of the level of B. cinerea resistance, we assessed the frequency by which B. cinerea inoculations inflict spreading lesions as a measure of the ability of B. cinerea to infect plant leaves in single- and double-stress situations. The experiment was designed in such a way that all measurements could be performed at the same time for all stress treatments (Fig. 1).

Fig. 2
figure 2

Phenotypic variation in B. cinerea resistance and effect of prior herbivory or drought stress in 346 A. thaliana accessions. aB. cinerea disease severity in average percentage of spreading lesions per accession (n = 6) at 3 days after pathogen inoculation. The accessions are ordered on the x-axis according to their level of disease severity under the single-stress condition (Bc; in blue). The B. cinerea disease severity values in the sequential double-stress treatments with P. rapae infestation (Pr + Bc) or drought (Dr + Bc) followed by B. cinerea inoculation are displayed with green and orange dots, respectively. In the panel on the right, photographs of representative plants of three contrasting accessions are depicted. The presented accessions are color coded with light blue, yellow and purple for Bur-0, Col-0 and Can-0, respectively, referring to their position in the graphs on the left. The six B. cinerea-inoculated leaves per plant are indicated with red dots. Leaves with damage caused by 24 h of P. rapae feeding are indicated with red arrows. b Relationship between the level of B. cinerea resistance in the 346 A. thaliana accessions in the single- versus the sequential double-stress conditions and between both sequential stresses. In the plots, a linear regression line was fitted with corresponding Pearson correlation (R2) values

Full size image

Phenotypic variation in plant responses to B. cinerea under sequential (a)biotic stress

Screening 346 naturally occurring A. thaliana accessions of the HapMap collection revealed that there is a large natural variation in plant susceptibility towards B. cinerea infection when the pathogen was applied as a single stress (Fig. 2a and Suppl. Table S2). Plant phenotypes ranged from highly susceptible, with all inoculated leaves heavily infected, to highly resistant phenotypes that hardly formed any spreading lesions. Herbivory by P. rapae caterpillars in the 24 h preceding the B. cinerea infection significantly reduced B. cinerea disease severity in 101 of the 346 A. thaliana accessions, whereas none of the accessions became significantly enhanced susceptible to the disease after herbivory (Student’s t test; p ≤ 0.05). Susceptibility to B. cinerea was also affected by drought stress in the 7 days preceding the infection, resulting in both increased susceptibility and increased resistance, depending on the accession. In 75 of the 346 tested accessions, drought stress significantly increased the level of B. cinerea resistance, while in 6 accessions the level of B. cinerea resistance was significantly reduced (Student’s t test; p ≤ 0.05). However, there was no correlation between the level of B. cinerea resistance and the level of drought tolerance in the tested set of 346 A. thaliana accessions (R2 = 0.029) (Suppl. Fig. S3). Figure 2b shows that the correlation between the single- and sequential double-stress datasets is moderate to strong in strength (R2 = 0.43–0.61). This suggests that the effects of both herbivory and drought stress on the level of B. cinerea resistance are partly dependent on the level of B. cinerea resistance observed in the single-stress condition.

Genome-wide association mapping of loci related to plant adaptive responses to B. cinerea under sequential (a)biotic stress

To analyze the genetic variation of the level of resistance to B. cinerea under conditions of single and sequential (a)biotic stresses, and to ultimately identify candidate genes that affect the level of B. cinerea resistance during multi-stress conditions, we performed a GWA study. For the GWA analysis, we used normally distributed and experimental design-corrected quantitative data of 346 A. thaliana accessions, representing the level of disease susceptibility in the single-stress condition and the two sequential double-stress conditions. Moreover, we calculated residuals for the relative effect of the sequential stress to the B. cinerea single stress (Suppl. Table S3). The residual data were used in the GWA analysis to search for SNPs that are associated with traits that influence the level of B. cinerea resistance, but are itself not part of the B. cinerea resistance response during the single-stress condition. To obtain an estimate of the breeding value for each trait, we calculated the narrow sense heritability (h2; Kruijer et al. 2015), which reflects the additive genetic effects that contribute to the observed phenotypes. Narrow sense heritability values ranged between 0.17 and 0.40 for the GWA input data, whereas for the residuals it ranged between 0.02 and 0.05 (Table 1). These values suggest that there is a moderate to low selection for the tested traits within the tested set of accessions.

Table 1 GWA input data summary
Full size table

To identify candidate genes with putative roles in B. cinerea resistance and in the effect of herbivory or drought stress on the level of B. cinerea resistance, we performed GWA mapping using the ~ 214k SNP set that is commonly used for GWA studies in A. thaliana (Kim et al. 2007; Atwell et al. 2010; Thoen et al. 2017). To this end, we used the B. cinerea disease severity scores and the residual data (Suppl. Table S3) to perform a GWA analysis using the factored spectrally transformed linear mixed model (FaST-LMM) algorithm (Lippert et al. 2011). SNP-trait associations of interest were selected by setting an arbitrary threshold with a LOD (− log10(p)) score of 4.0 (Fig. 3). Manhattan plots from the GWA analysis show several peaks with single (singleton) or multiple (string) SNP-trait associations (Fig. 3 and Table 2). In total, 117 SNP-trait associations were identified in 94 loci (14 strings and 80 singletons). Nine loci were shared between two traits, yielding a total of 85 unique loci with SNP-trait associations for B. cinerea resistance during single- and sequential double-stress conditions (Table 2 and Suppl. Table S4). Because the linkage disequilibrium (LD) of the natural population of A. thaliana accessions used for GWA mapping is estimated to be 10–50 kb (Nordborg et al. 2005; Kim et al. 2007), we considered all genes within 25 kb up- and downstream of each identified SNP to be candidates for the observed SNP-trait associations. When taking along a 50-kb window surrounding the 117 significant SNPs in the 85 loci, a total number of 233 candidate genes were identified. Of these 233 genes, 189 were found for only one trait, 43 genes were found for two traits, and 1 gene was found for three traits (Fig. 4 and Suppl. Table S5). This latter gene, AT1G20960 encoding an ATP-dependent RNA helicase, is shared between Pr + Bc sequential double stress, Pr + Bc sequential double-stress residuals and Dr + Bc sequential double stress.

Fig. 3
figure 3

GWA, QTL and fine mapping results of the effect of prior herbivory or drought stress on the level of B. cinerea resistance in 346 A. thaliana accessions. Manhattan plots (gray and dark gray) showing the − log10(p) values of the SNP marker-trait associations from the GWA mapping on A. thaliana chromosomes 1–5 (x-axis). Manhattan plots are given of the GWA analyses of B. cinerea resistance scores in 346 A. thaliana accessions exposed to B. cinerea alone (Bc), to herbivory by P. rapae followed by B. cinerea (Pr + Bc), and to drought stress followed by B. cinerea (Dr + Bc). The GWA output with the residual data (i.e., the relative effect of the prior stresses herbivory or drought on B. cinerea resistance in each accession) are given in Manhattan plots “Pr + Bc residuals” and “Dr + Bc residuals”. The dark blue line indicates the arbitrary LOD threshold of 4.0 (− log10(p) = 4.0). QTL mapping results with the A. thaliana MAGIC lines are indicated in green. Asterisks on the x-axis indicate significant QTLs (Kover et al. 2009; Kover and Mott 2012). Colored arrows indicate the nine GWA loci for which SNP-trait associations (LOD ≥ 4.0) were confirmed via fine mapping using the 50-kb window genome sequences of 164 of the 346 tested A. thaliana accessions

Full size image
Table 2 Number of SNP-trait associations, loci and candidate genes within a 50-kb window of the SNPs identified by GWA mapping
Full size table
Fig. 4
figure 4

Number of shared loci and candidate genes within a 50-kb window related to the SNP-trait associations obtained with the GWA analyses of B. cinerea resistance in single- and sequential double-stress conditions. Venn diagrams showing the number of shared loci with SNP-trait associations (a, LOD ≥ 4.0) and shared candidate genes within a 50-kb window of the SNP-trait associations (b, LOD ≥ 4.0) identified in the GWA analyses of B. cinerea resistance in 346 A. thaliana accessions exposed to B. cinerea alone (Bc; blue), to herbivory by P. rapae followed by B. cinerea (Pr + Bc; green), and to drought stress followed by B. cinerea (Dr + Bc; orange). The GWA output with the residual data (i.e., the relative effect of the prior stresses herbivory or drought on B. cinerea resistance in each accession) are given in “Pr + Bc residuals” (light green) and “Dr + Bc residuals” (yellow)

Full size image

Amongst the candidate genes are several genes previously associated with induced responsiveness of A. thaliana to either B. cinerea infection, herbivory by P. rapae, or drought stress. For instance, amongst the candidate genes associated with the B. cinerea single-stress response is LysM RLK1-INTERACTING KINASE 1 (LIK1), which encodes a LRR-RLK protein that is phosphorylated by the kinase CERK1 in response to chitin perception (Miya et al. 2007; Le et al. 2014). Amongst the candidate genes associated with the residuals of the P. rapae/B. cinerea sequential double-stress response is LOX2, which encodes one of the first enzymes of the JA biosynthesis pathway (Wasternack 2015) and is responsive to both herbivory by P. rapae and infection by B. cinerea (Coolen et al. 2016). Amongst the candidate genes associated with the response to the drought/B. cinerea double-stress treatment are YUCCA 7 (YUC7) and HISTIDINE KINASE 5 (HK5). YUC7 is implicated in drought tolerance, and ABA- and ET-mediated root growth inhibition, while HK5 is involved in stomatal closure and plant resistance against B. cinerea (Iwama et al. 2006; Desikan et al. 2008; Lee et al. 2012; Pham et al. 2012). These findings suggest that the followed GWA approach has the potential of identifying novel genes in plant adaptation to the single and double stresses applied in this study.

Support of GWA-identified loci by QTL mapping

To reduce the number of false positives and false negatives that are induced by the correction for population structure that is implemented within our GWA analysis, additional evidence is required to independently support the identified SNP-trait associations in the GWA mapping analyses (Nordborg and Weigel 2008; Bergelson and Roux 2010). Therefore, we performed quantitative trait loci (QTL) mapping using a set of MAGIC lines comprised of RILs descended from 19 heterogeneously intercrossed parental A. thaliana accessions (Kover et al. 2009; Kover and Mott 2012). Of these parental accessions, 12 were also present in the HapMap collection used in this study, representing accessions with levels of B. cinerea resistance that ranged from highly resistant to highly susceptible (Suppl. Fig. S4 and Suppl. Tables S6). In total, 431 RILs of the MAGIC population were tested for the level of resistance against B. cinerea under conditions of single or sequential double stress. The frequency distribution of the disease severity scores of the MAGIC lines under the different single- and sequential double-stress conditions (Suppl. Fig. S5; Suppl. Table S7) shows that the MAGIC lines displayed a large variation in susceptibility towards B. cinerea. From the frequency distribution of B. cinerea disease severity, it is clear that herbivory prior to B. cinerea inoculation significantly shifted the disease severity frequency distribution towards the lower disease severity classes (chi square, X2, p = 2.23−07), a phenomenon that was also observed in the 346 A. thaliana accessions of the HapMap collection (Fig. 2). Also the drought pre-treatment resulted in a shift in the B. cinerea disease severity frequency distribution, with significantly fewer MAGIC RILs in the lowest disease severity classes (X2, p = 0.05) and significantly more MAGIC RILs in the median disease severity class (X2, p = 0.05). QTL mapping of the obtained disease severity scores using available scripts in R (Kover et al. 2009; Kover and Mott 2012) revealed several significant QTLs that roughly map to positions on the A. thaliana genome as do some of the GWA-identified loci (Fig. 3, green line and Suppl. Table S8). No overlap was found with the high LOD score-associations found for the GWA mapping results on chromosome 1 (Fig. 3). The latter can be explained by the limited genetic variation in the MAGIC lines in comparison with the much larger genetic variation of the HapMap population. What is clear from the QTL mapping results is that there are two genomic regions, one on chromosome 2 and one on chromosome 5 (indicated by asterisks underneath Manhattan plots), that emerge after GWA mapping of the HapMap data sets and after QTL mapping of the MAGIC RILs dataset, and contain several significant marker-trait associations.

Furthermore, two of the significant QTL markers for the residuals of the P. rapae/B. cinerea sequential double-stress dataset are closely located to our GWA mapping results on chromosome 5 (Fig. 3), pointing to the locus AT5G55250 (Suppl. Fig. S6) containing the IAA CARBOXYLMETHYLTRANSFERASE 1 (IAMT1) gene, which is known to convert auxin indole-3-acetic acid to its methyl ester form MelIAA.

Fine mapping of GWA SNP-trait associations using genome sequences of 164 accessions

To further support the identified SNP-trait associations in the GWA mapping analyses, we made use of the genetic variation present in the available full genome sequences of the accessions used in this study (total 164 full genomes available, including genomes of 12 parental accessions of the MAGIC lines) to fine map the SNP-trait associations within the originally identified 85 unique loci of the GWA analysis (Fig. 3, Suppl. Table S4). For fine mapping, we used a Kruskal–Wallis test for significant SNP-trait associations with a minor allele frequency larger than 5% (MAF > 0.05), where significant associations exceed the false discovery rate (FDR)-corrected threshold for significance, and the 164 genomic sequences of each of the 85 GWA loci (spanning the 50-kb windows around the 117 SNPs with LOD ≥ 4.0). While the original GWA mapping procedure makes use of a selection of SNPs within the accessions relative to the reference genome of accession Col-0, the fine mapping approach makes use of all genetic variance within the used genomic regions of the 164 full genome sequences. In total, nine loci could be confirmed with significant SNP-trait associations (Fig. 5; Suppl. Fig. S7, and Suppl. Table S9), three for the B. cinerea single-stress data set (Fig. 3 blue arrows), two for the P. rapae/B. cinerea sequential double-stress data set (Fig. 3 green arrows), and four for the drought/B. cinerea sequential double-stress data set (Fig. 3, orange/yellow arrows).

Fig. 5
figure 5

Fine mapping results of confirmed GWA SNP-trait associations. Manhattan plots for the fine mapping results (gray dots) showing the − log10(p) values of the SNP marker-trait associations of the nine GWA loci for the Bc single stress (blue dots), Pr + Bc sequential double stress (green dots), Pr + Bc residual (light green dots), Dr + Bc sequential double stress (orange dots) and Dr + Bc residuals (yellow dots). On the x-axis the position in base pairs (bp) is shown. Genes indicated in blue, green, light green, orange and yellow show SNP-trait associations obtained with GWA mapping and the genes indicated in black correspond to the significant false discovery rate (FDR)-corrected) fine mapping associations

Full size image

For the B. cinerea single-stress data set, three significant loci were identified with fine mapping. The first significant association was found at the interval of AT2G33210 and AT2G33220 (Fig. 5), which contains two genes predicted to be transcribed in opposite directions from the point of the fine mapping association. The first gene is HEAT SHOCK PROTEIN 60-2 (HSP60-2; AT2G33210; Table 3), which has been described as being responsive to B. cinerea infection (Windram et al. 2012; Coolen et al. 2016). The second gene is AT2G33220, which codes for an oxidation–reduction process-related protein with homology to the mammalian cell death-related protein GRIM-19 (Huang et al. 2004). The second significant association was found in the genomic region of AT3G04130, AT3G04140 and AT3G04150, close to the GWA SNP-trait association located in AT3G04110. The first gene (AT3G04130) codes for a tetratricopeptide repeat (TPR)-like superfamily protein, the second gene (AT3G04140) for an ankyrin-repeat family protein, and the third gene (AT3G04150) for a RmlC-like cupin superfamily protein, all with unknown functions. The third significant association was found at position AT5G52040, containing a gene encoding ARGININE/SERINE-RICH SPLICING FACTOR 41 (RS41). RS41 is known to bind primary microRNA (pri-miRNA) in vivo and is involved in miRNA biogenesis (Chen et al. 2015). An A. thaliana double mutant rs40, rs41-containing mutations in RS41 and its closely related homologue RS40, was shown to have reduced DICER LIKE-1 (DCL1) expression. Mutation of DCL1 in the miRNA biogenesis mutant dcl1-7 was shown to result in enhanced susceptibility to B. cinerea (Weiberg et al. 2013).

Table 3 Fine mapping-confirmed candidate genes associated with B. cinerea resistance during single and sequential double stress
Full size table

For P. rapae/B. cinerea sequential double-stress data, a significant association was found with position AT4G16590, representing CELLULOSE SYNTHASE-LIKE A1 (CSLA01; Table 3). Impairment of cellulose synthases was previously shown to enhance disease resistance against necrotrophic fungi such as Plectosphaerella cucumerina, confirming its negative effect on necrotrophic fungi (Hernández-Blanco et al. 2007). Furthermore, a significant association was found for the P. rapae/B. cinerea sequential double-stress residual data with locus AT1G23760 containing a gene coding for POLYGALACTURONASE 3 (PG3). PG3 is a pectin-degrading enzyme that is involved in many stages of plant development (Cao 2012). However, pectin degradation is also part of the infection strategy used by B. cinerea and deletion of the B. cinerea PG2 gene was shown to result in delayed primary lesion formation and a strong reduction of virulence on tomato and broad bean (Kars et al. 2005). Furthermore, inhibition of polygalacturonase proteins in A. thaliana was shown to increase resistance to the flower fungal pathogen Fusarium graminearum (Ferrari et al. 2012).

The drought/B. cinerea sequential double-stress data yielded a significant association in our fine mapping results with position AT1G23730, representing BETA CARBONIC ANHYDRASE 3 (BCA3; Table 3). BCA3 is implicated in photosynthetic carbon fixation although it is not only localized in the cytosol of foliar tissue, but also in root tissue (Fabre et al. 2007). Furthermore, the expression of BCA3 was found upregulated 24 h after inoculation in the drought/B. cinerea sequential double-stress combination (Coolen et al. 2016). For the drought/B. cinerea sequential double-stress residual data, three significant associations were identified at the first locus, corresponding to a gene with unknown function (AT1G23340), a gene encoding a plant invertase/pectin methylesterase inhibitor (AT1G23350) and KNOTTED 1-LIKE HOMEOBOX GENE 6 (KNAT6, AT1G23380). The plant invertase/pectin methylesterase inhibitor was previously described as an important player in defense against B. cinerea infection, as overexpression was shown to restrict infection by B. cinerea (Lionetti et al. 2007). Homeodomain transcription factor KNAT6 was previously shown to be involved in plant flowering and activation of stress pathways that promote JA biosynthesis (Khan et al. 2015). Another strong and significant association for the drought/B. cinerea sequential double-stress residual data was confirmed at AT5G11610 on chromosome 5, representing a gene coding for an exostosin family protein with unknown function. Finally, a significant association was found at an F-box-associated ubiquitination effector family gene (AT5G41720) with unknown function.

To identify differences in the amino acid sequences of the 14 fine mapping-confirmed candidate genes that could explain the natural variation amongst the accessions in the level of B. cinerea resistance (Table 3), we analyzed the nucleotide sequences of these specific genes in the 164 full genomes in more detail. For 3 of the 14 fine mapping-confirmed genes, we identified differences in the predicted amino acid sequences that were significantly associated with the level of B. cinerea resistance (Suppl. Fig. S8). These include HSP60-2 (AT2G33210) and an ankyrin-repeat family protein (AT3G04140) for the single-stress data, and an exostosin family protein (AT5G11610) for the drought/B. cinerea residual data. In the predicted amino acid sequences of the other 11 candidate genes, no significant associations between amino acid sequence differences and phenotype could be detected, suggesting that the genotype–phenotype association for these genes may be caused by differences at the transcriptional level.

Expression patterns of the fine mapping-confirmed candidate genes

By looking at the expression patterns of A. thaliana (Col-0) in the first 24 h after B. cinerea inoculation in the single- and sequential double-stress treatments (Coolen et al. 2016), the fine mapping-confirmed genes are shown to be responsive to the applied stresses (Fig. 6). This may further indicate a role of these genes in the stress phenotypes that we observed. For instance, HSP60-2 (AT2G33210) is induced at 12 h after B. cinerea inoculation and continues to increase in expression until the last measured time point, indicating its responsiveness to B. cinerea. Furthermore, both P. rapae and drought-stress pre-treatments repress the induction of HSP60-2 until 18 h after B. cinerea inoculation, indicating that both pre-treatments change the single-stress B. cinerea expression pattern, making this an interesting candidate gene to study the possible consequence of the repression of HSP60-2 during early B. cinerea infection when preceded by another stress. Another example is the ankyrin-repeat family protein gene AT3G04140 that is repressed during the first 24 h of the B. cinerea infection and is induced until 12 h after inoculation in the drought/B. cinerea sequential double-stress combination, pointing to a possible involvement of this gene in the drought/B. cinerea sequential double stress. Also BCA3, from the drought/B. cinerea sequential double-stress dataset, shows to be differentially expressed compared to the single-stress B. cinerea treatment. BCA3 (AT1G23730) was shown to be repressed during the first 18 h after B. cinerea inoculation in the drought/B. cinerea sequential double-stress treatment, whereas this is not the case for single-stress B. cinerea. This is also the case for the exostosin family protein with unknown function, AT5G11610, which is induced in the B. cinerea single-stress treatment from 18 h onwards and is repressed in the drought/B. cinerea sequential double-stress data for the whole first 24 h after B. cinerea inoculation. And last, PG3 (AT1G23760) that is known to be involved in pectin degradation shows to be strongly repressed in the P. rapae/B. cinerea sequential double-stress combination at 24 h compared to the repression seen in the B. cinerea single stress. Since the inhibition of pectin degradation in A. thaliana was shown to affect the fungal pathogen F. graminearum (Ferrari et al. 2012), it would be interesting to study if the strong repression as observed in the P. rapae/B. cinerea sequential double-stress data also affects B. cinerea infection, leading to more resistant plants.

Fig. 6
figure 6

Expression patterns of fine mapping-confirmed candidate genes. Heatmap showing gene expression patterns of the 14 fine mapping-confirmed candidate genes (Table 3) with putative roles in the plant response to B. cinerea infection during single- or sequential double-stress conditions. Gene expression profiles are based on RNA-seq data obtained from Coolen et al. (2016) and represent the log2 fold change in gene expression, with addition of a pseudocount of 0.0001 during the first 24 h after inoculation with B. cinerea in the single- and sequential double-stress treatments. Colored dots indicate associations of genes with single or sequential double-stress treatment. Bc (blue dot), B. cinerea single stress; Pr + Bc (green dot), P. rapae/B. cinerea sequential double stress; Dr + Bc (orange dot), drought/B. cinerea sequential double stress. The candidate genes were clustered using hclust. Blue–red color keys show the change in gene expression level: − 2.0 > log2 fold change > 2.0

Full size image

All together, the combination of GWA mapping, QTL mapping and usage of fine mapping with full genome sequences confirmed nine loci with potential roles in the plant adaptive response to infection by the necrotrophic fungus B. cinerea during sequential biotic and abiotic stresses. Several of these genes have previously been implicated in plant defenses against necrotrophic pathogens (Table 3), confirming that the followed approach has the potential of identifying novel players in the plant’s adaptive response to the single and double stresses. Further research on the candidate genes will be required to reveal their role in plant stress signaling and adaptation during B. cinerea infection under sequential biotic and abiotic stresses.

Discussion

Plant responses to sequential double (a)biotic stresses

To identify genes that play a crucial role in the plant’s adaptive response to B. cinerea infection during sequential double stresses, we performed B. cinerea disease resistance bioassays with 346 A. thaliana accessions under single-stress and sequential double-stress conditions. Our screening readout for spreading B. cinerea lesions was obtained 3 days after pathogen inoculation, a so-called ‘endpoint’ measurement. Using this approach, one may expect to find genes that influence the outcome of the defense response in both early and late stages of the infection process. By inoculating the leaves with B. cinerea immediately after another stress (herbivory or drought stress), signals from the previous stresses are still present and can thus have an impact on the response of the secondary stress (Atkinson et al. 2013; Coolen et al. 2016; Davila Olivas et al. 2016). In a previous study focused on the transcriptional dynamics of B. cinerea-inoculated A. thaliana leaves during multi-stress conditions, we showed that both herbivory and drought stresses had an impact on B. cinerea-induced transcriptional reprogramming, especially within the first 12 h after B. cinerea inoculation (Coolen et al. 2016). Phenotypic analysis of the impact of prior herbivory or drought stress on the level of B. cinerea resistance in 346 A. thaliana accessions showed that a prior stress indeed significantly changed the level of B. cinerea resistance in many of the tested A. thaliana accessions (Fig. 2 and Suppl. Tabel S2). These findings suggest that there is natural genetic variation that can potentially be mined for yet unknown key players in the plants ability to adapt to changing environmental conditions, such as tested here with sequential double (a)biotic stresses in combination with B. cinerea.

Although timing and severity of stresses are expected to greatly influence the outcome of the plants adaptive response, a first stress encounter can leave a detectable and permanent imprint, e.g. at the level of gene transcription (Coolen et al. 2016), or at the level of second stress resistance. Prior herbivory, which in our experimental setup is often only a bite of less than a millimeter in size, inflicts major changes in the whole-plant transcriptional profile (Coolen et al. 2016). In the majority of the accessions, this prior herbivory resulted in a reduced susceptibility to B. cinerea (Fig. 2). This phenomenon was previously described to be the result of wound-induced resistance, and could partially be explained by increased and more rapid accumulation of camalexin (Chassot et al. 2008). These findings, together with the strong correlation that we observed between the level of resistance in the 346 A. thaliana accessions under conditions of single stress (Bc) and sequential double stress (Pr + Bc; Fig. 2b), suggest that the plant’s response towards P. rapae and B. cinerea make use of similar signaling components that act synergistically in mounting defense against B. cinerea. Both P. rapae and B. cinerea abundantly trigger JA-dependent defenses (Reymond et al. 2000; Windram et al. 2012; Coolen et al. 2016), making it likely that components of the JA response contribute to the synergistic effect of prior herbivory on the level of B. cinerea resistance in the majority of the A. thaliana accession.

For the drought sequential double stress, a moderate correlation between the level of B. cinerea resistance under single- and double-stress conditions was observed in the 346 A. thaliana accessions (Fig. 2b). This points to the contribution of accession-specific components on the effect of prior drought stress on the level of B. cinerea resistance. These accession-specific components may be related to water usage and drought-stress tolerance amongst the accessions (Suppl. Fig. S2). Drought is one of the most frequently experienced abiotic environmental stresses by plants. They are well able to sense a reduction in the availability of water in their rhizosphere and in response transmit a chemical signal to the shoots to reduce stomatal conduction and leaf growth (Schachtman and Goodger 2008). The plant hormone ABA is currently the best candidate for such a chemical signal (Leckie et al. 1998). ABA not only regulates the adaptive plant response to drought, it also modulates the relative strength of the MYC-branch of the JA response (predominantly effective against insect herbivory) and the ERF-branch of the JA response (predominantly effective against necrotrophic pathogens) (Verhage et al. 2011; Pieterse et al. 2012). Hence, differences in drought resistance amongst A. thaliana accessions may impact the level of B. cinerea resistance via components of the ABA signaling pathway.

Candidate genes for plant resistance to B. cinerea during sequential (a)biotic stresses

Using the phenotypic variation in plant resistance to B. cinerea during sequential double (a)biotic stresses, GWA mapping pinpointed genomic regions that are potentially responsible for the observed phenotypic variance amongst the accessions (Fig. 2). Using a combination of mapping approaches raises the probability that a real contributing factor can be identified, especially when traits are genetically complex (Nordborg and Weigel 2008). Therefore, to pinpoint truly interesting SNP-trait associations, we used a combined approach of GWA mapping, QTL mapping and fine mapping using full genome sequences. QTL mapping identified a number of map positions that were linked to the GWA mapping positions (Fig. 3). However, several GWA mapping positions, including ones associated with candidate genes previously implicated in B. cinerea resistance (e.g. LIK1 and LOX2), were not identified by QTL mapping. This can be explained by the fact that there is limited natural variation within these genes or within the natural genetic variation of the MAGIC lines (RILs of 19 A. thaliana accessions) used for QTL mapping, compared to the 346 A. thaliana accessions used for GWA mapping. Nevertheless, QTL mapping did result in several significant loci of which one QTL island for the P. rapae/B. cinerea sequential double-stress residual data mapped close to the GWA locus pointing to the primary auxin indole-3-acetic acid (IAA) methyltransferase gene IAMT1 (Fig. 3; Suppl. Fig. S6). IAMT is known to be involved in the regulation of plant development and auxin homeostasis through IAA methylation (Spiess et al. 2014). Although this locus could not be confirmed with fine mapping, both GWA mapping and QTL mapping supported involvement of this locus in the effect of P. rapae herbivory on B. cinerea resistance. Interestingly, IAMT1 was shown to be induced by P. rapae feeding, indicating its involvement in plant responses to herbivory (Coolen et al. 2016). Furthermore, IAA was previously implicated in the defense response of Nicotiana attenuata to Manduca sexta herbivory, where herbivory-induced IAA was shown to regulate a subset of systemic JA-dependent secondary metabolites involved in defense against herbivores (Machado et al. 2016).

Fine mapping of the 85 GWA-mapped loci resulted in the confirmation of nine loci with significantly SNP-trait associations. For B. cinerea as a single stress, three loci were identified pointing to candidate genes HSP60-2 or GRIM-16, a tetratricopeptide repeat (TPR)-like superfamily gene, an ankyrin repeat family gene, a RmlC-like cupin superfamily gene, and RS41 (Table 3). Of these genes, HSP60-2 has been shown to be stress inducible and is induced 24 h after B. cinerea inoculation (Windram et al. 2012; Coolen et al. 2016). We identified a difference in the HSP60-2 protein sequence amongst the accessions that is significantly linked to the natural variation in the level of B. cinerea resistance (Suppl. Fig. S8). Also for the ankyrin-repeat family protein gene (AT3G04140), which is repressed during the first 24 h after B. cinerea inoculation, an amino acid difference was found to be associated with natural variation for B. cinerea resistance. The miRNA biogenesis gene RS41 has previously been implicated in B. cinerea resistance as the double mutant rs40, rs41 displays reduced DCL-1 expression and the A. thaliana dcl1-7 mutant shows enhanced susceptibility to B. cinerea (Weiberg et al. 2013; Chen et al. 2015). A. thaliana DCL-1 processes hairpin mRNA into miRNA that targets specific mRNA transcripts for destruction. Since miRNAs were shown to be important for plant defense against B. cinerea, impairment of miRNA production is likely to result in enhanced disease susceptibility to B. cinerea (Jin and Wu 2015).

For the P. rapae/B. cinerea sequential double-stress condition, the cell wall biosynthesis gene CSLA01 emerged as a candidate affecting the level of B. cinerea resistance after prior exposure to herbivory (Table 3). This gene was previously found to be upregulated in response to P. rapae feeding (Coolen et al. 2016). Mutation of cellulose synthase genes CESA4, CESA7 and CESA8 was shown to enhance disease resistance against necrotrophic fungi such as P. cucumerina. Conversely, B. cinerea has been shown to downregulate the activity of plant cellulose synthases to increase virulence (Ramírez et al. 2014), highlighting the role of cellulose biosynthesis and breakdown in the interaction between A. thaliana and necrotrophs. Additionally, the polygalacturonase gene PG3 was found to be associated with a GWA- and fine mapping-verified locus linked to herbivory-mediated changes in the level of B. cinerea resistance. As a pectin-degrading enzyme PG3 is involved in many stages of plant development (Cao 2012) and since pectin degradation is also part of the infection strategy of necrotrophs, B. cinerea might make use of the plant’s PGs for plant invasion. Indeed, deletion of the B. cinerea PG2 gene was shown to result in delayed primary lesion formation and a strong reduction of virulence on tomato and broad bean (Kars et al. 2005), and inhibition of PG proteins in A. thaliana was shown to increase resistance to F. graminearum (Ferrari et al. 2012).

For the drought/B. cinerea sequential double-stress combination, four loci were confirmed with fine mapping pointing to the candidate genes BCA3, a plant invertase/pectin methylesterase inhibitor gene, KNAT6, an exostosin gene, and an F-box-associated ubiquitination effector family gene (Table 3). BCA3 was found to be induced at 24 h after inoculation in the drought/B. cinerea sequential double-stress combination, suggesting a possible involvement in drought-mediated plant responses to B. cinerea (Coolen et al. 2016). KNAT6 is known to be involved in plant flowering activated stress pathways that promote JA biosynthesis (Khan et al. 2015), potentially influencing JA-induced defenses against B. cinerea. The exostosin family protein gene (AT5G11610), which is induced from 18 h after B. cinerea inoculation onward, was repressed during the first 24 h after B. cinerea inoculation when the fungal infection was preceded by drought stress. The predicted protein sequence contains an amino acid difference that is associated with the natural phenotypic variation in B. cinerea resistance after prior drought, suggesting a role for this gene in B. cinerea resistance after encountering drought stress. Furthermore, the plant invertase/pectin methylesterase inhibitor was previously described in having an important role in B. cinerea resistance, as overexpression was shown to restrict infection by B. cinerea (Lionetti et al. 2007). Moderate drought stress, as we applied in this study, was previously shown to induce expression of plant invertase/pectin methylesterase inhibitors, potentially affecting B. cinerea in our setup (Clauw et al. 2015).

Besides the nine GWA loci that were confirmed by fine mapping, many of the other GWA loci showed a similar pattern of SNP-trait associations in the GWA and fine mapping approach, but because of the lack of power did not exceed the false discovery rate (FDR)-corrected threshold for significance (Fig. S7). Below, we will name a few that are worthwhile mentioning. For instance, RESISTANCE TO PSEUDOMONAS SYRINGAE PV MACULICOLA 1 (RPM1) INTERACTING PROTEIN 13 (RIN13, AT2G20310), which emerged from the drought/B. cinerea sequential double-stress data. RIN13 has been described to positively enhance RPM1 (NBS-LRR resistance protein)-mediated resistance to P. syringae while restricting the RPM1-induced hypersensitive response (HR) (Boyes et al. 1998; Al-Daoude et al. 2005). This could benefit B. cinerea as the hypersensitive response induced by P. syringae pv. tomato DC3000 (avrRPM1) was shown to promote B. cinerea disease development in A. thaliana (Govrin and Levine 2000), suggesting that RIN13 could be involved in HR-mediated disease development triggered by B. cinerea. Another candidate gene that emerged from the drought/B. cinerea sequential double-stress data is the putative TIR-NBS-LRR class disease resistance gene AT5G41750. Also for the other stress combinations, this locus displays several SNP-trait associations, albeit non-significant (Fig. S9). AT5G41750 was found to be induced by B. cinerea infection and shows similar expression patterns when preceded by drought or P. rapae feeding, suggesting that this gene is responsive to B. cinerea independent of the prior stresses (Coolen et al. 2016). Future research will be directed towards elucidating the roles of the identified candidate genes in the modulation of the defense response to B. cinerea.

Conclusions

In sum, this study shows that adaptive plant responses to a single stress (here B. cinerea infection) can be positively or negatively affected by prior exposure to another stress (here P. rapae herbivory or drought stress). We show that plants can have considerable natural genetic variation in their response to combinatorial stresses. By mining this natural genetic variation, we were able to pinpoint a number of candidate genes, some of which with known or expected effects on the level of B. cinerea resistance, and some with so far unknown functions in the combinatorial stress response. Our results show that the ability of plants to cope with multiple stresses at the same time is genetically tractable, which can be highly instrumental in future plant breeding programs that are directed to breed crops with multi-stress resistance traits.

Author contribution statement

SC, SCMVW and CMJP conceived and designed research. SC and JAVP conducted experiments. SC analyzed data. SC, SCMVW and CMJP wrote the manuscript. All authors read and approved the manuscript.