Comparative Transcriptomic Analysis of Virulence Factors in Leptosphaeria maculans during Compatible and Incompatible Interactions with Canola

Leptosphaeria maculans is a hemibiotrophic fungus that causes blackleg of canola (Brassica napus), one of the most devastating diseases of this crop. In the present study, transcriptome profiling of L. maculans was performed in an effort to understand and define the pathogenicity genes that govern both the biotrophic and the necrotrophic phase of the fungus, as well as those that separate a compatible from an incompatible interaction. For this purpose, comparative RNA-seq analyses were performed on L. maculans isolate D5 at four different time points following inoculation on susceptible cultivar Topas-DH16516 or resistant introgression line Topas-Rlm2. Analysis of 1.6 billion Illumina reads readily identified differentially expressed genes that were over represented by candidate secretory effector proteins, CAZymes, and other pathogenicity genes. Comparisons between the compatible and incompatible interactions led to the identification of 28 effector proteins whose chronology and level of expression suggested a role in the establishment and maintenance of biotrophy with the plant. These included all known Avr genes of isolate D5 along with eight newly characterized effectors. In addition, another 15 effector proteins were found to be exclusively expressed during the necrotrophic phase of the fungus, which supports the concept that L. maculans has a separate and distinct arsenal contributing to each phase. As for CAZymes, they were often highly expressed at 3 dpi but with no difference in expression between the compatible and incompatible interactions, indicating that other factors were necessary to determine the outcome of the interaction. However, their significantly higher expression at 11 dpi in the compatible interaction confirmed that they contributed to the necrotrophic phase of the fungus. A notable exception was LysM genes whose high expression was singularly observed on the susceptible host at 7 dpi. In the case of TFs, their higher expression at 7 and 11 dpi on susceptible Topas support an important role in regulating the genes involved in the different pathogenic phases of L. maculans. In conclusion, comparison of the transcriptome of L. maculans during compatible and incompatible interactions has led to the identification of key pathogenicity genes that regulate not only the fate of the interaction but also lifestyle transitions of the fungus.

Leptosphaeria maculans is a hemibiotrophic fungus that causes blackleg of canola (Brassica napus), one of the most devastating diseases of this crop. In the present study, transcriptome profiling of L. maculans was performed in an effort to understand and define the pathogenicity genes that govern both the biotrophic and the necrotrophic phase of the fungus, as well as those that separate a compatible from an incompatible interaction. For this purpose, comparative RNA-seq analyses were performed on L. maculans isolate D5 at four different time points following inoculation on susceptible cultivar Topas-DH16516 or resistant introgression line Topas-Rlm2. Analysis of 1.6 billion Illumina reads readily identified differentially expressed genes that were over represented by candidate secretory effector proteins, CAZymes, and other pathogenicity genes. Comparisons between the compatible and incompatible interactions led to the identification of 28 effector proteins whose chronology and level of expression suggested a role in the establishment and maintenance of biotrophy with the plant. These included all known Avr genes of isolate D5 along with eight newly characterized effectors. In addition, another 15 effector proteins were found to be exclusively expressed during the necrotrophic phase of the fungus, which supports the concept that L. maculans has a separate and distinct arsenal contributing to each phase. As for CAZymes, they were often highly expressed at 3 dpi but with no difference in expression between the compatible and incompatible interactions, indicating that other factors were necessary to determine the outcome of the interaction. However, their significantly higher expression at 11 dpi in the compatible interaction confirmed that they contributed to the necrotrophic phase of the fungus. A notable exception was LysM genes whose high expression was singularly observed on the susceptible host at 7 dpi. In the case of TFs, their higher expression at 7 and 11 dpi on susceptible Topas support an important role in regulating

