A complex network of additive and epistatic quantitative trait loci underlies natural variation of Arabidopsis thaliana quantitative disease resistance to Ralstonia solanacearum under heat stress

Abstract Plant immunity is often negatively impacted by heat stress. However, the underlying molecular mechanisms remain poorly characterized. Based on a genome‐wide association mapping approach, this study aims to identify in Arabidopsis thaliana the genetic bases of robust resistance mechanisms to the devastating pathogen Ralstonia solanacearum under heat stress. A local mapping population was phenotyped against the R. solanacearum GMI1000 strain at 27 and 30 °C. To obtain a precise description of the genetic architecture underlying natural variation of quantitative disease resistance (QDR), we applied a genome‐wide local score analysis. Alongside an extensive genetic variation found in this local population at both temperatures, we observed a playful dynamics of quantitative trait loci along the infection stages. In addition, a complex genetic network of interacting loci could be detected at 30 °C. As a first step to investigate the underlying molecular mechanisms, the atypical meiotic cyclin SOLO DANCERS gene was validated by a reverse genetic approach as involved in QDR to R. solanacearum at 30 °C. In the context of climate change, the complex genetic architecture underlying QDR under heat stress in a local mapping population revealed candidate genes with diverse molecular functions.

Unravelling the genetic and molecular mechanisms allowing plants to face pathogen attacks under elevated temperature therefore represents a promising strategy for sustainable disease resistance.
Being sessile, plants have developed a wide range of immune responses to face simultaneous and/or sequential stresses caused by various bioaggressors (Roux and Bergelson, 2016). Plant immunity relies on a surveillance system involving plasma membrane-anchored pattern recognition receptors (PRRs) that perceive microbial elicitors, called pathogen-or microbe-associated molecular patterns (PAMPs or MAMPs). PRR-triggered immunity (PTI) is efficient against a broad spectrum of pathogens (Cook et al., 2015). Adapted pathogens such as phytopathogenic bacteria trigger susceptibility thanks to secreted virulence factors called effectors that can inhibit PTI and promote pathogen invasion (effector-triggered susceptibility, ETS). The specific recognition of pathogen effectors by plant intracellular nod-like receptors (NLRs) triggers a more robust immune response called effector-triggered immunity (ETI), often associated with a cell death response or hypersensitive response (HR) that restricts pathogen invasion to the infection site. In general, ETI is specific to a single pathogenic species, and even to a single pathogenic strain. This specificity causes a strong selective pressure on virulent strains to bypass ETI, making in most cases this form of immunity not durable in crop field conditions (Roux et al., 2014).
Bacterial wilt, caused by the gram-negative bacterium Ralstonia solanacearum, is one of the most devastating bacterial diseases in the world. Indeed, this soilborne pathogen affects more than 200 species, including members of Solanaceae and Brassicaceae, and is responsible for dramatic yield losses not only in tropical and subtropical areas, but also in warm temperate regions (Elphinstone, 2005). In the model plant Arabidopsis thaliana, a broad-spectrum resistance response to R. solanacearum is conferred by the RPS4/RRS1-R locus that encodes a pair of NLR receptors cooperating molecularly to form homodimers (Deslandes et al., 2002;Birker et al., 2009;Narusaka et al., 2009;Williams et al., 2014). In addition, the leucine-rich repeat receptor-like kinase ERECTA was identified as underlying one of the three quantitative trait loci (QTLs) detected against the R. solanacearum strain 14.25 (Godiard et al., 2003).
The genetic architecture and the molecular mechanisms of plant responses to R. solanacearum in changing abiotic environments, and more particularly under elevated temperature conditions, remain elusive.
Recently, a genome-wide association study (GWAS) performed in A. thaliana and aimed at exploring the genetic bases associated with the natural variation of plant response to strain GMI1000 at 30°C led to the identification of the Strictosidine Synthase-Like protein 4 (SSL4) gene, although the underlying molecular mechanisms are still unknown (Aoun et al., 2017).
This study was based on 176 accessions of A. thaliana from a worldwide collection. While being informative, a limitation of this mapping population-based approach resides in an increased effect of the demographic history on genotype-phenotype association at large geographical scales.
Statistical methods controlling for confounding by population structure can reduce the rate of false-positive associations, but to the detriment of a loss of detection power (i.e., markers linked to causative genes that are lost after correcting for population structure; Roux, 2010, Brachi et al., 2010). In addition, because different QTLs and/or different alleles at the same QTL can be responsible for the same phenotypic values, the power of GWAS can be strongly reduced by the effects of genetic and allelic heterogeneity due to the increased probability of the presence of rare alleles at large geographical scales (Bergelson and Roux, 2010). To limit these drawbacks, GWA mapping can be combined with traditional linkage mapping (based on the use of experimental populations such as recombinant inbred lines, RILs), which is prone to identifying rare alleles and not subjected to the effect of population structure (Bergelson and Roux, 2010). Combining GWA mapping and traditional linkage mapping has been demonstrated to reduce the rates of false positives and negatives when applied to flowering time data in A. thaliana (Brachi et al., 2010), but remains time-consuming due to the need to phenotype thousands of experimental lines. To limit the drawbacks of GWA mapping performed at a worldwide scale, an alternative approach is to work at a small geographical scale (Bergelson and Roux, 2010). As reported in a GWAS performed on flowering in A. thaliana from a worldwide to a local scale (by using two highly polymorphic French mapping populations), a great reduction of confounding by population structure was observed at the smaller geographical scales (Brachi et al., 2013). In addition, the genetic architecture was highly specific on the considered geographical scale (Brachi et al., 2013).
In the present study, we therefore investigated the genetic bases of QDR to R. solanacearum under elevated temperature by performing a GWAS at a small geographical scale using the TOU-A local mapping population. This local mapping population offers several advantages, including (a) the detection of more than 1.9 million single nucleotide polymorphims (SNPs), only 5.6 times less than observed in a panel of 1,135 accessions collected at the worldwide scale (Frachon et al., 2017); (b) an extensive genetic variation for a large range of phenotypic traits, including QDR to the bacterial vascular pathogen Xanthomonas campestris pv. campestris; (c) a linkage disequilibrium (LD) decay below 3 kb allowing fine-mapping of genomic regions associated with phenotypic variation; (d) a strongly reduced confounding effect by population structure; and (e) an adaptation to local warming in fewer than eight generations (Brachi et al., 2013;Huard-Chauveau et al., 2013;Baron et al., 2015;Debieu et al., 2016;Frachon et al., 2017).
Interestingly, this work revealed a genetic architecture of natural variation of QDR to R. solanacearum that totally differs from the one previously described at the worldwide scale (Aoun et al., 2017).
In particular, at 30°C, we observed a playful dynamics of 12 QTLs, showing an increase or a decrease in significance, along the disease symptom progression, with most QTLs displaying complex epistatic relationships. Using a reverse genetic approach, we identified SOLO DANCERS (SDS) encoding a cyclin-like protein as the gene underlying one of the two additive QTLs detected at 30 °C.

