# Nitrogen Supply and Host-Plant Genotype Modulate the Transcriptomic Profile of Plasmodiophora brassicae

Kévin Gazengel, Yoann Aigu, Christine Lariagon, Mathilde Humeau, Antoine

Gravot, Maria M. Manzanares-Dauleux, Stéphanie Daval

► To cite this version:

Kévin Gazengel, Yoann Aigu, Christine Lariagon, Mathilde Humeau, Antoine Gravot, et al.. Nitrogen Supply and Host-Plant Genotype Modulate the Transcriptomic Profile of Plasmodiophora brassicae. Frontiers in Microbiology, 2021, 12, pp.701067. 10.3389/fmicb.2021.701067. hal-03336692

HAL Id: hal-03336692

https://hal.science/hal-03336692

Submitted on 7 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Nitrogen Supply and Host-Plant Genotype Modulate the Transcriptomic Profile of *Plasmodiophora brassicae*

Kévin Gazengel, Yoann Aigu, Christine Lariagon, Mathilde Humeau, Antoine Gravot, Maria J. Manzanares-Dauleux and Stéphanie Daval*

IGEPP, INRAE, Institut Agro, Université Rennes 1, Le Rheu, France

Nitrogen fertilization can affect the susceptibility of *Brassica napus* to the telluric pathogen *Plasmodiophora brassicae*. Our previous works highlighted that the influence of nitrogen can strongly vary regarding plant cultivar/pathogen strain combinations, but the underlying mechanisms are unknown. The present work aims to explore how nitrogen supply can affect the molecular physiology of *P. brassicae* through its life epidemiological cycle. A time-course transcriptome experiment was conducted to study the interaction, under two conditions of nitrogen supply, between isolate eH and two *B. napus* genotypes (Yudal and HD-018), harboring (or not harboring) low nitrogen-conditional resistance toward this isolate (respectively). *P. brassicae* transcriptional patterns were modulated by nitrogen supply, these modulations being dependent on both host-plant genotype and kinetic time. Functional analysis allowed the identification of *P. brassicae* genes expressed during the secondary phase of infection, which may play a role in the reduction of Yudal disease symptoms in low-nitrogen conditions. Candidate genes included pathogenicity-related genes ("NUDIX," "carboxypeptidase," and "NEP-proteins") and genes associated to obligate biotrophic functions of *P. brassicae*. This work illustrates the importance of considering pathogen's physiological responses to get a better understanding of the influence of abiotic factors on clubroot resistance/susceptibility.

**Keywords:** clubroot, *Brassica napus*, host-pathogen interaction, abiotic stress, RNA-seq, root, pathogenicity-related genes

## INTRODUCTION

Plants are exposed to a wide range of phytopathogens, such as fungi, bacteria, viruses, nematodes, or protists, which cause infectious diseases. The outcome of the plant/pathogen interaction can be influenced by multiple environmental factors, both biotic (e.g., plant-associated microbiota) and abiotic (e.g., temperature, light, water content, salinity, pH, or nitrogen). Nitrogen (N) is a central element in various metabolic pathways: synthesis of enzymes, coenzymes, amino acids, and proteins; photosynthesis; and respiration (Scheible et al., 2004; Xin et al., 2014). Then, the plant's N status is a critical determinant not only for growth, development, and productivity, but also for plant response (susceptibility and resistance/tolerance) to both abiotic constraints [such as water

**OPEN ACCESS**

**Edited by:**
Hossein Borhan,
Agriculture and Agri-Food Canada
(AAFC), Canada

**Reviewed by:**
Edel Pérez-López,
Laval University, Canada
Tyler Helmann,
United States Department
of Agriculture (USDA), United States

**\*Correspondence:**
Stéphanie Daval
stephanie.daval@inrae.fr
orcid.org/0000-0001-8848-4190

**Specialty section:**
This article was submitted to
Microbe and Virus Interactions with
Plants,
a section of the journal
Frontiers in Microbiology

**Received:** 27 April 2021
**Accepted:** 09 June 2021
**Published:** 08 July 2021

**Citation:**
Gazengel K, Aigu Y, Lariagon C, Humeau M, Gravot A, Manzanares-Dauleux MJ and Daval S (2021) Nitrogen Supply and Host-Plant Genotype Modulate the Transcriptomic Profile of Plasmodiophora brassicae. *Front. Microbiol.* 12:701067.
doi: 10.3389/fmicb.2021.701067
(Zhong et al., 2017)] and plant pathogens (Mur et al., 2016). In general, it has been considered that the effect of N on plant diseases was dependent on the lifestyle of the pathogens (Solomon et al., 2003): the development of biotrophic pathogens would be enhanced by N, while the opposite effect would be observed for necrotrophic pathogens. In fact, the effect of N on plant disease, described in many examples, can be opposite depending on the pathogen and on the plant host. Thus, an increase in severity of infection by N fertilization was reported, for example, in powdery mildew (Zhang et al., 2019), stripe rust (Devadas et al., 2014), wheat blast (Ballini et al., 2013), and rice blast (Huang et al., 2017). On the contrary, disease decreases due to N nutrition were described in take-all (Brennan, 1992) and in leaf spot diseases (Krupinsky et al., 2007). Moreover, for the same pathogen, the situation may not be as clear-cut. For example, high N fertilization rates led to higher plant susceptibility to *Botrytis cinerea* in grape (Rhouma et al., 1998), legumes (Davidson and Krysinska-Kaczmarek, 2007), sweet basil (Yermiyahu et al., 2006), and begonia (Pitchay et al., 2007), but to lesser symptoms in tomato in which infection also varied according to the pathogen isolates (Lecompte et al., 2010). In a review that summarized data from 114 studies, an increase, decrease, or no effect of N on plant disease was reported in 62, 42, and 10 cases, respectively (Sun et al., 2020). In addition, the different forms of N supply (ammonium $NH_4^+$ or nitrate $NO_3^-$) can have various effects on plant disease severity because of the use of different assimilation and metabolism pathways (Bolton and Thomma, 2008; Mur et al., 2016).

Plant N nutrition influences both constitutive and induced defenses, as an element in key signal molecules such as NO and polyamines (Takahashi, 2016, 2020). Different N conditions can modulate plant resistance based on pathogen-associated molecular patterns triggered immunity (PTI) and effector-triggered immunity (ETI) (Mur et al., 2016). N-biochemical-mediated defenses are achieved through plant metabolites (e.g., phytoalexins, antimicrobial proteins, amino acids, and organic acids) and defense-related enzymes (Ngadze et al., 2012; Wang et al., 2018). Plant N availability can also impact cellular structure and composition affecting the thickness of the plant's physical barrier (Zhang et al., 2017). N input modulates plant defense but also impacts disease susceptibility by increasing the availability of N compounds for exploitation by pathogens (Tavernier et al., 2007), shaping the plant architecture and providing an environment favorable for pathogen growth and spore production (Mur et al., 2016). Although several studies reported the importance of N supply on pathogenesis and effector delivery (Soulie et al., 2020; Sun et al., 2020), conflicting results have been observed. Thus, in *Magnarpothe oryzae*, the induction or repression of some effectors was found under N starvation conditions (Donofrio et al., 2006; Huang et al., 2017).

The abundant literature highlights the various and contradictory effects of N fertilization on the outcome of the plant/pathogen interaction and underlines the need to better understand the mechanisms underlying this process in different pathosystems.

In the case of clubroot, a soilborne disease caused by the obligate parasite *Plasmodiophora brassicae* and affecting the

Brassicaceae, the study of N effect is of prime importance, notably in oilseed rape (*Brassica napus*) since this crop requires the use of relatively high amounts of N fertilizer (Rathke et al., 2005). Resistance to plant pathogens and adaptation to low-input agricultural practices (especially low N input) are thus two major traits targeted in oilseed rape breeding efforts (Bouchet et al., 2014, 2016). In this pathosystem, several studies have suggested that a high N supply tends to reduce the damage caused by *P. brassicae* infection (Dixon, 2009; Gossen et al., 2014). On the other hand, a study on a set of 92 diverse oilseed rape accessions and 108 lines derived from a cross between "Darmor-bzh" (clubroot-resistant) and "Yudal" (clubroot-susceptible) showed that some genotypes classified as susceptible to clubroot became resistant under N constraint (Laperche et al., 2017). In this previous work, we have shown that N supply can modulate the effect of two QTLs controlling symptom development (Laperche et al., 2017) and spore production (Aigu et al., 2018). Surprisingly, no studies have been published about the effect of N directly on the clubroot pathogen *P. brassicae*, while knowing that N input could have an influence on many factors involved in the different steps of its epidemiological cycle. First, the resting spores, survival form in the soil, germinate when root exudates are sensed, leading to biflagellate primary zoospores that infect the root hairs. Then, zoospores multiply to form the primary plasmodia. During this primary phase of infection, metabolisms of cell wall chitin digestion, starch, citrate cycle, pentose phosphate pathway, pyruvate, trehalose, carbohydrates, and lipids are highly active (Schwelm et al., 2015a,b; Bi et al., 2016) making an interaction between carbon metabolism and N supply possible. In the secondary infection phase, secondary zoospores are released from the root hairs and directly invade host cortical cells, then forming the secondary plasmodia (Kageyama and Asano, 2009). During the second phase of infection, genes involved in basal and lipid metabolism (Bi et al., 2016), in signal transduction pathways (Bi et al., 2019), and in activation of cell differentiation, growth, reproduction, as well as genes coding for proto-oncogenes (Bi et al., 2018) are highly expressed. At this stage, *P. brassicae* forms galls in root because of hypertrophy and hyperplasia of infected roots by the modification of hormone levels (Siemens et al., 2006), leading to the obstruction of nutrient and water transport and then a reduction of plant growth. During both primary and secondary phases, it has been suggested that *P. brassicae* could secrete effectors, but the validation of their real roles in infection is still in progress (Kombrink and Thomma, 2013; Schwelm et al., 2015b; Rolfe et al., 2016; Zhang et al., 2016; Perez-Lopez et al., 2018; Daval et al., 2019, 2020), except for a predicted secreted methyltransferase interfering in the plant salicylic acid-induced defense (Ludwig-Muller et al., 2015).

Understanding how *P. brassicae* responds to N input in the different stages of its life cycle to influence infection is, then, of fundamental importance. In this context, the aim of the present study is to explore how N supply may influence the functions of the pathogen and its pathogenesis. For this, we followed clubroot disease progression through a time-course experiment [14, 27, and 42 days post-inoculation (dpi)] in two *B. napus* genotypes (Yudal and HD-018) under two N
supply conditions (1 and 8 mM). The time-course allowed to analyze the N effect at both stages of the *P. brassicae* life cycle. The host genotypes were chosen because of their different responses to clubroot according to N nutrition: both were susceptible to clubroot in normal N supply (8 mM), and only the genotype Yudal showed a drastic reduction of the disease symptoms under N constraint (1 mM) (Laperche et al., 2017). To identify underlying molecular pathogen mechanisms potentially involved in the modulation of the disease by N input, we assessed (i) the phenotype of the two host-plant genotypes harboring different responses to N starvation and (ii) the transcriptome of the pathogen using RNA-sequencing (RNA-seq).

# MATERIALS AND METHODS

## Plant Material and Pathogen Inoculation

Two oilseed rape genotypes, Yudal (called Y) and HD-018 (called HD), were used in this study: the *B. napus* cv. Yudal is a spring oilseed rape and HD-018 is a double haploid line from the cross between Darmor-bzh and Yudal (Laperche et al., 2017). The eH isolate of *P. brassicae* belonging to pathotype P1 (Some et al., 1996; Fahling et al., 2003; Daval et al., 2019) was used for inoculation.

