RNA-Seq and secondary metabolite analyses reveal a putative defence-transcriptome in Norway spruce (Picea abies) against needle bladder rust (Chrysomyxa rhododendri) infection

Background Norway spruce trees in subalpine forests frequently face infections by the needle rust fungus Chrysomyxa rhododendri, which causes significant growth decline and increased mortality of young trees. Yet, it is unknown whether trees actively respond to fungal attack by activating molecular defence responses and/or respective gene expression. Results Here, we report results from an infection experiment, in which the transcriptomes (via RNA-Seq analysis) and phenolic profiles (via UHPLC-MS) of control and infected trees were compared over a period of 39 days. Gene expression between infected and uninfected ramets significantly differed after 21 days of infection and revealed already known, but also novel candidate genes involved in spruce molecular defence against pathogens. Conclusions Combined RNA-Seq and biochemical data suggest that Norway spruce response to infection by C. rhododendri is restricted locally and primarily activated between 9 and 21 days after infestation, involving a potential isolation of the fungus by a hypersensitive response (HR) associated with an activation of phenolic pathways. Identified key regulatory genes represent a solid basis for further specific analyses in spruce varieties with varying susceptibility, to better characterise resistant clones and to elucidate the resistance mechanism.


Background
Biotic stress constitutes a major threat for forest trees, negatively affecting tree growth and survival, and compromising important ecological, economical, and social forest functions [46]. Altered climatic conditions due to climate change are expected to favour insect and microbial pathogen attacks [79], and to depreciate plant defence responses [43]. Some of the economically most relevant biotic agents are fungi, and numerous pathosystems of forest trees and either ascomycetes or basidiomycetes are known (e.g. [1,20,83]). Given that many of these pathosystems have co-evolved over very long time periods, evolutionary arm races between trees and fungi have resulted in significant phenotypic variation among trees growing on the same site [47]. As a consequence, some trees often remain healthy, whereas the vast majority develops symptoms of fungal infection. This phenotypic variation has been thoroughly studied at the DNA level by identifying the responsible major QTLs (Quantitative Trait Locus) (e.g. [24,57]). In contrast, the defence transcriptome working upon infection is poorly understood. This is particularly the case in conifers, owing to the fact that genomic resources, such as annotated reference genomes and -transcriptomes, are rarely available [40]. Furthermore, studying gene expression in trees was so far restricted to the use of microarrays, which suffer from a number of technological weaknesses (e.g. cross-hybridization of related sequences) and only allowed for interrogating the expression pattern of already annotated genes [95]. Such draw backs have recently been overcome by completion of numerous reference genomes, such as for P. abies [74] and emerging next-generation-sequencing technologies (NGS), such as RNA-Seq which allows for the analysis of differential gene expression across the entire genome [48,61,89]. RNA-Seq has nowadays become a standard tool for investigating differences in gene expression under contrasting abiotic and biotic treatments, and thus may also be used to compare the expression of fungus-infected and non-infected individuals [91].
Among the various defence mechanisms against fungi described in plants, acquired or induced resistance (IR), although important [2], is not well understood and in long-living plants, such as trees, largely unexplored [22]. However, two specific molecular defence mechanisms within the IR spectrum, have been described in tree species (as reviewed in [22]): first, an inducible chemical defence that mainly comprises low molecular weight compounds (i.e. secondary metabolites) belonging to different biosynthetic pathways. Amongst them are terpenoids, phenolic compounds, and alkaloids, which are either working as phytoalexins [31] or phytoanticipins [86], depending on whether they are synthesized upon or already before an infection. In conifers, a large range of these metabolites is synthesized and accumulated in high concentrations in the needles [76,92]. Second, an inducible defence system defined by a group of low molecular weight proteins (usually < 100 kDa) can be rapidly produced upon infection. These proteins are also called pathogenesis-related proteins (PR proteins) and are mainly involved in degrading and destroying the cell wall of the fungus, either directly or indirectly by enhancing the release of additional elicitors from the fungal cell wall (e.g. [9]).
Here, we present a study on a pathosystem, which comprises the economically important European conifer Norway spruce (Picea abies) and the rust fungus Chrysomyxa rhododendri [17]. The parasite undergoes a host shift between the main host rhododendron (Rhododendron ferrugineum or R. hirsutum) and Norway spruce and its occurrence is therefore restricted to higher elevation areas, where both hosts co-occur. Attacks have been reported since the end of the nineteenth century and continue to be regularly recorded by the forest health monitoring programmes carried out in different parts of the Alps [25,96]. The infection of Norway spruce naturally occurs in spring, when freshly released basidiospores are wind-dispersed and arrive at current-year needles of the trees, where they quickly germinate and enter the young needles [8]. The infection causes a characteristic yellow needle discoloration during summer and severe needle loss in autumn [25], and repeated infections lead to significant reductions in timber yield. Infections also cause severe problems for rejuvenation at high elevation sites [25], which is critical considering the important protective function of mountain forests. Infected trees show several anatomical, morphological and physiological impairments, such as a decrease in chlorophyll content and net photosynthesis in affected needles, lower biomass production, and reduced radial and height growth [7,25,65]. Recent studies also indicated changes of the phenolic needle profile following infection and revealed phenolic concentrations of needles to correspond to their infection susceptibility [27,28].
In the present study, we investigated the activation of defence mechanisms at both expression and metabolite level by combining next-generation-sequencing (RNA-Seq) with ultra high performance liquid chromatography-mass spectrometry (UHPLC-MS). In a greenhouse experiment, Norway spruce cuttings were infected with C. rhododendri and responses compared with uninfected controls over 39 days (for detailed study design see Fig. 1). Use of clonal ramets, derived from one tree, enabled us to focus on the analysis of changes in expression patterns, minimising the influence of genotypic variation. We hypothesized that C. rhododendri infection causes differential gene expression due to an inducible defence reaction. Differences in gene expression were expected to coincide with specific changes in concentration of phenolic compounds. A comparison of infected symptomatic needles and non-symptomatic needles on plants attacked by C. rhododendri should reveal if observed defence reactions are localized or systemic.