INTRODUCTION
Blackleg disease (stem canker) caused by Leptosphaeria maculans (Desm.) Ces. & De Not. is one of the major constraints to canola (Brassica napus L.) production worldwide (Fitt et al., 2006). Infection by the fungus is known to cause more than 50% yield losses in canola (Kutcher et al., 2013). The major difficulty for combating the pathogen lies in the understanding of its complex lifestyle, which includes alternative biotrophic, and necrotrophic, phases, along with a symptomless endophytic phase (Howlett et al., 2001;Van de Wouw et al., 2016). Management of blackleg disease includes crop rotations, seed treatment and fungicide applications, and preferably, disease-resistant cultivars, arguably the most effective approach (Delourme et al., 2006).
Canola shows two types of resistance against L. maculans: qualitative and quantitative. Single or few major genes that are known to be involved in a gene for gene interaction govern qualitative resistance. Major genes provide resistance particularly at the seedling stage whereas quantitative resistance involves many small effect genes that are mostly expressed during the adult plant stage (Raman et al., 2012;Huang et al., 2014). To date, several resistance genes with major effect have been identified in Brassica species, but only two, LepR3 and Rlm2, have been cloned and well characterized (Delourme et al., 2006;Long et al., 2011;Larkan et al., 2013Larkan et al., , 2015Van de Wouw et al., 2014). On the other hand, better progress has been achieved with L. maculans where 14 avirulence genes have been identified, and seven of them, namely AvrLm1, AvrLm4-7, AvrLm6, AvrLm11, AvrLmJ1, AvrLm2, and AvrLm3 have been cloned (Gout et al., 2006;Fudal et al., 2007;Parlange et al., 2009;Balesdent et al., 2013;Van de Wouw et al., 2014Ghanbarnia et al., 2015;Plissonneau et al., 2016). Interestingly, some of these avirulence genes have been found to be clustered, with clusters AvrLm1-2-6 and AvrLm3-4-7-9-AvrLepR1 being the notable examples (Balesdent et al., 2002;Ghanbarnia et al., 2012). For the most part, avirulence genes, including L. maculans Avrs, are small-secreted proteins (SSPs) with several cysteine residues, and often referred to as effectors (Stukenbrock and McDonald, 2009;Rouxel et al., 2011).
Effectors are key elements in fungal virulence against plants and particularly important during the biotrophic phase of infection (Kloppholz et al., 2011). L. maculans, being a hemibiotroph, will initially rely on effectors to suppress plant defenses, and then will subsequently use effectors to kill plant cells. In L. maculans, most putative or candidate effectors are localized in transposon-rich repetitive DNA and are affected by a repeat-induced point mutation (Rouxel et al., 2011). The putative effector genes are mostly over-expressed during primary leaf infection (Soyer et al., 2014). Such information about the genomic localization, gene organization, and expression dynamics is helpful to understand the host-pathogen interaction and more particularly for the identification of bona fide effectors. Similarly, transcription factors (TFs), and carbohydrate active enzymes (CAZymes) are known to play a pivotal role in hostpathogen interactions, and are, along with effectors, prime targets for studying virulence factors in fungi (Guo et al., 2011;Lombard et al., 2014;Lowe et al., 2014;Malinovsky et al., 2014).
Transcription factors are essential players in the signal transduction pathways. In L. maculans, TF LmStuA is found to be required for normal growth, perithecium formation, pathogenicity on oilseed rape leaves, and expression of effectors (Soyer et al., 2015). The silencing of LmStuA triggers drastic effects on the morphogenesis and pathogenicity of L. maculans, indicating that it may affect a large number of genes and pathways (Soyer et al., 2015). Similarly, several CAZymes in L. maculans genome have been predicted to have a functional role in pathogenesis . CAZymes are important to break down the polysaccharides of plant cell walls, to establish infection, and also, to facilitate access to nutrients during the necrotrophic and saprophytic growth phases. For instance, global transcriptomic analyses of the hemibiotroph Colletotrichum higginsianum revealed that genes encoding secreted proteins without a functional annotation are expressed predominantly during the initial biotrophic phase, whereas expression of secreted lytic enzymes (including CAZymes) was higher in the subsequent necrotrophic phase (O'Connell et al., 2012). A similar finding was observed in Leptosphaeria biglobosa, a necrotroph expressing more cell wall degrading genes than L. maculans . However, L. maculans expressed many genes in the carbohydrate binding module (CBM) class of CAZymes, particularly CBM50 genes, during early infection, and cell wall degrading enzymes at later stages of growth . This suggests that expression of secreted proteins without functional annotation is a general feature of biotrophy, whereas expression of cell wall degrading enzymes is generally associated with necrotrophy. Other important necrotrophy-related genes code for sirodesmin PL (Sir), a phytotoxin that belongs to the class of epipolythiodioxopiperazine (ETP). The production of sirodesmin by L. maculans is thought to be suppressed by brassinin, a phytoalexin of canola (Pedras et al., 1993). In L. maculans, a cluster of 23 genes including polyketide synthase (PKS), non-ribosomal peptide synthase (NRPs) genes, and 18 Sir genes have been identified (Gardiner et al., 2004).
The infection process is highly dependent on host recognition and molecular cross-talk between the host and the pathogen where pathogenicity-related genes play an important role.
However, any given host-pathogen interaction is a very complex phenomenon, which makes it difficult to understand the factors dictating compatibility or incompatibility. In Arabidopsis, Huibers et al. (2009) were able to discriminate genes induced  during compatible from incompatible interactions with the  downy mildew pathogen Hyaloperonospora arabidopsidis (downy  mildew). Similarly, other studies have been conducted to compare gene expression profiling under compatible and incompatible interactions (Wang et al., 2010;Sestili et al., 2011;Li et al., 2015).
Expression profiling during the development of disease is an effective approach to better understand the pathogenesis process. Current improvements in sequencing technologies have provided new opportunities to evaluate gene expression under different conditions by sequencing the entire transcriptome (Wang et al., 2009). Compared to other hybridization based transcriptome profiling platforms like microarrays, RNA-seq provides access to simultaneous transcript discovery and abundance estimation, identification of differentially-expressed genes (DEGs) and associated molecular cellular pathways, and alternative splicing variants (Wang et al., 2009;Trapnell et al., 2012). However, RNA-seq analysis requires scalable, fast, and statistically relevant software tools that can handle complex and large sequence data. Fortunately, considerable efforts have been devoted to the development of specialized software tools to perform effective RNA-seq analysis (Garber et al., 2011;Trapnell et al., 2012;Seyednasrollah et al., 2015;Sonah et al., 2016).
The L. maculans-canola compatible interaction has been previously addressed using the RNA-seq approach Haddadi et al., 2016). These studies made significant efforts toward the understanding of susceptible reaction in canola after L. maculans infection and provides the first indepth look into the transcriptomic profile of this interaction. In this study, our objective was to compare the pathogen responses against susceptible and resistant host genotypes in order to obtain a precise definition of the virulence factors expressed by L. maculans during its biotrophic and necrotrophic phase. For this purpose, we performed the RNA-seq transcriptome profiling of L. maculans inoculated on susceptible and resistant canola lines at four developmental stages over five biological replications. A particular emphasis was placed on the identification of DEGs including effectors, CAZymes, and other pathogenesis-related genes during blackleg disease development in canola.

Plant Material and L. maculans Inoculation
Canola (B. napus) breeding lines, Topas DH16516 (Topas-wild), a double haploid line susceptible to L. maculans and Topas-Rlm2, an introgression line resistant to L. maculans isolates carrying AvrLm2 (Larkan et al., 2016a,b) were used as plant material. Seven-day old seedlings of both resistant and susceptible canola lines grown under controlled environment were point inoculated with pycnidiospores suspension of L. maculans isolate D5. L. maculans isolate D5 contains known avirulence effectors AvrLm1, AvrLm4-7, AvrLm2, and AvrLmJ1 but lacks AvrLm6 (Raman et al., 2012). Five biological replicates each with leaves from eight plants per sample were collected at 3, 5, 7, and 11 days post inoculation (dpi). The L. maculans isolate D5 was cultured in V8 R agar plates at 25 • C. For RNA extraction, two fungal agar plugs were inoculated in V8 R liquid media and incubated at 22 • C. After 7 days, the mycelium was harvested for RNA extraction.

