Exploring Mechanisms of Quantitative Resistance to Leptosphaeria maculans (Blackleg) in the Cotyledons of Canola (Brassica napus) Based on Transcriptomic and Microscopic Analyses

Using resistant cultivars is a common approach to managing blackleg of canola/rapeseed caused by Leptosphaeria maculans (Lm). Quantitative resistance (QR), as opposed to major-gene resistance, is of interest because it is generally more durable, due to its multi-genetic basis. However, the mechanisms and genes underlying QR are mostly unknown. In this study, potential QR modes of action in “74-44 BL” was explored. This Canadian canola cultivar showed moderate but consistent race-nonspecific resistance at the cotyledon and adult-plant stages. A susceptible cultivar, “Westar”, was used as a control. After inoculation, the lesions developed more slowly on the cotyledons of 74-44 BL than those of Westar. We used RNA sequencing (-RNA-seq) to identify genes and their functions, putatively related to this resistance, and found that genes involved in programmed cell death (PCD), reactive oxygen species (ROS), signal transduction or intracellular endomembrane transport were most differentially expressed. ROS production was assessed in relation to Lm hyphal growth and lesion size; it occurred beyond the tissue colonized by Lm in 74-44 BL and appeared to trigger rapid cell death, limiting cotyledon colonization by Lm. In contrast, Lm grew more rapidly in Westar, often catching up with the ring of ROS and surpassing lesion boundaries. It appears that QR in 74-44 BL cotyledons is associated with limited colonization by Lm possibly mediated via ROS. The RNA-seq data also showed a link between ROS, signal transduction, and endomembrane vesicle trafficking, as well as PCD in the resistance. These results provide a starting point for a better understanding of the mechanisms behind QR against Lm in canola.


Introduction
Canola or rapeseed (Brassica napus L.) is an economically important oilseed crop cultivated worldwide. Blackleg, caused by Leptosphaeria maculans (Lm) Ces. and de Not., is a serious disease of canola, especially in Australia, Europe, and Canada [1]. Genetic resistance is a cornerstone of blackleg management and is usually classified as either qualitative or quantitative. The former is controlled by single resistance (R) genes, while the latter is often, though not always, polygenic [2,3]. While many of the R genes have been identified [4][5][6][7][8][9], quantitative resistance (QR) is not well understood. QR can be attributed to multiple genomic regions in B. napus [10], with many of the same loci found in multiple canola cultivars [11]. QR to blackleg in canola is believed to be expressed primarily in adult plants.   Three libraries (replicates) were produced for each of the four treatments, i.e., Westar and 74-44 BL with mock and Lm inoculation, respectively. A total of twelve libraries were used for RNA-seq. Approximately 14.1-17.8 million paired-end reads were obtained from each library (Supplementary  Table S1). When annotated against Lm genome, a higher percentage of reads was mapped for the inoculated Westar, as compared to the inoculated 74-44 BL (Figure 2A). Principle-component analysis (PCA) indicated that the treatments grouped tightly together in terms of their alignment to the Lm genome ( Figure 2B). Using the criteria of adjusted p value ≤ 0.05 and log 2 fold change ≥ 2, only 16 differentially expressed genes (DEG) of Lm were identified between inoculated Westar and 74-44 BL, with three DEGs up-regulated and thirteen down-regulated in the 74-44 BL, relative to those in Westar (Supplementary Tables S2 and S3). Three libraries (replicates) were produced for each of the four treatments, i.e., Westar and 74-44 BL with mock and Lm inoculation, respectively. A total of twelve libraries were used for RNA-seq. Approximately 14.1-17.8 million paired-end reads were obtained from each library (Supplementary  Table S1). When annotated against Lm genome, a higher percentage of reads was mapped for the inoculated Westar, as compared to the inoculated 74-44 BL (Figure 2A). Principle-component analysis (PCA) indicated that the treatments grouped tightly together in terms of their alignment to the Lm genome ( Figure 2B). Using the criteria of adjusted p value ≤ 0.05 and log2 fold change ≥ 2, only 16 differentially expressed genes (DEG) of Lm were identified between inoculated Westar and 74-44 BL, with three DEGs up-regulated and thirteen down-regulated in the 74-44 BL, relative to those in Westar (Supplementary Tables S2 and S3). The Lm genes generally were expressed at low levels; only eight of them showed a base-mean expression value over 10,000 (14,346 to 40,534), three between 1000 and 9999, 85 between 100 and 999,  The Lm genes generally were expressed at low levels; only eight of them showed a base-mean expression value over 10,000 (14,346 to 40,534), three between 1000 and 9999, 85 between 100 and 999, and 152 between 50 and 99. There were a total of 12,119 Lm genes with non-zero expression values, the vast majority of which had base-means under 50.

RNA-Seq Analyses
The three upregulated Lm DEGs in inoculated 74-44 BL, had sequence similarities to genes encoding a short-chain dehydrogenase/reductase, a pyoverdine biosynthesis, and a hypothetical protein, respectively. Pyoverdine is a siderophore biosynthesized by Pseudomonads [26]. Zwiers et al. [27] found a gene encoding an ABC-transporter with a pyoverdine biosynthesis motif in the fungus Mycosphaerella graminicola. ABC-transporters can play a role in the virulence of fungal pathogens [28,29].
It is unclear what role, if any, these Lm DEGs played during the infection of 74-44BL.