Seeds from the two *B. napus* genotypes were sown in pots filled with “Falienor 922016F3” medium (65% Irish peat, 20% black peat, 15% perlite, and 2% clay). Pots were maintained in a phytotron at 22°C (day) and 19°C (night) with a 16-h photoperiod. Two Hoagland nutritive solutions were used, one containing 1 mM of nitrogen (named N1) and the other 8 mM (named N8) (**Supplementary Table 1A**), in which nitrogen was present only as nitrate (without ammonium), the preferred nitrogen source of oilseed rape. Pots were watered periodically at the bottom of the pots twice a week from sowing to 21 days after sowing, and then three times a week. This frequency was based on the requirement of the plants under N1 condition. Since transpiration and water requirement were higher in the N8-treated plants, those plants were also additionally supplied with water when needed, so as to avoid water stress. Plants were inoculated with a resting spore suspension of *P. brassicae* eH isolate. Non-inoculated (healthy) plants were also grown to make sure that no symptoms appeared and that the plants were growing well. For plant inoculation, clubs propagated on the universal susceptible host Chinese cabbage (*B. rapa* spp. *pekinensis* cv. Granaat) were collected, homogenized in a blender with sterile water, and separated by filtration through layers of cheesecloth. The resting spores were then separated by filtration through 500-, 100-, and 55-µm sieves to remove plant cell debris. The spore concentration was determined with a Malassez cell and adjusted to 1 × 10⁷ spores ml⁻¹. Inoculation was done as described previously (Manzanares-Dauleux et al., 2000): 7-day-old seedlings were inoculated by pipetting 1 ml of the spore suspension onto the soil surface at the base of each seedling. Experiments were conducted in four replicates. Inoculated root samples for nucleic acid extraction were from 24 pooled plants at 14 dpi (primary phase of infection), 12 pooled plants at 27 dpi

(secondary phase of infection), and 6 pooled plants at 42 dpi, for each combination of N level and genotype, and for each of the four replicates (total = 16 samples per time-course point, i.e., 48 in all). An additional six plants per modality were grown for aerial and root biomasses weighing.

## Phenotyping: Plant Characterization and Disease Assessment

At each sampling date and for each replicate, the aerial and root parts of the plants dedicated to these measures were cut, dried, and weighed. After measurement of the disease (symptoms and pathogen quantification as described below), the roots of plants grown for nucleic acid extraction were cut, washed in sterile water, sliced into small pieces, frozen in liquid nitrogen, lyophilized, and grounded with a FastPrep-24 (MP biomedicals). The dry root biomass was measured and the powder was kept at -80°C until nucleic acid extraction (DNA for pathogen quantification and RNA for RNA-seq analyses).

Disease was assessed at each sampling date after inoculation with *P. brassicae*. First, clubroot symptoms were evaluated by a disease index calculated with the scale previously described (Manzanares-Dauleux et al., 2000). Secondly, 1 µl of DNA extracted from root samples (see next paragraph) was used for quantitative PCR on the 18S gene to quantify the *P. brassicae* amount as previously described (Daval et al., 2020).

For statistical analyses between modalities of aerial and root biomasses, disease data, and 18S gene levels, linear models were used (LMM function “lmer,” package “lme4”) (Bates et al., 2015). Estimated Marginal Means (EMMeans) were calculated using the “emmeans” function of the “emmeans” package (Searle et al., 1980). Pairwise comparisons of EMMeans were performed (α = 5%), using the “cld” function of the “emmeans” package.

## Nucleic Acids Isolation From Roots

The lyophilized and ground roots from pooled plants of each genotype and each treatment (24, 12, and 6 plants at 14, 27, and 42 dpi, respectively) were used for nucleic acid extraction.

DNA was extracted and used for *P. brassicae* quantification by qPCR (**Supplementary Table 1B**) as previously described (Daval et al., 2020).

Total RNA was extracted from 20 mg of lyophilized powder with the Trizol protocol (Invitrogen). RNA purity and quality were assessed with a Bioanalyzer 2100 (Agilent) and quantified with a Nanodrop (Agilent).

## Library Construction and Illumina Sequencing

RNA-seq analysis was performed on RNA extracted from root tissues of both genotypes Y and HD, infected with resting spores of *P. brassicae* (eH isolate), and grown in the two different N conditions (N1 and N8) for four replicates, at 14, 27, and 42 dpi.

The preparation of Shotgun Random mRNA was conducted by Genoscreen (Lille, France). Library sequencing was conducted from both ends on an Illumina HiSeq4000 (Genoscreen, Lille, France) using 2 × 150 bp. The purified mRNA was fragmented and converted into double-stranded cDNA with random
priming. Following end-repair, indexed adapters were ligated. The cDNA fragments of ~350 bp were purified with AMPure XP beads and amplified by PCR to obtain the libraries for sequencing. The libraries were multiplexed (six libraries per lane) and sequenced, resulting in an expected average of 80 paired-end millions of reads per library (and then per sample). Raw reads of the 48 samples are available at the European Nucleotide Archive database system under the project accession number PRJEB44381 (samples ERS6263817-ERS6263820, ERS6263825-ERS6263828, ERS6263833-ERS6263836, ERS6263841-ERS6263844, ERS6263849-ERS6263852, ERS6263857-ER S6263860, ERS6263865-ERS6263868, ERS6263873-ERS6263876, and ERS6263881-ERS6263896).

## Read Mapping to the Reference Genome and Transcript Counting

The read quality was undertaken for the quality scores of Q28 and for the read length of 50 nucleotides using PrinSeq. The high-quality reads were aligned to the reference genome of eH *P. brassicae* (Daval et al., 2019) files using STAR 2.5.2a_modified. Non-default parameters were minimum intron length 10, maximum intron length 25,000, and maximum gap between two mates 300. Thanks to genome annotation files, the mapped sequencing reads were assigned to genomic features using featureCounts v1.5.0-p1 and counted. After filtering of the read counts below the threshold value (at least 0.5 counts per million in four samples), the count reads were then normalized with the TMM method. As the numbers of reads in the libraries were very different between each sampling time (due to the different infection rates and progression of the pathogen during the kinetic), the normalization was performed separately for each time. So, analyses of *P. brassicae* transcriptomes were specific to each sampling time, preventing the data comparison between the time points.

## Identification of Differentially Expressed Genes (DEGs) and Gene Co-expression Clustering Analysis

Analyses were performed using the R package AskoR developed in the laboratory<fnref n="1" />, allowing the integration of several packages widely used in gene expression in an automatic way. In this package, Edge R package was used for differential expression analysis and a false discovery rate (FDR) cutoff of 0.05 was applied to account for multiple testing correction. The DEGs with FDR ≤ 0.05 from specific comparison lists were selected for analysis. AskoR generated heat maps and Venn diagrams using the heatmap3 and the VennDiagram package, respectively. The clustering of gene expression profiles was performed, through AskoR, using CoSeq package (Rau and Maugis-Rabusseau, 2018; Godichon-Baggioni et al., 2019) by making groups of co-expressed genes on the EdgeR normalized counts. The UpSetR package was also integrated into AskoR for the visualization of intersecting sets of DEGs (Conway et al., 2017).

## Gene Ontology Term (GO-term) Enrichment

Gene Ontology enrichment analysis of DEGs was achieved with AskoR using TopGO R package (Alexa and Rahnenfuhrer, 2007) with weight01 algorithm, Fisher's exact statistic test, and a nodeSize parameter set to 5 (to remove enriched GO-term with less than five genes in the genome). For each GO category (Molecular Function, Cellular Component, and Biological Process), the top 15 enriched GO-terms (p-value < 0.05), enrichment ratios (>2), and number of genes under each enriched GO-term were represented using the ggplot2 R package.

## Quantitative Real-Time PCR (RT-qPCR) Validation

The expression levels of 16 DEGs were determined by RT-qPCR to confirm the results of RNA-seq analysis. Total RNA from roots infected by *P. brassicae* was reverse transcribed with a set of two external RNA quality controls (Sabater-Munoz et al., 2006) as previously described (Daval et al., 2011, 2013). Seven hundred and fifty ng of total RNA were mixed with known quantities of the two external controls. Reverse transcription was carried out in 30 µl containing 375 ng of random primers, 1 × ImPromII reaction buffer, 3 mM MgCl₂, 125 µM of each dNTP, 30 U of RNasin Ribonuclease Inhibitor, and 1.5 µl of ImProm-IITM (Promega). The following parameters were applied: 5 min at 25°C, 1 h at 42°C, and 15 min at 70°C. Reactions without RNA or without reverse transcriptase were performed as negative controls. The oligonucleotides designed with the Primer 3 software are described in **Supplementary Table 1B**. Quantitative PCR reactions (20 µl) containing 1 µl of cDNA, 0.4 µM of each primer, and 1 × SybrGreen I Master (Roche) were performed on the LightCycler1480 Real-Time PCR System (Roche). The quantitative PCR profile consisted of an initial denaturation at 95°C for 5 min, followed by 45 cycles of 95°C for 15 s, and hybridization–elongation temperature for 40 s (Supplementary Table 1B). A dissociation stage was applied at the end of the PCR to assess that each amplicon generated was specific. Moreover, each specific amplicon was sequenced (Genoscreen, Lille, France) to confirm it corresponded to the expected sequence. The expression levels of transcripts were normalized using external RNA controls from four independent biological replicates, each with three technical PCR replicates. For this normalization, the expression levels of each quantified transcript were calculated from the $C_T$ values (Daval et al., 2011). Using the geNorm software, an accurate gene expression normalization factor of qPCR data was calculated based on the geometric average of the two external controls $C_T$. For each candidate gene and each sample, the ratio between the relative expression value and this normalization factor gave the corrected expression level of the gene. Data were analyzed using the ANOVA procedure of the R statistical analysis software. The correlation between the two sets of data corresponding to gene expressions obtained by RNA-seq and by RT-qPCR was determined using the Spearman rank correlation test.

<footnote n="1">¹https://github.com/askomics/askoR</footnote>
# RESULTS

## Modulation of Plant Clubroot Symptoms According to the Nitrogen Supply and the Host-Plant Genotype

The dry aerial parts of inoculated plants (Figure 1A) increased all over the entire experiment, especially in the N8 condition. For each N condition and each time-course point, no significant differences were measured between Y and HD genotypes. From 14 dpi, the plants grown in N8 displayed significantly higher aerial biomasses than those in the N1 condition, to reach in both genotypes an aerial dry weight 6–7 times greater in N8 than in N1 at the end of the experiment. Non-inoculated plants displayed comparable growth profiles to those of inoculated plants (data not shown) in function of time and N supply. The roots (Figure 1B) exhibited similar profiles to the aerial parts, except at 42 dpi in both N conditions where HD had a slightly higher weight than Y. Thus, there was no important effect of genotype but a high effect of N condition on the growth profile.

At each sampling time, disease severity was scored by determining the disease index (DI) and the DNA pathogen content (Figure 1C). The non-inoculated plants showed no symptoms and did not contain any clubroot pathogen according to our detection threshold (data not shown). In inoculated plants, the DI value (not yet measurable at 14 dpi) increased along the time points, but this increase was smaller for Y at N1 compared to the three other conditions. The DI was twice smaller ($p < 0.05$) in Y in N1 compared to HD in both N fertilization conditions and to Y in N8. The DNA pathogen content, well detectable from 27 dpi, followed almost the same pattern and increased over the time-course experiment. Although the reduction of symptoms was more pronounced between Y N1 and the three other conditions from 27 dpi, the DNA content was significantly lower only in Y N1 compared to HD N1 at 27 dpi, but significantly lower in Y N1 than all the three other conditions at 42 dpi.

## Overview, Mapping, and Validation of RNA-seq Data

The impact of N on *P. brassicae* transcriptome was investigated by RNA-seq analysis of inoculated roots. We obtained about 80 million reads by sample. The percentage of reads mapped to the eH reference genome increased with the experimental time course and was concordant with the observed increase in the pathogen amount during the infection (Supplementary Figure 1A). The standardized data per million counts were plotted on a PCA to estimate the variability of the biological experiments (Supplementary Figure 1B). Sixty-four percent of the variance were represented on the plot. The x-axis (55.48%) separated clearly the sampling time. The four replicates of each experimental condition were overall clustered together, which made it possible to validate our RNA-seq experiment. Pathogen gene expression's profiles, described by heatmaps at each sampling time (Supplementary Figure 1C), were clustered at 14 dpi first by the N condition and then by the host genotype, and at 27 and 42 dpi first by the host genotype and then by the N supply condition.

