Identification of quantitative trait loci controlling gene expression during the innate immunity response of soybean.

Microbe-associated molecular pattern-triggered immunity (MTI) is an important component of the plant innate immunity response to invading pathogens. However, most of our knowledge of MTI comes from studies of model systems with relatively little work done with crop plants. In this work, we report on variation in both the microbe-associated molecular pattern-triggered oxidative burst and gene expression across four soybean (Glycine max) genotypes. Variation in MTI correlated with the level of pathogen resistance for each genotype. A quantitative trait locus analysis on these traits identified four loci that appeared to regulate gene expression during MTI in soybean. Likewise, we observed that both MTI variation and pathogen resistance were quantitatively inherited. The approach utilized in this study may have utility for identifying key resistance loci useful for developing improved soybean cultivars.

Plants, like animals, are constantly under threat from different pathogens. To counter this threat, plants have developed a sophisticated system to detect pathogens and trigger a strong defense response (Jones and Dangl, 2006;Boller and Felix, 2009;Segonzac and Zipfel, 2011). The plant immune system is generally considered to have two separate, but interacting components (Jones and Dangl, 2006;Zipfel and Robatzek, 2010). In one component, using systems analogous to those used in animals, the plant recognizes the invading pathogen through detection of conserved structural motifs, termed microbe-or pathogen-associated molecular patterns (MAMPs or PAMPs, respectively). MAMPs are conserved molecules that may be essential for the pathogen's life cycle. Well-studied examples include lipopolysaccharides, chitin, flagellin (flg), and translation elongation factor Tu (EF-Tu; Boller and Felix, 2009;Segonzac and Zipfel, 2011;Thomma et al., 2011). MAMPs are recognized at the plant cell surface via pattern recognition receptors (PRRs). Once MAMPs are detected, a defense response, termed MAMP-triggered immunity (MTI), is mounted that leads to plant resistance (Zipfel, 2008;Katagiri and Tsuda, 2010). MTI generally produces broad resistance against a variety of nonadapted pathogens; however it is usually weak and can be overcome by the adapted pathogen. Perhaps due to this fact, plants have also developed a second means to recognize pathogens that results in stronger but more specific resistance. In this case, the plant recognizes specific effector proteins produced by the pathogen and often directly inserted into the plant's cytoplasm (Göhre and Robatzek, 2008;Katagiri and Tsuda, 2010). Recognition of the effector triggers this second line of plant defense against pathogens, hence the name effectortriggered immunity (ETI; Göhre and Robatzek, 2008;Segonzac and Zipfel, 2011). ETI is often associated with an accelerated and amplified MTI response and a hypersensitive cell death response at the infection site (Jones and Dangl, 2006). However, since effector proteins are species, race, pathotype, or strain specific and often require a specific resistance protein for their recognition by the host, ETI is quite specific, often protecting plants from only specific genotypes of any given pathogen species (Segonzac and Zipfel, 2011).
MTI is characterized by a variety of plant responses, which include changes in ion flux across the plasma membrane, changes in cytoplasmic calcium levels, production of reactive oxygen species (ROS), callose deposition, modification of phytohormone concentrations, and induction or repression of different plant genes that condition antibiotic molecules, such as the phytoalexin glyceollin produced by soybean (Glycine max; Zipfel, 2009;Boudsocq et al., 2010;Lygin et al., 2010;Mersmann et al., 2010;Luna et al., 2011). These responses can be detected within minutes to hours after MAMP treatment, and the magnitude of the response may be plant-species dependent (Segonzac and Zipfel, 2011). The MTI response can be regulated at the transcriptional, posttranscriptional, and posttranslational levels through chromatin modification, transcription factors, noncoding RNAs, and changes in protein glycosylation patterns (Alvarez et al., 2010;Zhang et al., 2011). We mention these examples to emphasize the complexity of the MTI response.
Several reports have demonstrated the ability of MTI to protect plants from a variety of bacterial and fungal pathogens (Gómez-Gómez et al., 1999;Segonzac and Zipfel, 2011). MTI is usually expressed as partial resistance and segregates as a quantitative trait that is inherited in an oligogenic or a polygenic manner with genetic effects being largely additive. In contrast, specific resistance (ETI) typically confers complete resistance and is usually pathotype specific and conditioned by single dominant resistance genes. However, ETI can be overcome by adaptation of the pathogen population through selection toward pathotypes that produce effectors not recognized by this system. In comparison to narrowly effective ETI-based resistance, broadly effective MTI-based resistance may be more durable with a broader spectrum of effectiveness against multiple pathogens or pathotypes of a pathogen. Indeed, in some host-parasite systems, such as Sclerotinia stem rot (Sclerotinia sclerotiorum) disease on soybean, partial resistance may be the only type of resistance available since no effective ETI system is known (Poland et al., 2008;Vuong et al., 2008). In a few cases, a direct link has been established between MTI and specific quantitative trait loci (QTLs) for partial resistance. For example, two major effect QTLs associated with resistance to Pseudomonas syringae pv phaseolicola were reported in Arabidopsis (Arabidopsis thaliana), one of which is flagellin-sensitive 2 (FLS2) that encodes a PRR receptor (Forsyth et al., 2010;Ahmad et al., 2011).
Soybean is a chief source of protein and oil for human and animal consumption, and is grown on about 6% of the world's arable land (Hartman et al., 2011). Soybean production in the United States increased from 75 million metric tons in 2000 to 91 million metric tons in 2010 with a concomitant increase in the value of the crop from $13 to $39 billion, respectively (USDA, 2010). Data compiled for the 2010 growing season indicates that the total yield loss from a variety of diseases totaled approximately 7% of U.S. (http://aes.missouri.edu/delta/research/ soyloss.stm) production, which equates to approximately $5 billion.
The MTI response has been widely studied, mainly in the model plant Arabidopsis and a few crops, such as tomato (Solanum lycopersicum) and rice (Oryza sativa). Few studies have examined MTI in soybean. There is significant potential in harnessing soybean MTI to improve cultivars to withstand a variety of pathogens and, therefore, enhance soybean production. In this study, we evaluated the MTI responses in different soybean genotypes by challenging with two known MAMPs: (1) the conserved 22-amino acid peptide from bacterial flg22, representing bacterial pathogens, and (2) crab shell chitin, representing fungal pathogens that contain chitin in their cell walls. We found clear variation between soybean genotypes in their MAMP-triggered oxidative burst, as well as in the expression of MAMP-responsive genes. Based on this analysis and the availability of a well-characterized population of soybean lines, we chose two parental soybean genotypes, one with a strong MTI response (LD00-2817P; Diers et al., 2010) and one with a weak MTI response (LDX01-1-65; Diers et al., 2005). These two genotypes also show significant differences in pathogen susceptibility when infected after MAMP treatment. Analysis of the population of F3 lines from the cross between the two genotypes identified a single QTL associated with MAMP-triggered oxidative burst. Similarly, the analysis of gene expression across the F3 lines led to the identification of three QTLs associated with MAMP-triggered gene expression (i.e. expression QTLs [eQTLs]). These QTLs associated with MTI expression, as well as the general approach utilized in this study, may be useful in the genetic improvement of biotic stress resistance in soybean.