Expression of the B. napus Genes
A heatmap based on the level of relative gene expression showed similar patterns among three replicates within a treatment ( Figure 3). Similarities, as well as differences, were observed between inoculated Westar and 74-44 BL, and also between two mock inoculations.
Upregulated B. napus Genes in Inoculated 74-44BL A total of 908 genes in the B. napus genome were differentially upregulated in inoculated 74-44 BL, relative to those in inoculated Westar ( Figure 4A). The majority of these DEGs had adjusted p values and log 2 -fold changes in expression relatively close to the cut-off levels ( Figure 5A). Two of the DEGs showed base-mean expression levels over 10,000, six between 5000 and 9999, and 65 between 4999 and 1000. Five DEGs putatively encoding peptidases were among those with the highest scores (Supplementary Table S4). The three highest-scored DEGs, BnaA01g17570D, BnaA01g04000D and BnaA09g52180D, all encode peptidases. BnaA01g17570D, with InterPro domains, points to a cysteine peptidase belonging to family C1 (papain-like). BnaA01g04000D is a putative legumain peptidase C13, also known as a vacuole processing enzyme (VPE). BnaA09g52180D is also related to a cysteine peptidase. The DEG BnaC02g00130D potentially encodes a protease involved in RuBisCO degradation. Additionally, many chlorophyll A-B binding proteins showed high basemeans and were more highly expressed in inoculated 74-44 BL (Supplementary Table S4), such as an ATPase of AAA-type with a protein BLAST similarity to RuBisCO activase. This protein is potentially involved in endoplasmic reticulum (ER) to Golgi membrane budding.
Genes that potentially encode glycoside hydrolases, including a beta-galactosidase (BnaA04g04110D) and a pullulanase-type alpha-1,6-glucosidase (BnaA10g25820D), were upregulated in inoculated 74-44 BL, with small p values (Supplementary Table S4). A gene putatively encoding lactate/malate dehydrogenase (BnaC02g00740D) was also upregulated, albeit with a less significant p value but a higher base-mean expression than BnaA04g04110D or BnaA10g25820D (Supplementary  Table S4, Figure 4). BnaA03g11710D, with a thiazole biosynthetic enzyme InterPro domain, is a DEG with a protein sequence similar to a ribulose-1,5-biphoshate synthetase.
Gene ontology (GO) term enrichment analysis of these 908 DEGs was consistent with the results described above (Supplementary Table S13, Figure 4); many of the GO terms with the lowest FDR were related to photosynthesis and light responses. Similarly, the KEGG metabolic-pathway annotation also included highly expressed DEGs linked to photosynthesis with low adjusted p values and extremely high log 2 -fold changes in expression (Supplementary Table S14). Furthermore, three GO terms were linked to hydrogen peroxide (HP). While none of the enriched GO terms suggested peptidase activities, the GO term with the second lowest FDR was associated with cysteine biosynthesis (Supplementary  Table S13). This was consistent with the putative cysteine peptidase activity found with the DEG BnaA01g17570D (Supplementary Table S4). Plants 2020, 9,  A total of 640 genes were more highly expressed in inoculated Westar, as compared to inoculated 74-44 BL ( Figure 4B). These DEGs generally displayed lower expression values relative to the 908 DEGs upregulated in inoculated 74-44BL, ranging from a base-mean of 3,410 to 1.25. Only 11 DEGs showed base-means over 1000, 28 between 500 and 999, 73 between 100 and 499, and the remaining 527 were under 100. However, these upregulated DEG clusters were noticeably more above the cut-off levels (log 2 -fold change, adjusted p value) than those in inoculated 74-44 BL ( Figure 5A).
The DEG with the highest basemean, BnaC09g20030D, showed similarity to the genes related to a Bax inhibitor-1 (Supplementary Table S13), while the DEG BnaCnng58090D, with a basemean of 2354, was linked to a development/cell death domain (DCD). Other significant DEGs included BnaC08g42820D, which was associated with a heat shock protein 70, as well as BnaA04g06220D and BnaA09g26960D, which were linked to Section 23/Sec24 and Section 61/SecY, respectively. Section 23 and Section 24 are part of the coat protein II (COPII) complex that is involved in endoplasmic reticulum (ER) to Golgi vesicle transport [30]. The DEGs BnaA08g26550D, BnaA06g05280D, BnaC06g24690D, BnaA07g09950D, and BnaCnng06680D, with basemeans ranging from 972 to 3100, appeared to be related to small GTPases (Supplementary Table S13, Figure 4A).
GO terms related to ER, ER stress, vesicle transport, and the cellular endomembrane system were enriched in Lm-inoculated Westar. However, none of the enriched GO terms were associated with PCD. One enriched GO term was related to the HP response (Supplementary Table S8). BnaCnng58090D did not appear to be associated with any GO terms. KEGG annotation revealed less information than the DEGs upregulated in 74-44 BL. The KEGG data did not identify a connection between the resistance and the endomembrane system or PCD (Supplementary Table S9). However, this was unsurprising since KEGG focuses on metabolic pathways. Multiple DEGs were annotated by KEGG in relation to purine metabolism.

Upregulated B. napus Genes in Mock-Inoculated 74-44 BL
There were 774 genes upregulated in mock-inoculated 74-44 BL, relative to mock-inoculated Westar, including InterPro domains (BnaC08g39910D), which suggest xylogucan fucosyl-transferase activity, putative C2 calcium-dependent membrane targeting genes (BnaA09g48570D, BnaA09g05080D), and those with WD40 repeats (BnaA05g22060D, BnaC03g03670D). In addition, BnaC09g37530D showed putative glycoside hydrolase and raffinose synthase activity. A potential lectin (BnaA01g28810D) and isocitrate dehydrogenase (BnaC03g01130D) were also upregulated (Supplementary Table S14). Only one GO term, with the molecular function of histone acetyltransferase activity, was enriched among these DEGs (Supplementary Table S14), with two and four DEGs showing log 2 -fold change values < −5 or log 10 adjusted p values < −20 ( Figure 5B). There were only three DEGs with scores (based on expression, log 2 -fold change in expression, and adjusted p value) < 1 × 10 3 , with the KEGG annotation relating to fructose and mannose or starch and sucrose metabolism, or to fatty acid biosynthesis and elongations (Supplementary Table S12).