| RE SULTS
2.1 | Impact of temperature on genetic variation for QDR to R. solanacearum among local A. thaliana accessions In this study, we tested 192 whole-genome sequenced local accessions of A. thaliana in response to the R. solanacearum GMI1000 reference strain, under growth chamber conditions. No germination was observed for six accessions that were therefore discarded from the study ( Table S1). The remaining 186 local accessions from the TOU-A population were challenged with GMI1000 at 27 and 30 °C by cutting the roots. The accessions were on average more susceptible at 30 °C than at 27 °C ( Figure 1a,b). For each temperature treatment, we observed a large genetic variation at most infection stages, that is, 5, 6, and 7 days after inoculation (dai) at 27 °C and 4, 5, 6, and 7 dai at

| Playful dynamics of QTLs at 27 and 30 °C
To increase the probability to discover QTLs with additive effects conferring QDR to R. solanacearum along the infection stages, we combined a genome-wide association mapping approach with a local score analysis (with tuning parameter ξ = 2) (Bonhomme et al., 2019).
In agreement with weak cross-temperature genetic correlation, no single top SNP was common to both temperatures, indicating a contrasted genetic architecture for natural variation of response to R. solanacearum GMI1000 between 27 and 30 °C. Next, we focused on the 14 most highly significant additive QTLs (i.e., top QTLs with a Lindley process >10, Figures 2 and 3 Table 2).