A validation of the RNA-seq data was performed by RT-qPCR assay whose amplification conditions are described in Supplementary Table 1B. For this, 16 genes, which were significantly differentially expressed at least in one of the six main comparisons (described in Table 1) and which provided a large range of expression levels, were used to compare log2FoldChange values between RT-qPCR and RNA-seq experiments (Supplementary Figure 2). We found a moderate correlation ($r^2 = 0.48$) between log2FoldChanges values of both types of experiments but a significant relationship according to Spearman correlation test ($p$-value = 6.11 $10^{-14}$), confirming the validity of the RNA-seq data and a reliable expression result generated by RNA-seq.

## Global Overview of the Differential Gene Expression Profiles

Table 1 shows the number of DEGs in *P. brassicae* in the main comparisons between conditions of interest, combining host genotype and N supply. For both host genotypes, N input had a low impact on the *P. brassicae* transcriptome at 14 dpi with only 19 and 37 DEGs when HD and Y were host, respectively. Interestingly, for the later time-course points, *P. brassicae* displayed a different profile of DEGs between N conditions, depending on the host genotype with a high number of DEGs at 27 dpi when infecting HD (600 DEGs between N1 and N8) and a high number of DEGs at 42 dpi when infecting Y (318 DEGs between N1 and N8). When both factors (N and host genotypes) were considered, a gradual increase over time of the DEG number between Y N8 and HD N1 occurred with 128, 316, and 544 at 14, 27 and 42 dpi, respectively. Concerning the comparison Y N1 versus HD N8, it is interesting to note that the number of DEGs was high, especially at 27 dpi, with 1,985 DEGs corresponding to about 20% of the total genes expressed.

To further decipher the transcriptomic fingerprints of *P. brassicae* as a function of plant genotype and nitrogen supply, we plotted the intersections between the lists of DEGs for each of the six main comparisons of interest on an Upset graph (Figure 2A). Among the lists of DEGs whose total numbers are given in Table 1, the higher number of specific genes differentially expressed (not encountered in other comparisons) was found in the HD N8/Y N1 comparison, with 146, 997, and 379 specific DEGs at 14, 27, and 42 dpi, respectively. This was followed by the HD N1/Y N8 comparison that showed the largest number of specific DEGs, with 93 and 360 genes at 14 and 42 dpi, respectively. The proportion of DEGs in at least two comparisons relative to the total number of DEGs was 24% at 14 dpi (79 genes out of 336), 45% at 27 dpi (1,101 genes out of 2,449), and 31% (392 genes out of 1,245) at 42 dpi. Thus, on the other hand, at each time-course point, the DEGs were mostly specific for a given contrast.

A clustering of DEGs considering the genes differentially expressed in at least one comparison allowed description of the different expression profiles at each time point (Figure 2B). Profiles I and II were found at the three time-course points of the experimentation. They displayed similar expression in HD N1 and Y N8, but they differed in the expression of HD N8 compared
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Influence of nitrogen supply on plant biomass and clubroot development. The dry aerial parts (A) and roots (B) were weighed for both inoculated genotypes (Y, Yudal and HD, HD-018) at different days post-inoculation (14, 27, or 42 dpi), and for each N condition supply (N1, 1 mM and N8, 8 mM nitrogen). The clubroot symptoms were estimated according to the disease index and the quantification of *P. brassicae* DNA by qPCR expressed as a ratio of the 18S DNA quantity relative to the total DNA (C). Data are means of four biological replicates (6–24 pooled plants per time, per plant genotype, and per replicate), and error bars represent standard errors of the means. Means with different letters are statistically significantly different according to the analysis of variance test (*P* < 0.05).</figcaption></figure>

<table><caption>**TABLE 1** | Number of *Plasmodiophora brassicae* DEGs for each comparison of interest.</caption><thead><tr><th rowspan="2">Infection stage</th><th colspan="6">Comparison between conditions</th></tr><tr><th>Y N1<br>compared to<br>Y N8</th><th>HD N1<br>compared to<br>HD N8</th><th>Y N1<br>compared to<br>HD N1</th><th>Y N8<br>compared to<br>HD N8</th><th>Y N1<br>compared to<br>HD N8</th><th>Y N8<br>compared to<br>HD N1</th></tr></thead><tbody><tr><td>14 dpi</td><td>37</td><td>19</td><td>55</td><td>1</td><td>194</td><td>128</td></tr><tr><td>27 dpi</td><td>99</td><td>600</td><td>742</td><td>450</td><td>1985</td><td>316</td></tr><tr><td>42 dpi</td><td>318</td><td>5</td><td>119</td><td>140</td><td>613</td><td>544</td></tr></tbody></table>

DEGs, differentially expressed genes; dpi, days post-inoculation; Y, Yudal genotype; HD, HD-018 genotype; N1, 1 mM nitrogen supply; N8, 8 mM nitrogen supply.

to that of Y N1: this expression was higher in HD N8 than Y N1 in profile I but lower in profile II. Profiles III and IV were found at 14 and 42 dpi but not at 27 dpi. Profile III was characterized by the weakest expression in HD N1 and the strongest in Y N8, and profile IV was characterized by the opposite, and both profiles displayed intermediate and equal expression in HD N8 and Y N1. Gene expression in HD N8 and Y N1 therefore seemed to be decisive in determining the DEGs clustering. Profile V was specific of the time 27 dpi, with the lowest expression for HD N1 and HD N8, and the highest for Y N1 and Y N8. For this descriptive global overview of expression profiles, the lists

of DEGs were not specifically extracted, because in total, they represented the complete list of DEGs of the experiment.

## Effect of the Nitrogen Supply Condition on the *P. brassicae* Transcriptome Independently of the Infected Host-Plant Genotype

The description of the DEG numbers between nitrogen supplies is depicted in **Table 2**. The number of DEGs between N1 and N8 for both genotypes was low at 14 dpi (52 genes) with almost
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Global view of *Plasmodiophora brassicae* transcriptome profiles at the three time-course points. (A) Plot of intersections between sets of DEGs. The intersected lists of DEGs and their sizes (number of DEGs) are presented in the horizontal bars on the left. The connected lines among the lists represent their intersections. The vertical bars and associated numbers correspond to specific overlap of DEGs sets. The significance of gene expression changes was inferred based on an adjusted *p*-value < 0.05. (B) Co-expression clusters. CoSeq package was used to cluster expression data. The number of genes assigned to a particular cluster is indicated with the cluster name. Y-axis represents, for each gene, the number of centered reduced counts (previously filtered and normalized with EdgeR). DEGs, differentially expressed genes; Y, Yudal; HD, HD-018; N1, 1 mM N supply; N8, 8 mM N supply.</figcaption></figure>

3/4 being underexpressed in N1 compared to N8. The largest number of DEGs was found at 27 dpi (645 genes in total), with approximately equal number of genes overexpressed at N1 (361 genes) and at N8 (284 genes). At the final time-course point, of the 321 DEGs in total, the majority were underexpressed at N1 (224 genes).

A GO-term enrichment analysis was achieved from the 52, 645, and 321 DEGs at each sampling time (Figure 3) and the genes responsible for the enrichment of the main GO-terms, as well as their expression levels in each condition, are described in **Supplementary Table 2**. The GO-terms enriched at 14 dpi were involved in general cell divisions: protein localization, a crucial mechanism controlling access to and availability of all types of molecular interaction partners, and actin binding, contributing to many cellular structure, functions, and motility. In the same way, at 27 dpi, among the 645 DEGs linked to nitrogen supply, 19 genes were involved in mitosis-linked process (“actin binding,” “histone deacetylase activity,” “histone binding,” “isomerase activity,” and “acetyltransferase activity”) and 8 were involved in methyltransferase activity, a key component of genetic regulation. Interestingly, the molecular function “fucosyltransferase activity” was enriched at

27 dpi, and enzymes were classified on the basis of sequence similarities in the CAZy (carbohydrate-active enzymes) and were important for host-microbe interactions. Finally, at 42 dpi, many mechanisms involved in cell divisions and cell-signaling systems were enriched (“nucleosome assembly,” “DNA biosynthetic process,” “chromosome organization,” “DNA recombination,” “DNA replication,” or “microtubule-based movement”). The enriched molecular function “unfolded protein binding” was found, corresponding to a key cellular mechanism to restore protein homeostasis.

DEGs between N1 and N8 were globally specific for a given time (Figure 4): 28 genes were differentially expressed exclusively at 14 dpi, 566 genes (out of 645) exclusively at 27 dpi, and 256 genes (out of 321) exclusively at 42 dpi. At each time-course point, a small number of genes were differentially expressed between N1 and N8, concurrently in Y and HD: 4 genes at 14 dpi (i.e., 7.7%), 54 genes at 27 dpi (i.e., 8.4%), and 2 genes at 42 dpi (i.e., 0.6%) (Table 2). These 60 genes showed the same expression profile between N1 and N8 for both host genotypes.