HP Production in Cotyledons
RNA-seq results suggested that ROS, such as HP, might play a role in QR of 74-44 BL against Lm. To validate this finding, cotyledon tissues were stained with 3,3-diaminobenzidine (DAB), to quantify the area of ROS production surrounding the inoculation site.
At 7 dpi, the lesion size, extent of Lm hyphal colonization, and area with ROS reaction varied on inoculated cotyledons, depending on the cultivar and parameter measured. In Westar, the area colonized by hyphae (as visualized by GFP fluorescence) and the area stained positive for HP were both larger than the necrotic lesions ( Figure 6A). In contrast, the lesion size and area colonized by the Lm hyphae were significantly smaller than the area stained for ROS in 74-44 BL. As with the results in Figure 1, the lesion size did not differ between the two cultivars at 7 dpi, while the area colonized by the Lm hyphae was substantially greater in Westar. The area with ROS staining did not differ between the cultivars ( Figure 6B).

HP Production in Cotyledons
RNA-seq results suggested that ROS, such as HP, might play a role in QR of 74-44 BL against Lm. To validate this finding, cotyledon tissues were stained with 3,3-diaminobenzidine (DAB), to quantify the area of ROS production surrounding the inoculation site.
At 7 dpi, the lesion size, extent of Lm hyphal colonization, and area with ROS reaction varied on inoculated cotyledons, depending on the cultivar and parameter measured. In Westar, the area colonized by hyphae (as visualized by GFP fluorescence) and the area stained positive for HP were both larger than the necrotic lesions ( Figure 6A). In contrast, the lesion size and area colonized by the Lm hyphae were significantly smaller than the area stained for ROS in 74-44 BL. As with the results in Figure 1, the lesion size did not differ between the two cultivars at 7 dpi, while the area colonized by the Lm hyphae was substantially greater in Westar. The area with ROS staining did not differ between the cultivars ( Figure 6B).
The HP reaction was also examined in a time-course post-inoculation, in order to obtain more detailed data than that presented in Figure 6. Most parameters measured tended to increase with time, but the responses varied between cultivars. In Westar, Lm colonization developed more rapidly than it did in 74-44 BL, with significant development beyond the lesions (p ≤ 0.05, Tukey adjusted) but generally within the boundary of DAB staining (Figure 7). In 74-44 BL, Lm hyphal growth was more restricted, relative to the DAB-stained areas at 5-11 dpi, and to both DAB staining and the lesion at 11 dpi (Tukey adjusted, p ≤ 0.05). The HP reaction was also examined in a time-course post-inoculation, in order to obtain more detailed data than that presented in Figure 6. Most parameters measured tended to increase with time, but the responses varied between cultivars. In Westar, Lm colonization developed more rapidly than it did in 74-44 BL, with significant development beyond the lesions (p ≤ 0.05, Tukey adjusted) but generally within the boundary of DAB staining (Figure 7). In 74-44 BL, Lm hyphal growth was more restricted, relative to the DAB-stained areas at 5-11 dpi, and to both DAB staining and the lesion at 11 dpi (Tukey adjusted, p ≤ 0.05).

Genomic DNA Degradation as an Indicator of PCD
As RNA-seq showed that PCD could also play a role in QR in 74-44 BL, degradation of genomic DNA was examined as a proxy for PCD. However, no apparent difference in genomic DNA degradation was observed between the inoculated Westar and 74-44 BL, based on either agarose gel

Genomic DNA Degradation as an Indicator of PCD
As RNA-seq showed that PCD could also play a role in QR in 74-44 BL, degradation of genomic DNA was examined as a proxy for PCD. However, no apparent difference in genomic DNA degradation was observed between the inoculated Westar and 74-44 BL, based on either agarose gel electrophoresis or Experion 12K (Figure 8).

Impact of Protease Inhibitors on Lm Infection of Cotyledons
Results from RNA-seq also led to the hypothesis that proteases could contribute to QR in 74-44 BL. Several protease inhibitors, when applied in liquid to the surface of cotyledons of Westar or 74-44 BL, failed to show any significant impact on either the lesion size or the extent of Lm hyphal colonization, although the latter was consistently greater in Westar (Figure 9).