Sample Preparation and Illumina Sequencing
Total RNA for all samples was extracted using a combined Trizol/Qiagen RNeasy mini kit (Qiagen, Mississauga, ON, Canada). The RNA quality and quantity was accessed using NanoDrop ND-1000 spectrophotometer (NanoDrop technologies) and further verified by agarose gel electrophoresis. Four micrograms of total RNA from all samples were used to make individual barcoded cDNA library using TruSeq RNA library preparation kit v2 with some modifications. Individual library was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies). The products resulted in a smear with an average fragment size of approximately 260 bp. A total of six individual libraries were pooled using the uniform amount of each library and the quality of the final library pools was also assessed. Single end, 100 bp sequencing was performed using an Illumina HiSeq 2000 sequencer at the Génome Québec Innovation Centre (McGill University, Quebec, Canada).

RNA-seq Data Analysis
RNA-seq reads were quality-checked with fastqc, which performs various quality checks for the raw reads. Read processing was performed by using Trimmomatic software (Bolger et al., 2014). Reads were trimmed from both ends until the average of all 5 bp sliding windows reached a Phred score of 25 or higher and all the sequences shorter than 35 bases were discarded. Processed reads were aligned to the L. maculans genome and transcriptome with Tophat2 (Trapnell et al., 2009). Most of the parameters in Tophat was set as default. A mismatch of two bases were allowed for the alignment. The minimum and maximum intron length was set to 50 and 500,000 respectively. Reads aligned to multiple sites were removed prior to further analysis.
Novel transcripts that did not overlap with any annotated transcripts were identified using Cufflink tools. The gene expression level for each annotated as well as non-annotated novel transcripts were estimated as the number of Fragments (reads) per kilobase of transcript per million mapped reads (FPKM), considering only uniquely mapped reads in exonic region by using the Cufflink software. The differentially expressed genes (DEGs) were identified using four different tools including Cuffdiff (Trapnell et al., 2012), EdgeR (Robinson et al., 2010), DESeq2 (Love et al., 2014), and CLC Genomics Workbench. We used FDR < 0.0001 and the absolute value of log 2 (Foldchange) >1.5 relative to axenic culture as the threshold for the identification of DEGs. We used quartile normalization that excludes the top 25% of expressed genes to improve detection of less abundant genes. With the Cufflink software, we used -M option to mask rRNA, -b, and -u option for bias correction and option to normalize for aligned tags instead of total tags. HTseq tool was used to count reads prior to DEGs identification. The DESeq2 package was used to estimate sample quality (PCA) and the expression level of transcripts. The regularized rlog transformation and variance stabilizing transformation were used for data visualization. Time course analysis for all the four time points were carried out to find the genes that reacted in a time-specific manner using DESeq2 package.