Variation in the Oxidative Burst across Soybean Genotypes
The oxidative burst triggered by applications of flg22, chitin, or a mixture of both MAMPs was measured in leaf discs from four different genotypes, originally chosen due to their use as parents of recombinant populations previously genotyped with molecular markers (Kim et al., 2011). This analysis showed significant differences in MAMP-triggered oxidative burst among the four genotypes (Fig. 1, A-C). The data also showed that the genotypes LD00-2817P and Ripley gave the strongest oxidative burst and the genotypes LDX01-1-65 and EF59 gave the weakest responses, which was also observed when the total ROS produced during the oxidative burst (hereafter called as total ROS production) was compared (Fig. 1, D-F). A pairwise comparison of the total ROS produced revealed that LDX01-1-65 treated with MAMP mixture (flg22 + chitin) produced significantly more ROS than EF59; no significant difference in the total ROS was observed between LD00-2817P and Ripley with this treatment, even though Ripley showed a stronger oxidative burst than LD00-2817P over the first 10 min of treatment ( Fig. 1, C-F). Likewise, LD00-2817P and Ripley showed a similar oxidative burst and total ROS production when treated with either flg22 or chitin alone (Fig. 1, A and B). However, the oxidative burst and total ROS production triggered by chitin across genotypes, except for EF59, was significantly lower than that triggered by the mixture of MAMPs (Fig. 1, B and C). These results clearly indicate significant genotypic variation in the MTI response among soybean genotypes, which suggests that QTLs associated with this trait may be identified after further genetic analyses.