Discussion
QR to Lm was expected to be more prominent in adult canola/rapeseed plants than in seedlings [2]. However, this study demonstrated a quantitative reduction of Lm infection in cotyledons of 74-44 BL, in the absence of major-gene involvement, by using an isolate without an Avr gene corresponding to any of the R genes in this resistant cultivar. Specifically, more restricted host tissue colonization and lesion development were observed on 74-44 BL, relative to Westar. Earlier studies [12,32] also found consistent QR with 74-44 BL in stems, as well as in cotyledons, against multiple Lm races, confirming the race-nonspecific nature of resistance. This study also found that the Lm hyphal growth often extended beyond the borders of visible lesions in Westar, while the growth was consistently restricted within the lesion in 74-44 BL. Huang et al. [25,33] also measured QR to Lm growth in young B. napus plants and in some cases, found that this restricted early growth correlated with reduced blackleg in more mature plants [33]. They also found partial overlap in quantitative trait loci (QTL) contributing to Lm resistance at both young and adult plant stages. RNA-seq was used to explore genes potentially involved in the seedling-stage QR in 74-44 BL as a first step to understanding the molecular mechanisms.
Many of the highly upregulated genes in Westar relate to the control of PCD, endomembrane vesicle trafficking between the ER and Golgi, as well as molecular chaperones, cation transporters, protein glycosylases, and degradation enzymes (Supplementary Table S13). The GO terms enriched in these DEGs also support a role with endomembrane vesicle transport (Supplementary Table S13). The upregulation of BnaCnng58090D (Supplementary Table S13), a gene with sequence similarity to a development/cell death (DCD) domain, could potentially stimulate a hypersensitive reaction, a form of PCD in plants [34], in response to infection. Other upregulated genes with putative roles in endomembrane transport to/from the ER are potentially related to ER stress, which can also trigger DCD-mediated PCD [35]. However, the enrichment analysis did not uncover any GO terms related to PCD (Supplementary Tables S6 and S7). Furthermore, the upregulation of BnaC09g20030D, with sequence similarity to a Bax inhibitor-1, might negate the effect of BnaCnng58090D by inhibiting PCD [32]. The hypersensitivity that DCD would otherwise be induced in response to infection might be prevented by the upregulation of BnaC09g20030D in Westar.
One of the DEGs upregulated uniquely in inoculated Westar is a heat-shock protein 70, which might be linked to the upregulation of Bax inhibitor-1. Qi et al. [36] noted that the overexpression of a heat-shock protein 70 inhibited PCD induced by HP. This might explain the results of this study, i.e., that the plant cell death (lesion) generally fell behind HP production and Lm hyphal growth in Westar. The lack of observed differences in genomic DNA degradation failed to support this hypothesis. Although fragmentation of genomic DNA can be associated with PCD, including that mimicking apoptosis (in animal cells) and that involved in normal plant developmental processes [37], the results of Ruberti et al. [38] showed complex roles of Bax inhibitor-1 in plant PCD, which is currently not well understood. Therefore, the results of indifference in genomic-DNA degradation as a proxy for PCD cannot be clearly interpreted at this point.
PCD in general, and Bax inhibitor-1 in particular, play a role in plant resistance to pathogens. For example, Babaeizad et al. [39] found that the overexpression of Bax inhibitor-1 in barley led to increased susceptibility to the biotrophic fungal pathogen Blumeria graminis f.sp. hordei (powdery mildew). This finding was consistent with the upregulation of a DEG related to Bax inhibitor-1 in Westar, corresponding to a greater biotrophic growth of the Lm hyphae asymptomatically in live cells beyond necrotic tissues (lesions) on cotyledons. Increased PCD is associated with the resistance to biotrophic infection and susceptibility to necrotrophic colonization by fungal pathogens. Scotton et al. [40] suggested that constitutive overexpression of Bax inhibitor-1 result in elevated resistance to necrotrophic fungal pathogens. Lm is considered a hemibiotroph [19], with an initial biotrophic infection, followed by extensive necrotrophic colonization of host leaf tissues.
In inoculated 74-44 BL, many genes related to peptidases were more highly expressed than in inoculated Westar; these included those related to papain cysteine peptidases (BnaA01g17570D, BnaA09g52180D) [41,42] and legumain peptidase C13 (BnaA01g04000D), also known as a vacuole processing enzyme (VPE). VPEs, as suggested by the name, are located in plant vacuoles [43,44], as illustrated in Figure 10. In addition, the gene BnaC02g00130D, which putatively encodes a protease involved in the degradation of RuBisCO, was upregulated. These genes might also be involved in PCD (reviewed by Zamyatnin [45]), which causes the plant cell vacuole to rupture, releasing proteases that degrade cellular components [46]. Protease-mediated PCD is essential for plant hypersensitive responses (reviewed by Sueldo and van der Hoorn [47]), which limit infection during the biotrophic phase. These proteases merit further investigation for potential involvement in QR. illustrated in Figure 10. In addition, the gene BnaC02g00130D, which putatively encodes a protease involved in the degradation of RuBisCO, was upregulated. These genes might also be involved in PCD (reviewed by Zamyatnin [45]), which causes the plant cell vacuole to rupture, releasing proteases that degrade cellular components [46]. Protease-mediated PCD is essential for plant hypersensitive responses (reviewed by Sueldo and van der Hoorn [47]), which limit infection during the biotrophic phase. These proteases merit further investigation for potential involvement in QR. The protease inhibitor experiments were intended to verify the role for some of the peptidases in limiting Lm hyphal growth in inoculated 74-44 BL, which could possibly occur through a role in PCD. However, there was a lack of significant protease-or peptidase-related GO terms in the enrichment analysis (Supplementary Table S14), which would call into question the significance of proteases in the resistance. The lack of differences between inhibitor treatments, coupled with the insignificant GO terms, did not support proteases as important players in QR with 74-44 BL, despite strong upregulation of several peptidase-related genes.
Several genes related to chlorophyll A-B binding proteins (which are a source of ROS [48]) were also upregulated in inoculated 74-44 BL and thus might be involved in QR. ROS, including HP, can act as pro-PCD signals (reviewed by Galvez-Valdivieso and Mullineaux [49]). This was supported by the findings from the current study where HP production occurred beyond the lesion development and Lm hyphal fronts in inoculated 74-44 BL (Figures 6 and 7). The HP production (ROS) appeared to result in rapid plant cell death that restricted Lm hyphal growth. In contrast, Lm hyphal growth kept in pace with the circle caused by HP beyond the lesion, showing little restriction in Westar. It appeared that ROS in 74-44 BL resulted in a near-hypersensitive reaction rapidly surrounding the Lm hyphae, limiting the biotrophic phase of infection beyond the lesion. Additional mechanisms aside from ROS might also play a role. BnaA03g11710D was another DEG in the photosynthetic process upregulated in inoculated 74-44 BL, which putatively encodes the thiazole biosynthetic enzyme or ribulose-1,5-bisphosphate synthetase. Thiazole is a precursor of vitamin B1 (thiamine) which can activate plant defenses [50]. Ahn et al. [51] and Boubakri et al. [52] found a relationship between vitamin B1-induced disease resistance and HP. Consistently, KEGG analysis also identified BnaA03g11710D as being involved in thiamine metabolism (Supplementary Table S14). These results showed a link between upregulation of BnaA03g11710D and increased HP production in inoculated 74-44 BL, further strengthening the case for this DEG being involved in QR. Additionally, three GO The protease inhibitor experiments were intended to verify the role for some of the peptidases in limiting Lm hyphal growth in inoculated 74-44 BL, which could possibly occur through a role in PCD. However, there was a lack of significant protease-or peptidase-related GO terms in the enrichment analysis (Supplementary Table S14), which would call into question the significance of proteases in the resistance. The lack of differences between inhibitor treatments, coupled with the insignificant GO terms, did not support proteases as important players in QR with 74-44 BL, despite strong upregulation of several peptidase-related genes.
Several genes related to chlorophyll A-B binding proteins (which are a source of ROS [48]) were also upregulated in inoculated 74-44 BL and thus might be involved in QR. ROS, including HP, can act as pro-PCD signals (reviewed by Galvez-Valdivieso and Mullineaux [49]). This was supported by the findings from the current study where HP production occurred beyond the lesion development and Lm hyphal fronts in inoculated 74-44 BL (Figures 6 and 7). The HP production (ROS) appeared to result in rapid plant cell death that restricted Lm hyphal growth. In contrast, Lm hyphal growth kept in pace with the circle caused by HP beyond the lesion, showing little restriction in Westar. It appeared that ROS in 74-44 BL resulted in a near-hypersensitive reaction rapidly surrounding the Lm hyphae, limiting the biotrophic phase of infection beyond the lesion. Additional mechanisms aside from ROS might also play a role. BnaA03g11710D was another DEG in the photosynthetic process upregulated in inoculated 74-44 BL, which putatively encodes the thiazole biosynthetic enzyme or ribulose-1,5-bisphosphate synthetase. Thiazole is a precursor of vitamin B1 (thiamine) which can activate plant defenses [50]. Ahn et al. [51] and Boubakri et al. [52] found a relationship between vitamin B1-induced disease resistance and HP. Consistently, KEGG analysis also identified BnaA03g11710D as being involved in thiamine metabolism (Supplementary Table S14). These results showed a link between upregulation of BnaA03g11710D and increased HP production in inoculated 74-44 BL, further strengthening the case for this DEG being involved in QR. Additionally, three GO terms linked to HP were also identified (Supplementary Table S14), adding support to ROS involvement in QR.
The DEGs unique to mock-inoculated Westar and 74-44 BL provide some insight into the differences between the two cultivars, independent of Lm infection. For example, a WD40 repeat-containing protein and a putative C2 calcium-dependent membrane targeting protein, both of which could play a role in signaling [53,54], were upregulated in the 74-44 BL, while a gene potentially involved in intracellular signaling was upregulated in Westar (Supplementary Tables S6 and S7). These might suggest that signaling networks differ somewhat between these two cultivars, in the absence of Lm infection. The higher expression of certain genes putatively encoding a late embryogenesis abundance protein, a disease resistance protein and thaumatin (which could also be involved in plant disease resistance [55]) in Westar, without exposure to Lm, might point to inherent differences between the cultivars (Supplementary Table S13). Clearly, these putative defense-related genes, at least at the levels they are expressed in Westar, appear insufficient to prevent infection of Westar by Lm.
The DEGs involved in endomembrane trafficking, such as small GTPases, Section 23/sec24, Section 61 and WD40 repeats, can be involved in ER stress and unfolded protein responses (UPR); if not resolved, both processes can lead to PCD (reviewed by Williams et al. [56]). Vesicle trafficking between the Golgi and ER, and vice versa, affects traffic to the vacuole, where VPE-mediated PCD takes place ( Figure 10). VPE is one of the proteases upregulated in inoculated 74-44 BL. UPR is triggered by improperly folded proteins in ER; if the stress to ER is severe, UPR could also induce PCD (reviewed by Cui et al. [51]). As differences were not observed in genomic DNA degradation, the UPR might only be a signaling mechanism related to QR, rather than PCD. Further research into the potential roles of UPR, endomembrane dynamics, and ROS in plant defenses is merited for a better understanding of QR against the blackleg of canola.