To go further in the description of mechanisms potentially modulated by N supply, we have particularly studied the common sets of DEGs between N1 and N8 (among the complete list of
<table><caption>**TABLE 2** | Number of *Plasmodiophora brassicae* DEGs depending on the nitrogen supply.</caption><thead><tr><th></th><th>Total</th><th>Underexpressed in N1 compared to N8</th><th>Overexpressed in N1 compared to N8</th><th>Common DEGs in both host genotypes</th></tr></thead><tbody><tr><td>14 dpi</td><td>52</td><td>37</td><td>15</td><td>4</td></tr><tr><td>27 dpi</td><td>645</td><td>361</td><td>284</td><td>54</td></tr><tr><td>42 dpi</td><td>321</td><td>224</td><td>97</td><td>2</td></tr></tbody></table>

DEGs, differentially expressed genes; dpi, days post-inoculation; N1, 1 mM nitrogen supply; N8, 8 mM nitrogen supply.

<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | GO enrichment on DEGs lists according to the nitrogen effect at each time-course point. DEGs, differentially expressed genes; dpi, days post-inoculation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | Number of *Plasmodiophora brassicae* DEGs between N supply conditions. The Venn diagrams show the number of significantly DEGs (*P* < 0.05) between N1 and N8 for both genotypes between sampling dates. DEGs, differentially expressed genes; dpi, days post-inoculation; N1, 1 mM nitrogen supply, N8, 8 mM nitrogen supply.</figcaption></figure>

DEGs between N1 and N8 described in **Supplementary Table 3** whose expression profiles quantitatively varied in function of N conditions but independently of host genotype.

At 14 dpi (**Supplementary Table 3A**), among the four DEGs in both HD and Y, one corresponded to a hypothetical protein, one was overexpressed in N1 compared to N8 and encoded a “superoxide dismutase,” and two genes underexpressed in N1 compared to N8 encoded “pentatricopeptide repeat-containing,” whose putative function could be involved in growth.

At 27 dpi (**Supplementary Table 3B**), among the 54 common genes, 32 were related to hypothetical proteins. From the 22 remaining genes identified as potential known proteins, 18 were underexpressed in N1 compared to N8, as those encoded for “AGC kinase” involved in development, “ankyrin repeat” (two genes), “serine carboxypeptidase,” or “NUDIX hydrolase,” a domain present in “Nucleoside diphosphate-linked moiety X (NUDIX)” effectors, a group of secreted proteins identified in a broad range of pathogens, “4-pyridoxolactonase” (EC3.1.1.27) involved in the vitamin B6 degradation, several “heat shock proteins” and “disulfide isomerase” that can function as intra-cellular chaperones for other proteins and as part of cellular stress response, “ion transporter,” “cytochrome C oxidase,” or “methyltransferase” involved in many reactions such as cellular respiration, and “NEP1-interacting-like-2-isoform X1” (Pldbra_eH_r1s029g11071) known for their cytotoxicity in plants. The four genes overexpressed in N1 compared to N8 referred to “NHL-repeat binding protein” found in a large number of eukaryotic and prokaryotic proteins (like serine/threonine protein kinase in some pathogenic bacteria), a ubiquitous “transcription factor,” a “glutamine synthetase” that plays an essential role in the metabolism of nitrogen, and a “histone H2B” involved in mitosis.
Finally, at 42 dpi (Supplementary Table 3C), only two genes were commonly differentially expressed when both host genotypes were infected, one encoded a hypothetical and one underexpressed in N1 compared to N8 and encoded for "sporangia induced dynein heavy" playing a role in spore production.

## Nitrogen Supply Effect on the *P. brassicae* Transcriptome According to the Infected Host-Plant Genotype

The number of DEGs between N1 and N8 according to the host genotypes is presented in **Figure 5A** for each sampling date. At 14 dpi, Y has twice as many genes as HD, and even the total DEGs remained low. At 27 dpi, the majority of DEGs between N conditions were found when genotype HD was infected, with 600 genes out of 645 (93%), whereas in contrast at 42 dpi, the majority of DEGs were found when the infected genotype was Y, with 318 out of 321 (99%). The number of *P. brassicae* DEGs common to both infected genotypes was low, showing a nitrogen regulation that is dependent on the host-plant genotype. The functions of the genes identified at 27 or 42 dpi and described in the section above can thus finally be linked to the HD and Y genotypes, respectively.

A more in-depth analysis was performed on DEGs at 27 dpi when HD was infected and at 42 dpi when Y was infected, taking into account the over- or underexpression depending on the N condition (Figure 5B). In both categories of *P. brassicae* genes (those overexpressed and those underexpressed in N1 compared to N8), few genes were common between DEGs at 27 dpi in infected HD and at 42 dpi in infected Y. This showed that the difference in the number of DEGs was not due to a simple time lag in infection between the two genotypes, but to a difference in the response of *P. brassicae* to nitrogen supply as a function of host genotype and interaction time.

## Identification of *P. brassicae* Transcriptome Features Specific to Low-Nitrogen-Induced Resistance in Yudal

To better understand the potential underlying mechanisms related to the difference of the low infection level of Y in the N1 supply compared to the high symptom level in the other conditions (Y in N8 and HD in both N supplies), the specific DEGs between Y N1 and the other three conditions were studied at each interaction time (Figure 6A). The largest number of DEGs was found at 27 dpi with 2227 DEGs compared to 237 at 14 dpi and 781 at 42 dpi. For each time-course point, the greatest number of DEGs was found between Y N1 and HD N8, corresponding to both the N effect and the genotype effect, this contrast representing 78-89% of the total number of DEGs. Genes differentially expressed simultaneously between Y N1 and the three other conditions were only 7 at 14 dpi, 55 at 27 dpi, and 40 at 42 dpi, i.e., 102 in total, described in Supplementary Table 4. In addition, Figure 6B shows that these genes were mostly specific

for each interaction time point since only 2 DEGs were common to 14 and 42 dpi, and 1 gene was common to 27 and 42 dpi. No gene common to all time-course points was differentially expressed between Y N1 and the three other conditions.

Among the 102 DEGs of **Figure 6A**, 52 were not identified ("hypothetical protein," Supplementary Table 4), and 50 were identified by a GO-Term (Table 3 and Supplementary Table 4). Among these 50 genes, 9 were overexpressed and 41 were underexpressed in Y N1 condition.

At 14 dpi, the four identified genes were underexpressed in Y N1 and encoded an "ATP GTP-binding protein kinase" and a "kinase domain" involved in cellular signalization, a "chromodomain-helicase-DNA-binding 7 isoform X1" playing a role in DNA replication, and a "vacuolar sorting-associated 13," a protein of unknown function.

At 27 dpi, many DEGs between Y N1 and the three other conditions (most of them underexpressed in Y N1) were involved in general metabolic pathways, such as: (i) development, cell growth and proliferation ("regulator of chromosome condensation," "AGC kinase," "centrosomal POC5," "30S ribosomal S16," "S-adenosylmethionine decarboxylase proenzyme," "intraflagellar transport 122 homolog," and "Y-family DNA polymerase"), (ii) energy metabolism or cellular respiration ("adenosine triphosphatase," "transcription factor IIH," and "S-adenosylmethionine-dependent methyltransferase"), (iii) chaperones as part of cellular stress response ("heat shock Ssa2," "disulfide isomerase" "luminal-binding 5," and "FKBP-type peptidyl-prolyl cis-trans isomerase"), and (iv) signaling pathway ("iq domain-containing g," and "tyrosine-kinase Src42A"). Such gene expression profiles related to pathogen development were concordant with the lower clubroot symptoms observed in Y under N1 supply. Interestingly, genes potentially involved in the pathogenicity of *P. brassicae* (such as putative effectors) or causing cytotoxicity in the plant were differentially underexpressed in Y N1 ("ankyrin repeat," "serine carboxypeptidase," "NUDIX hydrolase," and "NEP1-interacting-like-2-isoform X1"). In contrast, one gene encoding a potential pathogenicity facilitator ("glycosyltransferase") was overexpressed in Y N1. Finally, it is essential to underline that a gene encoding a "glutamine synthetase," a key enzyme for nitrogen metabolism involved in ammonium assimilation into amino acids, was overexpressed in *P. brassicae* in Y N1 condition compared to the three other conditions.

At 42 dpi, in the same way as at 27 dpi, many genes underexpressed in Y N1 were linked to mechanisms of growth and development of the pathogen ("centrosomal of 120 kDa-like," "centrosomal of 290 kDa-isoform X5," "DNA topoisomerase 1 beta," "chromosome segregation," "PIN2 TERF1-interacting telomerase inhibitor-1," "suppressor SRP40," "histone H1," "Katanin p60 ATPase-containing subunit A-like," "splicing factor 3a subunit 3," "nucleolar NOP5," "Rad52 22 double-strand break repair," "motor domain-containing," "pre-rRNA-processing esf1," "eukaryotic rRNA processing," "nucleolar 56-like," "protein BUD31 homolog," "viral A-type inclusion," and "like subfamily c-member-21"), to signaling pathway ("CMGC CDKL kinase"), or related to stress response ("calreticulin
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**FIGURE 5** | Number of *Plasmodiophora brassicae* DEGs between N1 and N8 according to the host genotype. The Venn diagrams show the number of significantly DEGs (*P* < 0.05) between N1 and N8 for each host *B. napus* genotype at each sampling date (A), and for both genotypes between sampling dates (B). DEGs, differentially expressed genes; dpi, days post-inoculation; Y, Yudal genotype; HD, HD-018 genotype; N1, 1 mM nitrogen supply; N8, 8 mM nitrogen supply.

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**FIGURE 6** | Number of *Plasmodiophora brassicae* DEGs between Y N1 and each of the other three conditions. (A) The Venn diagrams show the number of significantly DEGs (*P* < 0.05) between Y N1 and Y N8, or HD N1, HD N8 at each sampling date. (B) The Venn diagram focuses on the DEGs common between Y N1 and the other three conditions described at each of the three interaction time-course points (i.e., 7 genes at 14 dpi, 55 genes at 27 dpi, and 40 genes at 42 dpi). DEGs, differentially expressed genes; dpi, days post-inoculation; Y, Yudal genotype; HD, HD-018 genotype; N1, 1 mM nitrogen; N8, 8 mM nitrogen supply.
<table><caption>TABLE 3 | *Plasmodiophora brassicae* genes simultaneously differentially expressed in Y N1 compared to the three other conditions.</caption><thead><tr><th rowspan="2">Interaction time</th><th rowspan="2">P. brassicae gene</th><th colspan="4">P. brassicae gene expression level</th><th rowspan="2">Description</th></tr><tr><th>in Y_N1</th><th>in Y_N8</th><th>in HD_N1</th><th>in HD_N8</th></tr></thead><tbody><tr><td rowspan="4">14 dpi</td><td>Pldbra_eH_r1s002g00869</td><td>0.22</td><td>12.38</td><td>9.44</td><td>16.12</td><td>ATP GTP-binding protein kinase</td></tr><tr><td>Pldbra_eH_r1s008g04491</td><td>106.00</td><td>185.05</td><td>182.15</td><td>268.57</td><td>Kinase domain</td></tr><tr><td>Pldbra_eH_r1s017g08107</td><td>273.64</td><td>430.81</td><td>400.69</td><td>434.53</td><td>Chromodomain-helicase-DNA-binding 7 isoform X1</td></tr><tr><td>Pldbra_eH_r1s026g10379</td><td>32.25</td><td>83.65</td><td>84.29</td><td>84.58</td><td>Vacuolar sorting-associated 13</td></tr><tr><td rowspan="22">27 dpi</td><td>Pldbra_eH_r1s002g00900</td><td>2.17</td><td>0</td><td>0.36</td><td>0</td><td>Regulator of chromosome condensation (RCC1)</td></tr><tr><td>Pldbra_eH_r1s003g02060</td><td>1824.32</td><td>2679.49</td><td>2782.38</td><td>4241.63</td><td>family</td></tr><tr><td>Pldbra_eH_r1s004g02133</td><td>6.28</td><td>12.74</td><td>13.67</td><td>29.82</td><td>AGC kinase</td></tr><tr><td>Pldbra_eH_r1s005g02782</td><td>76.76</td><td>124.45</td><td>110.88</td><td>175.18</td><td>Ankyrin repeat</td></tr><tr><td>Pldbra_eH_r1s005g02887</td><td>58.66</td><td>107.09</td><td>37.02</td><td>104.68</td><td>Adenosine triphosphatase</td></tr><tr><td>Pldbra_eH_r1s006g03454</td><td>22.13</td><td>35.08</td><td>42.93</td><td>52.57</td><td>Heat shock Ssa2</td></tr><tr><td>Pldbra_eH_r1s006g03599</td><td>2839.60</td><td>4625.62</td><td>3867.11</td><td>7217.39</td><td>Ankyrin repeat</td></tr><tr><td>Pldbra_eH_r1s006g03625</td><td>79.35</td><td>108.00</td><td>102.69</td><td>117.83</td><td>Luminal-binding 5</td></tr><tr><td>Pldbra_eH_r1s008g04647</td><td>2.33</td><td>0.12</td><td>0.08</td><td>0</td><td>30S ribosomal S16</td></tr><tr><td>Pldbra_eH_r1s011g06084</td><td>81.67</td><td>109.84</td><td>116.44</td><td>130.85</td><td>Centrosomal POC5 (macronuclear)</td></tr><tr><td>Pldbra_eH_r1s012g06550</td><td>91.93</td><td>64.61</td><td>63.40</td><td>56.29</td><td>S-adenosylmethionine decarboxylase</td></tr><tr><td>Pldbra_eH_r1s014g07249</td><td>118.64</td><td>163.23</td><td>204.29</td><td>268.61</td><td>proenzyme-like</td></tr><tr><td>Pldbra_eH_r1s015g07600</td><td>12.42</td><td>27.25</td><td>23.66</td><td>34.14</td><td>General transcription factor IIH</td></tr><tr><td>Pldbra_eH_r1s017g08043</td><td>620.67</td><td>914.36</td><td>841.54</td><td>1154.60</td><td>Iq domain-containing g</td></tr><tr><td>Pldbra_eH_r1s018g08464</td><td>312.47</td><td>428.92</td><td>446.17</td><td>589.06</td><td>Intraflagellar transport 122 homolog</td></tr><tr><td>Pldbra_eH_r1s018g08619</td><td>7387.43</td><td>10149.25</td><td>10588.46</td><td>14572.31</td><td>Serine carboxypeptidase family</td></tr><tr><td>Pldbra_eH_r1s019g08834</td><td>156.34</td><td>108.15</td><td>97.69</td><td>57.93</td><td>Ankyrin repeat-containing domain</td></tr><tr><td>Pldbra_eH_r1s022g09660</td><td>178.76</td><td>293.35</td><td>250.30</td><td>373.98</td><td>S-adenosylmethionine-dependent</td></tr><tr><td>Pldbra_eH_r1s027g10649</td><td>691.14</td><td>1348.12</td><td>1151.69</td><td>2022.38</td><td>methyltransferase At5g38100</td></tr><tr><td>Pldbra_eH_r1s028g10823</td><td>40.83</td><td>23.77</td><td>23.22</td><td>18.69</td><td>Glutamine synthetase</td></tr><tr><td>Pldbra_eH_r1s029g11027</td><td>4.42</td><td>10.58</td><td>16.86</td><td>17.56</td><td>FKBP-type peptidyl-prolyl cis-trans isomerase</td></tr><tr><td>Pldbra_eH_r1s029g11071</td><td>436.48</td><td>574.60</td><td>592.03</td><td>894.01</td><td>Disulfide isomerase</td></tr><tr><td>Pldbra_eH_r1s033g11505</td><td>4.57</td><td>0.59</td><td>0.97</td><td>0.66</td><td>Y-family DNA polymerase</td></tr><tr><td>Pldbra_eH_r1s037g11877</td><td>998.20</td><td>1806.81</td><td>2029.42</td><td>3731.75</td><td>Tyrosine-kinase Src42A</td></tr><tr><td rowspan="24">42 dpi</td><td>Pldbra_eH_r1s001g00113</td><td>73.64</td><td>102.51</td><td>97.98</td><td>97.23</td><td>Centrosomal of 120 kDa-like</td></tr><tr><td>Pldbra_eH_r1s001g00127</td><td>129.25</td><td>171.67</td><td>162.73</td><td>167.09</td><td>Like subfamily c-member-21</td></tr><tr><td>Pldbra_eH_r1s001g00331</td><td>243.77</td><td>412.32</td><td>411.74</td><td>407.37</td><td>Centrosomal of 290kDa-isoform X5</td></tr><tr><td>Pldbra_eH_r1s001g00618</td><td>187.01</td><td>279.84</td><td>263.76</td><td>258.39</td><td>DNA topoisomerase 1 beta</td></tr><tr><td>Pldbra_eH_r1s001g00659</td><td>230.32</td><td>325.73</td><td>331.44</td><td>360.67</td><td>Chromosome segregation</td></tr><tr><td>Pldbra_eH_r1s002g00895</td><td>377.54</td><td>272.52</td><td>259.27</td><td>251.73</td><td>L-isoaspartate (D-aspartate)</td></tr><tr><td>Pldbra_eH_r1s002g01266</td><td>140.81</td><td>228.54</td><td>215.10</td><td>231.26</td><td>O-methyltransferase-like</td></tr><tr><td>Pldbra_eH_r1s002g01288</td><td>189.29</td><td>351.64</td><td>308.38</td><td>326.88</td><td>PIN2 TERF1-interacting telomerase inhibitor-1</td></tr><tr><td>Pldbra_eH_r1s004g02319</td><td>358.78</td><td>517.01</td><td>526.94</td><td>477.75</td><td>Suppressor SRP40</td></tr><tr><td>Pldbra_eH_r1s005g02870</td><td>28.51</td><td>73.95</td><td>63.75</td><td>67.87</td><td>Calreticulin precursor</td></tr><tr><td>Pldbra_eH_r1s008g04609</td><td>26.03</td><td>16.48</td><td>16.08</td><td>13.09</td><td>Viral A-type inclusion</td></tr><tr><td>Pldbra_eH_r1s009g05132</td><td>0.09</td><td>1.21</td><td>1.07</td><td>1.28</td><td>PAS domain S-box</td></tr><tr><td>Pldbra_eH_r1s009g05217</td><td>59.34</td><td>111.47</td><td>136.50</td><td>121.52</td><td>CMGC CDKL kinase</td></tr><tr><td>Pldbra_eH_r1s011g06038</td><td>70.99</td><td>93.60</td><td>101.08</td><td>93.24</td><td>Histone H1</td></tr><tr><td>Pldbra_eH_r1s012g06458</td><td>162.17</td><td>215.60</td><td>238.30</td><td>212.49</td><td>Katanin p60 ATPase-containing subunit A-like</td></tr><tr><td>Pldbra_eH_r1s013g06611</td><td>257.85</td><td>405.35</td><td>359.86</td><td>391.28</td><td>Splicing factor 3a subunit 3</td></tr><tr><td>Pldbra_eH_r1s014g07146</td><td>43.88</td><td>96.57</td><td>73.44</td><td>88.94</td><td>Nucleolar NOP5</td></tr><tr><td>Pldbra_eH_r1s015g07419</td><td>42.39</td><td>93.37</td><td>91.10</td><td>96.21</td><td>Rad52 22 double-strand break repair</td></tr><tr><td>Pldbra_eH_r1s015g07500</td><td>111.71</td><td>216.12</td><td>204.77</td><td>202.75</td><td>Motor domain-containing</td></tr><tr><td>Pldbra_eH_r1s015g07561</td><td>197.19</td><td>323.11</td><td>293.71</td><td>293.49</td><td>Pre-rRNA-processing esf1</td></tr><tr><td>Pldbra_eH_r1s021g09305</td><td>369.19</td><td>639.38</td><td>549.37</td><td>617.65</td><td>Eukaryotic rRNA processing</td></tr><tr><td>Pldbra_eH_r1s053g12569</td><td>126.38</td><td>96.04</td><td>90.53</td><td>87.83</td><td>Nucleolar 56-like</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Protein BUD31 homolog</td></tr></tbody></table>

Genes underexpressed or overexpressed in Y N1 compared to the three other conditions are in normal font or in italics, respectively.

precursor"). On the other hand, one gene putatively associated to repair/degradation of damaged proteins ["L-isoaspartate(D-aspartate) O-methyltransferase-like"] and one gene connected to

the stress response ("PAS domain S-box") were overexpressed in Y N1 condition compared to the three other conditions displaying higher infection level (Y N8, HD N1, and HD N8).
DISCUSSION

Nitrogen fertilization plays an important role in the outcome
of the plant–pathogen interaction and thus in the development
of the plant diseases. To explain this, it is wise to hypothesize
that nitrogen can play a role on the defense mechanisms of the
plant and/or on the development of the pathogen. Low-nitrogen
conditional resistance was previously reported in Y (Laperche
et al., 2017; Aigu et al., 2018). The genetic analysis of this trait in a
Darmor × Yudal progeny had highlighted a quantitative genetic
control, and the existence of at least two low-nitrogen conditional
resistance loci. To better understand those phenomena and the
clubroot modulation by N supply, it was also important to get a
better idea how plant nitrogen supply can affect the physiology
of *P. brassicae*, and if this influence could be different according
to the host plant genotype. Thus, in the present work, we used
a transcriptomic approach to investigate the pathogen molecular
mechanisms that can underline the observed different levels of
infection severity according to the presence of this low-nitrogen
conditional resistance related to host-plant genotype. In this
study, we aimed at focusing on the pathogen transcriptome
including in an original way simultaneously different factors to
be studied: nitrogen, plant genotypes, and interaction kinetics.

For this, *B. napus* cultures were conducted under two
fertilization regimes, in which the nitrogen was supplied as
nitrate, preferentially taken up by oilseed rape (Arkoun et al.,
2012), and not as ammonium. As potassium and calcium can
have effects on clubroot disease development (Dixon, 2009), the
concentrations of these elements have been adjusted between
the two nutrient solutions, making sure that the comparisons
between N1 and N8 only focused on nitrate content. Finally,
as the experiments were performed in potting soils and not in
natural soils, biological interactions between nitrogen supply and
potential soil microorganisms were limited, allowing accurate
and controlled monitoring of the nitrogen conditions.

First, we confirmed that the low-nitrogen conditional
resistance in Y was absent in HD, consistent with the data
obtained in Laperche et al. (2017) and Aigu et al. (2018).
Furthermore, in the present work, we showed that the low-
nitrogen conditional resistance in Y resulted not only in lower
symptoms but also in lower pathogen content in root tissues.
Both genotypes had similar root and aerial growth profiles
for a given nitrogen condition, except a very slightly less root
growth in Y compared to HD in both nitrogen conditions.
Moderate differences in root growth were already described
in other genotypes of *B. napus* that can display different root
developmental patterns (Daval et al., 2020). Both Y and HD
displayed significantly lower growth at N1 than at N8, showing
that N1 was a nitrogen constraint condition, as expected. The
control non-inoculated plants displayed similar profiles of aerial
and root weights, with no differences between genotypes, and
significantly lower growth at N1 than at N8 for both genotypes
(data not shown). Decreases in symptoms and pathogen content
specific to Y N1 (and not to HD N1) were therefore not simply
due to the impact of low nitrogen fertilization on plant growth.

*P. brassicae* transcriptome was influenced by host-plant
nitrogen supply, but this modulation was very much dependent

on host-plant genotype. Although numerous other studies have
reported the influence of abiotic factors on pathogen gene
expression (Fagard et al., 2014; Zarattini et al., 2021), probably
few of them have reported the need for cautious multi-genotype
evaluation before drawing general conclusions. The regulation of
only 102 of *P. brassicae* genes was clearly associated to the low-
nitrogen conditional resistance (i.e., to low rate pathogen growth
in infected roots).

The modulation of clubroot disease could be explained by the
influence of plant N nutrition on many factors involved in the
different steps of the pathogen epidemiological cycle. The time-
course study of host-pathogen interactions thus allowed us to
highlight how relationships between *P. brassicae* and N could
be established throughout the pathogen life cycle. Depending on
the pathogen development stage (and thus the considered time-
course point), we will discuss mechanisms by which nitrogen
supply could modulate infection. These mechanisms could
include (i) the availability of N compounds for exploitation
by the pathogen, the quantity and nature of host N-based
substrates acquired by the pathogen, particularly in relation to
the obligate biotrophic character of *P. brassicae*; (ii) the metabolic
pathways known to be specifically expressed at the primary or the
secondary phase of the cycle; and (iii) the pathways and molecules
involved in the pathogenicity, virulence, or aggressiveness of the
pathogen. We thereafter discuss how this specific set of genes
suggests a reduced expression of putative effectors and reduced
metabolic and cell growth activities.

At 14 dpi, the N effect on *P. brassicae* transcriptome was
moderate, particularly regarding the DEGs between N1 and N8
that were common to both genotypes, or the DEGs between Y
N1 and the three other conditions (only seven genes for this
last comparison). The gene coding for a "superoxide dismutase"
was overexpressed in N1 compared to N8 whatever the host-
plant genotype. The role of this enzyme has not been rigorously
investigated, but it has been speculated to be important in
the pathogenesis of bacterial and fungal infections (Cox et al.,
2003). It is as if at this early stage of interaction, *P. brassicae*
would attempt to infect more actively *via* this pathway to
counteract the N restriction. Other DEGs at 14 dpi ("ATP
GTP-binding protein kinase," "kinase domain," "chromodomain-
helicase-DNA-binding CHD," and "pentatricopeptide repeat-
containing") were underexpressed in N1 and were mainly
involved in basal cellular processes, including metabolism,
transcription, DNA repair, cell cycle progression, cytoskeletal
rearrangement, and differentiation. The CHD proteins have
thus been identified in a vast array of organisms as diverse
as protists, plants, amphibians, and mammals, and a common
theme for the main role of CHD proteins appears to be
linked to their chromatin-remodeling activity (Hall and Georgel,
2007). Pentatricopeptide repeat proteins, found in all eukaryotes
including protists, are RNA binding proteins with functions in
organelle RNA metabolism. In a protist, they were shown to
be involved in growth (Manna et al., 2013). At this early stage
of interaction corresponding to the primary phase of infection,
the carbon metabolism was particularly active (Schwelm et al.,
2015a,b; Bi et al., 2016) but seemed to be weakly modulated by
the N supply conditions.
All in all, at 14 dpi, N restriction had weak effects on *P. brassicae*, apart from a decrease in the expression of some genes involved in the growth of the pathogen. This is consistent with the fact that at this early stage, the resources allocated by the plant to the pathogen were not yet limiting because the plants showed equivalent levels of development regardless of the N treatment: differences in N resources were not yet very differentiated between treatments at this stage. The two metabolisms, nitrogen and carbon, which are intertwined and co-regulated (Nunes-Nesi et al., 2010), were therefore still little impacted by nitrate input.

At the 27 dpi interaction time, interesting regulations of *P. brassicae* gene expression seemed to occur depending on the N supply and the host plant genotype.

First, genes linked to growth and development were particularly underexpressed in Y N1. For instance, the “AGC-kinase” gene was described as required for development, invasion, and survival of the protist *Plasmodium falciparum* (Arecibia et al., 2013). In the fungus *Flammulina velutipes*, the “S-adenosylmethionine-dependent-methyltransferases,” which can transfer a methyl group to many kinds of biological molecules, play crucial roles in many important biological processes, such as small-molecule biosynthesis, and are involved in the fungal stipe elongation (Huang et al., 2019).

Other underexpressed genes in Y N1 were linked to a pathogen with a lower capacity of exploiting available nutrients and N compounds from infected plant cells. This was, for example, the genes coding for “ion transporters” or for “cytochrome-c-oxidase” linked to cellular respiration, or for “4-pyridoxolactonase,” involved in the degradation of vitamin B6 whose levels respond to an interaction with nitrogen metabolism (Mooney et al., 2009; Colinas and Fitzpatrick, 2016). This decrease in the capacity of *P. brassicae* to use nitrogen-related resources would partly explain the lower growth of the pathogen in the N1 condition.

Furthermore, in an interesting way, at 27 dpi (corresponding to the secondary phase of infection), many *P. brassicae* genes potentially involved in pathogenicity were underexpressed in N1 compared to N8. The molecular function “fucosyltransferase activity” was enriched, with an underexpression at N1 in both genotypes compared to N8 for the “alpha-1;3-fucosyltransferase” gene. The fucose-containing oligosaccharides on the cell surface of some pathogenic bacteria are thought to be important for host-microbe interactions and to play a major role in the pathogenicity of bacterial pathogens (Kajiwara et al., 2012). In the human pathogen *Toxoplasma gondii*, an obligate intracellular protist, a fucosyltransferase was shown to affect protein expression and virulence (Bandini et al., 2020). Other pathogenicity-related genes were particularly underexpressed in Y N1 compared to the three other conditions. This was the case for genes encoding “ankyrin repeat,” which are protein domains frequent in all Rhizarian (Glöckner et al., 2014), and present in effector candidates of plant pathogens (Pan et al., 2008). Ankyrin domains have been identified as mediators of protein-protein interactions, a key process involved in the interplay between effector-receptor interactions (Voronin and Kiseleva, 2008). In *P. brassicae*, this domain was detected in potential

effectors and enriched in predicted proteins of its secretome (Schwelm et al., 2015a,b). A gene encoding “carboxypeptidase” was also underexpressed in Y N1. Many pathogenic fungi secrete serine carboxypeptidases, and this family of enzyme could be involved in the pathogenicity process (Monod et al., 2002; Lopez-Llorca et al., 2010). In *P. brassicae*, such proteases were also enriched in secretome (Schwelm et al., 2015b), suggesting a role in infection. In addition, genes of the “NUDIX” family were underexpressed in N1 condition for both plant genotypes, but especially when Y was the host-plant genotype. These pathogenicity-related genes, found in many phytopathogenic organisms (Dong and Wang, 2016), were previously described as overexpressed in conditions where clubroot symptoms were high (Daval et al., 2020). This suggests that this gene family may play an important role in the infectious process of *P. brassicae* and disease progression and strengthens the hypothesis that it represents relevant pathogenicity family genes. Another gene underexpressed in N1 encoded a “disulfide isomerase,” known as a molecular chaperone and component of signal-transduction pathways that is involved in fungi in the secretory pathway (Ngiam et al., 2000) and plays an important virulence role during plant infection by pathogens (Meng et al., 2015). Moreover, the “NEP1-interacting-like-2-isoform X1” was also underexpressed in N1, and even more markedly in Y N1. The “Necrosis and Ethylene-inducing Protein” (NEP)-like proteins are secreted by a wide range of plant-associated microorganisms (bacteria, fungi, and oomycetes); they are known for their cytotoxicity in dicot plants that leads to the induction of rapid tissue necrosis and plant immune responses, and they act as a microbe-associated molecular pattern (MAMP) (Oome et al., 2014). In fungi and oomycetes, NEP-proteins are especially present in species interacting with plants, and predominantly in species that display a hemibiotrophic, necrotrophic, or obligate biotrophic lifestyle on plants (Qutob et al., 2006; Lenarčič et al., 2019). Functional studies showed that NEP-proteins contribute to the virulence and pathogenicity of pathogenic microorganisms (Bae et al., 2006; Oome and Van Den Ackerveken, 2014; Levin et al., 2019). Such functional studies that characterize the gene coding for NEP and its role in infection are therefore now needed in *P. brassicae*.

On the contrary, few genes were overexpressed in Y N1. This was the case, for example, of the gene coding for a “glycosyltransferase.” In a fungal pathogen of wheat, a gene encoding a predicted “glycosyltransferase” was found to be involved in virulence by its role in extending surface hyphae (King et al., 2017). A similar role in *P. brassicae* is not known. Finally, and interestingly at 27 dpi, a gene encoding a “glutamine synthetase” was specifically overexpressed in Y N1 compared to three other conditions. “Glutamine synthetase” is an enzyme that plays an essential role in the metabolism of N by catalyzing the condensation of glutamate and ammonia to form glutamine as the principal nitrogen source for protein and nucleic acid biosynthesis. The glutamine synthetase activity, involved in the detection of nitrogen starvation (Sauer et al., 1999), is described as negatively affected by high N supply and dependent on the type of nitrogen source (Ferreira et al., 2015). In our study, the nutrient solution provided nitrogen as NO₃, while the substrate of “glutamine synthetase” is NH₃. This means that the
“nitrate reductase” and the “nitrite reductase,” whose upstream “glutamine synthetase” activities produce NH₃ from NO₃, must also be highly active in Y N1. By overexpressing in part this pathway under N stress conditions, *P. brassicae* could attempt to compensate for the nitrogen deficiency.

The time-course point 27 dpi thus showed that many genes were under-expressed in N1 compared to N8. Despite a significant reduction in symptoms in Y N1 compared to the other three conditions, the number of DEGs specifically differentially expressed between Y N1 and these other conditions was low (only 55), but with interesting putative functions. The main pathways affected by N starvation when Y was the plant genotype were involved in the growth and pathogenicity mechanisms of *P. brassicae*, which was consistent with the observed symptoms and levels of disease. Inhibition of *P. brassicae* development and pathogenicity appeared to occur at 27 dpi, although at this time, the differences in pathogen DNA amounts between Y N1 and the three other conditions were not yet marked (the symptoms were). Thus, this gene underexpression at 27 dpi in Y N1 would explain the marked differences in pathogen DNA amounts at 42 dpi. The mechanisms by which the Y genotype limited the infection by *P. brassicae* specifically in N1 will have to be elucidated by functional studies on the candidate genes found in this study such as “NUDIX,” “carboxypeptidase,” and “NEP-proteins.”

At 42 dpi, many genes involved in growth, development, signaling pathway, or stress response were still underexpressed in Y N1, concordant with the lower infection in this condition. For example, the gene encoding “sporangia induced dynein heavy” and underexpressed in N1 was shown as important in flagella formation and zoospore motility in the genus *Phytophthora* (Narayan et al., 2010). Among the genes overexpressed in N1, one encoded “L-isoaspartate (D-aspartate) O-methyltransferase,” described as crucial for maintaining cellular integrity and increasing cell survival under physiological stresses (Banerjee et al., 2015). Another gene involved in stress response was overexpressed in Y N1, the “PAS domain S-box,” which is a signaling domain widely distributed in proteins from members of prokaryotes and eukaryotes, and whose functions are to sense stimuli (Taylor and Zhulin, 1999). Then at 42 dpi, *P. brassicae* was in a phase during which it felt the stimulus of N starvation and had a weak growth in Y N1.

## CONCLUSION

The description of the molecular response of *P. brassicae* during its interaction with two *B. napus* genotypes under different nitrogen supplies shed light on candidate genes that behave as growth and pathogenicity-related, some of which were previously reported. The successful colonization and growth of the pathogen was related to its ability to draw resources from the plant, depending on the host genotype, in connection with its obligate biotrophic features. Moreover, the mechanisms of infection *via* potential secreted effectors were regulated by N supply, again depending on the host genotype. The results obtained herein need to be completed by the study of *B. napus* molecular mechanisms regulated by nitrogen, by studying more

precisely the mechanisms underlying the plant genotype-specific response. An increased understanding of the molecular responses to nitrogen is critical to reduce input costs, to minimize the potential environmental impacts of nitrogen fertilizer that has been used to optimize yield, and so to further improve agricultural sustainability.

## DATA AVAILABILITY STATEMENT

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ebi.ac.uk/ena, PRJEB44381.

## AUTHOR CONTRIBUTIONS

SD, KG, MM-D, YA, and AG were involved in the conceptualization of the project, study design, critical inputs, and finalization of the manuscript. CL and MH were involved in wet lab experiments. KG and SD were involved in bio-informatics analyses and data compilation, and wrote the manuscript. MM-D, YA, and AG advised on the manuscript. All the authors contributed to the article and approved the submitted version.

## FUNDING

This work was supported by grants from the French Association for the Promotion of Oilseed Crops Breeding (PROMOSOL).

## ACKNOWLEDGMENTS

We thank the Biological Resources Center *BrACySol* (INRAE Rennes, France) for providing the Brassica seeds.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2021.701067/full#supplementary-material

**Supplementary Figure 1** | Overview of the *Plasmodiophora brassicae* transcriptome samples.

**Supplementary Figure 2** | Validation of RNA-seq approach by comparison of RT-qPCR and RNA-seq data.

**Supplementary Table 1** | Nutrient solutions and primers used in this study.

**Supplementary Table 2** | Description of the gene subsets responsible for the driving GO-term enrichments between the N supply conditions for each time.

**Supplementary Table 3** | Description of the *Plasmodiophora brassicae* DEGs (Differentially Expressed Genes) between N supply conditions for both host plant genotypes.

**Supplementary Table 4** | Description of the *Plasmodiophora brassicae* DEGs (Differentially Expressed Genes) between Y N1 and each of the other three conditions.
REFERENCES

Aigu, Y., Laperche, A., Mendes, J., Lariagon, C., Guichard, S., Gravot, A., et al. (2018). Nitrogen supply exerts a major/minor switch between two Qtls controlling *Plasmodiophora brassicae* spore content in rapeseed. *Plant Pathol.* 67, 1574–1581. doi: 10.1111/ppa.12867

Alexa, A., and Rahnenfuhrer, J. (2007). *Gene set enrichment analysis with topGO*. Available online at: http://compdiag.molgen.mpg.de/ngfn/docs/2007/sep/topGO_Exercises.pdf (accessed August 8, 2007)

Arencibia, J. M., Pastor-Flores, D., Bauer, A. F., Schulze, J. O., and Biondi, R. M. (2013). AGC protein kinases: from structural mechanism of regulation to allosteric drug development for the treatment of human diseases. *Biochim. Biophys. Acta* 1834, 1302–1321. doi: 10.1016/j.bbpap.2013.03.010

Arkoun, M., Sarda, X., Jannin, L., Lainé, P., Etienne, P., Garcia-Mina, J.-M., et al. (2012). Hydroponics versus field lysimeter studies of urea, ammonium and nitrate uptake by oilseed rape (*Brassica napus* L.). *J. Exp. Bot.* 63, 5245–5258. doi: 10.1093/jxb/ers183

Bae, H., Kim, M. S., Sicher, R. C., Bae, H.-J., and Bailey, B. A. (2006). Necrosis- and Ethylene-inducing Peptide from *Fusarium oxysporum* induces a complex cascade of transcripts associated with signal transduction and cell death in Arabidopsis. *Plant Physiol.* 141, 1056–1067. doi: 10.1104/pp.106.07 6869

Ballini, E., Nguyen, T. T. T., and Morel, J.-B. (2013). Diversity and genetics of nitrogen-induced susceptibility to the blast fungus in rice and wheat. *Rice* 6, 1–13. doi: 10.1186/1939-8433-6-32

Bandini, G., Agop-Nersesian, C., Van Der Wel, H., Mandalasi, M., Kim, H. W., West, C. M., et al. (2020). The nucleocytosolic O-fucosyltransferase Spindly affects protein expression and virulence in *Toxoplasma gondii*. *J. Biol. Chem.* 2020:015883. doi: 10.1074/jbc.RA120.015883

Banerjee, S., Dutta, T., Lahiri, S., Sengupta, S., Gangopadhyay, A., Kumar Karri, S., et al. (2015). Enzymatic attributes of an l-isoaspartyl methyltransferase from *Candida utilis* and its role in cell survival. *Biochem. Biophys. Rep.* 4, 59–75. doi: 10.1016/j.bbp.2015.08.015

Bates, D. M., Bates, D., Maechler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. *J. Stat. Softw.* 67, 1–48. doi: 10.18637/jss.v067.i01

Bi, K., Chen, T., He, Z., Gao, Z., Zhao, Y., Fu, Y., et al. (2018). Proto-oncogenes in a eukaryotic unicellular organism play essential roles in plasmodial growth in host cells. *BMC Genomics* 19:881. doi: 10.1186/s12864-018-5307-4

Bi, K., Chen, T., He, Z., Gao, Z., Zhao, Y., Liu, H., et al. (2019). Comparative genomics reveals the unique evolutionary status of *Plasmodiophora brassicae* and the essential role of GPCR signaling pathways. *Phytopathol. Res.* 1, 18–16. doi: 10.1186/s42483-019-0018-6

Bi, K., He, Z., Gao, Z., Zhao, Y., Fu, Y., Cheng, J., et al. (2016). Integrated omics study of lipid droplets from *Plasmodiophora brassicae*. *Sci. Rep.* 6:36965. doi: 10.1038/srep36965

Bolton, M. D., and Thomma, B. P. H. J. (2008). The complexity of nitrogen metabolism and nitrogen-regulated gene expression in plant pathogenic fungi. *Physiol. Mol. Plant Pathol.* 72, 104–110. doi: 10.1016/j.pmpp.2008.07.001

Bouchet, A.-S., Laperche, A., Bissuel-Belaygue, C., Snowdon, R., Nesi, N., and Stahl, A. (2016). Nitrogen use efficiency in rapeseed. *Rev. Agron. Sustain. Dev.* 36:38. doi: 10.1007/s13593-016-0371-0

Bouchet, A.-S., Nesi, N., Bissuel, C., Bregeon, M., Lariepe, A., Navier, H., et al. (2014). Genetic control of yield and yield components in winter oilseed rape (Brassica napus L.) grown under nitrogen limitation. *Euphytica* 199, 183–205. doi: 10.1007/s10681-014-1130-4

Brennan, R. F. (1992). Effect of superphosphate and nitrogen on yield and take-all of wheat. *Fertil. Res.* 31, 43–49.

Colinas, M., and Fitzpatrick, T. B. (2016). Interaction between vitamin B6 metabolism, nitrogen metabolism and autoimmunity. *Plant Signal. Behav.* 11:e1161876. doi: 10.1080/15592324.2016.1161876

Conway, J. R., Lex, A., and Gehlenborg, N. (2017). UpSetR: an R package for the visualization of intersecting sets and their properties. *Bioinformatics* 33, 2938–2940. doi: 10.1093/bioinformatics/btx364

Cox, G. M., Harrison, T. S., Mcdade, H. C., Taborda, C. P., Heinrich, G., Casadevall, A., et al. (2003). Superoxide dismutase influences the virulence of *Cryptococcus neoformans* by affecting growth within macrophages. *Infect. Immun.* 71, 173– 180. doi: 10.1128/iai.71.1.173-180.2003

Daval, S., Belcour, A., Gazengel, K., Legrand, L., Gouzy, J., Cottret, L., et al. (2019). Computational analysis of the *Plasmodiophora brassicae* genome: mitochondrial sequence description and metabolic pathway database design. *Genomics* 111, 1629–1640. doi: 10.1016/j.ygeno.2018.11.013

Daval, S., Gazengel, K., Belcour, A., Linglin, J., Guillerm-Erckelboudt, A.-Y., Sarniguet, A., et al. (2020). Soil microbiota influences clubroot disease by modulating *Plasmodiophora brassicae* and *Brassica napus* transcriptomes. *Microb. Biotechnol.* 13, 1648–1672. doi: 10.1111/1751-7915.13634

Daval, S., Lebreton, L., Gazengel, K., Boutin, M., Guillerm-Erckelboudt, A.-Y., and Sarniguet, A. (2011). The biocontrol bacterium *Pseudomonas fluorescens* Pf29Arp strain affects the pathogenesis-related gene expression of the take-all fungus *Gaeumannomyces graminis* var. *tritici* on wheat roots. *Mol. Plant Pathol.* 12, 839–854. doi: 10.1111/j.1364-3703.2011.00715.x

Daval, S., Lebreton, L., Gracianne, C., Guillerm-Erckelboudt, A. Y., Boutin, M., Marchi, M., et al. (2013). Strain-specific variation in a soilborne phytopathogenic fungus for the expression of genes involved in pH signal transduction pathway, pathogenesis and saprophytic survival in response to environmental pH changes. *Fungal Genet. Biol.* 61, 80–89. doi: 10.1016/j.fgb. 2013.09.009

Davidson, J. A., and Krysinska-Kaczmarek, M. (2007). Effects of inoculum concentration, temperature, plant age and interrupted wetness on infection of lentil (*Lens culinaris*) by *Botrytis* spp. conidia. *Australas. Plant Pathol.* 36, 389–396.

Devadas, R., Simpfendorfer, S., Backhouse, D., and Lamb, D. W. (2014). Effect of stripe rust on the yield response of wheat to nitrogen. *Crop J.* 2, 201–206. doi: 10.1016/j.cj.2014.05.002

Dixon, G. R. (2009). *Plasmodiophora brassicae* in its environment. *J. Plant Growth Regul.* 28, 212–228. doi: 10.1007/s00344-009-9098-3

Dong, S., and Wang, Y. (2016). Nudix Effectors: A Common Weapon in the Arsenal of Plant Pathogens. *PLoS Pathog.* 12:e1005704. doi: 10.1371/journal. ppat.1005704

Donofrio, N. M., Oh, Y., Lundy, R., Pan, H., Brown, D. E., Jeong, J. S., et al. (2006). Global gene expression during nitrogen starvation in the rice blast fungus, *Magnaporthe grisea*. *Fungal Genet. Biol.* 43, 605–617. doi: 10.1016/j.fgb.2006. 03.005

Fagard, M., Launay, A., Clément, G., Courtial, J., Dellagi, A., Farjad, M., et al. (2014). Nitrogen metabolism meets phytopathology. *J. Exp. Bot.* 65, 5643–5656. doi: 10.1093/jxb/eru323

Fahling, M., Graf, H., and Siemens, J. (2003). Pathotype separation of *Plasmodiophora brassicae* by the host plant. *J. Phytopathol.* 151, 425–430. doi: 10.1046/j.1439-0434.2003.00744.x

Ferreira, E., Novais, R., Santos, F. A., Ribeiro, C., and Barros, N. (2015). Nitrate reductase (NR) and glutamine synthetase (GS) can be used as indicators of nitrogen status in eucalyptus clones. *Aus. J. Crop Sci.* 9, 561–569.

Glöckner, G., Hülsmann, N., Schleicher, M., Noegel, A. A., Eichinger, L., Gallinger, C., et al. (2014). The genome of the foraminiferan *Reticulomyxa filosa*. *Curr. Biol.* 24, 11–18. doi: 10.1016/j.cub.2013.11.027

Godichon-Baggioni, A., Maugis-Rabusseau, C., and Rau, A. (2019). Clustering transformed compositional data using K-means, with applications in gene expression and bicycle sharing system data. *J. Appl. Stat.* 46, 47–65. doi: 10.1080/ 02664763.2018.1454894

Gossen, B. D., Deora, A., Peng, G., Hwang, S.-F., and Mcdonald, M. R. (2014). Effect of environmental parameters on clubroot development and the risk of pathogen spread. *Can. J. Plant Pathol.* 36, 37–48. doi: 10.1080/07060661.2013. 859635

Hall, J. A., and Georgel, P. T. (2007). CHD proteins: a diverse family with strong ties. *Biochem. Cell Biol.* 85, 463–476. doi: 10.1139/O07-063

Huang, H., Nguyen Thi, Thu, T., He, X., Gravot, A., Bernillon, S., et al. (2017). Increase of fungal pathogenicity and role of plant glutamine in nitrogen-induced susceptibility (NIS) to rice blast. *Front. Plant Sci.* 8:00265. doi: 10.3389/ fpls.2017.00265

Huang, Q., Mukhtar, I., Zhang, Y., Wei, Z., Han, X., Huang, R., et al. (2019). Identification and characterization of two new S-adenosylmethionine-dependent methyltransferase encoding genes suggested their involvement in stipe elongation of *Flammulina velutipes*. *Mycobiology* 47, 441–448. doi: 10. 1080/12298093.2019.1658332

Kageyama, K., and Asano, T. (2009). Life cycle of *Plasmodiophora brassicae*. *J. Plant Growth Regul.* 28, 203–211. doi: 10.1007/s00344-009-9101-z
Kajiwara, H., Toda, M., Mine, T., Nakada, H., and Yamamoto, T. (2012). Isolation of fucosyltransferase-producing bacteria from marine environments. *Microbes Environ.* 27, 515–518. doi: 10.1264/jsme2.me12058

King, R., Urban, M., Lauder, R. P., Hawkins, N., Evans, M., Plummer, A., et al. (2017). A conserved fungal glycosyltransferase facilitates pathogenesis of plants by enabling hyphal growth on solid surfaces. *PLoS Pathog.* 13:e1006672-e1006672. doi: 10.1371/journal.ppat.1006672

Kombrink, A., and Thomma, B. P. (2013). LysM effectors: secreted proteins supporting fungal life. *PLoS Pathog.* 9:e1003769. doi: 10.1371/journal.ppat.1003769

Krupinsky, J. M., Halvorson, A. D., Tanaka, D. L., and Merrill, S. D. (2007). Nitrogen and tillage effects on wheat leaf spot diseases in the Northern Great plains. *Agron. J.* 99, 562–569. doi: 10.2134/agronj2006.0263

Laperche, A., Aigu, Y., Jubault, M., Ollier, M., Guichard, S., Glory, P., et al. (2017). Clubroot resistance QTL are modulated by nitrogen input in *Brassica napus*. *Theor. Appl. Genet.* 130, 669–684. doi: 10.1007/s00122-016-2842-8

Lecompte, F., Abro, M. A., and Nicot, P. C. (2010). Contrasted responses of *Botrytis cinerea* isolates developing on tomato plants grown under different nitrogen nutrition regimes. *Plant Pathol.* 59, 891–899. doi: 10.1111/j.1365-3059.2010.02320.x

Lenarčić, T., Pirc, K., Hodnik, V., Albert, I., Borišek, J., Magistrato, A., et al. (2019). Molecular basis for functional diversity among microbial Nep1-like proteins. *PLoS Pathog.* 15:e1007951. doi: 10.1371/journal.ppat.1007951

Levin, E., Raphael, G., Ma, J., Ballester, A. R., Feygenberg, O., Norelli, J., et al. (2019). Identification and functional analysis of NLP-encoding genes from the postharvest pathogen *Penicillium expansum*. *Microorganisms* 7:7060175. doi: 10.3390/microorganisms7060175

Lopez-Llorca, L. V., Gómez-Vidal, S., Monfort, E., Larriba, E., Casado-Vela, J., Elortza, F., et al. (2010). Expression of serine proteases in egg-parasitic nematophagous fungi during barley root colonization. *Fungal Genet. Biol.* 47, 342–351. doi: 10.1016/j.fgb.2010.01.004

Ludwig-Muller, J., Julke, S., Geiss, K., Richter, F., Mithofer, A., Sola, I., et al. (2015). A novel methyltransferase from the intracellular pathogen *Plasmodiophora brassicae* methylates salicylic acid. *Mol. Plant Pathol.* 16, 349–364. doi: 10.1111/mpp.12185

Manna, S., Brewster, J., and Barth, C. (2013). Identification of Pentatricopeptide Repeat Proteins in the model organism *Dictyostelium discoideum*. *Int. J. Genomics* 2013:586498. doi: 10.1155/2013/586498

Manzanares-Dauleux, M. J., Divaret, I., Baron, F., and Thomas, G. (2000). Evaluation of French *Brassica oleracea* landraces for resistance to *Plasmodiophora brassicae*. *Euphytica* 113, 211–218. doi: 10.1023/a:1003997421340

Meng, Y., Zhang, Q., Zhang, M., Gu, B., Huang, G., Wang, Q., et al. (2015). The protein disulfide isomerase 1 of *Phytophthora parasitica* (PpPDI1) is associated with the haustoria-like structures and contributes to plant infection. *Front. Plant Sci.* 6:00632. doi: 10.3389/fpls.2015.00632

Monod, M., Capoccia, S., Léchenne, B., Zaugg, C., Holdom, M., and Jousson, O. (2002). Secreted proteases from pathogenic fungi. *Int. J. Med. Microbiol.* 292, 405–419. doi: 10.1078/1438-4221-00223

Mooney, S., Leuendorf, J. E., Hendrickson, C., and Hellmann, H. (2009). Vitamin B6: a long known compound of surprising complexity. *Molecules* 14, 329–351. doi: 10.3390/molecules14010329

Mur, L. A. J., Simpson, C., Kumari, A., Gupta, A. K., and Gupta, K. J. (2016). Moving nitrogen to the centre of plant defence against pathogens. *Ann. Bot.* 119, 703–709. doi: 10.1093/aob/mcw179

Narayan, R. D., Blackman, L. M., Shan, W., and Hardham, A. R. (2010). *Phytophthora nicotianae* transformants lacking dynein light chain 1 produce non-flagellate zoospores. *Fungal Genet. Biol.* 47, 663–671. doi: 10.1016/j.fgb.2010.04.008

Ngadze, E., Icishahayo, D., Coutinho, T. A., and Van Der Waals, J. E. (2012). Role of polyphenol oxidase, peroxidase, phenylalanine ammonia lyase, chlorogenic acid, and total soluble phenols in resistance of potatoes to soft rot. *Plant Dis.* 96, 186–192. doi: 10.1094/pdis-02-11-0149

Ngiam, C., Jeenes, D. J., Punt, P. J., Van Den Hondel, C. A. M. J. J., and Archer, D. B. (2000). Characterization of a foldase, protein disulfide isomerase A, in the protein secretory pathway of *Aspergillus niger*. *Appl. Environ. Microbiol.* 66, 775–782. doi: 10.1128/aem.66.2.775-782.2000

Nunes-Nesi, A., Fernie, A. R., and Stitt, M. (2010). Metabolic and signaling aspects underpinning the regulation of plant carbon nitrogen interactions. *Mol. Plant* 3, 973–996. doi: 10.1093/mp/qq049

Oome, S., and Van Den Ackerveken, G. (2014). Comparative and functional analysis of the widely occurring family of Nep1-Like Proteins. *Mol. Plant-Microbe Interact.* 27, 1081–1094. doi: 10.1094/mpmi-04-14-0118-r

Oome, S., Raaymakers, T. M., Cabral, A., Samwel, S., Böhm, H., Albert, I., et al. (2014). Nep1-like proteins from three kingdoms of life act as a microbe-associated molecular pattern in *Arabidopsis*. *Proc. Natl. Acad. Sci. U S A.* 111, 16955–16960. doi: 10.1073/pnas.1410031111

Pan, X., Lührmann, A., Satoh, A., Laskowski-Arce, M. A., and Roy, C. R. (2008). Ankyrin repeat proteins comprise a diverse family of bacterial type IV effectors. *Science* 320, 1651–1654. doi: 10.1126/science.1158160

Perez-Lopez, E., Waldner, M., Hossain, M., Kusalik, A. J., Wei, Y., Bonham-Smith, P. C., et al. (2018). Identification of *Plasmodiophora brassicae* effectors - A challenging goal. *Virulence* 9, 1344–1353. doi: 10.1080/21505594.2018.1504560

Pitchay, D., Frantz, J., Locke, J. C., and Fernandez, G. (2007). Impact of applied nitrogen concentration on growth of Elatior Begonia and New Guinea Impatiens and susceptibility of Begonia to *Botrytis cinerea*. *J. Am. Soc. Hort. Sci.* 132:193. doi: 10.21273/JASHS.132.2.193

Qutob, D., Kemmerling, B., Brunner, F., Küfner, I., Engelhardt, S., Gust, A. A., et al. (2006). Phytotoxicity and innate immune responses induced by Nep1-like proteins. *Plant Cell* 18, 3721–3744. doi: 10.1105/tpc.106.044180

Rathke, G. W., Christen, O., and Diepenbrock, W. (2005). Effects of nitrogen source and rate on productivity and quality of winter oilseed rape (*Brassica napus* L.) grown in different crop rotations. *Field Crops Res.* 94, 103–113. doi: 10.1016/j.fcr.2004.11.010

Rau, A., and Maugis-Rabusseau, C. (2018). Transformation and model choice for co-expression analysis of RNA-seq data. *Brief. Bioinform.* 19, 425–436.

R'houma, A., Chérif, M., and Boubaker, A. (1998). Effect of nitrogen fertilization, green pruning and fungicide treatments on Botrytis bunch rot of grapes. *J. Plant Pathol.* 80, 115–124.

Rolfe, S. A., Strelkov, S. E., Links, M. G., Clarke, W. E., Robinson, S. J., Djavaheri, M., et al. (2016). The compact genome of the plant pathogen *Plasmodiophora brassicae* is adapted to intracellular interactions with host *Brassica spp*. *BMC Genomics* 17:272. doi: 10.1186/s12864-016-2597-2

Sabater-Munoz, B., Legeai, F., Rispe, C., Bonhomme, J., Dearden, P., Dossat, C., et al. (2006). Large-scale gene discovery in the pea aphid *Acyrthosiphon pisum* (Hemiptera). *Genome Biol.* 7:11. doi: 10.1186/gb-2006-7-3-r21

Sauer, J., Görl, M., and Forchhammer, K. (1999). Nitrogen starvation in *Synechococcus* PCC 7942: involvement of glutamine synthetase and NtcA in phycobiliprotein degradation and survival. *Arch. Microbiol.* 172, 247–255. doi: 10.1007/s002030050767

Scheible, W. R., Morcuende, R., Czechowski, T., Fritz, C., Osuna, D., Palacios-Rojas, N., et al. (2004). Genome-wide reprogramming of primary and secondary metabolism, protein synthesis, cellular growth processes, and the regulatory infrastructure of Arabidopsis in response to nitrogen. *Plant Physiol.* 136, 2483–2499. doi: 10.1104/pp.104.047019

Schwelm, A., Dixelius, C., and Ludwig-Müller, J. (2015a). New kid on the block – the clubroot pathogen genome moves the plasmodiophorids into the genomic era. *Eur. J. Plant Pathol.* 145, 531–542. doi: 10.1007/s10658-015-0839-9

Schwelm, A., Fogelqvist, J., Knaust, A., Julke, S., Lilja, T., Bonilla-Rosso, G., et al. (2015b). The *Plasmodiophora brassicae* genome reveals insights in its life cycle and ancestry of chitin synthases. *Sci. Rep.* 5:11153. doi: 10.1038/srep11153

Searle, S. R., Speed, F. M., and Milliken, G. A. (1980). Population Marginal Means in the linear model: an alternative to Least Squares Means. *Am. Stat.* 34, 216–221. doi: 10.1080/00031305.1980.10483031

Siemens, J., Keller, I., Sarx, J., Kunz, S., Schuller, A., Nagel, W., et al. (2006). Transcriptome analysis of *Arabidopsis* clubrooms Indicate a key role for cytokinins in disease development. *Mol. Plant-Microbe Interact.* 19, 480–494. doi: 10.1094/mpmi-19-0480

Solomon, P. S., Tan, K. C., and Oliver, R. P. (2003). The nutrient supply of pathogenic fungi; a fertile field for study. *Mol. Plant Pathol.* 4, 203–210. doi: 10.1046/j.1364-3703.2003.00161.x

Some, A., Manzanares, M. J., Laurens, F., Baron, F., Thomas, G., and Rouxel, F. (1996). Variation for virulence on *Brassica napus* L amongst *Plasmodiophora brassicae* collections from France and derived single-spore isolates. *Plant Pathol.* 45, 432–439. doi: 10.1046/j.1365-3059.1996.d01-155.x
Soulie, M.-C., Koka, S. M., Floch, K., Vancostenoble, B., Barbe, D., Daviere, A., et al. (2020). Plant nitrogen supply affects the Botrytis cinerea infection process and modulates known and novel virulence factors. *Mol. Plant Pathol.* 21, 1436–1450. doi: 10.1111/mpp.12984

Sun, Y., Wang, M., Mur, L. A. J., Shen, Q., and Guo, S. (2020). Unravelling the roles of nitrogen nutrition in plant disease defences. *Int. J. Mol. Sci.* 21:ijms21020572. doi: 10.3390/ijms21020572

Takahashi, T. (2020). Plant Polyamines. *Plants* 9:511. doi: 10.3390/plants9040511

Takahashi, Y. (2016). The role of Polyamines in plant disease resistance. *Environ. Control Biol.* 54, 17–21. doi: 10.2525/ecb.54.17

Tavernier, V., Cadiou, S., Pageau, K., Laugé, R., Reisdorf-Cren, M., Langin, T., et al. (2007). The plant nitrogen mobilization promoted by *Colletotrichum lindemuthianum* in *Phaseolus* leaves depends on fungus pathogenicity. *J. Exp. Bot.* 58, 3351–3360. doi: 10.1093/jxb/erm182

Taylor, B. L., and Zhulin, I. B. (1999). PAS domains: internal sensors of oxygen, redox potential, and light. *Microbiol. Mol. Biol. Rev.* 63, 479–506.

Voronin, D. A., and Kiseleva, E. V. (2008). Functional role of proteins containing ankyrin repeats. *Cell Tiss. Biol.* 2, 1–12. doi: 10.1134/S1990519X0801001X

Wang, M., Gu, Z., Wang, R., Guo, J., Ling, N., Firbank, L. G., et al. (2018). Plant primary metabolism regulated by nitrogen contributes to plant-pathogen interactions. *Plant Cell Physiol.* 60, 329–342. doi: 10.1093/pcp/pcy211

Xin, C., Qing-Wei, Y., Jia-Lin, S., Shuang, X., Fu-Chun, X., and Ya-Jun, C. (2014). Research progress on nitrogen use and plant growth. *J. Northeast Agric. Univ.* 21, 68–74. doi: 10.1016/S1006-8104(14)60036-2

Yermiyahu, U., Shamai, I., Peleg, R., Dudai, N., and Shtienberg, D. (2006). Reduction of Botrytis cinerea sporulation in sweet basil by altering the concentrations of nitrogen and calcium in the irrigation solution. *Plant Pathol.* 55, 544–552. doi: 10.1111/j.1365-3059.2006.01388.x

Zarattini, M., Farjad, M., Launay, A., Cannella, D., Soulié, M. C., Bernacchia, G., et al. (2021). Every cloud has a silver lining: how abiotic stresses affect

gene expression in plant-pathogen interactions. *J. Exp. Bot.* 72, 1020–1033. doi: 10.1093/jxb/eraa531

Zhang, C., Dong, Y., Tang, L., Zheng, Y., Makowski, D., Yu, Y., et al. (2019). Intercropping cereals with faba bean reduces plant disease incidence regardless of fertilizer input; a meta-analysis. *Eur. J. Plant Pathol.* 154, 931–942. doi: 10.1007/s10658-019-01711-4

Zhang, D., Burroughs, A. M., Vidal, N. D., Iyer, L. M., and Aravind, L. (2016). Transposons to toxins: the provenance, architecture and diversification of a widespread class of eukaryotic effectors. *Nucleic Acids Res.* 44, 3513–3533. doi: 10.1093/nar/gkw221

Zhang, W., Wu, L., Ding, Y., Yao, X., Wu, X., Weng, F., et al. (2017). Nitrogen fertilizer application affects lodging resistance by altering secondary cell wall synthesis in japonica rice (Oryza sativa). *J. Plant Res.* 130, 859–871. doi: 10.1007/s10265-017-0943-3

Zhong, C., Cao, X., Hu, J., Zhu, L., Zhang, J., Huang, J., et al. (2017). Nitrogen metabolism in adaptation of photosynthesis to water stress in rice grown under different nitrogen levels. *Front. Plant Sci.* 8:01079. doi: 10.3389/fpls.2017.01079

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2021 Gazengel, Aigu, Lariagon, Humeau, Gravot, Manzanares-Dauleux and Daval. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.