| SOLO DANCERS is the gene underlying QTL4 involved in QDR to R. solanacearum GMI1000 at 30 °C
Next we investigated the molecular mechanisms underlying plant response to R. solanacearum at 30 °C in the TOU-A population.
For this, we focused on the additive QTL with the highest allelic TA B L E 1 Natural variation among TOU-A natural accessions for disease index at 27 and 30°C   At 27 °C, the sds-2 mutant response was not significantly different from Col-0 except at 6 dai ( Figure 6a and Table S2). At 30 °C, the wilting of sds-2 was significantly delayed compared to that of Col-0 from 3 to 5 dai ( Figure 6b and Table S2). As expected, sds-3 mutant and Ws-4 plants remained symptomless at 27 °C, from 3 to 7 dai ( Figure 6c and Table S2). By contrast, at 30 °C, wilting of sds-3 plants was strongly reduced during all infection stages ( Figure 6d and Table S2). Altogether, these data suggest that SDS plays a role in wilting disease development on infection with the R. solanacearum GMI1000 strain at 30 °C.

| Worldwide versus local genetic variation in A. thaliana facing R. solanacearum
In comparison with a temperature of 27 °C, local Arabidopsis accessions exposed at 30 °C were on average more susceptible to The solid coloured curve indicates the Lindley process (local score method with ξ = 2) calculated from left to right (y axis on the right). The two coloured dashed vertical lines indicate the QTL intervals detected, without taking into account the right part of the curve (Fariello et al., 2017;Bonhomme et al., 2019) the GMI1000 strain as indicated with a faster wilting disease progression than that observed at 27 °C. This observation is in line with the drastic impact of heat stress on Arabidopsis response to the GMI1000 strain previously monitored in a worldwide collection of A. thaliana (Aoun et al., 2017). This is also in accordance with a growing number of studies performed on crops infected by impact of heat stress on host and pathogen is still a matter of debate.

Alteration of plant immunity under heat stress has been reported
in several studies. For instance, a temperature of 28 °C inhibits "spontaneous lesion" phenotypes or autoimmune responses linked to autoactive alleles of genes encoding NLR resistance proteins (Zhu et al., 2010;Negeri et al., 2013). While heat stress increases

A. thaliana susceptibility to infection with Pseudomonas syringae
pv. tomato DC3000, it also promotes the plant-dependent bacterial multiplication (Huot et al., 2017). In this study, R. solanacearum grows faster at 30 °C than at 28 °Cunder in vitro conditions ( Figure   S3a) and a significant increase in bacterial multiplication was observed in planta at 30 °C ( Figure S3b), suggesting an effect of heat stress on bacterial multiplication and its pathogenicity.
In previous studies, the level of genetic variation for diverse phenotypic traits such as flowering time and QDR to the bacterial pathogen X. campestris was similar between the TOU-A population and a set of worldwide accessions (Brachi et al., 2013;Huard-Chauveau et al., 2013;Debieu et al., 2016). Here, the level of phenotypic and genetic variation for QDR to R. solanacearum in the TOU-A popula-  Brachi et al., 2013;Frachon et al., 2018;Frachon et al., 2019), it would be interesting to investigate the level of genetic variation of QDR to R. solanacearum GMI1000 (or to other strains) within those populations. This would provide valuable information on the local dynamics of QDR to a bacterial pathogen in A. thaliana at the metapopulation level (Ding et al., 2007;Vetter et al., 2012;Karasov et al., 2014;Roux and Bergelson, 2016).

| Complex genetic architecture of QDR to R. solanacearum at 27 and 30 °C
Combining GWA mapping related mixed models and genome-wide local score analysis increases the probability of discovering minor QTLs with additive effects (Fariello et al., 2017;Bonhomme et al., 2019). This is also well exemplified here with a more detailed charac-

At2g30820
Aspartyl/glutamyl-tRNA(Asn/Gln) amidotransferase subunit No genomic regions containing top SNPs were shared between the local and worldwide scales (Data S1). Interestingly, similar results were obtained in a study investigating the genetic determinants of flowering time scored on both local and worldwide mapping populations of A. thaliana in two environmental conditions simulating two seasonal germination cohorts (Brachi et al., 2013). The genetic architecture of flowering was highly dependent on both the geographical scale and the considered season (Brachi et al., 2013). Together, the data obtained from various phenotypic traits reinforce the need to account for the geographical scale of phenotypic variation when choosing accession panels for GWAS (Bergelson and Roux, 2010).