Materials and Methods
This manuscript includes the following five experiments on the cotyledons of Westar (susceptible) and 74-44 BL, with a level of QR against the different Lm races [11,12]: (1) RNA-seq and corresponding infection assessment, (2) a time-course evaluation on necrotic lesion development corresponding to Lm colonization post inoculation, (3) staining for HPwith DAB as an indicator for the production of ROS, (4) a protease-inhibitor study to validate the role of peptidases in resistance to Lm, and (5) an assessment of the level of fragmentation of genomic DNA as a proxy for PCD.

Fungal and Plant Material
Inoculum was prepared from L. maculans isolates 12CC09 carrying AvrLm6,7 and 12CC09-GFP, grown on V8 agar for about 10 days, until pycnidia were visible. The isolate, 12CC09, was selected for this study, as it lacked AvrLm1, AvrLm3, and AvrLmS, which would interact with the R genes Rlm1, Rlm3, or RlmS in 74-44 BL [11]. Pycnidiospores were harvested in sterile water, filtered through a Falcon™ Cell Strainer (70 µm pore size), diluted to 2 × 10 7 spores/mL, and stored at −20 • C until use. One week after planting, the cotyledons were wounded on each lobe with bent-tipped tweezers before being inoculated with 10 µL droplets of water or pycnidiospore suspension.
The isolate 12CC09-GFP was generated by transforming the isolate 12CC09 with a binary vector containing the GFP gene via Agrobacterium-mediated transformation. Aliquots of 1 mL Lm 12CC09 suspension (1 × 10 7 spores/mL) were mixed with 1 mL of log-phase Agrobacterium tumefaciens Agl1 cells carrying the GFP vector. The cell mixture was pelleted by centrifugation, resuspended with the liquid AB-MES medium [57], and spread on sterile 0.45 µm nitrocellulose membranes (GE Healthcare, Ottawa, ON, Canada) overlaid on a solid AB-MES medium. After 3-day incubation in darkness at room temperature, the membranes were cut into 0.5-cm strips, placed on a selection medium (solid V8 containing 50 µg/mL hygromycin and 50 µg/mL carbenicillium), and incubated at 24 • C for 7-10 days. Pycnidial ooze from single colonies was transferred individually to the selection medium and the resulting cultures (potential transformants) were assessed for GFP expression with a Zeiss Stereo-Lumar epifluorescence microscope. The isolate 12CC09-GFP expressed strong and consistent GFP signal in a pure culture and in canola leaf tissues ( Figure 1C,D).
74-44 BL is DEKALB ® hybrid with multi-genic Lm resistance and R genes Rlm1, Rlm3, and RlmS (Saskatchewan Seed Guide, 2019). The differential reaction of the 74-44 BL cotyledons to Lm isolates with and without AvrLm1 and AvrLm3 (Supplementary Figure S1) confirmed the presence of Rlm1/LepR3 and Rlm3. 74-44 BL also carries a level of QR against multiple Lm races in cotyledons [12], with lower lesion scores [32], more restricted Lm colonization [12], as well as QR displayed in adult plants in the absence of major-gene interactions [32].
Plants were grown in Sunshine #3 soil-less mix (Sun Gro Horticulture Canada Ltd., Vancouver, BC, Canada) to which 12.5 g L −1 Osmocote Plus 16-9-12 (N-P-K; Scotts Miracle-Gro Canada, Mississauga, ON, Canada) was added. For all experiments, except those involving the time series that did not involve DAB staining, canola plants were grown in 72-well flats. The flats were placed in a growth chamber set to 22 • C and 16 • C, during the 16 h of light (approximately 280-575 µmol m −2 s −1 ) and 8 h of darkness, respectively. Plants intended for the time-series microscopic examination were grown either as described above or in the greenhouse, in 10 cm square pots, exposed to a mix of natural and fluorescent (430W Philips high pressure sodium lamps) light, and inoculated with water, 12CC09, or 12CC09-GFP. Isolate 12CC09 was included as a control to determine if the fluorescence observed could be attributed to GFP.
Plants were divided into Lm-inoculated and mock-inoculated. Within each inoculation treatment, plants were split between cultivars. The RNA-seq and time-series microscopy experiments were repeated 3 times, as were the experiments involving staining for HP (ROS) with DAB. The protease inhibitor experiments were carried out five times.
For the RNA-seq experiments, within each replicate, there were six seedlings per treatment (Westar or 74-44 BL, mock, or 12CC09-GFP-inoculated), divided at random into two blocks of three plants. At 7 dpi, the cotyledon samples were taken from three of these seedlings and pooled as a replicate for RNA extraction and subsequent RNA-seq. The remaining three seedlings were kept until 14 dpi and rated for infection severity on the 0-9 scale [58,59]. The experiment was conducted three times, resulting in three independent replicates per treatment for the RNA-seq experiment.