Development of an Affymetrix Soybean Whole-Genome Transcript Array to Evaluate Genome-Wide Gene Expression Changes during the MTI Response
The first-generation Affymetrix Soybean Genome array was derived from a collection of over 350,000 ESTs and was predicted to allow transcriptional analysis of 35,611 distinct genes, or approximately 50% of the soybean genome (http://media.affymetrix.com/ support/technical/datasheets/soybean_datasheet.pdf). With publication of the full soybean genome sequence (Schmutz et al., 2010), it became apparent that the first-generation Affymetrix soybean array does not represent all genes in the soybean genome, and that unambiguous identification of duplicated homologous genes was not possible in every case using this array (Libault et al., 2010).
To overcome these issues and investigate the expression of all genes in the soybean genome, a custom Affymetrix soybean whole-genome transcript (sense orientation) array (termed here Soybean WT array) was developed. The design of this array paid special attention to the ability to distinguish soybean paralogs such that the Soybean WT array allows quantification of the expression of every exon from all 66,000 predicted soybean transcripts (Schmutz et al., 2010). The specific design of the array also allows for the detection of (1) transcripts with undefined 3# ends, (2) nonpolyadenylated mRNAs, (3) truncated transcripts, (4) alternative polyadenylation sites, and (5) alternative splicing events. A full description of the array design can be found in Supplemental Materials and Methods S1 and Supplemental Figure S5. To demonstrate the utility of this array, it was utilized to identify those soybean genes whose expression responded to MAMP treatment (see below and Supplemental Results S1).

Identification of MAMP-Responsive Soybean Genes
To identify MAMP-responsive genes the transcriptional profile of four soybean genotypes was analyzed after 30 min of MAMP treatment (this time point was chosen since it gave the strongest response to MAMP treatment in Arabidopsis; Navarro et al., 2004). We hybridized fragmented cDNA from MAMP-or mocktreated leaves to the Affymetrix Soybean WT array (24 GeneChip hybridization = four genotypes 3 two treatments [MAMP or mock] 3 three replicates). We used a linear mixed model (false discovery rate [FDR] , 0.05) with an additional cutoff of a 2-fold ratio in pairwise comparisons (i.e. LD00-2817P Treatment versus LD00-2817P Control ) to identify genes differentially expressed in each genotype. We observed considerable variation in the transcriptional profiles of the different genotypes ( Fig. 2; Supplemental Tables S3-S5). On average, we detected 986 genes (both up-regulated and downregulated) per genotype, with a minimum of 586 (LDX01-1-65) and a maximum of 1,733 (LD00-2817P; Fig. 2). Across all four genotypes, we detected a total of 4,249 differentially regulated genes. Only 437 common , 50 mg/mL chitin (B and E), or a mixture of 1 mM flg22 + 50 mg/mL chitin (C and F) in LD00-2817P (LD), LDX01-1-65 (LX), Ripley (Rip), and EF59 (EF) leaf discs measured in relative luminescence units (RLUs). A to C, Time kinetics of the MAMPtriggered oxidative burst. Each data point represents the average of five biological replicates with three technical replicates. Error bars represent 6 SE of the average. D to F, Total ROS produced during the oxidative burst over 30 min of MAMP treatment. Results are the average 6 SE of the ratio between MAMP/mock from five biological replicates with three technical replicates for each. Within each diagram, values sharing the common letters are not significantly different at P # 0.05 by ANOVA.
To confirm these DNA microarray results, the expressions of 40 randomly selected genes were analyzed via quantitative reverse transcription (qRT)-PCR (Fig. 3). The patterns of expression obtained using the Soybean WT GeneChip was confirmed for 35 genes (Fig. 3). This corresponds to about 90% validation of the Soybean WT GeneChip results by qRT-PCR.
Unique genes differentially expressed in each genotype were functionally classified using the MapMan gene functional classification system (Thimm et al., 2004;Usadel et al., 2009). Across all four genotypes, six functional categories were overrepresented: regulation, protein modification, regulation of transcription, hormones, enzyme families, and transport (Fig. 4). Likewise, most of the gene expression variation was detected in these functional categories (Fig. 4). Interestingly, we observed that most of the genes belonging to the regulation category encode different kinds of receptor-like kinases, including different Leu-rich repeat (LRR) and LysM receptors and calcium-dependent protein kinases (Supplemental Fig. S1). Likewise, the analysis revealed that LD00-2817P, Ripley, and EF59 express several genes involved in auxin (indole-3-acetic acid), abscisic acid, brassinosteroid (BR), and ethylene synthesis, whereas LDX01-1-65 expresses genes involved only in indole-3-acetic acid and BR production ( Fig. 4). Together, this transcriptional analysis is consistent with the differential response of these four soybean genotypes to MAMP treatment.

MTI Can Induce Pathogen Resistance in Soybean
MAMP induction of gene expression may not always be translated into pathogen resistance. Therefore, it was important to show that the observed genotypic differences reflected measures of pathogen susceptibility. To address this issue, the mapping parents LD00-2817P and LDX01-1-65, crossed originally to develop a recombinant inbred population to map QTLs controlling soybean cyst nematode resistance (Heterodera glycines Ichinohe; Kim et al., 2011), which show a clear divergent MTI response (in this work), were selected for further analysis. Detached leaves from each of these two soybean-mapping parents were treated for 24 h with flg22, chitin, or a combination of both. After incubation, individual leaves were challenged separately with either P. syringae pv glycenea (Psg) or S. sclerotiorum. Three days after inoculation with Psg, no significant differences (FDR , 0.05) in colony-forming units (CFUs) were found between the two mock-treated parents, consistent with preliminary experiments in whole plants in the absence of any MAMP treatment (   . qRT-PCR validation of the MAMP-responsive genes in four contrasting soybean genotypes treated over 30 min with MAMPs (flg22 + chitin). A total of 40 randomly selected genes were used for qRT-PCR validation. Log 2 fold change values (MAMP/mock) from the qRT-PCR data were plotted against log 2 (MAMP/mock) hybridization intensity ratio values from the Soybean WT array (Array). r, Pearson correlation coefficient. Data are the average from three biological replicates with consistent results in each one. qPCR and Soybean WT values from the 40 selected genes are provided in the Supplemental Table S12. log CFUs values were reduced 30% upon MAMP treatment of genotype LD00-2817P, while LDX01-1-65 plants showed a reduction of 7% (Fig. 5, A and C). Interestingly, we observed that flg22 treatment did not reduce colonization of Psg on LDX01-1-65 leaves (Fig.  5C). Somewhat similar results were found when plants were challenged with S. sclerotiorum. For example, flg22 pretreatment resulted in a significant decrease of S. sclerotiorum colonization in LD00-2817P, but not in LDX1-1-65 (Fig. 5, B and D). In contrast, chitin and the mixture of both MAMPs resulted in a decrease of S. sclerotiorum colonization in LDX01-1-65 but not LD00-2817P. Further statistical analysis on these data revealed that the difference between these two parents was statistically significant (FDR , 0.05), suggesting that interaction between genotype and treatment was relevant for soybean pathogen resistance triggered by MTI. Together, these results indicate that there was variable response by the two soybean genotypes to the MAMP treatments, and in general, the parent showing the stronger MTI response (i.e. LD00-2817P) exhibited stronger Psg resistance after MAMP treatment.

Soybean MTI Segregates as a Multigenic, Quantitative Trait
To identify specific loci controlling the soybean MTI response, we analyzed a population of 97 F3 lines developed from a cross between LD00-2817P and LDX01-1-65, which gave divergent MTI response based on both ROS and gene expression analysis (see above). In the first experiment, we analyzed the individual MTI response across the 97 lines by quantifying the total ROS production. In a second experiment, we analyzed MAMP-triggered gene expression across the population for 25 selected genes that showed higher expression in MAMP-treated LD00-2817P when compared to MAMP-treated LDX01-1-65 (see "Materials and Methods"). Levels of total ROS production among F3 lines were normally distributed (P , 0.05); whereas, the expression of only some of the MAMP-responsive genes showed a normal distribution.
QTL mapping was done by first forming a linkage map with the 390 single-nucleotide polymorphism (SNP) markers that were polymorphic between the parents LD00-2817P and LDX01-1-65. Using this map, QTLs associated with MAMP-induced ROS production and MAMP-induced gene expression were positioned through interval mapping and a logarithm of odds (LOD) threshold corresponding to an experiment-wide threshold of 0.05 defined for each phenotype. One QTL for total ROS production was mapped and the LOD peak was at the SNP marker ss107920734 on chromosome 3 (Fig. 6A). This QTL was significant at a LOD of 2.83 and explained 12.6% of the phenotypic variation for ROS production. The allele from LD00-2817P conferred greater ROS production. Interestingly, genes encoding NADPH oxidase respiratory burst oxidase homolog (RBOH; main ROS producer), MAP kinase kinase, and enzymes involved in the biosynthesis of secondary metabolites (Supplemental Table S14) were localized within 1 cM of this marker. eQTLs may act either in cis, when mapped within 5 cM of the gene used for the expression analysis, or in trans, when mapped to a chromosome different from the gene used to measure expression (Swanson-Wagner et al., 2009). Also, it is possible that the expression of a gene may be controlled by multiple QTLs (Potokina et al., 2008;Swanson-Wagner et al., 2009). According to these criteria, two cis-eQTLs and a third trans-eQTL were mapped in the LD00-2817P by LDX01-1-65 population. From the expression of the gene Glyma01g43420 (localized on chromosome 1; WRKY 12 transcription factor), one trans-eQTL was mapped near marker ss107931019 on chromosome 13 and explained 12.8% of the phenotypic variation, and had a LOD of 2.87 (Fig. 6D). The allele from LD00-2817P was associated with low expression levels of this MTI-responsive gene. One of the cis-eQTLs, from the expression of the gene Glyma11g19650 (Argonaute5), mapped near ss107919127 on chromosome 11 (Fig. 6C) and explained 16% of the gene expression variation, and had a LOD of 4.12. Similar to the trans-eQTL, the allele from LD00-2817P was associated with lower gene expression. Interestingly, analysis of the expression of the gene Glyma04g28420 (S-receptor kinase-like protein) identified a major effect cis-eQTL region on chromosome 4 and close to marker ss107928610 (Fig. 6B) and explained 49.4% of the phenotypic variation, and had a LOD of 14.19 (P , 0.05). The allele from LD00-2817P was associated with a higher level of gene expression, as compared to the response of the LDX01-1-65 allele. Most of the genes localized in the trans-eQTL and in the major effect cis-eQTL (LOD 14.19) encode for different receptor-like kinases (including LRR receptors), transcription factors (MYB, WRKY, and AP2-EREBP), protein modification enzymes, and transporters (zinc and calcium transporters). The cis-eQTL, in contrast, localized to the chromosome 11 encodes for either transporters (phosphorus and nitrate) or enzymes involved in secondary metabolism. However, all QTLs identified in this work map to regions where genes involved in the biosynthesis of hormones are located (Supplemental Table S14).

Broad Pathogen Resistance Triggered by MTI Is Inherited
Individual F3 lines were tested for their response to either Psg or S. sclerotiorum infection to determine whether the response of the two parental lines to pathogen attack is heritable. We selected two progeny lines that had a weak MTI response (based on the expression of the 25 MAMP-responsive genes used for the QTL analysis; Supplemental Fig. S3; Supplemental Table S16) and two other lines that had a strong MTI response. Representative data for one set of paired lines (one line with weak MTI and other with strong MTI) is shown in Figure 7, and information about the second pair of lines is shown in the Supplemental Figure S4. A significantly lower level of Psg colonization after MAMP treatment was found in the lines with the stronger MTI response (Fig. 7A). In contrast, in the lines selected based on their lower MAMP-triggered gene expression, no significant differences were found between the MAMP and mock treatments, indicating a weak MTI response (Fig. 7C). In the case of S. sclerotiorum-infected plants, the lines with a stronger MTI response had significantly less infection than the parental lines, about 50% less (Figs. 7B and 5,B and D). In contrast, no significant differences were observed after S. sclerotiorum inoculation in the line showing a weak MTI based on MAMP-triggered gene expression (Fig.  7D). These results suggest variation in MTI response is positively associated with the resistance level to these plant pathogens and is heritable.

DISCUSSION
Recognition of various MAMPs by their corresponding PRR receptors elicits a variety of biochemical and transcriptional responses that lead to plant-pathogen resistance (Katagiri and Tsuda, 2010;Segonzac and Zipfel, 2011). For example, MAMP-triggered oxidative burst, which is a hallmark of the early events of MTI, has been implicated in lignin production and in signal transduction during programmed cell death (Asai and Yoshioka, 2008). MAMP-triggered oxidative burst appears to play an important role in plant disease resistance. For example, Arabidopsis, tobacco (Nicotiana tabacum), and rice plants with NADPH oxidase (the main ROS producer during the oxidative burst) either silenced or inhibited by diphenyleneiodonium exhibited a significant increase in susceptibility to different pathogens (Torres et al., 2002;Yoshioka et al., 2003;Chi et al., 2009). In contrast, an enhanced oxidative burst through the application of different MAMPs leads to increased resistance to pathogen infection (Mersmann et al., 2010). Here, we showed that MAMP treatment of selected soybean genotypes results in a significant oxidative burst that varies by genotype. There is a rough correlation in the level of MAMP-triggered oxidative burst and the level of pathogen resistance after MAMP treatment. The results of this study provide support that the oxidative burst in soybean plays a role in resistance to two different pathogens.  To analyze genotypic variation in the MAMP response of different soybean genotypes (selected based on their MAMP-triggered oxidative burst and the availability of mapping populations), transcriptome analyses were performed on four soybean genotypes either mock treated or MAMP treated for 30 min, using a newly developed Affymetrix soybean whole-genome transcript array. The results showed significant variation in their transcriptional responses. MAMP-responsive genes map primarily into classes that consist of receptor kinases, protein modification enzymes, transcription factors, proteins involved in calcium regulation, and phytohormones. Similar results were observed in the transcriptional response of Arabidopsis and rice plants treated with different MAMPs (Fujiwara et al., 2004;Navarro et al., 2004;Zipfel et al., 2004;Wan et al., 2008). Consistent with the variation in MAMP-triggered oxidative burst among the four soybean genotypes, the transcriptional analysis also reflects genotypic variation. For example, the genotype LDX01-1-65, which was classified based on its oxidative burst as having a weaker MTI and also lower resistance to cyst nematode than LD00-2817P (Kim et al., 2011), expressed relatively few genes encoding for receptor-like kinases and transcription factors and, in general, showed a lower fold change (#2) in gene expression after MAMP addition. In contrast, those genotypes exhibiting a strong MAMPtriggered oxidative burst (e.g. LD00-2817P) also showed a strong transcriptional response to MAMP treatment. Analyses of Arabidopsis plants lacking genes belonging to these two functional categories, like PPRs (e.g. flg, chitin receptors), WRKY transcription factors, and different protein kinases, demonstrated that many are key players in the regulation of MTI and also important for plant disease resistance (Boudsocq et al., 2010;Segonzac and Zipfel, 2011).
It is well established that hormone levels also can strongly affect MTI. For example, variation in the levels of salicylic acid, jasmonic acid, ethylene, and BRs can influence plant disease resistance (Vert, 2008;Verhage et al., 2010). For example, ethylene and BR are important for FLS2 (flg receptor) accumulation and activation of BAK1 (a coactivator triggered by flg22 and EF-Tu), respectively. A defect in perception or accumulation of these hormones can negatively affect MTI and plant resistance (Boutrot et al., 2010;Oh et al., 2010;Jaillais et al., 2011). Consistent with these findings, MAMP treatment of soybean genotypes that exhibit a weak MTI resulted in little change in the expression of genes involved in hormone biosynthesis.
MTI represents a potential method to enhance soybean resistance to a broad array of pathogens and pests. For example, 24 h of treatment with flg22, chitin, or EF-Tu significantly reduced Arabidopsis colonization by P. syringae DC3000, Alternaria brassicicola, or Piriforomospora indica Wan et al., 2008;Jacobs et al., 2011). Several authors have pointed to the potential to harness MTI to protect crop plants from biotic stress (Ingvarsson and Street, 2011). By way of example, Lacombe and colleagues (2010) recently dem-onstrated that transfer of the EFR1 PRR receptor of the MAMP Elf-18, from Arabidopsis to either tomato or tobacco plants resulted in a significant increase in resistance to bacterial pathogens. However, a more direct approach would be to harness the endogenous MTI pathways in crop plants to elevate resistance to biotic stress through selection and breeding for stronger MTI response. Our current study is a step toward this goal, because breeding for enhanced MTI requires that genotypic variation in the MAMP-triggered responses be heritable with controlling genes accurately mapped to enable marker-assisted selection.
In this study, MTI appeared to be a quantitative trait. Therefore, we sought to identify specific QTLs that regulate soybean MTI by using either total ROS production or gene expression as a phenotype. This approach identified four QTLs, one related to total ROS production and the three others related to gene expression of MAMP-responsive genes. However, none of these QTLs colocalized with previous QTLs associated to cyst nematode resistance identified with this same mapping population (Kim et al., 2011). Interestingly, most of the genes localized in the trans-eQTL, as well as in the major effect cis-eQTL (chromosomes 13 and 4), are predicted to encode different receptor-like kinases (including LRR receptors), transcription factors, or different transporters (like zinc transporters). Hence, they may regulate MTI either by direct effects on intracellular regulation or by modifying membrane-localized signaling events (e.g. ion flux). By contrast, the QTL associated with total ROS production contains genes that encode for the NADPH oxidase RBOH, as well as enzymes involved in the biosynthesis of secondary metabolites. Because RBOH can directly modulate ROS levels during the oxidative burst (Asai and Yoshioka, 2008), this is probably the gene underlying this QTL. The minor effect cis-eQTL contains genes related to transport and secondary metabolism only.
In conclusion, results of this study demonstrate that the MTI responses among the four soybean genotypes varied depending on MAMP treatment and pathogen. Results also indicated that MTI variation among soybean genotypes was quantitatively inherited. Likewise, we propose that it may be possible to directly select for lines with stronger MTI through the development of markers that identify key QTLs that control plant innate immunity.

Plant Material and MAMP Treatment
This study was conducted by testing plants of four soybean (Glycine max) genotypes with contrasting genes for resistance to soybean cyst nematode (LD00-2817P [Diers et al., 2010] and LDX01-1-65 [Diers et al., 2005]) or sudden disease syndrome caused by Fusarium virguliforme (Ripley [de Farias Neto et al., 2007] and EF59 [Lightfoot et al., 2005]), as well as the 97 F3:4 (here after called as F3) lines from a cross between LD00-2817P (Diers et al., 2010) and LDX01-1-65. The lines from the LD00-2817P by LDX01-1-65 cross are the same lines that were used to test combinations of soybean cyst nematode (Heterodera glycines Ichinohe) resistance genes from PI 437654, PI 88788, and Peking that were bred into LD00-2817 and from Glycine soja PI 468916 that were bred into LDX01-1-65 (Kim et al., 2011). Seeds were surface sterilized (30% bleach, rinsed several times with sterile, double-distillated water), planted in soil-less medium (SogeMix SM-2 general purpose growing medium), and germinated in growth chambers (190 mmol photosynthetically active radiation m 22 s 21 , 23°C, 60% relative humidity, and 16-h photoperiod). Trifoliolate leaves from 3week-old plants were detached and then vacuum infiltrated with doubledistilled water for 2 min. Water-infiltrated trifoliolate leaves from five different plants of each genotype were pooled and cut into approximately 1-cm 2 slices. Equal amount of leaf slices (approximately 30 slices) were transferred into two different petri dishes and then floated overnight on autoclaved doubledistilled water. Water was removed from both petri dishes and was replaced with 5 mL of MAMP solution (1 mM flg22; Genescript) and 50 mg of crab shell chitin (called in this study as chitin; Sigma-Aldrich), or 5 mL of mock solution (autoclaved double-distilled water plus equivalent amount of dimethyl sulfoxide [DMSO; Fisher Scientific]). DMSO was included since it was contained in the solution used to dissolve the flg22 peptide. After a 30-min treatment, mock-and MAMP-treated leaf slices were harvested into different tubes and immediately frozen in liquid nitrogen. Samples were stored at 280°C until use. All procedures described above were performed under dark conditions.

Oxidative Burst Measurement in Soybean Leaf Discs
ROS production by leaf tissue was assayed by hydrogen peroxide-dependent luminescence of luminol (Keppler et al., 1989). Under dark conditions, four 4mm leaf discs per treatment from soybean water-infiltrated trifoliolate leaves were floated for 16 h in 200 mL of autoclaved double-distilled water in a 48-well plate. Continuing under dark conditions, water was removed and replaced with 150 mL of MAMP reaction buffer (20 mM luminol [Sigma-Aldrich], 1 mg horseradish peroxidase [Sigma-Aldrich], 1 mM flg22 [Genescript], 50 mg chitin [Sigma-Aldrich], or combination of both MAMPs) or mock reaction buffer (20 mM luminol, 1 mg horseradish peroxidase, and 2 mL of DMSO [Fisher Scientific]). Luminescence was immediately recorded over 30 min in a CCD photoncounting camera (Photek 216). The luminescence intensity of individual wells was calculated using Photek IFS32 software (Photek 216). The luminescence intensity data from MAMP-treated leaf discs were normalized using intensity values from mock-treated plants. Parental lines had five biological replicates with three technical replicates each, whereas three biological and technical replicates were performed for the F3 lines.

RNA Extraction, DNase Treatment, and qRT-PCR
Total RNA was purified from 0.5 g of MAMP-or mock-treated leaves by using Trizol reagent (Invitrogen), according to the manufacturer's instructions, and subsequently purified using chloroform extraction. Genomic DNA (gDNA) was removed from the purified RNA by using Turbo DNase (Ambion) following the manufacturer's instructions. Two micrograms of gDNA-free RNA were used to synthesize cDNA as described by Libault et al. (2010). The lack of gDNA contamination was verified by qPCR by using two different primers able to amplify gDNA.
qRT-PCR was performed as reported by Libault et al. (2008) using three housekeeping genes (con6, con16, and con4) to normalize the expression levels of the analyzed genes. Primer design was performed as described by Libault et al. (2010). Primer sequences are reported in Supplemental Tables S1 and S2. Gene expression levels were calculated according to the following equation: where P EFF is the primer efficiency calculated using LinRegPCR (Ramakers et al., 2003), and DCT is CT Housekeeping minus CT Gene , where CT Housekeeping is the geometric median from the cycle threshold (CT) of three different housekeeping genes and CT Gene is the CT from each gene. The gene expression levels were used to calculate the ratio between MAMP-and mocktreated leaves, and then the ratios were Log2 transformed to obtain the fold change. All the experiments were performed in 384-well plates. Parental lines had three biological replicates, whereas two biological replicates were performed for the F3:4 lines.

Affymetrix Soybean Whole Sense Transcript (Wild-Type) Exon-Array Design
A custom whole-genome Affymetrix Soybean Whole Sense orientation transcript exon array (Soybean WT array) was designed from the published soybean genome sequence (see Supplemental Materials and Methods S1). This array includes all the soybean genes, including both high-and low-confidence gene models, reported in Schmutz et al. (2010). To design probes for this array, sequences from the first draft assembly of the soybean genome (Schmutz et al., 2010) were used and downloaded from the Department of Energy-Joint Genome Institute Web site (pythozome: http://phytozome.net). The probe selection algorithm was developed by Affymetrix, and probes were selected and designed to span each exon of the predicted gene models, whenever possible. Each probe set has 11 different 25mers probes printed onto the Soybean WT array. In total, the Soybean WT array contains over 1,000,000 probe sets, some of them designed for background normalization purposes. An association probe set file is provided in Supplemental Table S15.

Affymetrix Soybean Whole Sense Transcript (Wild-Type) Exon-Array Hybridization and Analysis
Total RNA from both MAMP-and mock-treated leaves was used for Affymetrix Soybean WT array analysis following the manufacturer's protocols (Affymetrix). Briefly, 500 ng of gDNA-free total RNA were used to produce fragmented and biotin-labeled cDNA using the Ambion wild-type expression kit (Applied Biosystems) and the Affymetrix wild-type terminal labeling and hybridization kit (Affymetrix). Soybean WT array hybridizations were performed at the University of Missouri (http://biotech.rnet. missouri.edu/dnacore/), following standard Affymetrix procedures. The arrays were scanned with a GeneChip 7G plus high-resolution scanner. Both the gene level and the exon level expression values were summarized by Robust Multi-Array Average algorithm by using Expressionist Refiner 6.1 (GeneData). Briefly, the signal intensity cel files from scanner were analyzed by a workflow with GC background subtraction, quantile normalization, and median polish. The gene expression values from different samples were further used for gene differential statistical analysis. The exon level expression values were used in the microarray detection of alternative splicing (MiDAS) P value calculation. This MiDAS P value determines the gene splicing. To ensure the reliable exons in above MiDAS test, the detection above background P values were calculated for each probeset by comparing perfect mach probes to members of the background probes with the same GC content. A cutoff of detection above background P value of 0.05 was used for filtering and determining the exon present/absent. The full microarray datasets generated from the 24 Genechip used in this study were deposited in the Gene Expression Omnibus (http://www.ncbi.nlm. nih.gov/geo/) under the accession number GSE32642.

Identification of Differentially Expressed Genes
A linear mixed-effects model was applied to the normalized log-scale (base 2) expression measures separately for each gene using the bioconductor software R/mmanova (1.16.0 version). Each linear mixed model includes the genotype effect, treatment effect, the interaction between the genotype and treatment effects, and the biological replicate effect, where the biological replicate effect is considered to be random effect. The varianceshrinkage Fs tests (Cui et al., 2005) were conducted using R/mmanova to identify differentially expressed genes between the treatment and control groups for each of the four genotypes. P values for the variance-shrinkage Fs tests were obtained through 100 permutations via the Matest statement in the R/mmanova software, and a q value (Storey, 2002) was then computed for each P value to produce lists of differentially expressed genes with estimated FDRs of 5%. Among these significantly differentially expressed genes, genes with fold change above 2 were considered for functional classification through PFAM, KOG, and PANTHER domain predictions. Additionally, MapMan gene functional classification system was used (Thimm et al., 2004;Usadel et al., 2009). For MapMan analysis, an updated soybean mapping file (Yang et al., 2010), which allows to survey all the functional categories included in the software MapMan, was used.

Selection of Candidate Genes for QTL Analysis
A linear mixed model, as described above, was applied to the normalized GeneChip data from LD00-2817P and LDX01-1-65. The permutation Fs test (as described in the above section) was applied to conduct three comparisons: the comparison between the MAMP treatment and mock treatment for the two genotypes, and the comparison of fold changes between the LD00-2817P and LDX01-1-65. The intersection union principle (Berger, 1982) was utilized to identify genes that are statistically significantly up-regulated (P value # 0.05) with fold change above 1.5 in both genotypes (LD00-2817P and LDX01-1-65) but has statistically significant larger fold change in the line with strong MTI (LD00-2817P) than in the line with weak MTI (LDX01-1-65; P value , 0.05) with at least 2 times higher fold change in LD00-2817P. A total of 45 genes were identified with this criterion. Based on this list and their biological relevance in MTI, 25 genes from this analysis were selected to analyze their expression by qRT-PCR, as described above. Ratios between MAMP/mock were calculated to determine whether the MAMP treatment induced the expected response in both parental lines. After the logarithm transformation (base 2), the qRT-PCR data from the 99 treated plants (two parents and 97 F3 lines) were tested for normality using both the Lilliefors test and the Jarque-Bera test in Matlab software and normality plots were used as visual aid to identify genes that are normally distributed across the 99 plants. Genes that were declared as nonsignificant for both tests, which indicates no violation against the normality assumption, and have a close-to-straight-line normality plot were selected for QTL studies.

QTL Mapping Analysis
The population of F3 lines was evaluated with the Illumina GoldenGate 1,536 Universal Soy Linkage Panel 1.0 (USLP 1.0; Hyten et al., 2010) as described by Kim et al. (2011). The population segregated for 390 SNP markers from the USLP 1.0 and a linkage map was formed using these markers with Joinmap 3.0 (Van Ooijen and Voorips, 2001). QTLs controlling both ROS and the expression level of each gene assayed in the population through qRT-PCR were mapped through interval mapping with the program MapQTL 4.0 (Van Ooijen and Maliepaard, 1996;Van Ooijen, 2000). Genome-wide significance thresholds were determined by permutation tests in MapQTL 4.0.

Fungal Infection Assay
Cultures of Sclerotinia sclerotiorum were started 48-h prior to inoculation by subculturing actively growing edges of fungal colonies from stock cultures onto potato dextrose agar (DIFCO). For inoculation, 4-mm-diameter agar plugs with growing mycelium were cut from the edges of colonies using a cork borer.
Forty-eight hours before inoculation, trifoliolate leaves from 3-week-old soybean plants were detached and floated for 16 h in 20 mL of water in a petri dish. Solution was removed and the trifoliate leaves were treated with 0.01% Silwet L-77 for 20 min and immediately rinsed several times with doubledistilled water. The trifoliolate leaves were then treated for 24 h with 1 mM flg22, 50 mg chitin, 1 mM flg22 plus 50 mg chitin, or DMSO (as mock). Trifoliolate leaves were transferred into a petri dish (one per trifoliolate) that contained moistened Whatman paper. One agar plug per leaflet of each trifoliolate (three per trifoliolates in total) was placed in each petri dish. Petri dishes were sealed with Parafilm and then placed in a growth chamber with the same environmental conditions described above. Agar plugs were removed from the leaves 24 h after inoculation. Two days after inoculation, individual leaves from each trifoliate leaf were harvested and immediately frozen in liquid nitrogen and stored until used. S. sclerotiorum infection levels were determined by the abundance of S. sclerotiorum actin (Kasza et al., 2004) mRNA via qRT-PCR and following the conditions described in the section "RNA Extraction, DNase Treatment, and qRT-PCR." The experiment was repeated four times, each at different dates and with new inoculum, to obtain four biological replicates. Expression level of the S. sclerotiorum marker gene was logarithm base 2 transformed and then was applied to an ANOVA model with the treatment effect. The genotype effect and their interaction were estimated using the SAS software PROC GLM procedure. Forty pairwise comparisons were conducted using t tests to confirm the statistical significance of the differences among the four treatments for each genotype and the difference among the two parental lines as well as the four selected F3:4 lines for each treatment. The 40 P values were transformed to q value (Storey, 2002) to control the FDR level at 0.05.

Bacterial Infection Assay
Stocks of Pseudomonas syringae pv glycenea (Psg) were maintained on NYG agar medium (5 g/L bacto peptone, 3 g/L yeast extract, 20% [w/v] glycerol, 1.5% bacto agar) plates supplemented with rifampicin (50 mg/mL) and as 30% glycerol stock at 280°C. Forty-eight hours before infection, Psg was grown into NYG agar rifampicin (50 mg/mL) plates at 30°C. Just prior to the plant infection, a small amount of bacteria were scrapped from the plate and dissolved in 1 mL of 10 mM MgCl 2 , and the optical density at 600 nm was determined. From this bacterial solution, serial dilutions in 10 mM MgCl 2 were made to reach a bacterial concentration of approximately 10 5 CFU/mL. The CFU levels were determined by plating of serial dilutions onto NYG agar plates.
Forty-eight hours before inoculation, fifteen 4-mm diameter leaf discs from trifoliate leaves of 3-week-old soybean plants were floated for 16 h in 2 mL of autoclaved double-distilled water in a 12-well plate. The solution was removed, and the leaf discs were treated for 24 h with 1 mM flg22, 50 mg chitin mixture, 1 mM flg22 plus 50 mg chitin mixture, or DMSO (as mock). All these treatments were supplemented with 0.01% Silwet L-77. Solutions were removed, and the leaf discs were treated for another 3 h with a Psg cell suspension at a concentration of 10 5 CFU/mL. The bacterial suspension was removed, and the leaf discs were rinsed several times with autoclaved doubledistilled water and then floated for 3 d on water in the same 12-well plate. The quantification of bacterial colonization at 0 and 3 d post inoculation was determined by pooling and grinding four leaf discs in 10 mM MgCl 2 and then performing different serial dilutions with 10 mM MgCl 2 until 10 24 and 10 26 cells mL 21 were reached. These dilutions were sprayed into NYG agar rifampicin (50 mg/mL) plates, and then incubated at 30°C for 2 d. The number of colonies was expressed as CFU/cm 2 . This experiment was repeated four times with the same conclusion. CFU data were logarithm base 10 transformed and then applied to an ANOVA model with the treatment effect. The genotype effect, the day effect, and their interactions were estimated using SAS proc glm procedure. Two hundred and fifty two pairwise comparisons were conducted using t tests to check the statistical significance of the differences among the four treatments for each genotype on each day, the differences among two parental lines, as well as the four selected F3 lines for each treatment on each day and the differences among the days for each combination of treatments and genotypes. The 252 P values were transformed to q value (Storey, 2002) to control the FDR level at 0.05.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S2. Sclerotinia infection in whole plants. Figure S3. Expression level of four genes from the 25 genes used for QTL analysis.

Supplemental
Supplemental Figure S4. Variation of the pathogen resistance levels in MAMP-treated contrasting F3 lines.
Supplemental Figure S5. Illustration of the association between probe id and the predicted gene model.
Supplemental Table S1. Primers used to validate the Soybean WT array data by qRT-PCR.
Supplemental Table S2. Primer used for eQTL analysis.
Supplemental Table S12. Verification of the Soybean WT array results by qRT-PCR of randomly selected genes.
Supplemental Table S13. Marker name, chromosome, and LOD values associated from the QTLs associated to soybean MTI.
Supplemental Table S14. Genes localized at 1 cM in each QTL.
Supplemental Table S15. Association probe set file can be downloaded from http://seedgenenetwork.net/annotate#soybeanWT. Table S16. Expression level of 25 different genes tested in the F3:4 mapping population.
Supplemental Materials and Methods S1. Design of the Soybean WT array.