Functional Annotation and Gene Ontology
Standard gene ontology (GO) was used to describe DE gene functionality, a hypergeometric test and the p < 0.05 of Pearson Chi-Square test between the gene numbers of the two input dataset were used to map the DT genes to GO terms based on the BGI WEGO (Web Gene Ontology Annotation Plot, http://wego.genomics.org.cn/cgi-bin/wego/index.pl). Single enrichment analysis (SEA), a function of AgriGO was used to examine GO term enrichment using Magnaporthe grisea as a background reference using Fisher statistical test and 0.05 significance level and other default parameters. Transcription factors were identified using Fungal Transcription Factor Database (FTFD) (Park et al., 2008). The FTFD pipeline sorts fungal TFs initially based on the relevant InterPro terms like DNA-binding motifs, and then the false-positive TFs are filtered with different criteria. CAZymes were identified using dbCAN server (Yin et al., 2012). The dbCAN hosted an analytical pipeline compiled with CDD (conserved domain database) search, family specific hidden Markov model and literature curation (Yin et al., 2012). For the classification of putative secreted peptidases, the sequences for the secreted proteins predicted by WoLF PSORT (cutoff score = 15) were submitted to MEROPS Batch Blast analysis (http://merops.sanger.ac.uk, Rawlings et al., 2014).

Identification of Candidate Secreted Effector Proteins (CSEPs)
Predicted protein sequences from the L. maculans genome were retrieved from the JGI MycoCosm (http://genome.jgi.doe.gov/ Lepmu1/Lepmu1.home.html) (Rouxel et al., 2011;Grigoriev et al., 2014). SignalP (cutoff probability = 0.8), TargetP (cutoff probability = 0.8), Psort (cutoff score = 15), BlastP (cutoff e = 1e −05 ), and TMHMM software tools along with Secretool pipeline were used to predict small-secretory proteins (Petersen et al., 2011;Cortázar et al., 2014). To prioritize candidate effector genes, the entire secretome was analyzed by EffectorP a machine learning method optimized for the prediction of fungal effectors (Sperschneider et al., 2015). Crinkler type effectors were also searched in the L. maculans genome with the FEMO software tool implemented in MEME suit by using conserved domain LFLAK (http://meme-suite.org/tools/fimo). Proteins with LFLAK-domain within the initial 100 AA were considered for further analysis. To verify the search parameters, Crinkler effector search was also performed in Phytophthora sojae genome and compared with earlier report by Haas et al. (2009). After an initial search with FIMO, candidate Crinkler effectors were analyzed with secretool pipeline and subsequently with effectorP software tool.

Disease Progress on Compatible and Incompatible Host
On plants of either the compatible host Topas-wild or the introgression line Topas-Rlm2 carrying a major resistance gene, disease symptoms were not observed until 7 dpi. After this asymptomatic early growth stage, symptoms became visible exclusively on Topas-wild plants expanding into clear lesions and zones of chlorotic tissues at 11 dpi. In the case of Topas-Rlm2 plants, there was no visible lesion beyond the site of inoculation (Figure 1). The disease progress on susceptible Topas-wild cotyledons caused by L. maculans was similar to that previously described on Westar .

Transcriptome Sequencing with the Time Course of Disease Progress
A total of 1.6 billion single-end reads consisting of approximately an average of 33 million reads for each cDNA library were obtained. The entire data set was submitted to SRA NCBI database and can be accessed with accession number SRP078092. The raw reads obtained were uniform across the different libraries. After performing quality assessment and read processing, about 0.1% of the reads with poor quality and shorter length were discarded. Mapping of the processed reads to the L. maculans genome showed a very high percentage of mapping for the axenic samples (Supplementary Table 1). The percentage of reads from in planta samples mapped to the L. maculans reference genome increased over time from <0.5% at 3 dpi to about 12% at 11 dpi in compatible host Topaswild (Figure 2, Supplementary Table 1). By contrast, mapped fungal reads remain stable in the incompatible interaction with Topas-Rlm2 and were well below 1% even at 11 dpi (Figure 2, Supplementary Table 1), a result well in line with the phenotype observed in Figure 1. Comparing the percent of mapped reads in the early asymptomatic growth stage in Topaswild plants, only a marginal increase was observed between 3 and 5 dpi. This percentage more than doubled between 5 and 7 dpi and registered a 10-fold increase from 7 to 11 dpi (Supplementary Table 1).
FIGURE 1 | Disease symptoms at 11 days post inoculation on leaves of compatible (Topas-wild) and incompatible (Topas-Rlm2) canola host inoculated with Leptosphaeria maculans D5 isolate or water.
Frontiers in Plant Science | www.frontiersin.org The principal components (PCs) analysis highlighted a clear differential effect of the treatments along with the uniformity of the five biological replications within a treatment (Figure 3). The first principal component explained 58% of the expression FIGURE 2 | Percentage of mapped reads to the Leptosphaeria maculans genome at 3, 5, 7, and 11 days post inoculation (dpi). RNA-seq was performed for L. maculans isolate D5 inoculated to Topas-wild (compatible) and Topas-Rlm2 (incompatible) canola genotypes.
variation supporting the large phenotypic differences between the conditions over the time period.

Comparison of Gene Expression Profiling during In vitro and In planta Growth of L. maculans
The number of DEGs in L. maculans between in vitro axenic samples and samples at different in planta growth stages showed an interesting pattern of gene expression turnover over time. The overall pattern of DEGs identified with different software tools was similar (Supplementary Figure 1). The highest number of DEGs for both the compatible (Topas-wild) and incompatible (Topas-Rlm2) interactions was recorded in the early events (3 dpi) and was similar in both cases (Figure 4). This number was reduced by around three-fold at 5 dpi in both interactions and remained fairly level over the next sampling times in Topas-Rlm2 plants. On the other hand, it increased steadily in Topas-wild plants to exceed by roughly three times the number of DEGs found in Topas-Rlm2 plants (Figure 4).
When looking at the top 20 upregulated DEGs at 3 dpi in the compatible interaction (Topas-wild), most of them were linked to uncharacterized proteins for which functional annotation was not available (Table 1). Since the uncharacterized proteins are unique to L. maculans, homology-based annotation failed to FIGURE 3 | Principal components (PCs) of RNA-seq expression profile obtained for Leptosphaeria maculans in planta samples during compatible and incompatible interactions. A total of five biological replications for 3, 5, 7, and 11 days post inoculation (dpi) were used for RNA-seq expression profiling of L. maculans inoculated to compatible host Topas-wild (blue) and incompatible host Topas-Rlm2 (Red).

FIGURE 4 | Differentially expressed
Leptosphaeria maculans genes identified at 3, 5, 7, and 11 dpi during compatible and incompatible interactions. Differentially expressed genes were identified by comparing in planta samples with axenic cultures. Analyses were performed with five biological replicates and a threshold value of FDR < 0.0001 and log 2 fold-change > 1.5. characterize them. Genes involved in chitin binding, fasciclin, and related adhesion glycoproteins known to play roles in early stages of infection were found to be highly upregulated (Table 1). Interestingly, of the top 20 upregulated genes of L. maculans identified during the compatible interaction with Topas-wild, 18 were also found associated with the incompatible interaction (Topas-Rlm2) ( Table 1, Supplementary Table 2). However, the number of common genes gradually decreased over time to 11, 6, and 4 at 5, 7, and 11 dpi, respectively ( Table 1, Supplementary  Table 2). Time series analysis showed high expression of Avr genes during the biotrophic phase at 7 dpi while most of the highly expressed genes during the necrotrophic phase at 11 dpi were associated with molecular functions involved in catalase activity, hydrolases, CAZymes, peptidases, and transporters (Table 1, Figure 5).

Expression of CSEPs in L. maculans during Compatible and Incompatible Interactions
A total of 552 classically secreted proteins were identified using computational pipeline comprising SignalP and TMHMM software tools along with Secretool pipeline. Following further analyses, 134 genes were prioritized as high confidence CSEPs based on the results obtained with EffectorP software. In the case of the compatible interaction with Topas-wild, an important increase in upregulated CSEPs was observed at 7 and 11 dpi, while very limited differences over time were noted in the incompatible interaction (Figure 6).
Out of the 134 CSEPs, 35 genes were not expressed in either the compatible or incompatible host. However, a total of 28 genes showed the highest expression level at 7 dpi in the compatible interaction, which would link them to the biotrophic phase of the fungus ( Table 2). This profile was similar to that of known Avr genes, which were mostly found to be highly expressed at 7 dpi in the compatible interaction, with the obvious exception of AvrLm6 that remained unchanged throughout the sampling periods (Figures 5, 7). The expression patterns of AvrLm1, AvrLm2, AvrLm4-7, AvrLm11, and AvrLmJ1 were similar in the compatible interaction as their expression increased over time to reach a peak at 7 dpi, and subsided at 11 dpi. As for the incompatible interaction, limited expression was observed for all known Avr genes with no distinctive pattern (Table 2, Figure 7). Of additional significance, an expression pattern similar to Avr genes was found for 23 other effectors in the compatible interaction, including eight specific to L. maculans and not previously reported as effectors ( Table 2). Of these, gene_6114 and gene_2728 are particularly interesting for their differentially higher expression at 7 dpi in Topas-wild. Finally, gene_9004 showed the highest expression at 7 dpi during the compatible interaction at a level of almost 24,000 FPKM.
In addition to the effectors listed in Table 2, 15 other genes identified as CSEPs showed higher levels of expression exclusively during the compatible interaction at 11 dpi when the fungus had entered its necrotrophic phase (Table 3). Furthermore, the crinkler effector search performed initially in P. sojae identified exactly the same set of proteins reported earlier by Haas et al. (2009) (Supplementary Table 3). This validated the method used here to identify crinkler effectors in L. maculans genome. The initial search identified 63 proteins with LFLAK-like domain present within the initial 100 AA (Supplementary Table 4). Further analysis with secretool revealed only five proteins and EffectorP confirmed only one candidate (Gene_2728, Lema_T080290.1) as a crinkler effector (Gene_2728, Lema_T080290.1). Gene_2728 showed the highest expression at 7 dpi during the compatible interaction (Supplementary Figure 2).

Transcription Factors Associated with Effector Expression
TFs regulate the expression of a number of genes simultaneously, and their upregulation during in vivo conditions was indicative of their importance in the outcome of the interaction. Both compatible and incompatible interactions were regulated by a large number of upregulated TFs in L. maculans in the early events (3 dpi) followed by a sharp decline at 5 dpi ( Figure 8). As with the phenotypes, a clear distinction between compatible and incompatible interactions appeared at 7 dpi where nearly 100 TFs were upregulated in the former case compared to 40 in the latter, and these differences were even more manifest at 11 dpi (Supplementary Table 5). Among TFs, the LmStuA TF (gene_1191), a member of APSES domain containing proteins was found to be strongly upregulated at 7 dpi in the compatible interaction with Topas-wild, thus suggesting its importance in the infection process. Mostly Zn2Cys6 TFs were upregulated during the compatible interaction (Supplementary Table 5).

Expression of CAZymes in L. maculans during Compatible and Incompatible Interactions
Plant-fungal interactions involve a variety of CAZymes that are used by the fungal pathogen to infect its host. RNA-seq data for L. maculans infecting compatible host Topas-wild and incompatible host Topas-Rlm2 showed a high number of   differentially expressed CAZymes at 3 dpi (Figure 9). At 5 dpi, this number dropped drastically. It remained fairly constant thereafter in the case of the incompatible interaction (Figure 9). By contrast, we observed a rapid increase in differentially expressed CAZymes from 5 to 7 dpi in the compatible interaction, and this number exceeded 200 at 11 dpi, supporting their role in the necrotrophic phase of the fungus. Among the differentially expressed CAZymes, glycosyl hydrolase (GH), carbohydrate esterase (CE) domain containing genes were more prevalent followed by those belonging to auxiliary activities (AAs) class (Supplementary Table 6). During the early stages (3 dpi) of L. maculans infection, cellulose and pectin-degrading enzymes such as PL3, CE4, GH43, and GH3 were upregulated in both interactions compared to axenic growth. At the later stage (11 dpi), CE1, PL1, GT34, and GH28 were only upregulated in the compatible interaction. Carbohydrate binding molecules (CBM) were another prevalent group among DEGs. In L. maculans, enzymes with LysM motifs are a well-studied class of CAZymes. Among the LysM domain containing genes present in L. maculans genome, Lm2LysM (Gene_4592) and Lm5LysM (Gene_11099) genes showed a higher level of expression at 7 dpi during the compatible interaction (Supplementary Figure 3). However, Lm4LysM (Gene_7646) gene was not expressed during either the compatible or incompatible interaction. The expression pattern for some of the CAZymes at different growth stages of L. maculans during compatible and incompatible interaction is shown in Supplementary Figure 4.

Expression of Important Peptidases, Secondary Metabolites, and Necrosis-Inducing Proteins (NIPs)
At the early stages of infection, a similar number of differentially expressed peptidases were observed in both interactions. However, as the fungus transitioned from the biotrophic phase at 7 dpi to the necrotrophic phase at 11 dpi, the difference in the number of differentially expressed peptidases between the compatible and incompatible interaction steadily increases to a ratio of over 60 to 1 at 11 dpi (Figure 10, Supplementary Table 7). Among the different classes of peptidases, serine proteases and carboxypeptidases were highly upregulated during L. maculans compatible interaction (Figure 11, Supplementary Table 7).
In L. maculans, a cluster of 23 genes including PKS, NRPs genes and 18 Sir genes have been identified (Gardiner et al., 2004). Most of these genes showed a similar pattern of expression at all the stages of infection with a higher expression at 11 dpi compared to early growth stages of L. maculans in planta (Supplementary Figure 5A). By contrast, RNA-seq analysis was performed for L. maculans isolate D5 inoculated to Topas-wild (compatible) and Topas-Rlm2 (incompatible) canola genotypes. *Analyses were performed with five biological replicates. dpi, days post inoculation; FPKM, Fragments per kilo-base of transcript per million mapped reads; 1, (Haddadi et al., 2016); 2, .
Sir O (oxidoreductase) and Sir T (thioredoxin reductase) showed higher level of expression at both 3 dpi and 11 dpi (Supplementary Figure 5A). However, comparative analyses indicated that there was no difference in expression at 3 dpi between the interactions (Supplementary Figure 5B). Most of the NRPs showed high level of expression at 3 dpi and 11 dpi but NRP gene_904 showed high expression at 7 dpi (FPKM 876) that receded at 11 dpi (FPKM 112.2). For the two necrosisand ethylene inducing proteins (Neps) in L. maculans, Nep1like proteins (NPP1, gene_11090) showed higher expression during the pathogenic phases of the fungus (Supplementary Table 8).

Functional Categorization and Gene Ontology (GO) Enrichment of the DEGs
Differential molecular responses in L. maculans under compatible and incompatible interactions were observed based on the functional categorization of DEGs (Figure 11,  Supplementary Figure 6). The DEGs were grouped into three main classes, molecular function, cellular component, and biological process, and several subclasses based on gene ontology terms assigned using WEGO tool. During early asymptomatic stage at 3 dpi, upregulated genes with similar functional categories were observed under both conditions (Supplementary Figure 6). However, the molecular responses of L. maculans between the compatible and incompatible hosts differed significantly over time, and were more prominent at 7 and 11 dpi (Figure 11; Supplementary  Figure 6). The most striking differences were found in the number of genes expressed in each sub-category and in the prevalence of functional categories linked to a transition between the biotrophic and necrotrophic phases such as carboxypeptidases, ligase, peroxidases, cell-wall biogenesis, and stimuli (Figure 11).  Hierarchical clustering and GO enrichment of DEGs during the compatible interaction was also performed using AgriGO tool. Most of the DEGs were found to belong to binding and catalytic functions (Supplementary Figure 7). Localization (GO:0051179, 3.21e-62), cellular process (GO:0009987, 2.37e-258), signal transduction (GO:0007165, 3.48e-13), response to stress (GO:0006950, 8.44e-11), and transcription regulator activity (GO:0030528, 3.3e-15), were the most significantly enriched GO terms at 3 and 5 dpi. At 7 dpi, localization (GO:0051179, 2.35e-61), transport (GO:0006810, 1.68e-61), catabolic process (GO:0009056, 1.03e-11), metabolic process (GO:0008152, 0), gene expression (GO:0010467, 3.93e-38), binding (GO:0005488, 6.84e-319), and transferase activity (GO:0016741, 8.4e-09) were linked to the highest number of genes, and at 11 dpi, hydrolytic activity, mostly monooxygenase activity (GO:0051179, 5.86e-09), transcription regulator activity (GO:0030528, 2.49e-21), transcription factor activity (GO:0051179, 3.14e-14), and regulation of biological process (GO:0050789, 7.09e-26) were the most common enriched GO terms (Supplementary Figure 7).

FIGURE 8 | Number of transcription factors (TFs) upregulated in
Leptosphaeria maculans at 3, 5, 7, and 11 days post inoculation (dpi) on compatible host Topas-wild and incompatible host Topas-Rlm2. Upregulated genes were identified by comparing in planta samples at different growth stages with axenic cultures. Analyses were performed with five biological replicates and a threshold value of FDR < 0.0001 and log 2 fold-change > 1.5.
FIGURE 9 | Number of carbohydrate active enzymes (CAZymes) upregulated at 3, 5, 7, and 11 days post inoculation (dpi) of Leptosphaeria maculans inoculated to compatible host Topas-wild and in-compatible host Topas-Rlm2. Upregulated genes were identified by comparing in planta samples with axenic cultures. Analyses were performed with five biological replicates and a threshold value of FDR < 0.0001 and log 2 fold-change > 1.5.

DISCUSSION
In the present study, we performed a comprehensive analysis of the L. maculans transcriptome profile during compatible and incompatible interactions with canola. Based on comparative analyses, key genes that dictate both the interaction between canola and L. maculans and the different pathogenic stages of the fungus were highlighted. Among the genes of particular significance, our results have identified candidate effectors, TFs, CAZymes, peptidases, and other pathogenesis-related genes that are specifically upregulated as L. maculans initiates a biotrophic interaction with the plant, and transitions to a necrotrophic phase. The differential expression of genes during compatible FIGURE 10 | Number of upregulated peptidases observed at 3, 5, 7, and 11 days post inoculation (dpi) of Leptosphaeria maculans inoculated to compatible host Topas-wild and in-compatible host Topas-Rlm2. Upregulated genes were identified by comparing in planta samples with axenic cultures. Analyses were performed with five biological replicates and a threshold value of FDR < 0.0001 and log 2 fold-change > 1.5. and incompatible interactions offers a precise insight into the mechanisms of pathogenesis in the L. maculans-canola interactions.
The effectiveness of transcriptome analyses in plant-pathogen interactions depends primarily on the approach of expression quantification, statistical methods avoiding possible errors and reliance on normalized comparisons, experimental design with sufficient replications and appropriate plant or pathogen material to address relevant biological questions (Williams et al., 2014). In this work, introgression lines Topas-wild and Topas-Rlm2 were used to study L. maculans molecular responses during disease development. Topas-wild is a common cultivar from Canada well known and exploited for its susceptibility to L. maculans isolate D5 (Larkan et al., 2013). On the other hand, Topas-Rlm2 is a recently developed cultivar that carries Rlm2, a major resistance gene that prevents infection from L. maculans isolate D5 . As such, this cultivar provided a unique opportunity to investigate the subtle elements that distinguish the ability of L. maculans to infect or not its host. The importance of obtaining good reproducible phenotypes for transcriptomic analyses cannot be overstated as it remains the reference basis for all analyses. On the basis of visual observations, this condition was clearly met as Topas-wild plants exhibited clear symptoms of infection over the course of the experiment that culminated with the presence of extensive necrotic tissues at 11 dpi while the infection never extended beyond the point of inoculation in Topas-Rlm2 plants. Other studies have reported similar symptom progression at varying time points Haddadi et al., 2016), thus suggesting that experimental conditions can influence the rapidity with which L. maculans can infect its host. For this reason, it is critical to ensure that sampling procedures include enough replications that will capture an accurate biological variability within a condition. In this work we have used five biological replications to enhance both the biological and statistical power to compare gene expression across stages and conditions. We have also collected the whole infected cotyledons for analysis to achieve a thorough FIGURE 11 | Functional categorization of upregulated genes at 7 days post inoculation (dpi) during compatible and incompatible canola-Leptosphaeria maculans interactions. The right y-axis indicates the number of genes in a category. The left y-axis indicates the percentage of a specific category of genes in the main category. Upregulated genes were identified by comparing in planta samples with axenic cultures. Analyses were performed with five biological replicates and a threshold value of FDR < 0.0001 and log 2 fold-change > 1.5.
understanding of pathogen transcriptome activities at cotyledon level during infection. This approach was validated by our PC analysis (see Figure 3) where the clustering of samples confirmed a uniformity within a given treatment and a variability among treatments thereby supporting that observed phenotypes were indeed associated with differential gene expression. This further supported subsequent statistical analyses of DEGs.

Expression of Known Avrs and CSEPs during Biotrophic and Necrotrophic Phases
Several Avr and R genes have been identified and/or proposed to play a role in the L. maculans-canola interaction (Balesdent et al., 2002;Ghanbarnia et al., 2012). Following our analyses, all known Avr genes observed by Lowe et al. (2014) were also found to be highly expressed at 7 dpi, and reduced in their expression at 11 dpi in the compatible interaction. Based on a comparison with Topas-Rlm2, it is apparent that the period 5-7 dpi harbors the biotrophic phase of L. maculans and 11 dpi is more consistent with the necrotrophic phase. As such the differential expression of Avrs at 7 dpi indicates that some of them are clearly involved primarily in the establishment of biotrophy, and possibly the transition to necrotrophy but are no longer relevant when L. maculans has entered its necrotrophic phase.
With the availability of full genome sequences and more advanced computational tools and pipelines, we were able to identify and characterize effectors in L. maculans and better address their potential role/functionality. From strict computational predictions, there are 552 classically secreted proteins in L. maculans, which represents an unrealistic number in terms of functional effectors (Sonah et al., 2016). Recently developed tools based on machine learning were used here to further prioritize 134 CSEPs that showed evidence of upregulation in RNA-seq data. At 3 dpi, which corresponded to an asymptomatic growth phase, a common set of effectors highly expressed in both the compatible and incompatible conditions was observed, which suggests they play a minimal role in the fate of the interaction. By contrast, the expression pattern of effectors varied drastically between the compatible and incompatible conditions at 7 dpi. This approach allowed to narrow down the list of possible functional effectors to 28 that were uniquely upregulated at 7 dpi only under compatible conditions and thus presumed to be important for the establishment and maintenance of the biotrophic phase. As a matter of fact, many of the identified effectors here were either known Avr, or genes for which a role in pathogenicity was suggested (Haddadi et al., 2016). In addition, eight new effectors are proposed on the basis of their features and expression, and should be interesting candidates in future functionality assays. The studies by Haddadi et al. (2016) and Lowe et al. (2014) used different strains of L. maculans, which explains the differences in the number of differentially expressed effectors. For instance, AvrLm6 was not expressed in our study since the L. maculans strain D5 does not have AvrLm6. This suggests that it should be important to properly assess resistant germplasm with specific L. maculans strains present in a given region.
Compared to biotrophic effectors, very few effectors responsible for the necrotrophic phase are known (Lo Presti et al., 2015). In the present study, on the basis of comparative expression in compatible and incompatible interactions, we have identified 15 effectors that were distinct from Avr genes (and associated effectors) in their chronology of expression. This suggests that these effectors are specific to the necrotrophic phase of the fungus and that L. maculans, as a hemibiotroph, has indeed evolved different mechanisms to support its biotrophic and necrotrophic phase. At the same time, high expression of a single crinkler-type effector at 7 dpi in L. maculans suggests its involvement in the transition from a biotrophic to a necrotrophic phase. However, the occurrence of a single crinkler in L. maculans compared to 74 in P. sojae suggests that L. maculans does not have the crinkler-mediated mechanism leading to the necrotrophic phase found in P. sojae and other oomycetes.
A biotrophic-necrotrophic effector system that has been well studied is that of Phytophthora infestans in its interaction with potato (Whisson et al., 2007). During the biotrophic phase, P. infestans secretes AVR3a from its haustoria to suppress celldeath; as the oomycete moves to a necrotrophic stage, AVR3a is downregulated (Whisson et al., 2007). Similarly, our results showed that all 28 genes, including known Avr genes, were highly expressed during the biotrophic phase at 7 dpi and their expression reduced during the necrotrophic phase at 11 dpi. The same genes showed limited or no expression in the incompatible interaction. This suggests that these genes produce functional effectors and are key factors responsible for crosstalk between L. maculans and its host. Another gene found here, Lm5LysM, is highly expressed at 7 dpi and shows homology with the SLP1 gene of M. grisea that is expressed at the interface between the fungal cell wall and host cell plasma membrane during biotrophic invasion. The LysM-type effectors have been previously associated with biotrophy and shown to be involved in plant-fungus interactions (Gust et al., 2012;Kombrink and Thomma, 2013;Lowe et al., 2014).

Expression Pattern of NLPs Confirming the Necrotrophic Stage in L. maculans
In the present study, the highest level of expression for Lm-NLP (gene_11090) gene at 11 dpi, as L. maculans entered into its necrotrophic phase, is an observation consistent with the results of Haddadi et al. (2016). Indeed, these authors used the NLP gene expression profile as a means to distinguish genes related to the biotrophic or necrotrophic phase in L. maculans; the induction of necrosis by Lm-NLP was confirmed with a transient assay in tobacco (Haddadi et al., 2016). Our results do confirm that the highest expression of NLPs is synchronized with the necrotrophic phase, while that of known Avrs with the biotrophic phase. Avrs and NLPs can thus be considered as valid markers to characterize CSEPs as biotrophic or necrotrophic effectors. Similar findings were observed in P. infestans where INF1 and Nep1-like effectors are secreted at later stages of infection that correspond to P. infestans transitioning from a biotrophic to a necrotrophic stage (Kanneganti et al., 2006).

Expression Dynamics of CAZymes and TFs
The increasing number of upregulated CAZymes from 5 to 11 dpi in the compatible interaction can be related to the biphasic life style of L. maculans as previously suggested by Lowe et al. (2014). This is particularly relevant when compared to the incompatible interaction, where this number remained relatively unchanged over the course of infection. Our results showed that GH, CE, and AA were the most prevalent groups of CAZymes among upregulated genes. During the early stages (3 dpi) of L. maculans infection, mostly cellulose and pectin-degrading enzymes were the most prevalent groups of CAZymes among upregulated genes. On the other hand, GH, CE and AA families were predominantly upregulated during the necrotrophic stage. Pathogenic fungi face the plant cell wall as a first barrier to establish infection. Plant cell walls are mainly composed of carbohydrates and glycoproteins. To breakdown this barrier, plant pathogenic fungi need to secrete a diverse range of carbohydrate-active enzymes (CAZymes). However, our results clearly showed that the release of CAZymes in the early stages is not a key determinant of the interaction since their number and expression were similar in both interactions.
CAZymes are also involved in nutrient uptake and those prominently expressed at later stages of infection are thought to be involved in uptake of amino acid and sugars from the host. Our results bring stronger support to the concept suggested by Lowe et al. (2014) that CAZymes play an important role in L. maculans-canola interaction, notably during the establishment of necrotrophy.
Biphasic expression turnover similar to that of CAZymes was also observed with TFs. As the infection progressed in the compatible interaction, a much higher number of differentially expressed TFs was noted especially from 5 to 7 dpi. This period is critical in the fate of the interaction since it is clearly synchronized with the establishment of the infection in the compatible interaction. Of particular importance, the LmStuA TF, a member of APSES domain-containing proteins, showed its highest expression at 7 dpi in L. maculans interaction with Topas-wild. The expression profile of APSES domaincontaining genes observed in this study is well aligned with an earlier report by Soyer et al. (2015). The StuA TF was found to be involved in morphogenesis, metabolites production and effector regulation (Baeza-Montañez et al., 2015;Soyer et al., 2015). Recently, StuA was suggested to play a key role in the L. maculans-canola interaction since its silencing led to a reduced expression of AvrLm1, AvrLm6, and AvrLm4-7, at 7dpi (Soyer et al., 2015). Moreover, the Zn2Cys6 TFs were upregulated in the compatible interaction supporting its role in regulating pathogenicity-related genes during the biotrophic phase. The Zn2Cys6 TFs have been reported to be involved in different regulatory functions. They are also unique to fungi and have been extensively studied in Saccharomyces cerevisiae and Aspergillus nidulans (Shimizu et al., 2003;Vienken et al., 2005).

Secondary Metabolites Involved in the Pathogenesis of L. maculans
Secondary metabolites like NRPs and phytotoxins such as sirodesmin were found to have a differential expression pattern between the compatible and incompatible interactions. Differences in gene expression pattern between the compatible and incompatible interaction were observed for NRPs, known to be involved in the production of phytotoxins, siderophores, and pigments, mostly at 11 dpi. Another important phytotoxin involved in the establishment of necrosis is sirodesmin, which is regulated by Sir genes (Gardiner et al., 2004;Haddadi et al., 2016). Most Sir genes were found to be highly expressed at 11 dpi compared to in planta L. maculans early growth stages in this study. A notable exception was SirO and SirT, which were highly expressed at 3 dpi and are thought to be involved in the production of sirodesmin early in the infection process. Indeed, Gardiner et al. (2004) have also reported significant amount of sirodesmin production at 4 dpi with L. maculans. On the other hand, the fact that SirT was equally expressed in the compatible and incompatible interactions at 3 dpi in our work would suggest that additional factors must complement their activity for infection to occur.
Based on our WEGO functional annotation, the respective vigorous and restricted growth of L. maculans on the compatible and incompatible hosts correlated very well with enrichment of genes in functional categories. Toward the later stages of the experiment, very few classes of functional categories were represented in the incompatible interaction, an observation that clearly confirms the inability of L. maculans to infect Topas-Rlm2. However, it is not possible to conduct robust statistical analyses because of the disparity in the number of genes between the compatible (572) and incompatible interactions (10). Nevertheless, these analyses give a good qualitative visualization of the genes and functions, and their relative importance in the case of compatibility. Interestingly, these analyses also highlighted in the compatible interaction how specific functional classes were indicative of the biotrophic or necrotrophic stage, and the transition from one to the other.
Validation of DEGs identified with transcriptome profiling is always a concern to build a confidence about the significance of results. The validation of DEGs by qPCR is expected for the microarray studies or transcriptomic studies with no or limited biological replications (Deshmukh et al., 2010;Xie et al., 2016). For instance, microarray studies mostly relied on qPCR to validate the results since there is a limit for the number of probes hybridized to the target in microarrays. Therefore, the highest level of gene expression is never truly represented for highly expressed genes. This also leads to a bias of data normalization and reduced linear range. By contrast, RNA-seq provides digital expression for the entire transcripts without any minimum or maximum count limit. In our study, we have used five replications, in line with the robust standards accepted for RNA-seq experiments (Fang and Cui, 2011). This is further supported by the fact that all previously identified Avr genes (AvrLm1, AvrLm4-7, AvrLm11, AvrLmJ1, AvrLm2) had an expression pattern consistent with expectations and previous reports Haddadi et al., 2016).

