Transcriptomic response in symptomless roots of clubroot infected kohlrabi (Brassica oleracea var. gongylodes) mirrors resistant plants

Background Clubroot disease caused by Plasmodiophora brassicae (Phytomyxea, Rhizaria) is one of the economically most important diseases of Brassica crops. The formation of hypertrophied roots accompanied by altered metabolism and hormone homeostasis is typical for infected plants. Not all roots of infected plants show the same phenotypic changes. While some roots remain uninfected, others develop galls of diverse size. The aim of this study was to analyse and compare the intra-plant heterogeneity of P. brassicae root galls and symptomless roots of the same host plants (Brassica oleracea var. gongylodes) collected from a commercial field in Austria using transcriptome analyses. Results Transcriptomes were markedly different between symptomless roots and gall tissue. Symptomless roots showed transcriptomic traits previously described for resistant plants. Genes involved in host cell wall synthesis and reinforcement were up-regulated in symptomless roots indicating elevated tolerance against P. brassicae. By contrast, genes involved in cell wall degradation and modification processes like expansion were up-regulated in root galls. Hormone metabolism differed between symptomless roots and galls. Brassinosteroid-synthesis was down-regulated in root galls, whereas jasmonic acid synthesis was down-regulated in symptomless roots. Cytokinin metabolism and signalling were up-regulated in symptomless roots with the exception of one CKX6 homolog, which was strongly down-regulated. Salicylic acid (SA) mediated defence response was up-regulated in symptomless roots, compared with root gall tissue. This is probably caused by a secreted benzoic acid/salicylic acid methyl transferase from the pathogen (PbBSMT), which was one of the highest expressed pathogen genes in gall tissue. The PbBSMT derived Methyl-SA potentially leads to increased pathogen tolerance in uninfected roots. Conclusions Infected and uninfected roots of clubroot infected plants showed transcriptomic differences similar to those previously described between clubroot resistant and susceptible hosts. The here described intra-plant heterogeneity suggests, that for a better understanding of clubroot disease targeted, spatial analyses of clubroot infected plants will be vital in understanding this economically important disease. Electronic supplementary material The online version of this article (10.1186/s12870-019-1902-z) contains supplementary material, which is available to authorized users.