RNA Extraction, Library Preparation, and Sequencing
Samples, measuring 5-10 mm × 5-10 mm, were collected from the area adjacent to and containing the lesion on each lobe of the cotyledons at 7 dpi ( Figure 1A). Samples were flash frozen in liquid nitrogen and stored at −80 • C until RNA extraction. Only one of the inoculated lobes was taken from each plant and pooled for the replicate (3 plants). A total of 12 pooled samples from four treatments, three replicates per treatment, were used for RNA extraction.
Cotyledon tissue was ground in liquid nitrogen by vortexing in 50 mL Nalgene Oak Ridge tubes containing two metal beads. RNA was extracted from 40-50 mg of the ground and frozen tissue using the QIAGEN RNeasy Plant Mini Kit on a QIAcube with a DNase I on-column digestion. The concentration and integrity of the resulting RNA was assessed via Nanodrop and Experion (Bio-Rad Canada, Mississauga, ON, Canada) automated electrophoresis, respectively.
Sequencing libraries were prepared using a Illumina ® TruSeq™ RNA Sample Preparation Kit. Each treatment consisted of three biological replicates, with each including samples pooled from three plants. A library was prepared for each replicate, with a total of 12 libraries from the four treatments. These libraries were sequenced on the Illumina HiSeq 2500 using one lane of V4 PE 125 bp at the Genome Quebec Innovation Center (McGill University, Montreal, QC, Canada).