Results
cDNA sequencing and mapping of reads to the Norway spruce reference genome A total of 19.9 million reads (median across samples) were obtained from cDNA sequencing of the libraries.
The number of input reads (after quality control trimming) varied between 18.3 million and 22.5 million depending on treatment and time point (Table 1). The number of reads that uniquely aligned to a location in the reference genome varied between 6.6 million and 17.6 million, corresponding to 36 and 78% of uniquely aligned reads, respectively. Average input read length varied slightly between 72 and 73 nt. Inter-replicate correlation between the three technical replicates per treatment and time point combination was high in all cases, Fig. 1 Study design. Design of plant treatment and sampling. Genetically identical cuttings were positioned for 39 days in a control and an infection tent (n = 3, respectively), the latter equipped for the first 21 days with rhododendron branches covered by mature C. rhododendri spore stocks. Fans improved spore distribution and microscopic slides collected spores for concentration monitoring. One current-year shoot per plant and time point was cut, the needles detached and stored at −80°C for RNA and phenolic analyses. Note that on day 39 infection symptoms were clearly visible on treated cuttings and needles were separated into healthy and symptomized ones  suggesting that cDNA sequencing produced reliable and robust results for downstream analyses (Additional file 1: Figure S1).

Differential gene expression patterns among treatments and time points
Principal component analysis of the overall differential gene expression profile revealed few differences between control and infected trees until 9 days post infection (dpi). This pattern changed drastically from 21 dpi on, when infected and non-infected trees showed significantly different gene expression profiles (Fig. 2). At that time, the first symptoms of infection were detectable, i.e., small blisters on the needle surface, visible under the microscope only. In contrast, at 39 dpi, infected trees showed the characteristic yellow needle discoloration (Additional file 2: Figure S2). Needles without symptoms sampled from the infected trees formed one cluster together with control samples, whereas needles with symptoms showed a different expression and were significantly separated along PC2.

Global analysis of differentially expressed genes (DEGs)
Overall 301, 72, 1851 and 2444 transcripts were differentially expressed in infected compared to control plants in needles at 4, 9, 21, and 39 dpi, respectively ( Fig. 3 and Additional file 3: Table S1A-D). For this and subsequent analysis, 39 dpi symptomatic needles were considered as only 11 transcripts were differentially expressed in 39 dpi non-symptomatic needles compared to control plants (Additional file 3: Table S1E). 3320 transcripts were specifically over-and under-expressed at 21 dpi and 39 dpi but not at other time points (Fig. 3), wherebỹ 30% of over-expressed and~15% of under-expressed transcripts were shared between both time points. During the infection period, the proportions of under-and over-expressed transcripts changed from~35% overexpressed transcripts at 4 and 9 dpi, to~60% at 21 and 39 dpi (Fig. 4).

Gene ontology (GO) enrichment analysis
Significantly enriched GO terms occurred in both, overand under-expressed gene sets at all analyzed time points (Fig. 5). During the onset of infection (4 dpi) metabolic processes, response to abiotic stimulus and response to stress (category biological processes), as well as catalytic, hydrolase and transferase activity (molecular functions), whereas no cellular component term were enriched for over-expressed gene sets. Regarding the under-expressed gene set, similar enriched terms plus the cellular component thylakoid were found. At 9 dpi, low number of transcripts were differentially expressed and only under-expressed DEGs were significantly enriched. Among biological processes, the most enriched terms were biosynthetic, metabolic and cellular processes, whereas among molecular functions, the most enriched terms were catalytic activity, binding and lipid binding (Fig. 5), while no cellular component term was enriched. At 21 dpi, numerous enriched GO terms occurred both in over-expressed (41 terms) and underexpressed transcripts (31 terms). Among the over- Fig. 3 Venn's diagrams of over-and under-expressed transcripts. Overlap between differentially expressed transcripts in Picea abies plants infected with Chrysomyxa rhododendri compared to control plants at the four time points after infestation (4 dpi, 9 dpi, 21 dpi, 39 dpi symptomatic needles) expressed gene set, the most enriched terms were related to protein metabolic processes, response to biotic stimulus, response to stress, metabolic processes, cellular processes, response to abiotic stimulus, signal transduction and biosynthetic process (biological processes), intracellular aspects, membrane and plasma membrane (cellular components), catalytic activity, transferase activity, binding and kinase activity (molecular function). For underexpressed transcripts, the distribution of GO categories differed, and again changes of the cellular component thylakoid were indicated. Finally, at 39 dpi, 42 and 35 GO terms were enriched in the over-and underexpressed transcript sets, respectively, and observed patterns were very similar to 21 dpi. Most striking differences were observed in the biological processes category, with enriched GO terms in over-expressed transcripts related to cell death and growth, suggesting pronounced changes after 21 dpi (Fig. 5).

