- Research article
- Open Access
RNA-Seq and secondary metabolite analyses reveal a putative defence-transcriptome in Norway spruce (Picea abies) against needle bladder rust (Chrysomyxa rhododendri) infection
BMC Genomics volume 21, Article number: 336 (2020)
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.
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.
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.
Biotic stress constitutes a major threat for forest trees, negatively affecting tree growth and survival, and compromising important ecological, economical, and social forest functions . Altered climatic conditions due to climate change are expected to favour insect and microbial pathogen attacks , and to depreciate plant defence responses . 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 . 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 . 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 . Such draw backs have recently been overcome by completion of numerous reference genomes, such as for P. abies  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 .
Among the various defence mechanisms against fungi described in plants, acquired or induced resistance (IR), although important , is not well understood and in long-living plants, such as trees, largely unexplored . However, two specific molecular defence mechanisms within the IR spectrum, have been described in tree species (as reviewed in ): 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  or phytoanticipins , 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. ).
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 . 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 . The infection causes a characteristic yellow needle discoloration during summer and severe needle loss in autumn , and repeated infections lead to significant reductions in timber yield. Infections also cause severe problems for rejuvenation at high elevation sites , 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.
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, 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), whereby ~ 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% over-expressed 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, over- and 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 under-expressed transcripts (31 terms). Among the over-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 under-expressed 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 under-expressed 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 (plant-pathogen 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).
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 . 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 . 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. . At 39 dpi several key enzymes related to endocytosis, phagosome  and peroxisome [38, 45] were over-expressed. 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 .
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 ). 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 Ca2+ 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 . 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 .
Furthermore, the MAPK signaling pathway and plant hormone signal transduction were differentially regulated: In both cases, at least 11 key enzymes were over-expressed 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 . Jagodzik et al.  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 chitin-containing fungal pathogens [42, 44, 51, 87]. In addition, plant defensins well known efficient antimicrobial peptides against fungi  were also over-expressed (MA_5320g0010, MA_1489g0010, MA_4353361g0010) in response to infection . 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  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 .
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 . 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 non-mevalonate 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 . 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 . 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 over-expressed transaldolase on the pentose phosphate pathway (MA_30531g0010), which produces the precursor D-erythrose 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 .
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 . 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 , (ii) metabolization of those compounds by the fungus (see e.g. ) or (iii) conversion, modification, and incorporation into the cell wall .
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  and Heterobasidion annosum . 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 non-symptomatic 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 . 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.
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 . 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 . 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 . 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 mm2. 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 NanoDrop 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  by using the open-source software STAR . We used the alternate protocol 1 as described in Dobin & Gingeras  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 , genome indices were produced by using a 189 GB RAM server according to the required Genome size/bytes ratio described in Dobin & Gingeras . 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 , which is implemented in the Bioconductor platform in R , 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 . 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 non-symptomatic 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 . Nowadays, ConGenIE database (http://congenie.org/enrichment)  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  using the KEGG Automatic Annotation Server (KAAS; https://www.genome.jp/tools/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 (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).
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 . 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 SuperScript 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 efficiency-corrected 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. non-symptomatic needles, were calculated with the comparative 2-ΔΔCT method . 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.
Identification and quantification of phenolic compounds
Phenolic analysis was conducted as described in Ganthaler et al.  and included the identification and quantification of eleven flavonoids (kaempferol, kaempferol 3-glucoside, kaempferol 7-glucoside, kaempferol 3-rutinoside, quercetin, quercetin 3-glucoside, quercitrin, naringenin, taxifolin, catechin, gallocatechin), five stilbenes (astringin, isorhapontin, piceid, piceatannol, resveratrol), three simple phenylpropanoids (picein, gallic acid, chlorogenic acid) and one precursor (shikimic acid). Briefly, needles were freeze-dried, homogenized, and extracted two times with 1 ml 95% (v/v) ethanol, containing 2 μmol L− 1 orientin, pinosylvin and naringin as internal quantification standards. Liquid chromatography-mass spectrometry (UHPLC-MS) was conducted using an ekspert ultraLC 100 UHPLC system with a reversed-phase UHPLC column (NUCLEODUR C18 Pyramid, EC 50/2, 50 × 2 mm, 1.8 μm, Macherey-Nagel, Düren, Germany) coupled to a QTRAP 4500 mass spectrometer (AB SCIEX, Framingham, MA, USA) operated in negative ion mode using multiple reaction monitoring (MRM). Based on retention time and MRM transition (for detailed parameters see ), and calibration curves of authentic samples of all substances, peaks were automatically detected and normalized relative to the internal standards, and concentrations were calculated using the software Analyst (version 1.6.3) and MultiQuant (version 2.1.1, both AB SCIEX, Framingham, MA, USA).
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 (two-tailed) 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).
Availability of data and materials
Raw sequence data have been submitted to the NCBI Short Read Archive (SRA) under accession number PRJNA579736.
Bi-directional best hit
Calcium-dependent protein kinase
Differentially expressed genes
Days post infection
False discovery rate
KEGG automatic annotation server
Kyoto encyclopedia of genes and genomes
Mitogen-activated protein kinase
Multiple reaction monitoring
Next generation sequencing
Pathogen-associated molecular patterns
Principal component analysis
Quantitative trait locus
Quantitative reverse transcription PCR
Ultra high performance liquid chromatography-mass spectrometry
Adams GC, Roux J, Wingfield MJ. Cytospora species (Ascomycota, Diaporthales, Valsaceae): introduced and native pathogens of trees in South Africa. Austral Plant Pathol. 2006;35:521–48.
Agrawal AA, Tuzun S, Bent E. Induced plant defenses against pathogens and herbivores. St. Paul: APS Press; 1999.
Ali MB, McNear DH. Induced transcriptional profiling of phenylpropanoid pathway genes increased flavonoid and lignin content in Arabidopsis leaves in response to microbial products. BMC Plant Biol. 2014;14:1.
Aliferis KA, Faubert D, Jabaji S. A metabolic profiling strategy for the dissection of plant defense against fungal pathogens. PLoS One. 2014;9(11):e111930.
Bahnweg G, Schubert R, Kehr RD, Müller-Starck G, Heller W, Langebartels C, et al. Controlled inoculation of Norway spruce (Picea abies) with Sirococcus conigenus: PCR-based quantification of the pathogen in host tissue and infection-related increase of phenolic metabolites. Trees. 2000;14:435–41.
Bardou P, Mariette J, Escudie F, Djemiel C, Klopp C. Jvenn: an interactive Venn diagram viewer. BMC Bioinformatics. 2014;15:293.
Bauer H, Plattner H, Volgger W. Photosynthesis in Norway spruce seedlings infected by the needle rust Chrysomyxa rhododendri. Tree Physiol. 2000;20:211–6.
Bennell AP. Rhododendron rust - taxonomic and horticultural implications. Notes R Bot Gard Edinb. 1985;43:25–52.
Bonello P, Gordon TR, Herms DA, Wood DL, Erbilgin N. Nature and ecological implications of pathogen-induced systemic resistance in conifers: a novel hypothesis. Physiol Mol Plant Pathol. 2006;68:95–104.
Breen S, Williams SJ, Outram M, Kobe B, Solomon PS. Emerging insights into the functions of pathogenesis-related protein 1. Trends Plant Sci. 2017;22:871–9.
Brignolas F, Lacroix B, Lieutier F, Sauvard D, Drouet A, Claudot A, et al. Induced responses in phenolic metabolism in two Norway spruce clones after wounding and inoculations with Ophiostoma polonicum, a bark beetle-associated fungus. Plant Physiol. 1995;9:821–7.
Carrasco A, Wegrzyn JL, Durán R, Fernández M, Donoso A, Rodriguez V, et al. Expression profiling in Pinus radiata infected with Fusarium circinatum. Tree Genet Genomes. 2017;13:46.
Cheng X-J, He B, Chen L, Xiao S-Q, Fu J, Chen Y, et al. Transcriptome analysis confers a complex disease resistance network in wild rice Oryza meyeriana against Xanthomonas oryzae pv. Oryzae. Nature. 2016;6:38215.
Chong J, Poutaraud A, Hugueney P. Metabolism and roles of stilbenes in plants. Plant Sci. 2009;177:143–55.
Coll NS, Epple P, Dangl JL. Programmed cell death in the plant immune system. Cell Death Differ. 2011;18:1247–56.
Danielsson M. Chemical and transcriptional responses of Norway spruce genotypes with different susceptibility to Heterobasidion spp. infection. BMC Plant Biol. 2011;11:154.
De Bary A. Aecidium abietinum. Bot Z. 1879;37:761–74 777–789, 801–811, 825–830, 840–847.
Devoto A, Muskett PR, Shirasu K. Role of ubiquitination in the regulation of plant defence against pathogens. Curr Opin Plant Biol. 2003;6(4):307–11.
Dobin A, Gingeras TR. Mapping RNA-seq reads with STAR. Curr Protoc Bioinformatics. 2015;51:11–4.
Ennos RA. Resilience of forests to pathogens: an evolutionary ecology perspective. Forestry. 2015;88:41–52.
Evensen PC, Solheim H, Hoiland K, Stenersen J. Induced resistance of Norway spruce, variation of phenolic compounds and their effects on fungal pathogen. For Pathol. 2000;30:97–108.
Eyles A, Bonello P, Ganley R, Mohammed C. Induced resistance to pests and pathogens in trees. New Phytol. 2010;185:893–908.
Fossdal CG, Nagy NE, Hietala AM, Kvaalen H, Slimestad R, Woodward S, et al. Indications of heightened constitutive or primed host response affecting the lignin pathway transcripts and phenolics in mature Norway spruce clones. Tree Physiol. 2012;32:1137–47.
Freeman JS, Potts BM, Vaillancourt RE. Few Mendelian genes underlie the quantitative response of a Forest tree, Eucalyptus globulus, to a natural fungal epidemic. Genetics. 2008;178:563–71.
Ganthaler A, Bauer H, Gruber A, Mayr M, Oberhuber W, Mayr S. Effects of the needle bladder rust (Chrysomyxa rhododendri) on Norway spruce: implications for subalpine forests. Eur J For Res. 2014;133:201–11.
Ganthaler A, Mayr S. Temporal variation in airborne spore concentration of Chrysomyxa rhododendri: correlation with weather conditions and consequences for Norway spruce infection. Forest Pathol. 2015;45:443–9.
Ganthaler A, Stoggl W, Mayr S, Kranner I, Schuler S, Wischnitzki E, et al. Association genetics of phenolic needle compounds in Norway spruce with variable susceptibility to needle bladder rust. Plant Mol Biol. 2017a;94:229–51.
Ganthaler A, Stöggl W, Kranner I, Mayr S. Foliar phenolic compounds in Norway spruce with varying susceptibility to Chrysomyxa rhododendri: analyses of seasonal and infection-induced accumulation patterns. Front Plant Sci. 2017b;8:1173.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5:R80.
Gudesblat GE, Torres PS, Vojnov AA. Stomata and pathogens: warfare at the gates. Plant Signal Behav. 2009;4:1114–6.
Hammerschmidt R, Nicholson RL. A survey of plant defense responses to pathogens. In: Agrawal AA, Tuzun S, Bent E, editors. Induced plant defenses against pathogens and herbivores. St. Paul: APS Press; 1999. p. 55–71.
Hammerbacher A, Ralph SG, Bohlmann J, Fenning TM, Gershenzon J, Schmidt A. Biosynthesis of the major tetrahydroxystilbenes in spruce, astringin and isorhapontin, proceeds via resveratrol and is enhanced by fungal infection. Plant Physiol. 2011;157:876–90.
Hammerbacher A, Schmidt A, Wadke N, Wright LP, Schneider B, Bohlmann J, et al. A common fungal associate of the spruce bark beetle metabolizes the stilbene defenses of Norway spruce. Plant Physiol. 2013;162:1324–36.
Hammerbacher A, Raguschke B, Wright LP, Gershenzon J. Gallocatechin biosynthesis via a flavonoid 3′,5′-hydroxylase is a defense response in Norway spruce against infection by the bark beetle-associated sap-staining fungus Endoconidiophora polonica. Phytochemistry. 2018;148:78–86.
Hammerbacher A, Kandasamy D, Ullah C, Schmidt A, Wright LP, Gershenzon J. Flavanone-3-hydroxylase plays an important role in the biosynthesis of spruce phenolic defenses against bark beetles and their fungal associates. Front Plant Sci. 2019;10:208.
Hietala AM, Kvaalen H, Schmidt A, Jøhnk N, Solheim H, Fossdal G. Temporal and spatial profiles of chitinase expression by Norway spruce in response to bark colonization by Heterobasidion annosum. Appl Environ Microbiol. 2004;70:3948–53.
Hoch G, Körner C. Growth and carbon relations of tree line forming conifers at constant vs. variable low temperatures. J Ecol. 2009;97:57–66.
Hu J, Baker A, Bartel B, Linka N, Mullen RT, Reumann S, et al. Plant peroxisomes: biogenesis and function. Plant Cell. 2012;24:2279–303.
Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45:D353–61.
López de Heredia U, Vázquez-Poletti JL. RNA-seq analysis in forest tree species: bioinformatic problems and solutions. Tree Genet Genomes. 2016;12:30.
Jagodzik P, Tajdel-Zielinska M, Ciesla A, Marczak M, Ludwikow A. Mitogen-activated protein kinase cascades in plant hormone signaling. Front Plant Sci. 2018;9:1387.
Jalil SU, Mishra M, Ansari MI. Current view on chitinase for plant defence. Trends Biosci. 2015;8(24):6733–43.
Jamiesen MA, Trowbridge AM, Raffa KF, Lindroth RL. Consequences of climate warming and altered precipitation patterns for plant-insect and multitrophic interactions. Plant Physiol. 2012;160:1719–22.
Jøhnk N, Hietala AM, Fossdal CG, Collinge DB, Newman MA. Defense-related genes expressed in Norway spruce roots after infection with the root rot pathogen Ceratobasidium bicorne (anamorph : Rhizoctonia sp.). Tree Physiol. 2005;25:1533–43.
Kao Y-T, Gonzalez KL, Bartel B. Peroxisome function, biogenesis, and dynamics in plants. Plant Physiol. 2018;176(1):162–77.
Kautz M, Meddens AJH, Hall R, Arneth A. Biotic disturbances in northern hemisphere forests - a synthesis of recent data, uncertainties and implications for forest monitoring and modelling. Glob Ecol Biogeogr. 2017;26:533–52.
Kjaer ED, McKinney LV, Nielsen LR, Hansen LR, Hansen JK. Adaptive potential of ash (Fraxinus excelsior) populations against the novel emerging pathogen Hymenoscyphus pseudoalbidus. Evol Appl. 2012;5:219–28.
Kovalchuk A, Zeng Z, Ghimire RP, Kivimäenpää M, Raffaello T, Liu M, et al. Dual RNA-seq analysis provides new insights into interactions between Norway spruce and necrotrophic pathogen Heterobasidion annosum s.l. BMC Plant Biol. 2019;19(1):2.
Krajnc AU, Novak M, Felicijan M, Kraševec N, Lešnik M, Zupanec N, et al. Antioxidative response patterns of Norway spruce bark to low-density Ceratocystis polonica inoculation. Trees. 2014;28:1145–60.
Krokene P, Christiansen E, Solheim H, Franceschi VR, Berryman AA. Induced resistance to pathogenic fungi in Norway spruce. Plant Physiol. 1999;121(2):565–70.
Kumar M, Brar A, Yadav M, Chawade A, Vivekanand V, Pareek N. Chitinases—potential candidates for enhanced plant resistance towards fungal pathogens. Agriculture. 2018;8:88.
Kusumoto N, Zhao T, Swedjemark G, Ashitani T, Takahashi K, Borg-Karlson AK. Antifungal properties of terpenoids in Picea abies against Heterobasidion parviporum. For Path. 2014;44:353–61.
Lacerda AF, Vasconcelos ÉAR, Pelegrini PB, Grossi de Sa MF. Antifungal defensins and their role in plant defense. Front Microbiol. 2014;5:116.
Lam E. Controlled cell death, plant survival and development. Nat Rev Mol Cell Biol. 2004;5(4):305–15.
Lecourieux D, Ranjeva R, Pugin A. Calcium in plant defence-signalling pathways. New Phytol. 2006;171(2):249–69.
Lim GH, Singhal R, Kachroo A, Kachroo P. Fatty acid- and lipid-mediated signaling in plant defense. Annu Rev Phytopathol. 2017;55:505–36.
Liu JJ, Schoettle AW, Sniezko RA, Sturrock RN, Zamany A, Williams H, et al. Genetic mapping of Pinus flexilis major gene (Cr4) for resistance to white pine blister rust using transcriptome-based SNP genotyping. BMC Genomics. 2016;17:753–64.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods. 2001;25(4):402–8.
Lorrain C, Gonçalves dos Santos KC, Germain H, Hecker A, Duplessis S. Advances in understanding obligate biotrophy in rust fungi. New Phytol. 2019;222:1190–206.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Lundén K, Danielsson M, Durling MB, Ihrmark K, Nemesio Gorriz M, Stenlid J, et al. Transcriptional responses associated with virulence and defence in the interaction between Heterobasidion annosum s.s. and Norway spruce. PLoS One. 2015;10(7):e0131182.
Lieutier F, Brignolas F, Sauvard D, Yart A, Galet C, Brunet M, et al. Intra- and inter-provenance variability in phloem phenols of Picea abies and relationship to a bark beetle-associated fungus. Tree Physiol. 2003;23:247–56.
Matić S, Bagnaresi P, Biselli C, Orru L, Carneiro GA, Siciliano I, et al. Comparative transcriptome profiling of resistant and susceptible rice genotypes in response to the seedborne pathogen Fusarium fujikuroi. BMC Genomics. 2016;17:608.
Mayr S, Siller C, Kriss M, Oberhuber W, Bauer H. Photosynthesis in rust-infected adult Norway spruce in the field. New Phytol. 2001;151:683–9.
Mayr S, Schwienbacher F, Beikircher B, Dämon B. Damage in needle tissues after infection with Chrysomyxa rhododendri increases cuticular conductance of Picea abies in winter. Protoplasma. 2010a;243:137–43.
Mayr S, Schwienbacher F, Dämon B. Increased winter transpiration of Norway spruce needles after infection by the rust Chrysomyxa rhododendri. Protoplasma. 2010b;243:137–43.
Melotto M, Zhang L, Oblessuc PR, He SY. Stomatal defense a decade later. Plant Physiol. 2017;174(2):561–71.
Meng X, Xu J, He Y, Yang KY, Mordorski B, Liu Y, et al. Phosphorylation of an ERF transcription factor by Arabidopsis MPK3/MPK6 regulates plant defense gene induction and fungal resistance. Plant Cell. 2013;25:1126–42.
Meyer J, Berger DK, Christensen SA, Murray SL. RNA-Seq analysis of resistant and susceptible sub-tropical maize lines reveals a role for kauralexins in resistance to grey leaf spot disease, caused by Cercospora zeina. BMC Plant Biol. 2017;17(1):197.
Mitou G, Budak H, Gozuacik D. Techniques to study autophagy in plants. Int J Plant Genomics. 2009;451357:14. https://www.hindawi.com/journals/ijpg/2009/451357/.
Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35(Suppl 2):W182–5.
Nagy NE, Sikora K, Krokene P, Hietala AM, Solheim H, Fossdal CG. Using laser micro-dissection and qRT-PCR to analyze cell type-specific gene expression in Norway spruce phloem. PeerJ. 2014;2:e362.
Nagy NE, Fossdal CG, Krokene P, Krekling T, Lönneborg A, Solheim S. Induced responses to pathogen infection in Norway spruce phloem: changes in polyphenolic parenchyma cells, chalcone synthase transcript levels and peroxidase activity. Tree Physiol. 2004;24:505–15.
Nystedt B, Street NR, Wetterbom A, Zuccol A, Lin Y-C, Scofield DG, et al. The Norway spruce genome sequence and conifer genome evolution. Nature. 2013;497:579–84.
Peng Y, van Wersch R, Zhang Y. Convergent and divergent signaling in PAMP-triggered immunity and effector-triggered immunity. Mol Plant Microbe Interact. 2018;31(4):403–9.
Porth I, Hamberger B, White R, Ritland K. Defense mechanisms against herbivory in Picea: sequence evolution and expression regulation of gene family members in the phenylpropanoid pathway. BMC Genomics. 2011;12:608.
Plumed-Ferrer C, Väkeväinen K, Komulainen H, Rautiainen M, Smeds A, Raitanen J, et al. The antimicrobial effects of wood-associated polyphenols on food pathogens and spoilage organisms. Int J Food Microbiol. 2013;164:99–107.
Ramachandran SR, Yin C, Kud J, Tanaka K, Xiao F, Hulbert SH. Effectors from wheat rust fungi suppress multiple plant defense responses. Phytopathol. 2017;107:75–83.
Roy BA, Guswell S, Harte J. Response of plant pathogens and herbivores to a warming climate. Ecology. 2004;85:2570–81.
Schwachtje J, Baldwin IT. Why does herbivore attack reconfigure primary metabolism? Plant Physiol. 2008;146:845–51.
Sharma B, Joshi D, Bhatt TK. Role of ubiquitin-mediated degradation system in plant biology. Front Plant Sci. 2016;7:806.
Smith AM, Stitt M. Coordination of carbon supply and plant growth. Plant Cell Environ. 2007;30:1126–49.
Sniezko RA. Resistance breeding against nonnative pathogens in forest trees—current successes in North America. Can J Plant Pathol. 2006;28(Suppl 1):270–9.
Sundell D, Mannapperuma C, Netotea S, Delhomme N, Lin Y-C, Sjödin A, et al. The plant genome integrative explorer resource: PlantGenIE.org. New Phytol. 2015 Dec;208(4):1149–56.
Üstün S, Sheikh A, Gimenez-Ibanez S, Jones A, Ntoukakis V, Börnke F. The proteasome acts as a hub for plant immunity and is targeted by pseudomonas type III effectors. Plant Physiol. 2016;172:1941–58.
VanEtten HD, Mansfield JW, Bailey JW, Farmer E. Two classes of plant antibiotics: phytoalexins versus “phytoanticipins”. Plant Cell. 1994;6:1191–2.
Veluthakkal R, Sundari BKR, Dasgupta MG. Three chitinases – stress and developmental-driven gene regulation. Forest Pathol. 2012;42:271.
Vogt T. Phenylpropanoid biosynthesis. Mol Plant. 2010;3:2–20.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57–63.
Wang F, Zhang F, Chen M, Liu Z, Zhang Z, Fu J, et al. Comparative transcriptomics reveals differential gene expression related to Colletotrichum gloeosporioides resistance in the octoploid strawberry. Front Plant Sci. 2017;8:779.
Westermann AJ, Gorski SA, Vogel J. Dual RNA-seq of pathogen and host. Nat Rev Microbiol. 2012;10:618–30.
Witzell J, Martín JA. Phenolic metabolites in the resistance of northern forest trees to pathogens – past experiences and future prospects. Can J For Res. 2008;38:2711–27.
Xie F, Xiao P, Chen D, Xu L. Zhang B: miRDeepFinder: a miRNA analysis tool for deep sequencing of plant small RNAs. Plant Mol Biol. 2012;80:75–84.
Xing P, Zhang X, Bao Y, Wang Y, Wang H, Li X. Comparative transcriptome analyses of resistant and susceptible near-isogenic wheat lines following inoculation with Blumeria graminis f. sp. tritici. Int J Genomics. 2017;2017:7305684.
Zhu QH, Stephen S, Kazan K, Jin G, Fan L, Taylor J, et al. Characterization of the defense transcriptome responsive to Fusarium oxysporum-infection in Arabidopsis using RNA-seq. Gene. 2013;512:259–66.
Zottele F, Salvadori C, Corradini S, Andreis D, Wolynski A, Maresi G. Chrysomyxa rhododendri in Trentino: a first analysis of monitoring data. Balt For. 2014;20:28–36.
We thank the Landesforstdirektion Tirol, Landesforstgärten Tirol, and the Botanical Garden of the University Innsbruck for their assistance and support.
This study was supported by the Austrian Science Fund (FWF) P29896-B22. FWF had no role in the design of the study, collection, analysis, and interpretation of data, neither in writing the manuscript.
Ethics approval and consent to participate
Not applicable for this manuscript.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Inter-replicate correlation plots. Inter-replicate correlation of (a) sample R1 (Baseline control T1), (b) sample R4 (Control T4), (c) sample R9 (Control T9), (d) sample R21 (Control T21), (e) sample R39 (Control T39), (f) sample R4 (Infected T4), (g) sample R9 (Infected T9), (h) sample R21 (Infected T21), (i) sample R39 (Infected T39 no symptoms), (j) sample R39 (Infected T39 clear symptoms). (k) Pearson correlation coefficients.
Symptoms of C. rhododendri infection at 21 and 39 dpi. At 21 dpi, first symptoms of infection were detectable under the microscope, i.e. small blisters on the needle surface become visible (see arrows). At 39 dpi, several current-year-needles of infected trees showed the characteristic yellow discoloration and first aecio spore stocks were formed.
Differentially expressed genes. Transcripts differentially expressed in infected compared to control plants in needles at 4 dpi (A), 9 dpi (B), 21 dpi (C). At 39 dpi, infected symptomatic needles were compared to control plants (D), infected non-symptomatic needles were compared to control plants (E), and infected symptomatic needles were compared to infected non-symptomatic needles (F). Over- and under-expressed genes are highlighted in red and green respectively. (XLS 1218 kb)
KEGG Orthology (KO) assignment and proportion. KO assignment to the Picea abies differentially expressed genes during infection by Chrysomyxa rhododendri and proportion of assignments at each time point. At 39 dpi symptomatic needles were compared to controls. (XLS 284 kb)
KEGG pathways of differentially expressed genes by Picea abies infected by Chrysomyxa rhododendri. Metabolic pathway and orthology-oriented functional annotations of the protein sequences were performed against the KEGG database using the KEGG Automatic Annotation Server. KEGG Orthology (KO) assignment was applied using the Bi-directional Best Hit (BBH) method. Numbers refer to key enzymes represented by Picea abies KAAS assigned orthologous. At 39 dpi symptomatic needles were compared to controls. (XLS 47 kb)
KEGG Orthology (KO) detailed output. This analysis was applied at each time point for differentially over- and under-expressed genes. Resulting hierarchical text (Htext) was converted to tabulated excel file and KO assignments resumed on a Norway spruce pathway hierarchy model for each time point. At 39 dpi symptomatic needles were compared to controls. (XLS 962 kb)
List of putative defense related genes differentially expressed by Picea abies infected by Chrysomyxa rhododendri. Genes arranged according to the related KEGG pathways (several genes are shared by more than one pathway). Selection was done based on literature, therefore also pathways and genes not included in this list could be relevant for defense response (see Additional file 6: Table S4). For 39 dpi, only symptomatic needles were compared with controls. (XLS 696 kb)
KEGG pathways related to plant defence. A) Plant-pathogen interaction, B) MAPK signaling pathway-Plant, C) Plant hormone signal transduction, D) Phenylalanine, tyrosine and tryptophan biosynthesis, E) Phenylpropanoid biosynthesis, F) Stilbenoid, diarylheptanoid and gingerol biosynthesis, G) Flavonoid biosynthesis, H) Flavone and flavonol biosynthesis, I) Terpenoid backbone biosynthesis, J) Cutin, suberine and wax biosynthesis. At 39 dpi symptomatic needles were compared to controls.
Thirty nine dpi symptomatic vs. non-symptomatic (S-NS) differential expression analysis. A) DEGs overlapped between 39 dpi symptomatic vs. control (S-C) and 39 dpi symptomatic vs. non-symptomatic (S-NS), B) Gene ontology (GO) term enrichment analysis for symptomatic vs. non-symptomatic (S-NS) DEGs, C) KEGG Orthology assignment for symptomatic vs. non-symptomatic (S-NS) DEGs, D) KEGG pathways for symptomatic vs. non-symptomatic (S-NS) DEGs, E) KEGG Orthology detailed output for symptomatic vs. non-symptomatic (S-NS) over-expressed genes, GF) KEGG Orthology detailed output for symptomatic vs. non-symptomatic (S-NS) under-expressed genes. (XLS 61 kb)
RT-qPCR validation of selected DEGs. Relative expression changes obtained by RT-qPCR (red bars; mean fold changes ± SE, n = 3) and RNA-seq (blue bars) compared to control plants at 4 dpi, 9 dpi, 21 dpi and 39 dpi (S-C). Symptomatic needles at 39 dpi were additionally compared to non-symptomatic needles (S-NS). Asterisk means significant (p < 0.05) difference compared to non-infected needles.
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.
Command S1. STAR command. The detailed STAR settings used for mapping and genome index generation.
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)
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)
RT-qPCR assay details. Primer sequences of RT-qPCR assays.
RT-qPCR data. Mean efficiency-corrected Cq values obtained by RT-qPCR for (A) experimental samples and (B) determination of reference gene stability. (XLS 40 kb)
About this article
Cite this article
Trujillo-Moya, C., Ganthaler, A., Stöggl, W. et al. RNA-Seq and secondary metabolite analyses reveal a putative defence-transcriptome in Norway spruce (Picea abies) against needle bladder rust (Chrysomyxa rhododendri) infection. BMC Genomics 21, 336 (2020). https://doi.org/10.1186/s12864-020-6587-z
- Forest tree
- Fungal infection
- Phenolic compounds
- PR proteins
- RNA sequencing