RNA-Seq Data Analysis
Adapter sequences were removed with Trimmomatic (version 0.32) [60]. Subsequently, both paired and unpaired reads were aligned to the B. napus and Lm reference genomes ( [61] and [62], respectively) via STAR (version 2.4.2a) [63]. Gene annotations for Lm and B. napus used the reference genomes from the same sources. Unpaired reads were retained to reduce the risk of introducing bias that might result from discarding these potentially shorter reads. As all samples were treated in the same manner, comparisons between treatments are still valid. Next, the gene models were defined using the GenomicFeatures package in R, and the reads were counted using the R package GenomicAlignments [64]. Differential expression analysis was conducted in R (version 3.3.1 or 3.3.2) using the DESeq2 package for non-normalized read counts (as a count matrix within an assay of a RangedSummarizedExperiment) as required by the package [65]. DESeq2 normalizes samples based on the library size. Due to the importance of the relative changes in expression, statistical significance of these changes were computed between treatments and the overall gene expression (basemean); genes were considered differentially expressed if they had a log-base-2 (log 2 ) fold of change in expression > 2 or < −2, we well as an adjusted p value ≤ 0.05. Differentially expressed genes (DEGs) were scored based on expression (basemean), adjusted p value (adjp) and log 2 -fold changes in expression (basemean): Venn diagrams (Figure 4) were used to identify DEGs that were unique to each combination of contrasting treatments-inoculated Westar versus inoculated 74-44 BL, mock-inoculated Westar versus mock-inoculated 74-44 BL, mock versus Lm-inoculated Westar, and mock versus Lm-inoculated 74-44 BL. DEGs were also subdivided into those with higher expression in the former of the two treatments, being contrasted (positive, Figure 4A) and those upregulated in the latter of the two contrasted treatments (negative, Figure 4B).
KEGG and gene ontology (GO) annotation for DEGs in different groups were extracted from an Ensemble Plant database, using the R biomaRt package. GO enrichment was performed using the enricher function of the R clusterProfiler package. GO annotation and enrichment was also performed using the Blast2Go-pro suite [66]. In this method of GO enrichment, all B. napus genes were searched against the non-redundant protein database from the National Center for Biotechnology Information using BLASTX algorithm with an E-value threshold of 1e −5 . All BLAST hits were mapped onto the GO database to retrieve the terms associated with each hit. Subsequently, all B. napus genes were searched against the InterPro database and annotated by merging the search results from both Blast2GO and InterPro. All expressed B. napus genes were used as a background for the enrichment if they were annotated in GO. The p value and false discovery rate (FDR) cut-off were both set at 0.05.
Volcano plots were performed using the R plot function, based on the results of the DEG analysis, using the R DESeq2 package. All expressed genes were clustered to plot a heatmap using the R pheatmap package, based on gene expression patterns. Z-score transformation was applied prior to heatmap plotting for normalization.

Colorimetric Detection of Hydrogen Peroxide
The area staining for the ROS HP in B. napus cotyledons infected by Lm isolate 12CC09-GFP was measured initially at 7 dpi. Images were collected with a Zeiss Stereo-Lumar epifluorescence microscope, equipped with a NeoLumar S 0.8× objective and an Axiocam 512 camera. Light was provided by a KL-2500 LCD (bright field) or a HBO100 mercury (fluorescent images) bulb. Subsequently, the detached cotyledons were placed in a solution of DAB at room temperature. After 40 (first two experiments) or 90 (third experiment) min, the cotyledons were vacuum infiltrated with DAB for approximately 2 to 3 h. Samples were then boiled in 95% ethanol for approximately 10 to 20 min at 70 • C, to remove the chlorophyll, making the DAB staining more visible. The cotyledons were stored in 95% ethanol, prior to measurement of the area stained brown for HP under a dissecting microscope.
The lesion size, the area colonized by GFP-tagged Lm hyphae, and the area of DAB staining was measured with the aid of ZEN 2 pro or ZEN 2.3 lite (blue edition, © Carl Zeiss Microscopy GmbH, 2011) software, which automatically accounted for magnification.

Time-Series Infection in Cotyledons
In an initial experiment, the inoculated cotyledons were detached from the plant for microscopic examination at 3, 7, 10, and 14 dpi. The experiment was repeated with narrower time intervals of sampling, in order to obtain more detailed time-series data at 3, 5, 7, 9, and 11 dpi. The later experiments also included colorimetric staining for HP. This was done to compliment and confirm colorimetric detection of HP at 7 dpi, as described above. In each experiment, destructive sampling was used at each time point. As the lesion area and the area colonized by the Lm hyphae at 3 dpi were frequently zero, data from this time-point was not included in the analyses.
Bright field images of the top surface of each cotyledon were used to measure the area of the lesion (mm 2 ), the distance from the inoculation point to the most distant edge of the lesion (mm), and the area stained HP (mm 2 ). The area colonized by hyphae (mm 2 ) and the distance from the edge of inoculation wound to the furthest hyphal tip (mm) (first set of time-series experiments only) were quantified using fluorescent images. The Zen Active Contour or Polygon Contour and Length tools were used to collect the area and distance data, respectively.