Metabolic pathway analysis
With the emphasis on biochemical pathways, Kyoto Encyclopedia of Genes and Genomes (KEGG) database pathways were used as an alternative approach to categorize gene functions. On average, 40% of DEGs could be assigned to key enzymes involved in a total of 118 biological pathways, most of them at 21 and 39 dpi (Additional file 4: Table S2). Due to the complexity of the P. abies genome, several genes were assigned to the   same key enzymes of pathways involved in metabolism, genetic information processing, environmental information processing, cellular processes and organismal systems ( Table 2; Additional file 5: Table S3, Additional file 6: Table S4). Particularly notable is the over-expression of several key enzymes in the carbohydrate metabolism pathways, glycolysis/gluconeogenesis and starch and sucrose metabolism, and the under-expression of many key enzymes in the energy metabolism pathways, photosynthesis and carbon fixation. Translation, protein folding, sorting and degradation were also differentially regulated upon rust infection. Several pathways were identified to be probably involved in the defense ( Table 2) of P. abies against C. rhododendri comprising a total set of 152 DEGs (Additional file 7: Table S5). Several play a role in signal transduction and environmental adaptation (plantpathogen interaction, MAPK signaling and plant hormone signal transduction), biosynthesis of plant secondary metabolites (phenylpropanoid, flavonoid, flavone and flavonol, stilbenoid, diarylheptanoid and gingerol, terpenoid backbone), lipid biosynthesis (cutin, suberine and wax), and amino acid metabolism (phenylalanine, tyrosine and tryptophan biosynthesis). At the early stage of infection (4 dpi), several key enzymes of these pathways were under-expressed ( Table 2, Fig. 6), and after 21 dpi, the majority of key enzymes were over-expressed (Additional file 8: Figure S3). Interestingly, many key genes related to translation machinery were under-expressed in later infection stages, and few genes from cutin, suberine and wax biosynthesis were under-expressed during the entire period investigated.

Symptomatic versus non-symptomatic needles
Only few transcripts were differentially expressed in 39 dpi non-symptomatic needles compared to control plants (Additional file 3: Table S1E) and we thus focused our analysis on transcripts differentially expressed in 39 dpi symptomatic needles compared to control plants. However, it is remarkable that of the differentially expressed transcripts found by comparing 39 dpi symptomatic needles to 39 dpi non-symptomatic needles (Additional file 3: Table S1F), 61.45% were shared with 39 dpi symptomatic needles vs. control DEGs (Additional file 9: Table S6A). In addition, Gene Ontology (GO) enrichment analysis and metabolic pathway analysis revealed similar patterns for both comparisons (Additional file 9: Table S6B-F).

RT-qPCR validation of selected DEGs
23 DEGs assigned to the significantly enriched pathways plant-pathogen interaction, MAPK signaling and flavonoid biosynthesis were validated by RT-qPCR (Fig. 7). Overall RNA-Seq and RT-qPCR results were in a good agreement. Expression changes detected by RNA-Seq could be verified for the majority of transcripts. Only the over-expression of serine/threonine-protein kinase CTR1 (MA_35694g0010) at 21 dpi and 39 dpi symptomatic vs. non-symptomatic was not detected by RT-qPCR (Additional file 10: Figure S4). The most prominent expression changes were detected by RT-qPCR at 21 dpi and 39 dpi symptomatic needles (Fig. 7). The highest levels of over-expression at 21 dpi were found for the basic endochitinase B gene (CHIB, MA_8921185g0010: 738-fold increase; and MA_10313114g0010: 33-fold) and the pathogenesis-related protein 1 (PR1, MA_ 53673g0010: 22-fold). Under-expression among the selected transcripts was only detected at the early time points of infection. No major changes were found in non-symptomatic needles at 39 dpi.