Background
Clubroot disease is one of the most important diseases of Brassica crops worldwide accounting for approximately 10% loss in Brassica vegetable, fodder, and oilseed crops [1]. Clubroot is caused by Plasmodiophora brassicae, an obligate biotrophic protist, taxonomically belonging to Phytomyxea within the eukaryotic supergroup Rhizaria [2,3]. This soil borne pathogen has a complex life cycle. Zoospores infect root hairs where primary plasmodia form. These plasmodia develop into secondary zoospores, which are released into the soil and re-infect the root cortex where secondary plasmodia develop [4]. The secondary plasmodia mature into resting spores, which are released into the soil. In infected host tissue division and elongation of cells is triggered upon infection, which leads to hypertrophies of infected roots resulting in the typical root galls or clubroots (Fig. 1).
Plasmodiophora brassicae can only be grown and studied in co-culture with its host. This has hampered both, targeted and large scale studies on the molecular basis of P. brassicae and the interactions with its host [5]. Because of the economic importance of clubroot disease, numerous studies analysed specific aspects of the biology, physiology and molecular biology of the plantpathogen interaction to better understand and control the disease. The first of these experimental studies were based on the Arabidopsis/Plasmodiophora pathosystem e.g. [6,7]. During the last years an increasing number of Brassica (host) genomes became available [8][9][10][11], as did several P. brassicae genomes [12][13][14][15] permitting new research approaches, including targeted transcriptome studies. There are analyses of (plant) transcriptomes of roots of clubroot infected plants compared with uninfected plants [14,[16][17][18], of host varieties susceptible and tolerant to clubroot [14,19], or the host response to different P. brassicae isolates [6,16,20,21].
Plants infected with P. brassicae show marked physiological changes including cell wall biosynthesis, plant hormone metabolism and plant defence related processes. Expansin genes, involved in plant cell expansion and elongation [22], were up-regulated in P. brassicae infected roots [6,16,23]. In P. brassicae infected roots enzymatic activity of xyloglucan endotransglycosylases/ hydrolases (XTHs) increases [24], while an early response to P. brassicae infection was up-regulation of the phenylpropanoid pathway that provides lignin precursors [18]. With progression of clubroot development, the lignification of clubroot tissue was reduced [25] and genes involved in lignification processes were downregulated [26]. On the other hand cell wall thickening and lignification was suggested to limit the spread of the pathogen in tolerant B. oleracea [27] and B. rapa [28].
The development of clubroot symptoms is accompanied by changes of plant hormone homeostasis [6,17,29]. During clubroot development, auxins mediate host cell divisions and elongation. Auxins increase over time during clubroot development and accumulate in P. brassicae infected tissues in a sink like manner [29]. In addition, genes belonging to the auxin conjugating GH3 protein family are regulated differentially during clubroot development [30] with one GH3 protein gene (PbGH3;  .1) identified in the P. brassicae genome [12]. Cytokinins (CKs) increase initially, but decrease again with the onset of gall formation [17]. At the same time P. brassicae plasmodia produce minute amounts of CKs [31]. Therefore, CKs play a crucial role in disease development not only through their regulation of cell division but also through their interference in the sugar metabolism and invertase production, which might be crucial for the nutrition of P. brassicae [32,33].
Stress-and defence related phytohormones like salicylic acid (SA), jasmonic acid (JA), brassinosteroids (BR), and ethylene (ET) and their regulatory pathways also change in response to pathogen infection [34]. The accumulation of SA plays a key role in plant defence against biotrophic pathogens, often resulting in a localized hypersensitive response and induction of pathogenesis-related (PR) genes. Systemic acquired resistance (SAR) is a form of induced resistance that is activated by SA throughout a plant after being exposed to elicitors from microbes or chemical stimuli [35]. High endogenous levels of SA and exogenous SA reduced the infection of the host by P. brassicae [20,36]. In tolerant hosts SA related genes are induced upon infection [16,21,37]. The SAR-deficient npr1-1 and SA-deficient isochorismate synthase 1 (ICS1) sid2 Arabidopsis mutants showed an increased susceptibility to P. brassicae, whereas the bik-1 mutant, with elevated SA levels, was more resistant [38]. Pathogenesis-related defence proteins were induced by SAR and expressed more highly in resistant than in susceptible B. rapa and Arabidopsis species [14,39,40]. P. brassicae might counteract the plant SA-mediated defence via a secreted methyltransferase (PbBSMT; AFK13134.1). This SABATH-like methyltransferase has been shown to convert SA to methyl-salicylate (MeSA) in vitro [41]. The proposed function in planta is the removal of SA in local infected tissue as MeSA is volatile. Arabidopsis mutants expressing the PbBMST gene showed a higher susceptibility towards P. brassicae [42].
The P. brassicae PbGH3 was also able to conjugate JA with amino acids in vitro [12]. In general, JA is associated with resistance against necrotrophic microbes [43,44]. In A. thaliana Col-0 several JA-responsive genes were induced in infected root tissues and JA accumulates in galls [6,45]. Jasmonate resistant 1 (jar1) mutants, impaired in JA-Ile accumulation, showed a higher susceptibility to P. brassicae [16,45]. Thus, JA responses contributed to a basal resistance against some strains of P. brassicae in A. thaliana Col-0 [45]. But in partially resistant Arabidopsis Bur-0 only weak JA responses compared with the susceptible Col-0 were found [46]. Generally, clubroot susceptible hosts show a high level of JA response, whereas it is reduced in resistant hosts [14,19]. Those differences might be due to occur if aliphatic or aromatic glucosinolate production is induced by JA in the particular host [47].
The aim of this study was to generate the first data set of root tissue specific transcriptomic response of individual plants during clubroot development. Usually, clubroot infected plants do not develop symptoms uniformly on all root parts with some roots showing strong symptoms and others not showing symptoms at all (Fig. 1). We collected samples of kohlrabi (Brassica oleracea var. gongylodes) infected with P. brassicae (PbraAT) from a field in Austria. We compared morphologically different clubroots and symptomless roots of the same infected plants. We analysed similarities and differences of their transcriptomic profile focussing on cell wall metabolism, hormone metabolism, and defence response.

Transcriptome analyses
In total 143 million good quality reads with an average length of 125 bp were obtained from all libraries (Additional file 2: Table S1). A total of 10,940 genes including isoforms were predicted for P. brassicae and 42,712 for B. oleracea. About 50% of the P. brassicae and 85% of the kohlrabi transcripts could be functionally annotated using eggNOG-mapper. Only 0.0005% of the reads of the symptomless root samples (SL) matched P. brassicae transcripts, which indicates that these roots were not infected by P. brassicae. In the white spindle galls (WG) and brownish spindle galls (BG) libraries 23 and 33% of the reads matched P. brassicae, respectively (Additional file 1: Figure S1).
The transcriptomes (SL, WG, BG) contained a total of 5204 differentially expressed genes (DEGs). Compared with SL, in WG 1619 DEGs were up-regulated and 2280 were down-regulated (Fig. 2, Additional file 2: Table S2), while in BG 942 DEGs were up-and 2571 were down-regulated. Of all DEGs, 790 were assigned to the COG (Cluster of Orthologous Groups) category "Information and Storage Processing", 1401 to "Metabolism", 1245 to "Cellular Process and Signalling" and 1768 to "Poorly Characterized" by eggNOG-mapper (Fig. 2, Additional file 2: Table S3). Only 19 plant genes were differentially expressed between BG and WG (Additional file 2: Tables S4).

Plant cell wall metabolism
In B. oleracea 161 of the 5204 DEGs within infected plants (SL vs WG vs BG) were involved in cell wall synthesis, modification, degradation, or phenylpropanoid metabolism. Cellulose, hemicellulose, pectin, and lignin synthesis genes involved in cell wall reinforcement were up-regulated in SL compared with gall tissues (Fig. 3), whereas genes involved in cell wall modification and degradation were down-regulated (Fig. 3). The changes in expression of cell wall genes were more prominent between SL and BG than between SL and WG. In SL a UDP-D-glucuronate 4-epimerase 6 (GAE6) homolog responsible for the synthesis of UDP-D-glucuronic acid, the main building block for pectins [48] was upregulated compared with gall tissue. Predicted expansin (EXP) and expansin-like (EXL) genes were mainly downregulated in SL compared with WG and BG (Fig. 3). Genes coding for XTHs were among the strongest down-regulated DEGs in SL, with XTH24 being the strongest down-regulated transcript of all DEGs. The phenylpropanoid pathway was up-regulated in SL compared with WG and BG (Fig. 3). This includes the phenylalanine ammonia-lyase 1 (PAL1) homolog, a key enzyme in lignin biosynthesis. Lignin is also part of the xylem and xylogenesis genes were up-regulated in SL compared with BG. Flavonoid metabolism was also induced in SL (Additional file 1: Figure S2).

Plant hormones
The CK and auxin metabolism was altered in WG and BG compared with SL (Fig. 4). So was a homolog of CKX6 (cytokinin oxidase/dehydrogenase 6) downregulated in SL compared with WG. Homologs of CKX5, CK receptors, and a CK-regulated UDPglucosyltransferase were up-regulated in SL. The CK synthesis genes CYP735A2 (cytochrome P450) and LOG1 (lonely guy 1) were up-regulated in SL compared with root galls (Fig. 4). Down-stream in CK-signalling we found DEGs within ARR (Arabidopsis response regulator; CK-signalling target) genes (Fig. 4). Type-B ARRs were differentially expressed in WG and BG whereas differentially expressed Type-A ARRs were exclusively found in BG. A putative AHK4 (Arabidopsis histidine kinase 4; CK-receptor) gene was up-regulated in SL compared with BG where it was not detected although expression values in SL were low. No DEGs of AHPs (Arabidopsis histidine phosphotransfer protein) were found in SL compared with root galls. Most auxin related DEGs, auxin response factors (ARFs and IAAs) and IAA amino acid conjugate synthetases (GH3) were up-regulated in SL compared with gall tissue (Additional file 1: Figure S3). However, an IAA7 (repressor of auxin inducible genes) and a homolog of GH3.2 were downregulated. Expression of PIN-FORMED 1 (PIN1) genes was reduced in SL ( Fig. 4), whereas SAUR (small auxin up-regulated RNA) and AIR12 (auxin-induced in root cultures protein 12-like) genes, were found in up-and down-regulated DEGs (Fig. 4). Myrosinases and nitrilases were up-regulated in SL compared with galls (Additional file 1: Figure S4). In BG two transcripts related to auxin synthesis and regulation were down-regulated compared with WG (Additional file 2: Table S4).

Plant defence
Generally, genes for disease resistance proteins were upregulated in SL compared with WG and BG (Fig. 5, Additional file 1: Figure S6). From pathogen recognition genes via signalling proteins and transcription factors to pathogenesis related (PR) proteins, the whole signal cascade of pathogen defence was affected. Of the predicted defence related DEGs within infected plants (SL vs WG vs BG) 60 were assigned as TIR-NBS-LRR (Toll/ interleukin-1 receptor nucleotide-binding site leucinerich repeat) class proteins. Within R-genes only PP2-A6 and PP2-A8 (phloem protein 2-A) could be functionally annotated and were up-regulated in SL compared with root galls (Fig. 5). Expressed genes for disease resistance protein RSP4 were not altered between the tissue types (Additional file 5). Signalling genes encoding for MLO (mildew resistance locus o) and MKS1 (MAP kinase substrate 1) were upregulated in SL (Fig. 5), while no differences in the expression of CHORD or SRFR1 (suppressor of RSP4-RDL1) genes was observed (Additional file 5). Two transcription factors involved in the biotic stress The BIK1 (botrytis-induced kinase 1), which interacts with BRI1 and BAK1 (BRI1-associated receptor kinase) to induce defence responses was up-regulated in SL compared with galls (Additional file 1: Figure S7).
In the SL samples JA-related genes such as LOX2 (lipoxygenase 2), AOC (allene oxide cyclase), and HPL (hyperoxide lyase) were down-regulated, while other LOX genes and the JA amido synthetase genes JAR1 were up-regulated in SL compared with galls (Fig. 4).
One down-regulated isoform of JAR1 was found between BG and WG. We found no glucosinolate biosynthesis genes in the DEGs in our study.
Genes for SA modification, like the SA methylating SABATH methyl transferase genes (PXMT1, GAMT1) and a SA-glucosidase (UGT74F2) were up-regulated in SL compared with galls (Fig. 4). The SA induced PR1 gene was induced in SL (Fig. 5, Additional file 1: Figure  S6). The PR-gene expression regulator NPR1 (nonexpressor of PR1) was not differentially expressed in our samples. Of genes that regulate PR1 expression via NPR1, WRKY70 was down-regulated in SL whereas NPR3 and TGA3 were up-regulated (Fig. 5, Additional file 1: Figure S5). Genes for the TAO1 disease resistant protein, which induces PR1 expression were up- Fig. 4 Kohlrabi phytohormone metabolism. Heatmaps of log 2 fold change values of DEGs. Genes involved in cytokinin, jasmonic acid, salicylic acid, and abscisic acid metabolism were up-regulated in SL. Genes coding for brassinosteroids clustered into two groups: later stages in BR biosynthesis (down-regulated) and early sterol biosynthesis (up-regulated). Genes involved in auxin metabolism were found within the up-and down-regulated DEGs in SL compared with root galls. One DEGs (JAR1) was found between WG and BG. Up-regulated genes are shaded in purple and down-regulated genes in green. Arabidopsis homologs are given. Genes are annotated and categorized according to MapMan/ Mercator, except for cytokinin metabolism for which genes were categorized as described by [17]. NA: not assigned by MapMan. Additional file information for the genes including KEGG, eggNog, and TAIR annotation are given in Additional file 4 regulated in SL, as well as the NDR1 (non-race specific disease resistance 1) gene, required for the establishment of hypersensitive response and SAR. Genes for the disease resistant protein RPS2, activated by NDR1, were also up-regulated in SL (Fig. 5, Additional file 1: Figure  S6). Additionally, other defence related genes coding for protease inhibitor genes, R-genes or some chitinases were down-regulated in SL (Fig. 5).

P. brassicae gene expression
The P. brassicae genes with the highest FPKM values belonged to growth and cellular process related COG categories such as translation, transcription, and signal transduction, but also to energy conversion and carbohydrate and lipid metabolism (Additional file 1: Figure S8). The P. brassicae PbBSMT gene was amongst the highest expressed pathogen genes (Additional file 2: Tables S5, S6). Other highly expressed genes were HSPs (heat shock proteins), a glutathione-S-transferase, an ankyrin repeat domain-containing protein, ribosomal genes, and genes of unknown function (Additional file 2: Tables S5,  S6). The PbGH3 gene was not expressed in our samples. The P. brassicae protease gene PRO1, proposed to be involved in spore germination [49], was expressed in both WG and BG.
Between WG and BG samples only five P. brassicae DEGs were identified, coding for a HSP, a chromosomal maintenance protein, a DNA-directed RNA-polymerase, a retrotransposon and a Scl Tal1 interrupting locus protein (Additional file 2: Table S7).
Cumulating all FPKM values revealed that most sequenced reads from P. brassicae RNA extracted from root galls (WG and BG) mapped to the COG categories "Post-translational modification, protein turnover, chaperon functions" and "Translation" (Additional file 1: Figure S9). Very few P. brassicae reads were found in the data obtained from SL (Additional file 1: Figures S1, S10), those were most likely from attached spores or contamination via soil particles.

Symptomless roots of clubroot infected plants show transcriptomic traits of clubroot resistant/tolerant plants
We found that symptomless roots and clubroots originating from the same plant showed differences in their transcriptomic profile similar to previously described differences of roots between resistant and susceptible plants [40,50]. Gene expression patterns of symptomless roots were similar to the patterns described for resistant hosts, while in clubroot tissue patterns were similar to those observed in susceptible plants.
Reinforcement of cell walls has previously been reported to hamper the development of P. brassicae in resistant B. oleracea [27] and B. rapa callus cultures [28]. Lignin biosynthesis genes were up-regulated in SL tissue compared with root gall tissue (Fig. 3). Induced lignification processes were observed in shoots of infected plants [23] and between resistant and susceptible B. oleracea cultivars [21]. PAL1, a key enzyme in lignin, SA (discussed below), and flavonoid biosynthesis [51][52][53] was up-regulated in symptomless roots compared with clubroot tissue (Fig. 3). Increased lignin biosynthesis and upregulation of PAL1 has been described for a clubroot resistant oilseed B. rapa line carrying the resistance gene Rcr1 [53], while callus cultures overexpressing PAL1 were resistant to infection by P. brassicae [28]. Root reinforcement, therefore, seems to be a part of the tolerance mechanisms of plants against P. brassicae, which has only a limited arsenal of plant cell wall degrading enzymes in its genome [12,13] and infects its hosts via mechanical force with a specialised extrusosome called "Stachel and Rohr" [4]. Once P. brassicae successfully infected its host, movement and spread within root tissues has been discussed to happen via the plasmodesmata [27,54,55]. Increased stability of the cell wall, as indicated by gene expression patterns in symptomless roots (Fig. 3) or described for resistant plants [27,28], requires higher mechanical force for successful (primary) infection of yet uninfected roots or movement between host cells. Therefore, it is likely that cell wall reinforcement of the host cells is a considerable obstacle for P. brassicae. The result of increased lignification in SL compared with galls could also be the result of repressed xylem development in P. brassicae infected hosts [56].
In symptomless roots of infected plants, defence related pathways regulated by hormones showed patterns usually linked to induced plant defence (Fig. 4). Clubroot tissue on the other hand showed suppression of SA-defence related processes. Genes of SA biosynthesis were up-regulated in SL, but were down-regulated in galls. Salicylic acid can be synthesised via isochorismate (ICS pathway) or from phenylalanine (via PAL pathway) [36,43]. The up-regulation of the PAL1 gene in symptomless roots could therefore also be linked to the PALdependent synthesis of SA, but the majority of SA in clubroot is produced via the ICS pathway [36]. Based on the higher expression of ICS1 (not significant) and WRKY28 (Additional file 1: Figure S5) the synthesis of SA was likely induced in SL. In clubroot resistant hosts SA-defence is usually induced [14,36,39], whereas the SA deficient sid2 (ICS1) mutant of Arabidopsis was more susceptible to P. brassicae infection [38].
High SA levels in plant tissues helped to reduce new P. brassicae infections [20], but SA alone is not sufficient to induce resistance against P. brassicae [36]. Because SA levels increase in clubroot tissue, P. brassicae is thought to secret a SABATH-type methyltransferase [PbBSMT; 42,43], which was one of the highest expressed genes of P. brassicae in this study (Additional file 2: Tables S5, S6). The PbBSMT has been shown to methylate salicylic acid, contributing to a local reduction of SA in the galls [41]. MeSA is the major transport form of SA in the plant and has a key role in inducing SAR [57,58]. So based on our data we hypothesise that P. brassicae reduces SA concentrations in the galls via PbBMST mediated methylation [41]. The produced MeSA could trigger SA-related defences in distant plant parts, which become resilient towards new pathogen infection.
Processes downstream from SA are mediated via NPR1 (not differentially expressed in our dataset). NPR1 induces PR-gene expression when bound to TGA3 [59]. In its interaction with WRKY70, NPR1 serves as a negative regulator of the SA biosynthesis gene ICS1 [60]. Thus the observed up-regulation of TGA3 in SL (Additional file 1: Figure S5) would lead to an induction of expression of PR-genes in SL. The reduced expression of WRKY70 in SL compared with the galls would lead to a higher SA production in the symptomless root tissue. Additionally, we found an up-regulation of NPR3 genes, coding for repressors of SA-defence genes [61], in SL compared with the galls (Fig. 5). When bound to SA, NPR3 would lose its function of repressing SA-defence genes [61]. The higher expression of NPR3 in the SL might be necessary to compensate for negative effects of the SA synthesis in symptomless roots.
Jasmonic acid contributed to a basal resistance against some strains of P. brassicae. Arabidopsis thaliana Col-0 and A. thaliana mutants impaired in JA-Ile accumulation showed a higher susceptibility to P. brassicae [16,45]. In A. thaliana Col-0 several JA-responsive genes were induced in infected root tissues and JA accumulates in galls [6,45]. But in partially resistant Arabidopsis Bur-0 only weak JA responses compared with the susceptible Col-0 were found [46]. Those differences might be due to if aliphatic or aromatic glucosinolate production is induced by JA in the particular host [47]. Generally, clubroot susceptible hosts show a high level of JA response, whereas it is reduced in resistant hosts [14,19]. In our samples, JAZ genes were upregulated in SL compared with the gall tissue (Fig. 4). JAZ proteins act as JA co-receptors and transcriptional repressors in JA signalling [62] reducing the JA synthesis in SL. This was previously observed in B. oleracea plants, where JAZ expression was up-regulated in resistant plants and JA synthesis was highly induced in susceptible plants [21]. Among JA-metabolism genes HPL, and LOX2 expression were up-regulated in the galls. LOX2 is essential for the formation of oxylipin volatiles [63]. The HPL protein competes with substrates essential for JA-synthesis, producing volatile and non-volatile oxylipins [64]. The higher expression of HPL and LOX2 in the galls might lead to the production of volatile aldehydes rather than a JA accumulation in the galls.
In SL BR-synthesis and BR-signalling genes were down-regulated (Fig. 4). BRs are necessary for the development of clubroot tissue [65], hence a reduction in BRs impairs P. brassicae growth and development in SL. The receptor-like cytoplasmic kinase BIK1, a negative regulator of BR-synthesis [66], was induced in SL (Additional file 1: Figure S7). Arabidopsis bik1 knockout mutants have an increased tolerance to P. brassicae lacking the typical pathogen phenotype [14,38]. For BIK1 the expression of the SL roots differed to resistant plants, as Arabidopsis bik1 knockout mutants have an increased tolerance to P. brassicae [14,38]. In Arabidopsis bik1 clubroot resistance was likely not due to the regulatory function of BIK1 for BR, but to increased PR1 expression in this mutant [38].

Transcriptomes of symptomless roots and clubroots of the same plant are markedly different
Clubroots and symptomless roots of the same P. brassicae infected plants showed markedly different gene expression patterns. The morphological changes of the clubroots go in hand with expression of genes that reduce cell wall stability and growth related processes. Genes for cell wall-loosening processes such as expansins and XTHs [67,68], were up-regulated in gall tissue (Fig. 3). Up-regulation of expansins was reported from Arabidopsis clubroot tissue [6,23] while XTH activity was reported in B. rapa clubroots [24]. Suppression of xyloglucan, xylan, hemicellulose, pectin, and lignin synthesis in gall tissues implies additional reduction of the cell wall stability and rigidity in gall tissues similar to previous findings [21,65]. GAE6 expression was reduced in the galls supporting clubroot development: Arabidopsis gae1, gae6, and gae1/gae6 mutants contained lower levels of pectin in their leaf cell walls making them more susceptible to Pseudomonas syringae and Botrytis cinerea [69]. A similar mechanism might benefit P. brassicae. Induction of the lignin pathway was an early response (48 h) in Arabidopsis [18]. In the brownish kohlrabi galls investigated in our study, lignification genes were downregulated suggesting P. brassicae suppresses host lignification in tissues where it is already established. In B. napus reduced lignification was also implied by the down-regulation of CCoAOMT (caffeoyl-CoA Omethyltransferase) [26]. Here, although not statistically significant, CCoAOMT was also lower expressed in root galls compared with SL, implying a decreased lignin biosynthesis in clubroots.
The marked hypertrophies of the plant roots go hand in hand with changes in the homeostasis of the growth hormones CK and auxins which appear to be host and time dependent [29,40]. With the exception of CKX6, cytokinin related genes were up-regulated in SL (Fig. 4). Fine tuning of the hormone balance appears to be essential in clubroot disease development [29]. Elevated CK levels are important for the onset of disease development via increasing cell divisions. However, at the onset of gall formation, CK metabolism genes including CK synthesizing and degrading enzymes are repressed [6,17]. The more active CK metabolism in the SL might interfere with clubroot development as CKX overexpressing Arabidopsis mutant showed reduced gall formation [6]. In our kohlrabi samples, CKX6 was strongly down-regulated in SL compared with WG but not compared with BG (Fig. 4). The high expression of CKX6 in WG indicates the presence of plasmodia, as the gene was found to be strongly up-regulated only in cells containing P. brassicae plasmodia [65]. Although AHK4 was described to be highly induced in infected A. thalinana plants 10 dpi [6], in a different study expression of AHK and AHP genes did not differ between infected and uninfected Arabidopsis at later time points [17]. Similar to these findings, in our much longer infected kohlrabi plants AHK and AHP genes were not differential expressed between SL and root galls. Expression levels were also very low indicating no major role of those genes in the analysed roots. The ARR5 gene was downregulated in infected roots and hypocotyl tissues 16 and 26 dpi in A. thaliana [17] but in a different study ARR5 was up-regulated 10 dpi in infected roots [6]. For Brassica hosts, the expression of ARR5 in clubroot infected roots has not been described. In our samples ARR5 levels increase in BG, but no difference is seen between SL and WG. The ARR5 gene does not appear to have a prominent role at the stage of the here investigated infections. Also evidence that P. brassicae interferes with the CK balance via the PbCKX of P. brassicae in gall tissues has not been seen in this study as the PbCKX gene was not expressed.
Myrosinases and nitrilases that can synthesize auxins from secondary metabolites or aromatic amino acids were up-regulated in SL compared with galls (Additional file 1: Figures S4), but were previously reported to be induced in Arabidopsis galls [6,70]. The auxin-induced GH3 gene family, which conjugates IAA to several amino acids, is involved in various responses of plants to abiotic and biotic stresses. The GH3.2 gene was shown to be specifically expressed in clubroots of Arabidopsis [30], and was also up-regulated in the kohlrabi galls. Expression of the P. brassicae PbGH3 gene was not detected in this study and appears not to play a role in the auxin (or JA) homeostasis at the stage of our gall samples.
Besides the described defence responses in SL that are similar to defence responses of resistance plants, we observed up-regulation of further defence related genes in SL (Fig. 5, Additional file 1: Figure S6). Homologues of the Arabidopsis Toll-IL-1 receptor disease resistance proteins TAO1, NDR1 and RPS2 genes were upregulated in SL. Those genes confer resistance to biotrophic bacterial pathogens in Arabidopsis when recognizing effectors [71,72]. In roots of beans, NDR1 also suppressed nematode parasitism by activating defence responses [73]. Thus, we found indications that those proteins might also be involved in a defence response towards P. brassicae. The identified homologs of differentially expressed phloem proteins PP2-A6 and PP2-A8 contain TIR domains and were therefore classified as Rgenes [74]. However, phloem proteins play different roles in plants [75] and their differential expression could also be a result of altered phloem development during clubroot disease [17,33,56].
As a result of gall formation and a reduced number of fine roots, clubroot infected plants face abiotic stress like lower water and nutrient supply. The differential expression of ABA related genes (Fig. 4) are therefore likely a response to the abiotic stress in the galls. Lower water supply could also be a consequence of reduced xylem production [56].

Conclusions
Clubroots and symptomless roots of the same P. brassicae infected plants showed very different gene expression patterns. The differences in the plant hormone metabolism might be responsible for the different outcomes in gall tissue and in symptomless roots as is the increased cell wall stability in symptomless roots. These results highlight, that interpreting clubroot transcriptomes or any other data originating from whole root systems might result in a dilution of biologically relevant signatures. This clearly calls for further studies analysing intra-and inter-tissue specific patterns of clubroot infected plants. As genes involved in resistance responses to P. brassicae were up-regulated in symptomless roots, this might aid the identification of novel traits for resistance breeding.

Sampling
Kohlrabi "Purple Vienna" plants with clubroots were collected from a P. brassicae infested field in Ranggen (Tyrol, Austria, 47°15′27″N 11°12′37″E) in August 2016 with the consent of the farmer. No other permissions were necessary to collect these samples. Root samples were classified by visible properties into symptomless roots (SL), smaller white spindle galls with waxy appearance (WG) and larger brownish spindle galls (BG) (Fig. 1). Samples were taken in triplicates from three individual clubroot infected plants. Galls and roots were thoroughly washed with tap water before samples were taken (categories SL, WG, and BG) and transferred to RNAlater™ (Ambion, Austin, TX, USA) until RNA extraction.

RNA extraction and sequencing
The outer layer of the root galls was trimmed-off and the trimmed galls and symptomless roots were snapfrozen in liquid nitrogen and transferred to 1.5 mL tubes containing RNase free zirconia beads (0.5 mm and 2 mm in diameter). Samples were homogenized using a Fas-tPrep (MP Biomedicals, Santa Ana, CA, USA) for 40 s at 6 m s − 1 followed by manual grinding with pestles after repeated snap-freezing. Total RNA was extracted using the Qiagen RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions, but with an additional 80% ethanol column wash prior elution. RNA quantity and quality were determined using Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA). Additional RNA quality assessment, polyA selection (SENSE mRNA-Seq Library Prep Kit; Lexogen, Vienna, Austria), library construction (9 libraries; 3x SL, 3x WG, and 3x BG) and sequencing was performed at the VBCF NGS Unit (Vienna, Austria). Sequencing was performed on the Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) with a strand specific paired end library (2 × 125 bp) using v4 chemistry.

