Biological and RNA samples representing all aspects of the life cycle of the fungus in interaction with its host
The RNA-Seq experiment was designed to generate samples corresponding to all stages of the interaction of L. maculans with its host plant, either alive (Fig. 1e–c) or dead (stem residues after harvest; Fig. 1d), in controlled conditions or in field experiments under natural inoculum pressure, over periods of time ranging from a few days to months or years (Fig. 1; For a detailed description of biological samples, see the “Methods” section). As a reference, we included a series of samples grown in axenic conditions promoting a particular fungal developmental stage (mycelium growth, pycnidium differentiation, pseudothecium differentiation, resting conidia, conidial germination) (Fig. 1e), to distinguish between these fundamental biological processes and those specific to interaction with the plant. In total, 102 biological samples, corresponding to 37 different sets of conditions, were subjected to RNA-Seq. Several quality control checks and parameter optimization were required before this heterogeneous dataset could be subjected to statistical analyses. In fine, 32 sets of conditions were retained for statistical analyses, whereas five sets of conditions were excluded because too few fungal reads were available (Additional files 1, 2 and 3: S1 Text, S1 and S2 Tables).
Quality control analyses highlight the consistency of the RNA-Seq dataset
Clustering based on pairwise Pearson’s correlation analyses of the 32 sets of conditions identified four groups of clustered conditions (Fig. 2). In each of these four groups, principal component analysis (PCA) showed a low level of variation between biological replicates, except for a few field samples and in vitro crosses (Additional files 4 and 5: S2 Text, S1A-D Fig). The four clusters (Fig. 2) corresponded to the following: (i) the early infection stage on cotyledons and leaves and colonization of the petiole (group 2), (ii) the late infection stage on stems (group 1), (iii) the saprophytic lifestyle (group 3), in controlled conditions or in the field, and (iv) all the in vitro growth conditions promoting the differentiation of pseudothecia and/or pycnidia (group 4).
Within each group, lifestyle provided a second level of sample discrimination, with samples from biotrophic or endophytic stages distinguished from necrotrophic-stage samples (Fig. 2). For example, within group 2, samples displaying necrotrophic behavior on petioles 14 days post infection (DPI) were more closely related to the necrotrophic samples collected from cotyledons 12 and 15 DPI than to the corresponding biotrophic samples (petioles 7 DPI) (Fig. 2). Similarly, in group 1, the residues collected immediately after harvest grouped with the last samples collected from the stem base, even though the residues did not originate from the same growing season, or from the same rapeseed variety, suggesting that the transition from stem necrosis to saprophytic life takes place after harvest (Fig. 2).
Some axenic growth conditions were also grouped with some in planta conditions. Within group 2, germinating conidia grouped with the earliest samples collected from cotyledons, consistent with the use of conidia to inoculate plant tissues in controlled conditions (Fig. 2). Samples grown in axenic culture on solid V8 medium under conditions promoting mycelial growth over sporulation clustered within group 3, indicating that mycelial growth is a major component of the saprophytic lifestyle in the wild. By contrast, mycelium grown in static liquid medium grouped with all conditions in which the fungus colonized stem tissues, possibly reflecting growth under low-oxygen conditions in plant vessels or intercellular spaces.
Analytical strategy for identifying genes involved in pathogen-plant interactions
For the identification of fungal genes specifically involved in at least one stage of the pathogen-plant interaction, we used the 10 sets of in vitro growth conditions as 10 sets of control conditions, to exclude all genes involved in basic fungal metabolism or life traits other than those associated with plant infection (Fig. 3). We choose very stringent criteria (LogFC > 4; p value < 0.01) to focus on genes specifically and robustly overexpressed in planta relative to in vitro control conditions. We analyzed the expression patterns of these genes, their functional annotation, their genomic localization, and their location with respect to histone modifications during axenic growth. We then focused on genes encoding SSPs from the new repertoire generated in this study (Additional file 6: S3 Text). Finally, the expression profiles of co-expressed SSP genes were used as a template for the identification of waves of genes displaying a high level of coregulation (R-squared > 0.80) and exemplifying specific pathogenic strategies (Fig. 3).
Genes overexpressed during interaction with the plant define eight biologically relevant expression clusters
Relative to axenic growth conditions, 1207 genes were found to be overexpressed in at least one set of conditions in planta; thus, 9.2% of L. maculans genes are specifically mobilized during the interaction of the fungus with its host, alive or dead. Eight expression profiles were identified, referred to hereafter as clusters 1 to 8, defining a complex landscape of expression profiles (Fig. 4). We analyzed the enrichment of each cluster in GO terms, focusing on CAZymes (carbohydrate-active enzymes) associated with degradation activities, to identify the functions and pathways associated with each expression profile (Additional files 7, 8, 9 and 10: S2-S5 Fig).
Clusters 1, 2, and 3 included genes expressed sequentially during cotyledon infection and were representative of three consecutive stages occurring during the 2-week period of cotyledon colonization by L. maculans. Cluster 1 encompassed 97 genes highly expressed at the first stages of cotyledon infection, when the conidia germinate and the hyphae penetrate the wounded plant tissues (2–5 DPI). It was enriched in GO terms corresponding to carbohydrate metabolic processes (p = 1 × 10− 2) and catalytic activities (p = 1 × 10− 3), including peptidase and hydrolase activities (all p ≤ 1 × 10− 2). An enrichment in cutinase activity (p = 1 × 10− 2), specific to this cluster, was also observed (Additional files 7, 8 and 10: S2-S3, S5 Figs).
The 148 genes of cluster 2 displayed peaks of expression during all stages of asymptomatic growth within plant tissues (cotyledon colonization 5–9 DPI, first stage of petiole colonization, first two dates of stem colonization in the field). The genes of this cluster were annotated with few GO terms or functions, but were strongly enriched in SSP genes (p = 2 × 10− 16; Additional file 8: S3 Fig). Cluster 3 encompassed 87 genes highly expressed during the necrotrophic stage on both cotyledons (9–15 DPI) and petioles, but not during the necrotrophic stem canker stage in the field. This cluster was enriched in carboxylic ester hydrolase activities (p = 3 × 10− 3) suggesting a role in degradation activities specifically targeting the plant cell wall.
Cluster 4 encompassed 238 genes highly expressed during the shift from biotrophy to necrotrophy in cotyledons (7–9 DPI) and in stem bases infected in the field (7–8 months post sowing, “MPS”). This cluster displayed the strongest enrichment in catalytic processes (p = 2 × 10− 10) and hydrolase activities (p = 1 × 10− 12), and a specific enrichment in proteolysis processes (p = 2 × 10− 2; Additional file 8: S3 Fig). It also displayed the largest number of CAZymes (63 genes; Additional file 10: S5 Fig) targeting different substrates, mostly plant cell-wall components, such as pectin, arabinose, and galactan (Additional file 10: S5 Fig).
Two clusters (cluster 5 and 6) displayed an increase in expression late in stem infection, but differed in terms of the timing of expression and the functions involved. Cluster 5 contained 128 genes highly expressed during the asymptomatic colonization of stems in the field (at 7–8 MPS) and on residues immediately after harvest. It was enriched in catalytic activities (p = 1 × 10− 3), mostly linked to saccharide degradation, suggesting involvement in activities relating to the uptake of sugar resources (carbon-oxygen lyase activity acting on polysaccharides; hydrolase activity; cellulose binding; Additional file 8: S3 Fig). Cluster 6 contained 141 genes highly expressed during colonization of the stem base in the field, with maximal expression at 10 MPS, corresponding to the time at which stem canker develops. This cluster was enriched in catalytic activities (p = 2 × 10− 4) but was also linked to cofactor binding (p = 3 × 10− 6), monooxygenase (p = 1 × 10− 2), and oxidoreduction activities (p = 5 × 10− 5), involved in detoxification processes (Additional files 8 and 9: S3 and S4 Figs). The higher proportion of xylanase genes in cluster 6 than in cluster 5 (Additional file 10: S5 Fig) suggests a change in the mode of nutrition during the necrotrophic stage on stems.
Finally, clusters 7 and 8 were specific to saprophytic behavior on stem residues. The 198 genes in cluster 7 displayed an increase in expression from the last sampling on necrotic stems until the last sampling on stem residues 1 year after harvest. This cluster was enriched in genes relating to monooxygenase activities (p = 1 × 10− 4) and contained the largest number of genes encoding LPMOs (lytic polysaccharide monooxygenases), cellulases, and xylanases. Cluster 8 encompassed 170 genes specifically expressed during all stages occurring on senescent or dead plant tissue, with the highest levels of expression during saprophytic growth on residues. It was also enriched in oxidoreductase processes (p = 1 × 10− 3), in response to oxidative stresses and oxidant detoxification associated with an enrichment in catalase (p = 1 × 10− 3) and peroxidase (p = 3 × 10− 3).
The different clusters of gene expression thus can be designated as follows: “Penetration and establishment” (cluster 1), “Biotrophy” (cluster 2), “Cotyledon and petiole necrotrophy” (cluster 3), “Biotrophy to necrotrophy transition” (cluster 4), “Stem biotrophy” (cluster 5), “Stem necrotrophy” (cluster 6), “Stem canker and saprophytism” (cluster 7), and “Saprophytism” (cluster 8).
Contribution of candidate pathogenicity genes to fungal life in interaction with the plant
Genes involved in fungal pathogenicity in other fungal species
We wondered whether the 9.2% of infection-specific genes were conserved in other phytopathogenic fungi and if we could relate their phase-specific expression pattern with the lifestyle of these fungi. Of the corresponding 1207 proteins, 235 (19.5%) were specific to L. maculans. Of the remaining 979 proteins sharing homologies with proteins predicted in other fungal species (Additional file 11: S6 Fig), 316 (32.2%) had matches with hypothetical or predicted proteins, while the remaining ones matched with proteins with an annotated function or domain. The cluster 2 comprised the lowest proportion of proteins with orthologs in other species (64.84%) compared to other clusters, likely due to its high content of species-specific effectors (Additional file 11: S6 Fig). Globally, most hits were with related phytopathogenic Dothideomycete species (genera Bipolaris, Alternaria, Pyrenophora, Parastagonospora, etc.), for which only limited RNA-Seq resources are available, preventing us from comparing the clusters of expression of L. maculans with theirs. In contrast, a non-negligible number of hits were found with species of the Colletotrichum or Fusarium genera, and to a lesser extent with Z. tritici and B. cinerea, species with either an hemibiotrophic or a necrotrophic pathogenicity behavior, for which extensive RNA-Seq data are available. For these species, we identified a very low number of homologous annotated proteins present in the eight expression profiles (Additional file 12: S7 Fig). The cluster 4 contained the highest number of homologous proteins annotated in Fusarium and Colletotrichum species. They corresponded to catalytic activities such as protease and polysaccharide degradation enzymes. These results suggest that the degrading enzymes involved in the transition from biotrophy to necrotrophy are conserved among the hemibiotrophic fungi.
Several genes previously shown to be involved in fungal pathogenicity were found in one of the eight expression clusters and were mostly associated with biotrophic stages of the interaction. One RALF (rapid alkalinization factor) gene known to favor fungal infection [39] was identified in cluster 1 (Lmb_jn3_05329). Two of the six salicylate hydroxylase genes, known to compromise plant salicylic acid signaling in other models [40], were found in clusters 1 (Lmb_jn3_12582) and 2 (Lmb_jn3_11892). Three of the four LysM-domain containing genes of the genome, known to protect chitin from degradation by plant enzymes or to scavenge chitin oligomers [41], were found in clusters 2 (two genes, Lmb_jn3_00461, Lmb_jn3_10047) and 6 (one gene, Lmb_jn3_04300), and LysM genes from cluster 2 matched with Colletotrichum spp. LysM genes expressed in planta. Three isochorismatase genes (Lmb_jn3_11641, Lmb_jn3_10045, Lmb_jn3_06552), interfering with plant salicylate and jasmonate defense signaling in other phytopathogens, were present in clusters 1, 4, and 6, respectively. Last, all the proteins involved in the biosynthesis of ABA in B. cinerea matched L. maculans proteins annotated in the ABA biosynthesis cluster, but the ABA genes were not expressed in planta in B. cinerea while they are expressed in the biotrophy cluster 2 in L. maculans (Additional file 12: S7 Fig). Finally, one of the two genes encoding NEP (necrosis and ethylene-inducing protein), frequently associated with necrotrophic behavior [42], was found to be overexpressed in planta, and belonged to cluster 4 (Lmb_jn3_10035).
Secondary metabolite gene clusters
The genome of L. maculans contains 22 genes predicted by SMURF tools to encode secondary metabolite enzymes (SM) and annotated as “polyketide synthase” (PKS, 11) “non-ribosomal peptide synthetase” (NRPS, 5), or “others” (dimethylallyltryptophan synthase, 1; NRPS-like, 4; PKS-like, 1). Seven of the 1207 genes in the eight expression clusters encoded SM enzymes from the PKS or NRPS family. Four of these genes were highly expressed during necrotrophic stages of the life cycle within the stem (cluster 6), whereas the other three were associated with clusters 2, 4, and 5 (Additional file 13: S3 Table).
Among the SM genes analyzed in previous studies [27, 43, 44], the PKS gene involved in abscisic acid (ABA) synthesis was linked to biotrophic behavior (cluster 2), whereas genes involved in the synthesis of the toxins sirodesmin PL and phomenoic acid were associated with stem necrotrophy, together with one other PKS and one NRPS-like gene (cluster 6) (Additional file 13: S3 Table). Secondary metabolite enzymes encode an SM backbone that may be modified by tailoring enzymes, forming complex secondary metabolite gene clusters (SMGC). Here, only three of the seven SM enzyme genes overexpressed in planta showed coregulation with all or part of surrounding genes of their SMGC: (i) the PKS responsible for ABA production (all genes of the ABA SMGC being included in cluster 2), (ii) the NRPS responsible for sirodesmin PL production (16 of the 25 genes of the SMGC included in cluster 6), and (iii) an unknown NRPS-like gene, co-expressed with nine other surrounding genes also included in cluster 6 (Additional file 14: S8 Fig).
Genes involved in L. maculans pathogenicity
Thirty-six of the 52 genes previously identified as candidate pathogenicity genes in L. maculans (Additional file 15: S4 Table) did not belong to the SM or avirulence effector gene (see below) categories. Fifteen of these genes were identified within a single cluster (Additional file 15: S4 Table): cluster 1 (2 genes), cluster 2 (5 genes), cluster 3 (1 gene), cluster 4 (4 genes), cluster 5 (2 genes), or cluster 7 (1 gene). Nine had previously been selected as putative pathogenicity genes on the basis of their overexpression during infection, for inactivation with CRISP-cas9, but the inactivation of these genes individually resulted in no visible pathogenicity defect at the cotyledon stage [45] (Additional file 15: S4 Table). Only one gene, encoding a 3-ketoacyl-thiolase (cluster 1), was shown to be involved in pathogenicity, but the mutant also displayed impaired in vitro growth and germination [28] (Additional file 15: S4 Table). The other five genes have not been validated by functional approaches and are, thus, still considered to be candidate pathogenicity-related genes. Four of these genes (three GH genes and one ABC transporter) were included in the late expression clusters 4, 5, and 7, suggesting that they may be bona fide pathogenicity genes (Additional file 15: S4 Table).
By contrast, another 21 putative pathogenicity genes did not belong to any of the expression clusters. Nineteen had previously been analyzed in functional studies (knockout or silenced mutants), and 10 were shown to induce a loss or reduction of pathogenicity following inactivation (Additional file 15: S4 Table). However, the inactivation of six of these genes also induced defects in growth, sporulation, or germination [30, 33, 36, 45, 46], suggesting a general effect on different stages of the fungal life cycle with no specific overexpression during interactions with the plant.
Validated or candidate SSP effectors of L. maculans
SSPs have been identified in L. maculans and were classified as “early” effectors (including all known avirulence proteins, inducing an avirulence phenotype when matching the corresponding resistance genes) or “late” candidate effectors by Gervais et al. [14]. The role of avirulence proteins as effectors suppressing plant defense responses or increasing/decreasing the size of cotyledon lesions has been investigated in a few cases following comparisons of near-isogenic L. maculans isolates, gene silencing or complementation (Additional file 15: S4 Table, [21, 22, 47, 48]). The most recently published version of the genome assembly of L. maculans isolate JN3 ([49]; Additional file 16: S4 Text) included a number of major changes. We therefore produced a new repertoire of genes encoding putative SSPs, encompassing 1070 genes (8.2% of the entire L. maculans gene set) (Additional files 6, 17, 18 and 19: S3 Text, S9 Fig, S5 and S6 Tables).
All expression clusters displayed a significant enrichment in genes encoding SSPs (Additional files 8: S3 Fig). We found that 272 of the 1207 genes upregulated during the L. maculans interaction with B. napus (22.5%) encoded SSPs. This enrichment was limited in clusters 1, 3, 6, 7, and 8, which contained only 14–17% of SSP-encoding genes (0.013 < p < 3.10− 5), but three other clusters displayed a much higher level of enrichment in genes encoding SSPs: cluster 4 (23.1%, p = 7 × 10− 16), cluster 5 (31.2%, p = 2 × 10− 16), and cluster 2 (45.3%, p = 2 × 10− 16). Interestingly, both clusters 2 and 5 are associated with biotrophic behavior of the fungus. All the currently known AvrLm genes belonged to the “Biotrophy” cluster (cluster 2) (Additional file 15: S4 Table). This cluster displayed sequential rounds of overexpression/repression, with three peaks of overexpression occurring during the three biotrophic stages (5–9 DPI on cotyledons, leaf infection in the field; 7 DPI on petioles; 6–7 MPS in field plant stems) (Fig. 5; Additional file 20: S10 Fig).
The ten “late” effectors, expressed during stem colonization in controlled conditions [14], were scattered over four different clusters (3, 4, 5, 6), along with a series of other genes encoding SSPs (cluster 3, 14 SSPs; cluster 4, 55 SSPs; cluster 5, 40 SSPs and cluster 6, 22 SSPs) (Fig. 5). The “late” effectors in clusters 5 and 6 were specific to stem colonization (and the first stages of saprotrophy on residues), but those in clusters 3 and 4 were also expressed during cotyledon infection. These “late” effector candidates were associated with either biotrophic behavior on stems (cluster 5) or necrotrophic behavior on cotyledons and stems (clusters 3 and 6).
Clusters 1, 7, and 8 were also found to be enriched in genes encoding SSPs (p value < 5 × 10− 3). They contained 16, 33, and 24 SSP genes, respectively. Candidate effectors may, therefore, be mobilized very early in colonization (cluster 1) and during saprophytic life (clusters 7 and 8) (Additional file 20: S10 Fig).
Additional RNA-Seq samples (ascospores germinating on unwounded cotyledons 24 and 48 h post inoculation; stem base samples 2, 3 and 5 MPS) were available. They provided too few L. maculans reads for statistical analyses, but these samples allowed us to confirm that genes belonging to cluster 2, and “early” effector genes in particular, were among the genes most strongly expressed both during ascospore germination on cotyledons and very early stages of stem colonization in winter (Additional files 21 and 22: S5 Text, S7 Table).
Histone methylation associated with gene regulation
All cloned avirulence genes of L. maculans identified to date are associated with H3K9me3. Furthermore, H3K9me3 and H3K27me3 domains are significantly enriched in genes encoding effectors, regardless of their expression [38]. During axenic culture, 7373 L. maculans genes are associated with H3K4me2, 104 with H3K9me3, 2020 with H3K27me3 and 101 with the two heterochromatin histone modifications [38]. We combined the genome-wide histone maps previously generated in vitro with our analysis of gene expression during infection, to determine whether (i) all genes, (ii) only (putative) effector genes, or (iii) only certain subsets of genes upregulated in planta were associated with a particular chromatin landscape. H3K4me2-domains encompass 56% of L. maculans genes, but only 20% of the genes upregulated in planta were associated with H3K4me2 (χ2 test = 662; p < 2 × 10− 16). By contrast, the genes upregulated in planta were significantly enriched in genes located in an H3K27me3 domain, as 44% of the upregulated genes were located in such a domain (χ2 test = 749, p < 2 × 10− 16). Strikingly, only 205 genes were located in H3K9me3 or H3K27me3 domains in vitro, but 85 of these genes were upregulated in planta. Thus, regardless of the wave of expression concerned, the genes upregulated in planta were enriched in genes associated with heterochromatin in vitro. Only cluster 2 was significantly enriched in genes associated with the H3K9me3 histone modification (Fig. 6). Interestingly, all clusters of genes upregulated in planta relative to axenic culture were enriched in H3K27me3-associated genes in vitro, whereas these clusters contained significantly fewer genes located in H3K4me2-domains in vitro than the rest of the genome (Fig. 6). The integration of RNA-Seq and ChIP-Seq data revealed that the genes upregulated in planta were not randomly located in the genome of L. maculans, with H3K4me2-domains significantly depleted of genes upregulated during rapeseed infection, and both types of heterochromatin domains significantly enriched in such genes.
The co-expression of effector genes defines waves of coregulated gene expression
Almost all the expression clusters were enriched in genes encoding SSPs. However within each cluster, the genes encoding SSPs could have different expression profiles, due to the clustering method (Additional file 20: S10 Fig). We aimed here to increase the robustness of the detection of SSP expression waves. We defined eight SSP reference expression profiles by calculating the mean levels of expression for SSP genes within the eight clusters. We selected genes displaying significant coregulation (SSP and non-SSP genes), by analyzing the linear regression between the expression patterns of the entire set of 13,047 genes and that of each SSP reference expression profile. Genes significantly co-expressed with the reference expression wave tested (R2 > 0.80, p < 0.05) were included in the new SSP expression waves and were considered to be “highly coregulated genes.”
For clusters 1, 3, and 8, no highly coregulated genes could be identified (Fig. 7). For the other five clusters, waves of highly coregulated genes were identified, each containing five to 39 highly coregulated SSP genes (including 90 SSP genes from the 178 previously assigned to these clusters), together with 186 highly coregulated genes encoding other types of proteins (Fig. 7). Two waves of expression, corresponding to clusters 5 and 6, displayed coregulation of a limited number of SSP genes (seven and five, respectively, i.e., 17% of the SSP genes of these clusters), together with 20 and 25 genes encoding other proteins, respectively. However, only 15 (cluster 5) and 13 (cluster 6) of these genes belonged to the clusters (Fig. 7), suggesting that additional coregulated genes are either expressed in axenic conditions, or do not display specific regulation during plant infection. For clusters 2, 4, and 7, a higher proportion of highly coregulated SSP genes was found, with 58% (exact confidence interval [45–70%]), 34% ([22–48%]), and 63% ([45–79%]) of the SSP genes of the clusters retained in the waves of highly coregulated genes, respectively (Fig. 7). Thus, clusters 2 and 7 displayed the highest proportion of highly coregulated SSP genes, suggesting that the regulation of expression is much tighter in these two clusters than in the others. The highly coregulated genes in the “Biotrophy” wave (contained in cluster 2) included 39 SSPs, 23 (59%) of which displayed no sequence matches to the NCBI nr protein database (Additional file 23: S8 Table) and could be considered L. maculans-specific effectors. These genes also included other genes mentioned above, such as the genes of the ABA cluster, providing additional support for their role in the biotrophic parts of the interaction. The highly coregulated genes in the “Stem canker and necrotrophy” wave included most of the SSP genes of the cluster, but only two of these genes (9.5%) were specific to L. maculans. This wave also included a large number of genes (24) encoding CAZymes.