Changes in phenolic profiles
Needles of control plants showed a rather constant concentration of shikimic acid during the experiment, while flavonoid concentrations constantly decreased until 21 dpi and stilbenes strongly accumulated towards the end (Fig. 8). In cuttings attacked by the fungus, significantly higher levels of shikimic acid were detected at 9 dpi, which decreased constantly to far below the level of control plants until 39 dpi. Flavonoid concentrations significantly increased in infected plants at 9 and 21 dpi, mainly due to an increase of kaempferol and quercetin on 9 dpi and of quercetin 3-glucoside on 21 dpi (Fig. 8, Additional file 11: Figure S5). Additionally, taxifolin significantly increased in infected symptomatic needles at 39 dpi. In contrast, stilbene accumulation was reduced in attacked cuttings, mainly due to lower concentration of the compound astringin. Healthy and symptomized needles originating from the same shoot of treated plants at 39 dpi differed strongly in their phenolic composition (Fig. 8). The healthy needles corresponded precisely to needles of control plants with respect to all analysed metabolites (Fig. 8, Additional file 11: Figure S5).

Discussion
Plant defence systems against pathogens are well studied in model organisms and crops, but largely unexplored in conifers. In the present study, we combined RNA-Seq and LC-MS to investigate the defense response of Norway spruce infected by C. rhododendri [25,27,28]. This combined approach has been proven to be very useful in the study of molecular defense mechanisms of host plants against pathogens [63,69,90]. RT-qPCR was used to validate the RNA-seq results. Expression changes detected by RNA-seq could be verified for 22 transcripts out of the 23 DEGs selected for validation (Additional file 10: Figure S4). Differences in the relative expression values might be explained by the different normalisation processes implemented by the two methods. The combination of both datasets, however, improves the identification of the most affected targets, especially for those P. abies transcripts that were assigned to identical gene functions.

General host response
An overall analysis of GO Enrichment terms showed that from the beginning of the infection process (4 dpi) onwards, numerous genes related to stress responses were differentially regulated (Fig. 5). In particular, enriched terms at 21 dpi indicated an activation of genes related to plant defense, plant secondary metabolism, signal transduction and kinases, which may play important roles in signaling, communication and coordination of defense reactions. Also, enriched terms indicated effects at the thylakoids, which is in agreement with reductions in chlorophyll content and photosynthesis observed in field studies [64]. Over-expressed transcripts at the end of the infection experiment (39 dpi) related to cell death suggest the occurrence of a hypersensitive response (HR), which prevents the spread of infections by the rapid death of cells surrounding the infected area [54]. Also many key enzymes from the ubiquitin mediated proteolysis and proteasome were over-expressed from 21 dpi on, which play a key role in plant immunity [18,81,85] and degradation of unneeded or damaged proteins in cells undergoing programmed cell death (Table 2). Accordingly, distinct necrotic areas in partly infected Norway spruce needles were described in Mayr et al. [66]. At 39 dpi several key enzymes related to endocytosis, phagosome [70] and peroxisome [38,45] were overexpressed. The onset of this strong defense response occurred relatively late compared to other pathosystems (reactions within hours or days; e.g. in [90,94]) which may be due to overall low metabolic activity in trees adapted to high elevation and extreme climatic conditions [37].

Plant-pathogen interaction and signal transduction
During infection, plant-pathogen interaction pathways were considerably differentially regulated ( Table 2, Fig. 6, Additional file 7: Table S5, Additional file 8: Figure S3A). Among identified DEGs, some are involved in recognition of pathogen-associated molecular patterns (PAMPs) by plasma membrane receptors (PAMP-triggered immunity PTI), while others are involved in the detection of pathogen effectors by intracellular or plasma membrane-localized immune receptors (effector-triggered immunity -ETI) (review by [75]). Interestingly, at 4 dpi two calmodulins (MA_10204459g0010; MA_10280648g0010) and two HSP90 (MA_10431770g0010 and MA_10434030g0010) were under-expressed, indicating blockage of the Ca 2+ signaling pathway and thus the defense response via HR, cell wall reinforcement and stomata closure, and probably enabled C. rhododendri to colonize other tissues [55,67]. Rust fungus can interfere with a variety of host defense pathways and suppress multiple plant defense responses [59,78]. In contrast, at 21 dpi, 9 key enzymes (22 genes) involved in both PTI and ETI and including several calmodulines as well as HSP90 were over-expressed. It is likely that Norway spruce, at this infection phase, activated defense reactions against the fungus via stomatal closure, hypersensitive response (HR) and cell wall reinforcement [30]. This response is also promoted by the over-expression of MA_ 9458g0010, a calcium-dependent protein kinase (CDPK) putatively involved in abscisic acid-activated signaling pathway and intracellular signal transduction upon phosphorylation by fungal elicitors. Finally, DNA expression promoted by this signaling pathway enhanced the expression of two (MA_501572g0010, MA_53673g0010) pathogenesis-related protein 1 (PR1), which is induced in response to a variety of pathogens and was shown to be a useful molecular marker for the salicylic acid-mediated disease resistance [10].
Furthermore, the MAPK signaling pathway and plant hormone signal transduction were differentially regulated: In both cases, at least 11 key enzymes were overexpressed from 21 dpi on, which again indicates that the fundamental plant defense reactions started between 9 and 21 dpi (Additional file 8: Figure S3B-C). In the MAPK signaling pathway, mitogen-activated protein kinase 6 (MPK6: MA_10437020g0010) is a key protein in several signaling routes and plays a critical role in plant disease resistance [68]. Jagodzik et al. [41] reported a complex crosstalk between the MAPK cascade and hormone signaling pathways to orchestrate plant defense. Accordingly, in infected Norway spruce, over-expressed genes included serine/threonine-protein kinase (CTR1: MA_35694g0010), mitogen-activated protein kinase kinase 9 (MKK9: MA_10194g0020), ethylene-insensitive protein 3 (EIN3: MA_20516g0010) and transcription factor MYC2 (MYC2: MA_10435905g0020), promoting defense responses via DNA expression of several genes as basic endochitinases B (MA_10313114g0010; MA_ 8921185g0010), an effective defense against chitincontaining fungal pathogens [42,44,51,87]. In addition, plant defensins well known efficient antimicrobial peptides against fungi [53] were also over-expressed (MA_ 5320g0010, MA_1489g0010, MA_4353361g0010) in response to infection [72]. In parallel, plant growth is affected by under-expression of several genes involved in the auxin signaling pathway, contributing to an optimized use of resources during this biotic stress.