CONCLUSIONS
Hemibiotrophic fungi such as L. maculans have a distinctive life style that involves complex molecular processes and the expression turnover of thousands of genes. The atlas of gene expression provided here will be helpful to understand the molecular crosstalk between canola and L. maculans as it relates to compatibility or incompatibility, and could be exploited toward the deployment of novel strategies to overcome blackleg disease. The comparison made between the compatible and incompatible interaction highlighted the role of specific CSEPs, CAZymes, TFs, and secondary metabolites involved in the infection process. The differential expression pattern observed for these classes of pathogenicity-related genes can serve as a valuable resource to differentiate compatible and incompatible interactions in the context of developing resistant canola germplasm. In addition, our time series analysis based on a comprehensive and comparative differential gene expression as L. maculans infects a susceptible or resistant host has brought to light elements that define the biotrophic and necrotrophic phases of the fungus, as well as some of the mechanisms involved in the transition between the two phases.

AUTHOR CONTRIBUTIONS
HS performed RNA extractions, library preparations, data analysis, and compilation; XZ and WF performed all bio-assays, maintained the fungal cultures, and performed RNA extractions, HB developed introgression line and provided seed material, RD and RB contributed to data analysis and overall conclusions, HS and RB wrote initial draft of the MS with final contributions from XZ, RD, HB, and WF. RB designed and directed the project.

FUNDING
This study is funded by the Agri-Innovation program Growing Forward 2, SaskCanola and Agriculture and Agri-Food Canada, and the Canada Research Chair in plant protection to RB.