Exploring the genetics of lesion and nodal resistance in pea (Pisum sativum L.) to Sclerotinia sclerotiorum using genome‐wide association studies and RNA‐Seq

Abstract The disease white mold caused by the fungus Sclerotinia sclerotiorum is a significant threat to pea production, and improved resistance to this disease is needed. Nodal resistance in plants is a phenomenon where a fungal infection is prevented from passing through a node, and the infection is limited to an internode region. Nodal resistance has been observed in some pathosystems such as the pea (Pisum sativum L.)‐S. sclerotiorum pathosystem. In addition to nodal resistance, different pea lines display different levels of stem lesion size restriction, referred to as lesion resistance. It is unclear whether the genetics of lesion resistance and nodal resistance are identical or different. This study applied genome‐wide association studies (GWAS) and RNA‐Seq to understand the genetic makeup of these two types of resistance. The time series RNA‐Seq experiment consisted of two pea lines (the susceptible ‘Lifter’ and the partially resistant PI 240515), two treatments (mock inoculated samples and S. sclerotiorum‐inoculated samples), and three time points (12, 24, and 48 hr post inoculation). Integrated results from GWAS and RNA‐Seq analyses identified different redox‐related transcripts for lesion and nodal resistances. A transcript encoding a glutathione S‐transferase was the only shared resistance variant for both phenotypes. There were more leucine rich‐repeat containing transcripts found for lesion resistance, while different candidate resistance transcripts such as a VQ motif‐containing protein and a myo‐inositol oxygenase were found for nodal resistance. This study demonstrated the robustness of combining GWAS and RNA‐Seq for identifying white mold resistance in pea, and results suggest different genetics underlying lesion and nodal resistance.


| INTRODUCTION
Sclerotinia sclerotiorum (Lib.) de Bary, the causal agent of white mold disease, is one of the most destructive plant pathogens worldwide.
S. sclerotiorum is capable of infecting more than 400 host plants and causes millions of dollars of crop yield losses each year (Bolton, Thomma, & Nelson, 2006). Several studies have reported different secondary metabolites, effectors, and pathogenicity factors of S. sclerotiorum that are involved in establishing the infection (Bolton et al., 2006;Mbengue et al., 2016;Wei & Clough, 2016). One of the wellknown virulence strategies is the production of oxalic acid, which creates a low pH and acidic environment for infection (Xu, Xiang, White, & Chen, 2015). Oxalic acid suppresses reactive oxygen species (ROS) produced by plants at the beginning of infection and generates a reducing status that favors colonization (Williams, Kabbage, Kim, Britt, & Dickman, 2011). Fine-tuned redox homoeostasis from the initial reducing status to the later oxidative status in plant tissues is important for S. sclerotiorum to switch from the initial biotrophic lifestyle to the later necrotrophic lifestyle (Kabbage, Yarden, & Dickman, 2015). Studies searching for plant resistance to S. sclerotiorum have found quantitative interactions (McCaghey et al., 2017), and potential resistance genes included those with functions to maintain ROS and redox stresses during S. sclerotiorum infection (Girard et al., 2017;Ranjan et al., 2017;Zhou, Sun, & Xing, 2013).
Pea (Pisum sativum L.) is an important legume crop in the United States (Tayeh et al., 2015), and white mold continuously causes substantial damage and yield reduction (Biddle, 2001).