Primary and secondary metabolism
Both the primary and secondary metabolism [4] of infected plants were highly affected, especially from 21 dpi on (Table 2, Fig. 6). The over-expression of genes involved in the biosynthesis of secondary metabolites (Additional file 8: Figure S3 D-I) is one of the frequently observed response patterns [3,12,13], and frequently comes at the cost of an under-expression of key enzymes involved in photosynthesis and carbon fixation. An extensive list of genes involved in photosynthesis were under-expressed, in accordance with the yellow discoloration of rust infected needles and the documented reduction in chlorophyll content and photosynthetic activity [7,64]. In contrast, many key enzymes of starch and sucrose metabolism were over-expressed by the end of the infection process and DEGs indicated increased glycolysis. Lipid metabolism was also affected and several genes involved in cutin, suberin and wax biosynthesis were under-expressed during infection (Additional file 8: Figure S3 J). Infected cuttings obviously changed their metabolism from autotrophy (photosynthesis) to a consumption of resources [80,82] and optimized limited resources by use of lipids on alternative pathways [56].
Many transcripts differentially expressed in infected needles are involved in the biosynthesis of secondary metabolites, among them phenolic compounds and terpenoids. The antifungal activity of these compounds has been extensively described [14,32,34,35,52,77] and they are considered to play a central role in conifer response to fungal pathogens [5,11,14,62,92], including spruce response to needle rust [28]. Accordingly, numerous pronounced effects were observed in the infection experiment: In the terpenoid backbone biosynthesis pathway genes belonging to the mevalonate pathway (cytoplasm) were over-expressed in comparison to the nonmevalonate pathway (plastid) (Additional file 8: Figure S3 I). Variations in phenolic profiles were detected from 4 dpi on, whereby concentrations of individual compounds showed highest changes at 9 dpi (increased concentration levels) and 39 dpi (mainly decreased concentration levels). Please note that constitutive concentration levels (in healthy spruce needles) also changed considerably during the experiment (Fig. 8, Additional file 11: Figure S5), as the phenolic profile changes during needle development [28]. A key molecule in the plant secondary metabolism is shikimic acid, as it is the precursor of phenylalanine, tryptophan, tyrosine (Additional file 8: Figure S3D), and subsequently for the main phenolic pathway [88]. Therefore, the increase and decline of shikimic acid at 4 and 9 dpi can be directly related to a suppressed or enhanced phenolic production downstream. In addition, at 9 dpi an overexpressed transaldolase on the pentose phosphate pathway (MA_30531g0010), which produces the precursor Derythrose 4-phosphatemay explain the increased concentrations in infected plants. It remains unclear, why later (39 dpi), shikimic acid was reduced in infected compared to control trees: production may have been limited due to reduced photosynthetic capacity and respective resources, or due to under-expression of the first key enzyme on the shikimate pathway (3-deoxy-7-phosphoheptulonate synthase; MA_54864g0010). Alternatively, shikimic acid may have been consumed downstream by the production of plant secondary metabolites, as indicated by over-expression of many phenylpropanoid and flavonoid biosynthesis genes from 21 dpi on [16].
Due to the complexity and interconnection of phenolic pathways and the different nature of involved enzymes it is hardly possible to directly link over-and under-expressed genes in the individual pathways directly to measured metabolite concentrations. However, chalcone synthase (CHS: MA_20764g0010 and MA_5735g0010), a key enzyme related to the flavonoid biosynthesis, was under-expressed at 4 dpi and over-expressed at 21, and corresponded well to flavonoids concentration at these time points (Fig. 8, Additional file 8: Figure S3G). Up-regulation of chalcone synthase was shown to be central in the defence mechanisms of Norway spruce [73]. In addition, over-expression of a flavonoid 3`-monooxygenase (CYP75B1: MA_ 76780g0010) at 21 dpi and 39 dpi corresponded to concentration changes in quercetin 3-glucoside. Nevertheless, Flavanone-3-hydroxylase and Flavonoid 3′,5′-hydrolase genes that plays an important role in the biosynthesis of spruce phenolic defenses against bark beetle-associated fungus Endoconidiophora polonica [34,35] were not differentially expressed here. Different pathogens and infected tissues may result in different responses even some similarities can be found. Here, high concentration of taxifolin in symptomatic needles is remarkable as it is a relevant antifungal flavonoid induced by various infections [11,21,28,35,49].
On the other hand the decline of picein and stilbene concentrations at 39 dpi, when needle tissues were densely pervaded by fungal hyphae, may be explained either (i) through the reduction on investments in the plant defence when tissue damages are already severe [7], (ii) metabolization of those compounds by the fungus (see e.g. [33]) or (iii) conversion, modification, and incorporation into the cell wall [23].
Like in the flavonoid biosynthesis, also in the phenylpropanoid biosynthesis pathway a large number of genes were over-expressed at 21 and 39 dpi, but partially also at 4 dpi (Additional file 7: Table S5, Additional file 8: Figure S3E). In particular, several phenylalanine ammonia lyases (PAL), and some peroxidases (PX3), but likewise also a class IV chitinase (CHI4; MA_ 10427514g0010 were differentially over-expressed as reported also for spruce infected by Ceratocystis polonica [72] and Heterobasidion annosum [36]. In any case, combined RNA-Seq and biochemical data suggest that the rust infection provoked complex and dynamic regulative changes in the phenolic pathways, leading to altered metabolite concentration in damaged spruce needles.

Local versus systemic defence response
The comparison of 39 dpi symptomatic and nonsymptomatic needles confirms the local activation of molecular defence mechanisms based on both inducible chemical defence (plant secondary metabolites) and low-molecular-weight proteins (Fig. 6). Interestingly, we observed no systemic induction of defence response in healthy needles of infected cuttings, neither on the RNA expression nor on the metabolic level (Figs. 2 and 8). Similar defence response was found under Ceratocystis polonica infection [50]. Systemic acquired resistance represents, besides constitutive and inducible defence, a third defence strategy, and was hypothesized to be a common phenomenon also in conifers [9,49,92]. However, as the complex life cycle of C. rhododendri limits infection events to a short period of the year, an acquired resistance may be of little advantage and thus not outweigh the respective costs for the plant.

Conclusions
Norway spruce shows a locally activated response to infection by C. rhododendri, which starts at 4 dpi but is pronounced between 9 and 21 dpi and mainly aims at isolating the fungus by a hypersensitive response (HR) in conjunction with an activation of phenolic pathways. This has been extensively described also for other pathosystems involving plants and fungi [15,90] and probably differs between resistant and susceptible spruce phenotypes [16]. The presented study represents a solid basis for further specific analyses of identified key regulatory genes in clones with varying resistance.

Plant material and experimental design
The experiment was conducted with five-year-old cuttings gained from the Norway spruce clone ASS-7, originating from a highly needle rust affected subalpine site in Tyrol, Austria [27]. The adult tree ASS-7 was shown to exhibit a slightly lower percentage of infected needles compared to surrounding trees in the field, but nevertheless showed clear symptoms of needle rust infection. We therefore expected this clone to show pronounced responses upon infection on the molecular level. By using genetically identical cuttings for both the infection treatment and controls we avoided genotype influences on gene expression patterns and enabled three clonal replicates within the groups. Cuttings were grown in pots under field conditions in the forest garden Bad Häring, Tyrol, Austria and, one month before start of the treatment, six equally grown individuals were moved to the Botanical Garden of the Department of Botany of the University of Innsbruck.
The experiment was performed in a greenhouse to control climatic and infection conditions, but simulated the field situation, where Norway spruce is constantly exposed to this pathogen spores during the infection period. Inoculation treatment was started when buds opened and needles started to unfold, which corresponds to the period when needles are susceptible to the fungus [26]. Cuttings were split into a treatment (n = 3) and control (n = 3) group and placed in two different custom made plastic tents sized 1x1x1m. Tents were equipped with fans for better fungal spore distribution and with microscopic slides for monitoring of spore densities (Fig. 1). In the infection tent, cuttings were exposed to C. rhododendri spores by placing rhododendron branches with mature spore stocks above the cuttings. Spore exposure was continued from day 1 to day 21, to mimic a situation that plants face in the natural environment. All plants were regularly watered and repeatedly randomly rearranged within the tent to ensure uniform spore exposure. Spore deposition within the exposure time equalled 2.04 ± 0.17 spores per mm 2 . On infected cuttings, 51 to 76% of current year needles developed clear infection symptoms, while control cuttings remained completely healthy. First infection symptoms (small blisters at the needle surface) appeared 21 days after start of the treatment, while the clear yellow discoloration was visible one week later (Additional file 2: Figure S2).
Needle samples were taken at 1, 4, 9, 21, and 39 days post infection (dpi) by cutting one current-year shoot per plant and distributing needles on two RNase-free tubes for (a) RNA-Seq analysis and (b) phenolic analysis. At 1 dpi, only control plants were sampled. At 39 dpi, samples from infected trees were further subdivided into symptomized and healthy needles (Fig. 1), as infection symptoms were clearly visible. Samples were immediately submersed in liquid nitrogen and stored at − 80°C until further analysis.
RNA isolation, sequencing and mapping to reference genome RNA was extracted from a total of 30 Norway spruce samples by homogenizing needle samples with mortar and pestle and subsequently mixing with lysis solution in order to preserve RNA integrity. The Protocol B Sigma-Aldrich Spectrum Plant Total RNA Kit (PN: STRN50-1KT) was used for extraction according to the manufacturer's guidelines (Sigma Aldrich, St. Louis, USA). RNA concentrations were measured with a Nano-Drop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham MA, USA) and approximately 1 mg of DNAse I treated RNA was diluted in a total reaction volume of 10 μl. Quality of DNAse I treated total RNA was determined by using an Agilent 2100 Bioanalyzer with RNA 6000 Nano kit as plant run mode class (Agilent technologies, Santa Clara CA, USA). 5 μl of DNAse I digested total RNA was processed with QuantSeq 3′ mRNA-Seq FWD Library Preparation Kit and each reaction was spiked in with SIRV-Set 3 variant controls by following the manufacturer's guidelines (Lexogen, Vienna, Austria). Quality of the libraries was determined with an Agilent 2100 Bioanalyzer DNA High Sensitivity Kit (Agilent technologies, Santa Clara, USA). Samples were pooled in equimolar ratio and the library pool was quantified using aQubit dsDNA HS assay kit (Thermo Fisher Scientific, Waltham MA, USA). Sequencing was performed on an Illumina NextSeq 500 system with SR75 High Output Kit at Lexogens (Lexogen, Vienna, Austria). The obtained reads were mapped to the Norway spruce reference genome [74] by using the open-source software STAR [19]. We used the alternate protocol 1 as described in Dobin & Gingeras [19] and only those contigs that had a gene annotation in the Norway spruce genome assembly were retained for mapping. Since the Norway spruce reference genome is approximately 20 GB in size [74], genome indices were produced by using a 189 GB RAM server according to the required Genome size/bytes ratio described in Dobin & Gingeras [19]. The detailed STAR settings used for mapping and genome index generation can be found in the Additional file 12: Command S1.

Differential gene expression analysis
For differential gene expression analysis the DESeq2 package [60], which is implemented in the Bioconductor platform in R [29], was used. Briefly, DESeq2 transforms read counts per gene obtained from RNA-Seq data into systematic contrasts of gene expression across experimental conditions. For this, DESeq2 uses so called shrinkage estimators which were shown to better estimate and compare fold changes in gene expression across treatments compared to other available methods [60]. The null hypothesis tested implies that the logarithmic fold change in gene expression between infected and control plants is exactly zero. We first used the DESeqDataSetFromHTSeqCount function for importing read count data and pre-filtered for rows with zero counts. Subsequently, differential gene expression analysis was performed by using the DESeq function which generates the log2fold changes with adjusted p-values for each treatment per time point combination (infected vs control at 4, 9, 21, 39 dpi symptomatic and nonsymptomatic needles). Differential gene expression analysis was performed also for 39 dpi symptomatic vs 39 dpi non-symptomatic needles.
Venn diagrams of over-and under-expressed transcripts were created by using jvenn [6]. Nowadays, Con-GenIE database (http://congenie.org/enrichment) [84] provides useful information of the ongoing P. abies sequencing project, even if genes are not yet annotated hindering the work of mining. Nevertheless, gene role or function can be investigated further by using as starting point its PFAM domains or Gene ontology (GO) terms.
GO enrichment analysis with Fisher exact tests on differentially expressed transcripts was performed with functional enrichment analysis tools available at ConGenIE, using p-values corrected with a false discovery rate of (FDR) < 0.05. GO-Slim was run to reduce complexity of GO terms for gene class analysis. Based on GO-Slim annotations, over-and under-expressed transcripts were classified into three ontological categories: cellular component, biological process, and molecular function.
Homology is a fruitful tool for mapping conserved genes (orthologous) and here we try to elucidate changes in well described and characterized pathways and genes. Metabolic pathway and orthology-oriented functional annotations of the protein sequences (.fasta files downloaded at ConGenIE) were performed against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [39] using the KEGG Automatic Annotation Server (KAAS; https://www.genome.jp/tools/kaas/) [71]. KEGG Orthology (KO) assignment was applied using the Bidirectional Best Hit (BBH) method and all datasets of dicot plants (excluding Camelina sativa, Lupinus angustifolius and Cucurbita moschata), monocot plants and a basal magnoliophyta (Amborella family) were selected as reference organisms (Additional file 13: Table S7). This analysis was applied at each time point for differentially over-and under-expressed genes (DEGs). KEGG Orthology (KO) output htext was converted to tabulated excel file (Additional file 6: Table S4). Norway spruce pathway hierarchy model was created by merging several plant reference pathway hierarchies available at KEGG (Additional file 14: Table S8) and this was used to standardize the KO assignment results (Additional file 6: Table S4).

RT-qPCR
Twenty three DEGs were selected for RT-qPCR validation. Primers were designed using the PrimerQuest tool (Integrated DNA Technologies, Coralville, IA, USA). RT-qPCR assay details are listed in Additional file 15: Figure  S6. The amplification efficiency of each assay was calculated based on the standard curve from a four-fold dilution series of cDNA (0.025-100 ng). Seven candidate reference genes (RGs) were tested on a selection of 21 P. abies test RNAs that reflect the experimental setup. The expression stability of the RGs was assessed using the RefFinder tool [93]. Two of the most stably expressed genes, actin (ACT) and ubiquitin (UBI), were selected for normalisation (Additional file 15: Figure S6). Reverse transcription (RT) of 1 μg total RNA was done with the Super-Script III First-Strand Synthesis SuperMix (Invitrogen, Carlsbad, CA, USA) with oligo (dT) primers according to the manufacturer's protocol. No RT controls (without enzyme mix) were included to monitor the amplification of contaminating DNA. RT-qPCR reactions were conducted in 20 μl reaction volumes including 1x HOT FIREPol EvaGreen qPCR Mix Plus ROX (Solis BioDyne, Tartu, Estonia), 200 nM of each primer and 25 ng cDNA. Each sample was analysed in duplicates on the AriaMx Real-time PCR System (Agilent, Santa Clara, CA, USA) with the following temperature protocol: initial activation at 95°C for 12 min, 40 cycles of 95°C for 15 s and 60°C for 1 min, followed by a melting curve analysis step (60°C -95°C). The efficiencycorrected target gene Cq values were normalised to the geometric mean of ACT and UBI, and relative expression changes compared to non-infected controls and for 39 dpi also for symptomatic vs. nonsymptomatic needles, were calculated with the comparative 2 -ΔΔCT method [58]. Mean Cq values obtained by RT-qPCR are enclosed in Additional file 16: Table S9. Statistical analyses (two-tailed unpaired t-test with Welch's correction) of the RT-qPCR data were done GraphPad Prism 5 software (GraphPad Software, San Diego, CA, USA). A p value of < 0.05 was considered as significant.
For phenolic concentrations, differences were tested (1) between time points, (2) between the control and treatment group, and (3) between healthy and symptomized needles of the treatment group at 39 dpi. Differences were analysed with the Bonferroni test (for data with homogeneity of variance) or Tamhane test (no homogeneity of variance) after testing for homogeneity (Levene test) and for Gaussian distribution (Kolmogorov-Smirnov test). All tests (twotailed) were performed pairwise at a probability level of 5% using SPSS (version 24; SPSS, IL, USA), and all values are given as mean ± standard error (SE). flavonol biosynthesis, I) Terpenoid backbone biosynthesis, J) Cutin, suberine and wax biosynthesis. At 39 dpi symptomatic needles were compared to controls.
Additional file 11: Figure S5. Concentration changes of all phenolic needle metabolites. Concentration changes of individual phenolic compounds during the experiment in needles of control (filled symbols, mean ± SE, n = 3) and spore exposed spruce cuttings (open symbols, mean ± SE, n = 3). Triangles on day 39 indicate concentrations in healthy needles of treated cuttings (mean ± SE, n = 3). Concentration values are given as μmol g − 1 dry weight; compounds with needle concentrations below the quantification threshold (naringenin, quercitrin, piceid, piceatannol, resveratrol, chlorogenic acid, gallic acid) are not shown. Significant differences between the control and treated group are indicated with asterisks and between healthy and symptomized needles of treated cuttings on day 39 with a circle.
Additional file 12: Command S1. STAR command. The detailed STAR settings used for mapping and genome index generation.
Additional file 13: Table S7. Sequence datasets of plant species applied for KEGG Orthology (KO) assignment using the KEGG Automatic Annotation Server (KAAS). KEGG Orthology (KO) assignment was applied using the Bi-directional Best Hit (BBH) method and all datasets of dicot plants (excluding Camelina sativa, Lupinus angustifolius and Cucurbita moschata), monocot plants and a basal magnoliophyta (Amborella family) were selected as reference organisms. (XLS 1108 kb) Additional file 14: Table S8. Picea abies pathway hierarchy model. This model was created by merging several plant reference pathway hierarchies available at KEGG and was used to standardize the KO assignment results (Additional file 6: Table S4). (XLS 25 kb) Additional file 15: Figure S6. RT-qPCR assay details. Primer sequences of RT-qPCR assays.