Bioinformatics
Raw reads were quality checked using FastQC (https:// www.bioinformatics.babraham.ac.uk/projects/fastqc/). Illumina adapters were removed and good quality reads were kept using Trimmomatic v0.36 (sliding window 5 bp; average quality score > 20) [76]. Only reads with a minimum length of 75 bp were processed further after a repeated FastQC check to confirm quality improvement. Transcripts were de novo assembled using Trinity v2.2 [77] with strand specific library type (RF) and jaccard clip options. Expression estimation was performed using Trinity embedded RSEM [78] keeping only transcripts with more than at least one fragment per kilobase per million (FPKM > 1) and an isoform percentage (IsoPct) > 1%.
The assembled transcripts were blasted using BlastN [79] against the coding sequences (CDS) of B. oleracea [10; http://brassicadb.org/brad/datasets/pub/Genomes/Brassica_oleracea/V1.1/] and a custom database containing the CDS of the P. brassicae isolates e3 [12] and PT3 [13] to identify, if the transcript derived from the pathogen or host (E-value < 10 − 5 ). Transcripts with blast hits in both reference databases were analysed manually to identify their origin according to sequence identity and E-value. Transcripts with no hit in any reference were blasted (BlastP) against National Center for Biotechnology Information (NCBI) non redundant protein database and manually assigned to the corresponding species or discarded for further analysis. Transcripts with a best hit to a Brassicaceae reference sequence were assigned to the host transcriptome. Transcripts with hits to P. brassicae sequences were assigned to the pathogen transcriptome. Open reading frames (ORFs) were predicted using TransDecoder v3.0.1 (https://github.com/Trans-Decoder/). Only the longest ORF per transcript was used for further analysis. Translated peptide sequences were annotated using the KEGG (Kyoto Encyclopedia of Genes and Genomes) Automatic Annotation Server (KAAS) [80] and eggNOG-mapper v0.99.3 [81]. Kohlrabi genes were additionally annotated using Mercator [82] with default settings. Mercator categories were used for MapMan v3.6.0RC1 [83] to bin predicted genes into groups. Putative secreted proteins of P. brassicae were predicted with Phobius v1.01 [84] and SignalP v4.1 [85] in combination with TMHMM v2.0 [86]. Carbohydrate active enzymes were predicted using dbCAN [87].

Additional files
Additional file 1: Figure S1. Kohlrabi (Brassica oleracea) and Plasmodiophora brassicae reads. Samples from the symptomless roots (SL) of infected plants contained less than 0.0005% P. brassicae reads. P. brassicae read partition increased from 23% in WG to 33% in BG. Figure  S2: Kohlrabi flavonoid metabolism. Clustered heatmaps of log 2 fold change values of DEGs. No DEGs were present comparing WG with BG. Arabidopsis homologs are given. NA: not assigned. Figure S3: Kohlrabi auxin response factors. Clustered heatmaps of log 2 fold change values of DEGs. One DEG (AXR3) was present comparing WG with BG. Arabidopsis homologs are given. Figure S4: Kohlrabi myrosinases and nitrilases. Clustered heatmaps of log 2 fold change values of DEGs. No DEGs were present comparing WG with BG. Arabidopsis homologs are given. NA: not assigned. Figure S5: Kohlrabi WRKY and bZIP transcription factors. Clustered heatmaps of log 2 fold change values of DEGs. No DEGs were present comparing WG with BG. Arabidopsis homologs are given. NA: not assigned. Figure S6: KEGG map for plant-pathogen interaction. Up-regulated genes in SL compared with galls are shaded in purple and down-regulated genes in green. Genes shaded in purple and green were up-and down-regulated in different homologs or isoforms. Figure S7: Kohlrabi kinases. Clustered heatmaps of log 2 fold change values of DEGs. No DEGs were present comparing WG with BG. Arabidopsis homologs are given. NA: not assigned. Figure  S8: Numbers of P. brassicae genes in clubroot infected kohlrabi roots per COG category. Bars indicate total genes found across all libraries. Unassigned genes (n = 5482) are not illustrated. Figure S9: Cumulated FPKM values of P. brassicae reads obtained from WG and BG samples. Unassigned genes are not illustrated. Figure S10: Cumulated FPKM values of P. brassicae reads obtained from SL samples. Unassigned genes are not illustrated. (DOCX 2519 kb) Additional file 2: Table S1. Total (raw) and good quality reads of the sequenced libraries. Only reads with a minimum length of 75 bp and a sliding window of 5 bp and an average quality score > 20 were kept. Illumina adapters were removed. Table S2: Kohlrabi DEGs compared across the three root tissue types.  Table S4: Complete list of kohlrabi DEGs between BG and WG. The Functional annotations are given as returned by eggNOG mapper. Up-regulated genes are highlighted in red, down-regulated genes in blue. Table S5: Twenty highest expressed genes of Plasmodiophora brassicae in WG. Annotation given as predicted with eggNOG mapper. Putative secreted proteins are tagged (•). FPKM: fragments per kilobase (of exons) per million reads. NA: not assigned. Table S6: Twenty highest expressed genes of Plasmodiophora brassicae in BG. Annotation given as predicted with eggNOG mapper. Putative secreted proteins are tagged (•). FPKM: fragments per kilobase (of exons) per million reads. NA: not assigned. Table S7: Plasmodiophora brassicae DEGs between BG and WG. Annotations are given as predicted with eggNOG mapper. Up-regulated genes in BG are highlighted in red, downregulated genes in BG in blue. (DOCX 31 kb) Additional file 3: Kohlrabi cell wall metabolism heatmaps including transcript identifiers and cluster dendrograms. Leading zeroes were added to the gene identifier for better visualization. Additional data for each individual heatmap is provided in individual spreadsheets. Log2 fold change values and additional annotation information are given. These spreadsheets include also non-differentially expressed genes not plotted in the heatmaps. n. s.: not significantly differentially expressed. (ODS 1488 kb) Additional file 4: Kohlrabi phytohormone metabolism heatmaps including transcript identifiers and cluster dendrograms. Leading zeroes were added to the gene identifier for better visualization. Additional data for each individual heatmap is provided in individual spreadsheets. Log2 fold change values and additional annotation information are given. These spreadsheets include also non-differentially expressed genes not plotted in the heatmaps. n. s.: not significantly differentially expressed. (ODS 1435 kb)