Sclerotinia sclerotiorum infection begins when ascospores of
S. sclerotiorum colonize blooms and invade through petioles into the stem. Severely infected plants will wilt and lodge. Resistance to white mold in pea has been observed via two different phenotypes. The first is lesion size where the length of stem lesion is measured after inoculation. The second phenotype is referred to as nodal resistance, which appears to be a unique mode of resistance and has been observed in some varieties of pea and soybean (Calla, Voung, Radwan, Hartman, & Clough, 2009;Porter, 2011;Porter, Hoheisel, & Coffman, 2009). Nodal resistance can be defined as the inhibition of lesion expansion at a node limiting pathogen colonization of plant stem tissue. Restriction of lesion expansion at the nodes has also been observed for stem-infecting fungi such as Diaporthe and Macrophomina species on soybean and cowpea (Hobbs, Schmitthenner, & Ellett, 1981;Muchero, Ehlers, Close, & Roberts, 2011). However, nodal resistance has been rarely documented and other than knowing lignin content is negatively correlated with nodal resistance in soybean (Peltier, Hatfield, & Grau, 2009), our understanding is limited.
Transcriptomics and differential expression (DE) analysis using RNA-Seq have become a standard approach to identifying resistance genes for white mold, and studies have applied this approach to oilseed rape (Brassica napus) and pea (Girard et al., 2017;Seifbarghi et al., 2017;Zhuang, McPhee, Coram, Peever, & Chilvers, 2012).
While most of these studies focused on the expression comparisons between a resistant and a susceptible variety, the genetic diversity of white mold resistance in B. napus might be underestimated using only this approach. Genome-wide association study (GWAS) is a robust approach to map white mold resistance and to capture the resistance diversity in a germplasm collection (Moellers et al., 2017;Wei et al., , 2017. The GWAS approach has been demonstrated in soybean (Glycine max Merr. L.) resistance to S. sclerotiorum where numerous single nucleotide polymorphisms (SNPs) associated with this quantitative resistance were discovered (Bastien, Sonah, & Belzile, 2014;Moellers et al., 2017;Wen et al., 2018;Wu, Zhao, Liu, et al., 2016). However, mapping results may discover SNPs that locate in intergenic genomic regions, and the interpretation of a confidence interval relies on the size of linkage disequilibrium (Bush & Moore, 2012). RNA-Seq and GWAS both have their advantages, and combining them provides a powerful tool to discover not only active genes that express in response to treatments, but also genetic diversity and SNPs associated with the treatment. This combined strategy has been applied to understand white mold resistance and yields in B. napus (Lu et al., 2017; and soybean (Wen et al., 2018), but not in pea. Because genes that can be found by both GWAS and RNA-Seq will have higher potential in contributing to white mold resistance, this study aimed to understand and compare the genetics of lesion and nodal resistance by applying both GWAS and RNA-Seq approaches in the pea-S. sclerotiorum pathosystem.

| GWAS: data source and analysis
Data used for GWAS were published in Porter et al. (2009). Briefly, there were 282 pea lines with a mean lesion resistance rating. The white mold fungus (S. sclerotiorum) Scl02-05 isolated from pea in Quincy, Washington, USA in 2003 was used for inoculations (Porter et al., 2009). A mean lesion resistance was measured in centimeters of lesion size, smaller values representing higher resistance. The data were collected after 72 hpi in a humid greenhouse and day/ night temperature ranges around 28°C/15°C. There were 266 pea accessions with nodal resistance ratings. Nodal resistance was measured using an ordinal scale from 0 to 5 after 2 weeks postinoculation, where 0 = dead plant; 1 = lesion expanded down the stem from the fourth inoculated node to the first node; 2 = lesion expanded from the fourth to the second node; 3 = lesion expanded from the fourth node to the third node; 4 = lesion did not expand beyond the initial inoculation point at the fourth node (Porter et al., 2009). There were four to eight replications to represent each accession. The USDA Pea Single Plant Plus Collection with SNP data was included in this study (Holdsworth et al., 2017). Association test was conducted in PLINK version 1.9 (Purcell et al., 2007).
Population stratification was controlled using a pairwise identity-bystate (IBS) clustering with a maximum clustering node of 2 and a p value cutoff of 0.05 for the pairwise population concordance test.
The IBS clustering matrix was included in a basic association test, and a minor allele frequency of 0.05 was applied. The empirical q value at 0.01 from an adaptive permutation test with default parameters was used to determine association significance. The genotyping-by-sequencing (GBS) raw reads containing significant SNPs were searched against the Trinity de novo transcriptome (assembled in the following sections) using BLASTN to acquire annotations at an E value cutoff of 10 À5 .

| Plant inoculations for RNA-Seq
A white mold-susceptible pea (P. sativum L.) cultivar 'Lifter' (PI 628276) and a white mold-partially resistant pea accession, PI 240515, were used in this study. These two lines were also used as parents in the development of a recombinant inbred line for investigating resistance, in a separate study.

| RNA extraction and sequencing
The

| De novo transcriptome assembly
Illumina raw reads were quality checked using FastQC version 0.11.5 (Andrews, 2010) and quality controlled using FASTX-toolkit version 0.0.14 (Gordon, 2014). Reads with 90 percent length above Phred score 30 were kept for analyses. Trimmomatic version 0.33 was used in default mode to remove adapters and to separate paired reads and single reads (Bolger, Lohse, & Usadel, 2014), and only paired reads were used for de novo assembly. All samples were pooled and aligned to the complete nearly gapless S. sclerotiorum genome sequence (Derbyshire et al., 2017) using the sensitive mode of Bowtie2 version 2.2.6 and Tophat2 version 2.1.0 (Kim et al., 2013). Reads unmapped to S. sclerotiorum genome were de novo assembled by Trinity version 2.4.0 using Kmer size 25, 29, and 32 (Grabherr et al., 2011;Haas et al., 2013).

| Differential expression, heatmap clustering, and gene ontology (GO) analyses
A k-mer index of 31 bp was built for the Trinity de novo transcriptome and paired-end reads were pseudo-aligned to the index using Kallisto version 0.43.0 with 1,000 bootstrap (Bray, Pimentel, Melsted, & Pachter, 2017). DE analysis was conducted using Sleuth version 3 in default mode using transcripts per million (TPM) normalization (Bray et al., 2017

| Expression verification using reverse transcription quantitative polymerase chain reaction (RT-qPCR)
A factorial experiment was set up with two pea varieties ('Lifter' and PI 240515), two treatments (mock and S. sclerotiorum inoculation), three time points (

| Phenotype data for lesion and nodal resistance
There were 282 and 266 pea germplasm lines screened for lesion and nodal resistance, respectively (Supporting Information   Table S1).  Table S2).

| Genome-wide association study
A total of 35,658 SNPs were included in the association analysis using PLINK (Purcell et al., 2007). There were 206 and 118 SNPs found to be significantly associated with lesion and nodal resistance, respectively (Supporting Information Table S3 and S4, respectively).
Without a standard genome for pea, the position and chromosome information for SNPs were deficient and it also made the annotation to these SNPs difficult. In order to understand the annotations of these significant SNPs, the original GBS raw reads harboring each SNP were retrieved (Holdsworth et al., 2017) and searched against our RNA-Seq de novo transcriptome using BLASTN.

| De novo transcriptome assembly
Using a stringent quality control threshold that keeps only raw reads with 90 percent of bases above Phred score of 30 (error rate 0.001%, one error per thousand bases), paired-end reads were mapped to the complete nearly gapless S. sclerotiorum genome using Tophat2 (Derbyshire et al., 2017;Kim et al., 2013). The de novo transcriptome using k-mer of 29 bp, which contains 96,588 transcripts including isoforms, resulted in the highest assembly quality (Table 1) similar to a previous de novo transcriptome of pea (Kerr, Gaiti, Beveridge, & Tanurdzic, 2017), and the re-mapped rate at 80% was satisfactory based on an empirical threshold of Trinity (Grabherr et al., 2011;Haas et al., 2013). Accordingly, the k-mer 29 de novo transcriptome was selected and a total of 60,598 transcripts were extracted from the longest representative isoform per gene model in the k-mer 29 de novo transcriptome (Table 1).   | 5 orthologous soybean gene (Table 3, Supporting Information   Table S4). In comparing the GWAS for lesion and nodal resistance, only one SNP (TP13557) can be found in both cases, and the de novo transcript containing this SNP was annotated as a putative glutathione S-transferase (GST) (orthologous to soybean gene Glyma.06G117800). Together with the weak phenotype correlation, the results suggest the genetics of lesion and nodal resistance may be different.

Seq samples
In     than 'Lifter' after S. sclerotiorum inoculation, including three transcripts that encode LRR receptor-like kinase (LRR-RLK). The first LRR-RLK is TRINITY_DN22904_c0_g1_i2, which had nearly zero expression in mock samples, but the expressions were induced higher in PI 240515 than 'Lifter' after S. sclerotiorum inoculation (Figure 6a).
The expressions of TRINITY_DN23231_c0_g2_i2 was higher in PI 240515 than 'Lifter' after S. sclerotiorum inoculation, but their expressions were down-regulated after S. sclerotiorum inoculation compared to mock samples (Figure 6b). On the other hand, the expression of TRINITY_DN4777_c0_g1_i1 was generally higher in 'Lifter' than PI 240515, and S. sclerotiorum inoculation caused up-regulation more in 'Lifter' than PI 240525 (Figure 6d) and the expression of TRINI-TY_DN18054_c0_g1_i1 and TRINITY_DN21848_c0_g1_i1 were higher in mock samples than in S. sclerotiorum-inoculated samples (Figure 6c,e).
It is worth noticing that the power for finding candidate resistance genes solely using expression difference is limited, and although the assumptions narrowed down the candidate pool, they may be subjective and other important genes might be neglected. Accordingly, we combined the power of GWAS and RNA-Seq to search candidate resistance genes with both genetic mapping and expression evidence.

| Integration of GWAS and RNA-Seq results
Integration of results from DE analyses and GWAS identified additional candidate resistance genes; however, most transcripts were down-regulated after S. sclerotiorum inoculation (Figure 3), and only a
about 6% of resistance (Wisser et al., 2011). Other studies also pointed out the importance of GST in potato, rice, and tobacco to Phytophthora infestans, Magnaporthe oryzae, and two Colletotrichum species, respectively (Dean, Goodwin, & Hsiang, 2005;Leonards-Schippers et al., 1994;Wisser, Sun, Hukbert, Kresovich, & Nelson, 2005). Moreover, most studies focusing on B. napus resistance to S. sclerotiorum also identified GST regardless of the methodologies, GWAS or RNA-Seq (Girard et al., 2017;Wu, Zhao, Liu, et al., 2016). A recent GWAS for soybean resistance to S. sclerotiorum also identified soybean GST  which was also noted to have high expression in a transcriptomic study (Calla et al., 2014). Consistently, our results also pointed out GST of pea may play a fundamental role in lesion and nodal resistance to S. sclerotiorum. GST has diverse molecular functions in a cell to balance redox homoeostasis, and glutathione is important for maintaining a reducing status for cell survival (Tew, 2007). Accordingly, GST may be involved in prohibiting the switch from biotrophic to necrotrophic stage of S. sclerotiorum. Another GST function is to detoxify phytotoxins and oxidative substances such as ROS (Wisser et al., 2011), and these functions may slow S. sclerotiorum infection and plant cell death. Other than GST, redox-related genes were up-regulated after   Tables S3 and S4; . Together, our results support the importance of redox homoeostasis for S. sclerotiorum resistance, and we identified many potential redox-related transcripts as well as others with roles in basal resistance to white mold.

| Lesion resistance
Several LRR-containing genes were found in transcriptomic studies of B. napus-S. sclerotiorum Wu, Zhao, Yang, et al., 2016). The results of our study also discovered several up-regulated LRR-containing transcripts for lesion resistance but not nodal resistance. Two transcripts with evidence from both GWAS and RNA-Seq are a putative CC-NBS-LRR protein and an LRR-RLK protein (Figure 7b,c). Both transcripts had higher expression in PI 240515 at 12 hpi but the expressions dropped over time to the expression level of 'Lifter.' Although it is well known that LRR-containing proteins contribute to R-gene based resistance in plants to biotrophic pathogens (Kushalappa, Yogendra, & Karre, 2016), it is unclear how much these LRR-containing transcripts are involved in lesion resistance to S. sclerotiorum. Moreover, these LRR-containing transcripts were not discovered in GWAS for nodal resistance, which indicated the possibility that the lesion resistance relies on LRR-containing proteins more than nodal resistance. Other than the LRR type of tandem repeats, several U-Box/ARM repeat-containing proteins were found for lesion resistance by GWAS (Table 2, Supporting Information Table S3). It has been shown that an U-Box/ ARM-containing ligase in rice negatively controls resistance (Li et al., 2012;Sharma & Pandey, 2015;Zeng et al., 2004), indicating higher expression in 'Lifter' might favor S. sclerotiorum infection ( Figure 7d). While plant cell wall synthesis enzymes such as cellulose synthase were identified, only a putative UDP-arabinopyranose mutase was up-regulated after S. sclerotiorum inoculation (Figure 7f). Similarly, two pleiotropic drug resistance ABC transporters were found (Table 2) but only one displayed up-regulation after S. sclerotiorum inoculation (Figure 7g). It has been shown that a pleiotropic drug resistance ABC transporter is involved in resistance to Botrytis cinerea, a closely-related fungal species to S. sclerotiorum (Stukkens et al., 2005). Accordingly, our results suggested diverse mechanisms were involved in lesion resistance to limit S. sclerotiorum expansion. Arabidopsis resistance to B. cinerea while overexpression increased susceptibility (Wang, Hu, Pan, & Yu, 2015). In addition, overexpression of Arabidopsis VQ5 and VQ20 demonstrated enhanced susceptibility to B. cinerea and Pseudomonas syringae (Cheng et al., 2012). These discoveries might help to explain the potential functions of this VQ motif-containing protein in pea-S. sclerotiorum interaction.

| Nodal resistance
Another transcript with earlier and higher expression in 'Lifter' was a putative myo-inositol oxygenase (Figure 7d). Myo-inositol is a product catalyzed from glucose-6-phosphate by the myo-inositol-1phosphate synthase, and it can be further metabolized into UDP-glucuronic acid by myo-inositol oxygenase (Kanter et al., 2005). One of the functions of UDP-glucuronic acid is being the precursor of plant Gene expression difference for myo-inositol metabolism was also reported in resistant and susceptible soybeans to S. sclerotiorum (Calla et al., 2009). Additionally, myo-inositol is also the precursor of galactinol, which has been suggested to induce resistance against syncytia development for the cyst nematode Heterodera schachtii (Siddeigue et al., 2014). When myo-inositol oxygenase processes myo-inositol into UDP-glucuronic acid, the metabolism bypasses and reduces the production of galactinol for inducing disease resistance signaling (Cho et al., 2010;Kim et al., 2008). Under this scenario,

| Summary
In this study, we localized significant SNPs identified from GWAS using a de novo transcriptome. We also incorporated expression analyses to understand the responses of candidate resistance genes to S. sclerotiorum infection. While we revealed SNPs exclusively for F I G U R E 8 Time course expressions of candidate nodal resistance transcripts identified from both differential expression analyses and genome-wide association studies. M represents mock samples, and WM represents Sclerotinia sclerotiorum-inoculated samples. (a) TRINITY_DN5298_c0_g1_i1 (b) TRINITY_DN25769_c0_g1_i1 (c) TRINITY_DN23515_c1_g1_i4 (d) TRINITY_DN21524_c0_g1_i1 (e) TRINITY_DN16214_c1_g2_i1 either lesion or nodal resistance, it is worth noticing that many SNPs appear to locate in the intergenic regions of the pea genome and many de novo transcripts might be too short to be annotated perfectly. The availability of a pea reference genome will not only improve the precision of SNP analysis but also provide a comprehensive understanding on coding sequences of pea. Nonetheless, our strategy integrating the advantages of GWAS and RNA-Seq indicated that aside from a single SNP located within a transcript encoding GST, there are likely different genetic mechanisms underlying lesion and nodal resistance.

ACKNOWLEDGMENTS
This study was supported by the National Sclerotinia Initiative, grant ID: 58-5442-9-239, USDA-ARS. We thank Patricia Santos and J. Alejandro Rojas for helpful discussions in preparation of this manuscript.

ACCESSI ON N UMBER
The Illumina sequences were deposited at SRA database under BioProject accession number PRJNA261444.

CONF LICT OF I NTEREST
The authors claimed no conflict of interest. ing, and all the co-authors refined and approved the manuscript.