Assessment of Genomic DNA Degradation as a Marker of Programmed Cell Death
Samples of canola cotyledons were collected as described for the RNA extraction. The samples were freeze dried, and ground to a fine powder in 2 mL tubes, with one 3 mm tungsten carbide bead per tube in a TissueLyser (Qiagen), at room temperature for 5 min at 25 hertz. Genomic DNA was extracted using the QIAGEN DNeasy Plant Mini Kit, according to the manufacturer's instructions. Extracted DNA was diluted to 50 ng/µL. The integrity of the resulting DNA was assessed using Experion DNA 12K analysis kit (Bio-Rad Canada) on an Experion automated electrophoresis system, according to the manufacturer's instructions.

Statistical Analysis
Statistical analyses were done using SAS (version 9.3). Data were assessed for homogeneity of variance and normality, respectively, using Bartlett's Test and the Shapiro-Wilk Test. Data from mock-inoculated plants, which consisted exclusively of zeros, were excluded from statistical analysis.
For ratings of infection severity in parallel with RNA-seq, a randomized complete block design (RCBD) was used, with 3 replicates per treatment and 3 subsamples (plants) per replicate ( Figure 1B). These data were pooled from all plants in a given experiment and y = log 10 (x + 10) transformed. Means were compared via a t-test, using Proc GLM. For the time-series experiments that did not include HP assessment ( Figure 1E), data were y = log 10 (x + 10) transformation. Each parameter (lesion area, area colonized by Lm hyphae, distance to the edge of lesion, or to furthest hyphal tips from the inoculation site) and time-point (7, 10 or 14 dpi) were analyzed individually using t-test to compare the means between cultivars. As much of these data did not meet all requirements for ANOVA, they were analyzed using the Wilcoxon two-sample (for cultivar) or Kruskal-Wallis multi-sample (for comparisons between parameters) test, using the Proc NPAR1WAY. When significant differences were found for a given factor (cultivar or parameter), the means of treatments were separated using Tukey's adjusted multiple comparison. For the protease inhibitor experiment (Figure 9), a two-factor (cultivar and protease inhibitor) analysis was performed, and Tukey's adjustment was also used for mean comparison.

Conclusions
QR observed in 74-44 BL cotyledons appeared to result from restricted Lm colonization due, at least in part, to ROS, which might trigger a rapid PCD via signal transduction and endomembrane vesicle trafficking. The upregulation of genes related to a putative Bax inhibitor-1 might reduce the PCD response in the susceptible Westar. The restriction of cotyledon colonization by Lm could be of significant benefit to blackleg resistance as it might limit or prevent the pathogen from spreading into the stem, where the most damaging form of the disease takes place. This study served as a prelude to further research into mechanisms involved in QR against blackleg. Such work, in conjunction with studies of more resistant canola cultivars, might facilitate judicious deployment of QR traits for improved blackleg management.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2223-7747/9/7/864/s1, Figure S1: Differential reaction of Westar and 74-44 BL cotyledons to inoculation with L. maculans isolates carrying AvrLm1 (S7) or AvrLm3 (19.4.24), or lacking AvrLm1 and AvrLm3 (12CC-09), Table S1: RNA-seq read mapping, Table S2. Genes of Leptosphaeria maculans that are more highly expressed in inoculated 74-44 BL than in inoculated Westar, Table S3. Genes of Leptosphaeria maculans that are more highly expressed in inoculated Westar than in inoculated 74-44 BL, Table S4. Genes of Brassica napus that are more highly expressed in Lm-inoculated cotyledons of 74-44 BL than those of Westar (score ≤ −2 × 10 4 ), Table S5. Gene Ontology (GO) term enrichment of differentially expressed genes (DEGs) in Brassica napus that are more highly expressed in Lm-inoculated 74-44 BL than in inoculated Westar, Table S6. KEGG annotation of differentially expressed genes (DEGs) in Brassica napus that are more highly expressed in Lm-inoculated 74-44 BL than in inoculated Westar (scores > 1.0 × 10 3 ), Table S7. Genes of Brassica napus that are more highly expressed in Lm-inoculated cotyledons of Westar than in those of 74-44 BL (scores > 2 × 10 4 ), Table S8. GO term enrichment of DEGs in Brassica napus that are more highly expressed in Lm-inoculated Westar than in inoculated 74-44 BL, Table S9. KEGG annotation of differentially expressed genes (DEGs) in Brassica napus that are more highly expressed in Lm-inoculated Westar than in inoculated 74-44 BL (scores > 1.0 × 10 3 ), Table S10. Genes of Brassica napus that are more highly expressed in mock-inoculated cotyledons of 74-44 BL than in those of mock Westar (scores < 2 × 10 4 ), Table S11. GO term enrichment of DEGs in Brassica napus that are more highly expressed in mock-inoculated 74-44 BL than in mock Westar, Table S12. KEGG annotation of differentially expressed genes (DEGs) in Brassica napus that are more highly expressed in mock-inoculated 74-44 BL than in mock-inoculated Westar (scores > 1.0 × 10 3 ), Table S13. Genes of Brassica napus that are more highly expressed in mock-inoculated cotyledons of Westar than in those of 74-44 BL (scores > 8.6 × 10 3 ), Table S14. GO term enrichment of DEGs in Brassica napus that are more highly expressed in mock-inoculated Westar than in mock 74-44 BL, Table S15. KEGG annotation of differentially expressed genes (DEGs) in Brassica napus that are more highly expressed in mock-inoculated Westar than in mock 74-44 BL (scores >1.0 × 10 3 ). Funding: This study was supported by the Canola Agronomic Research Program Project No. CARP 2015-12 entitled "Understanding the mechanisms for race-specific and non-specific resistance for effective use of cultivar resistance against blackleg of canola in Western Canada" administered by Canola Council of Canada.