2-Oxoglutarate (2OG) and Fe(II)-dependent oxygenase superfamily protein
Consequently, this would help to get a better view of the genetic architecture flexibility of phenotypic traits.
Theoretical predictions suggest that phenotypic changes in ontogenetic time (typically time-to-event or time-to-failure traits such as flowering time or death time) are often driven by the temporal regulation of QTLs (Johannes, 2007 Interestingly, it turns out that transcriptional responses of plants to combined stresses are unique and cannot be predicted from that of individual stress (Atkinson and Urwin, 2012;Suzuki et al., 2014;Pandey et al., 2015). In addition, combined stresses induce a major transcriptional reprogramming characterized by the regulation of the expression of a greater number of genes than observed with individual stresses (Rasmussen et al., 2013;Suzuki et al., 2014). For instance, more Arabidopsis genes were differentially expressed when nematode infection was combined with water stress compared to plants only subjected to nematode infection (Atkinson et al., 2013). Similar observations were made in other pathosystems, including A. thaliana exposed to the simultaneous application of virus and heat or virus and drought (Prasch and Sonnewald, 2013;Pandey et al., 2015). This may indicate that the more complex the environment is, the more the plants establish a response with a polygenetic architecture involving different genetic pathways.
Epistatic networks involving long-distance LD among physically unlinked loci were reported to represent the main fraction of phenotypic variance for herbivore resistance in A. thaliana at a worldwide scale (Brachi et al., 2015), yeast growth (Forsberg et al., 2017), or body weight and abdominal fat content in chicken (Carlborg et al., 2006;Li et al., 2013). While complex epistatic relationships among QTLs, including higher-order epistasis, may be therefore more frequent than anticipated (Carlborg and Haley, 2004;Roux et al., 2005;Pettersson et al., 2011), the functional validation of epistatic QTLs remains challenging but feasible if we consider a multi-CRISPR-Cas9 approach to create double, triple, quadruple, etc. mutants. Nonetheless, it would be interesting to determine whether such an epistatic network is Least-square means ± SE of the LS means from three independent inoculations (n = 72 plants per genetic line × temperature combination. Symbols *, **, and *** denote significant difference observed between each wild-type background and its corresponding mutant at p < .05, p < .01, and p <.001, respectively. Coloured stars indicate significant differences after a false-discovery rate correction. DI, disease index; dai, days after inoculation restricted to the TOU-A population by estimating LD between these five QTLs in other local highly polymorphic populations or at a larger geographical scale. However, we should stress that some limitations in this study preclude a full description and understanding of the epistatic network underlying QDR to R. solanacearum. First, in contrast to traditional mapping populations such as F2 populations or RILs, the number of accessions among the haplotypes expected between two (or more) epistatic QTLs was clearly unbalanced in the TOU-A population. While it may reflect the maintenance of haplotypes with extreme phenotypes (i.e., susceptible versus resistant) by selective processes (Brachi et al., 2015), it precludes testing with sufficient power the magnitude and type of epistasis. Secondly, epistatic relationships were only tested on QTLs with additive effects that were first identified by combining a GWA mapping approach with a local score analysis. Although computationally intensive, a complementary step will be to test the significance of all pairwise interactions among the 981,617 SNPs used in this study, which will require controlling the individual and joint effect of population structure on both SNPs tested in interaction (Wang et al., 2018).

| Various molecular functions are involved in QDR to the GMI1000 strain at 30°C
Consistent with the molecular functions of previously cloned QDR genes (Roux et al., 2014), the nature of most candidate genes under- This is consistent with previous data showing that the production of this metabolite by the R. solanacearum effector RipTPS plays an important role in pathogen virulence (Poueymiro et al., 2014).
For QTL2, the top SNPs fall in the At2g19050 gene encoding a GDSL-like lipase/acylhydrolase superfamily protein. Interestingly, overexpression of GLIP1 that also belongs to the Arabidopsis GDSL LIPASE-LIKE gene family was shown to confer enhanced resistance to several pathogens, including Alternaria brassicicola, Erwinia carotovora and P. syringae) (Kwon et al., 2009). Therefore, these proteins might also play a role in plant immunity against R. solanacearum.
The molecular functions of the candidate genes underlying the 12 major QTLs detected at 30 °C are even more diverse.
Interestingly, these functions may reflect different plant responses to face virulence strategies developed by the bacteria to colonize plant tissues and promote its multiplication within the xylem vessels. For instance, candidate genes underlying QTL7A and QTL10 are involved in the synthesis or signalling of hormones that may contribute positively to pathogen resistance and in plant response to combined biotic and abiotic stress. In particular, JA is known to interfere with GA signalling through the degradation of transcriptional repressors such as JAZ9 (the candidate gene underlying QTL7A) to balance plant defence response and growth (Yang et al., 2012).
From 4 to 7 dai, QTL4 was detected with increasing significance on chromosome I. The corresponding candidate gene, SDS, encodes an atypical meiotic cyclin-like protein related to A-and B-type cyclins, previously described as being required for DNA double-strand break (DSB) repair (Azumi et al., 2002;De Muyt et al., 2009). To our knowledge, SDS has never been associated with plant disease susceptibility. Interestingly, two allelic null sds mutants were found to be more resistant at both 27 and 30 °C to GMI1000, albeit the allelic effect was different between the two genetic backgrounds.
The functional validation of SDS as a susceptibility gene represents the first demonstration of its involvement in plant defence response to a bacterial pathogen under heat stress. It is noteworthy that (a) SDS acts together with CYCB3;1 in suppressing unscheduled cell wall synthesis (Bulankova et al., 2013); and (b) the two candidate genes underlying QTL12A and QTL12B encode, respectively, proteins involved in cell wall and lignin polymerization. Two cyclin-L type proteins, MOS12 (Modifier of SNC1, 12) and MOS4-associated complex (Modifier of SNC1, 4), have also been shown to participate in the alternative splicing of SNC1 and RPS4 genes, thereby enabling the fine-tuning of NLR gene expression (Xu et al., 2012). As several NLR genes have been described to be alternatively spliced without knowing the regulatory mechanism (Xu et al., 2012), it is tempting to hypothesize that SDS would participate in the regulation of NLR functions under combined R. solanacearum and elevated temperature conditions through the production of splicing variants. Because the top SNPs are located in the promoter region of SDS, the next step to decipher the underlying molecular mechanisms would be to investigate the natural variation of SDS expression in the TOU-A population and its link to the QDR.

| Bacterial strain, plant material, and growth conditions
The wild-type R. solanacearum GMI1000 strain was grown on complete bacto glucose (BG) agar as previously described (Plener et al., 2010

| Plant inoculation and phenotyping
Four-week-old plants were root-inoculated with the R. solanacearum GMI1000. The R. solanacearum GMI1000 reference strain was grown in complete BG agar, supplemented with 6 ml of glucose (20%) and 1 ml of triphenyltetrazolium chloride (1%), and incubated at 28 °C for 48 hr then left at room temperature for 24 hr. One day before inoculation, one colony was grown in liquid BG medium and grown overnight at 28 °C under shaking. Plants were inoculated with a bacterial suspension at OD 600 nm between 0.8 and 1.
Before inoculation, roots were cut with scissors 1 cm from the bottom of the Jiffy pot (Deslandes et al., 1998). This method gives the bacteria direct access to the xylem vessels. During inoculation, plants were soaked in a bacterial suspension at 10 7 bacteria/ml for 15 min. Inoculated plants were incubated in growth chambers at 27 or 30 °C (75% RH, 12 hr light, 100 μmol⋅m −2 ⋅s −1 ). Disease symptoms were scored daily from 3 to 9 dai using a disease index scale from 0 to 4 as previously described (Deslandes et al., 1998) with the scores from 0 to 4 corresponding to healthy and fully wilted plants, respectively.  where µ is the overall mean of the phenotypic data, "block" accounts for differences in microenvironmental conditions between the three experimental blocks, "accession" corresponds to the genetic differences among the TOU-A natural accessions, covCol is a covariate accounting for tray effects within blocks (phenotypic mean of the three Col-0 replicates per tray was used as a covariate), and "ε" is the residual term. The factor "block" was considered as a fixed factor and the factor "accession" as a random factor. The significance of the random effect was determined by likelihood ratio tests of model with and without this effect. Residuals were normally distributed so no transformation was applied on raw phenotypic data. For GWA mapping analyses, we used best linear unbiased predictors (BLUPs) obtained for each natural accession. Because A. thaliana is a highly selfing species (Platt et al., 2010), BLUPs correspond to genotypic values. Using a formula adapted from Gallais (1990), broad-sense heritabilities (H 2 ) at each time point of phenotyping were estimated from the variance component estimates of the "block" and "accession" terms obtained with the VARCOMP procedure in SAS v. 9.4 (SAS Institute Inc.).

| GWA mapping with local score analysis
To fine map the genomic regions with additive effects associated with natural disease index variation at each time of phenotyping for each temperature treatment, a mixed model implemented in the software EMMAX was adopted (Efficient Mixed-Model Association eXpedited; Kang et al., 2010). To control for the effect of population structure, we included as a covariate an identity-by-state kinship matrix K based on the 1,902,592 SNPs identified in the TOU-A population (Frachon et al., 2017). Because rare alleles increase the rate of false positives when included in mixed models, we considered a threshold of minor allele relative frequency (MARF) >7% and ended up with 981,617 SNPs (Brachi et al., 2010;Kang et al., 2010). As previously described in Frachon et al. (2017), a threshold of 7% corresponds to the MARF value above which the p value distribution obtained from the mixed model is not dependent on MARF values in the TOU-A population.
Thereafter, we implemented a local score approach on the set of p values provided by EMMAX. The local score allows detection of significant genomic segments by accumulating the statistical signals from contiguous markers such as SNPs (Fariello et al., 2017). In a given QTL region, the association signal, through the p values, will cumulate locally due to LD between SNPs, which will then increase the local score (Bonhomme et al., 2019). Briefly, a sequence of scores is calculated along the chromosome as X i = − log 10 (p i ) − ξ, where p i is the p value of marker i and ξ is a tuning parameter with an optimal value that can be fixed at 2 or 3 in a GWAS context (Bonhomme et al., 2019). Then, finding segments that accumulate strong signals is equivalent to finding peaks along a Lindley process defined as h i = max(0, h i−1 + X i ) along the chromosome, with h 0 = 0. Significant SNP-phenotype associations were identified by estimating a chromosome-wide significance threshold for each chromosome (Bonhomme et al., 2019). (1) disease index = + block i + accession j + covCol c +

| Detecting QTL epistasis
In order to detect epistatic interactions among our set of 14 top QTLs identified with additive effects by combining a GWA mapping approach with a local score analysis, we followed the procedure adopted in Brachi et al. (2015). We first identified within each QTL region the SNP with the highest association score estimated by EMMAX, hereafter named bait top SNP. For each of the 14 QTLs, we then computed LD estimates between the bait top SNP and all the other SNPs in the TOU-A population. LD between two biallelic (homozygous) SNPs was calculated using the absolute value of r statistic (correlation coeffi-

| Estimates of allelic effect
To display the allelic effect of the bait SNPs after controlling for the effects of population structure, BLUPs estimated from model (1) were adjusted by fitting them with a kinship matrix. Kinshipadjusted BLUPs were computed under the R environment 3.6.1 (R_Core_Team, 2019). In order to avoid pseudoreplication due to the presence of SNPs in stretches of LD, we first pruned the SNP data set with the snpgdspLD command using the following parameters: ld.threshold = 0.8, slide.max.bp = 500, maf = 0.07 ("gdsfmt" and "SNPRelate" packages), leaving 365,952 SNPs for the estimation of the kinship matrix. The kinship matrix was then estimated using the popkin function (allowing missing data in the SNP matrix) in the popkin package, with the subpopulation vector set to NULL. Because the resulting matrix was not positive semi-definite, the function make.positive.definite() from the package lqmm was used. Finally, the kinship-adjusted BLUPs were calculated with the function kin.
blup from package rrBLUP. Keeping the notations from model (1), the parameters were: accession as geno, disease index as pheno, the above-mentioned kinship matrix as K, GAUSS = F, indicating that the genotypes are not independent and follow G = K V G , block as fixed effect and covCol as covariate. The kinship-adjusted BLUPs were then extracted using the command $pred.

| Analysis of the SDS candidate gene
For each temperature treatment, an experiment of 288 plants was set up according to a RCBD of three temporal experimental blocks.
Each block was represented by one tray of 96 positions, corresponding to 24 replicates of each genotype, that is, the sds-2 and sds-3 mutants with their corresponding wild-type background Col-0 and Ws-4, respectively.
For each temperature treatment, we tested whether each null mutant differs from its corresponding wild-type background along the infection stages by using the following model (MIXED procedure in SASv. 9.4; SAS Institute Inc.): where µ is the overall mean of the phenotypic data, "block" accounts for differences in microenvironmental conditions between the three experimental blocks, "genotype" corresponds to the genetic differences between the T-DNA mutant and its corresponding wild-type background, "block × genotype" accounts for variation between genotype differences among blocks, and "ε" is the residual term. All factors were considered as fixed.

ACK N OWLED G EM ENTS
We are grateful to Raphaël Mercier for providing T-DNA insertion mu-

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.