Drought stress described by transcriptional responses of Picea abies under pathogen Heterobasidion parviporum attack

The major threats to the sustainable supply of forest tree products are adverse climate, pests and diseases. Climate change, exemplified by increased drought, poses a unique threat to global forest health. This is attributed to the unpredictable behavior of forest pathosystems, which can favor fungal pathogens over the host under persistent drought stress conditions in the future. Currently, the effects of drought on tree resistance against pathogens are hypothetical, thus research is needed to identify these correlations. Norway spruce (Picea abies) is one of the most economically important tree species in Europe, and is considered highly vulnerable to changes in climate. Dedicated experiments to investigate how disturbances will affect the Norway spruce Heterobasidion sp. pathosystem are important, in order to develop different strategies to limit the spread of H. annosum s.l. under the predicted climate change. Here, we report a transcriptional study to compare Norway spruce gene expressions to evaluate the effects of water availability and the infection of Heterobasidion parviporum. We performed inoculation studies of three-year-old saplings in a greenhouse (purchased from a nursery). Norway spruce saplings were treated in either high (+) or low (-) water groups: high water group received double the water amount than the low water group. RNA was extracted and sequenced. Similarly, we quantified gene expression levels of candidate genes in biotic stress and jasmonic acid (JA) signaling pathways using qRT-PCR, through which we discovered a unique preferential defense response of H. parviporum-infected Norway spruce under drought stress at the molecular level. Disturbances related to water availability, especially low water conditions can have negative effects on the tree host and benefit the infection ability of the pathogens in the host. From our RNA-seq analysis, 114 differentially expressed gene regions were identified between high (+) and low (-) water groups under pathogen attack. None of these gene pathways were identified to be differentially expressed from both non-treated and mock-control treatments between high (+) and low (-) water groups. Finally, only four genes were found to be associated with drought in all treatments.


Introduction
Norway spruce (Picea abies) is one of the most economically important forestry tree species in Europe. Currently, massive inputs of anthropogenic greenhouse gases (i.e. CO2, N2O, CH4) into the atmosphere have resulted in increasing atmospheric temperatures in an effect known as 'global warming' [1]. This effect has several consequences; one of which, climate change is a potential driver in influencing forest health. Under global warming, extreme climate events, e.g. drought, are expected to increase in frequency, duration and intensity [2]. In fact, this has started to take its toll on Norway spruce that is sensitive to abiotic disturbances. Under climate change, there is an increase in droughtassociated stress in this particular tree species, rendering them more susceptible to threats such as pests and pathogens, which can compromise their overall tree health. This phenomenon can be explained by the disease triangle model that shows the relationship between three factors: (1) a susceptible plant host, (2) a pathogen that causes the disease in the plant host, and (3) an environment that favors the pathogen [6,7]. An environment that weakens the defense capability of the plant host and supports the growth and spread of the pathogen can greatly influence the success rate for a pathogen to infect its targeted host. This means that changes in climate do not limit the new challenges of Norway spruce to adjust only to abiotic disturbances, as the projected change in climate will favor certain pathogens in forests [2,5,[8][9][10]. The change in tree resistance under global warming against fungal pathogens still remain hypothetical, thus research is needed to identify these correlations. Discriminating responses at the genomic level between biotic and abiotic stressors and understanding the potential trade-offs in plant gene expression will allow for more mechanistic predictions about future scenarios driven by complex shifts in climate change [11,12].
Heterobasidion parviporum is one of the native fungal pathogen species (the other two being H. annosum and H. abietinum) belonging to a root-rot fungus, H. annosum sensu lato (s.l.) complex, which is considered to be the most devastating Norway spruce pathogen in Europe [13]. Being necrotrophic in nature, it mainly derives its nutrients by killing its tree host and feeding off its contents. There are two modes in which H. parviporum can infect Norway spruce: (1) through deposition of basidiospores on fresh cut stumps, or (2) via root contacts [14]. Detection of H. parviporum in Norway spruce is often challenging, as the infected tree might not necessarily show signs of symptoms (stringy white rot, necrosis in sapwood, presence of fruiting body) over decades. This further promotes the spread of H. parviporum to neighboring trees via root contacts. In addition, H. parviporum mainly colonizes the heartwood of Norway spruce, causing decay and even tree mortality [15].
The induction of a disease by a pathogen in a tree host is often linked to the environmental conditions in the ecosystem as well. In this particular Norway spruce -H. parviporum forestry pathosystem, climate plays an important role in promoting the spread of H. parviporum in several ways: (1) short winters that provide a longer period of hyphae growth and spore transmission, (2) warmer summers that encourages fungal hyphae to grow faster and produce more spores, and (3) windy days that can promote the dispersion of basidiospores [2,16,17]. Increased frequency and intensity of drought-associated stress in Norway spruce is expected in the following years [18]. Considering how Norway spruce is a valuable tree species to the boreal regions, it is imperative to study the interactions between this tree species and its associated fungal pathogens. This way, appropriate and effective measures can be implemented to mitigate this inevitable climate change.
Rapid advancements in molecular research have provided insights into the Norway spruce -H. parviporum forestry pathosystem through various studies such as pathogen characterization and detection [19], constructing gene expression libraries [20] and screening of resistant trees based on genetic components [21]. Norway spruce resistance against H. annosum s.l. is a quantitative trait [21][22][23] and is associated by several known genes (defense responses) with variation [21,24]. Signaling molecules are a group of plant secondary metabolites with a critical role in defense processes [25]. The plant hormones, especially salicylic acid (SA) and jasmonic acid (JA), help mediate information about the attack beyond the point of invasion [25]. By altering the hormonal balance in plants, pathogens use the defensive machinery of plants to their advantage and either induce or suppress the processes relevant for cell death and accumulation of antimicrobial compounds [25,26]. Jasmonic acid (JA)-mediated signaling is the prioritized module in the Norway spruce defense-signaling network against H. parviporum [27][28][29]. Similarly, beside JA, the ethylene (ET) signaling pathway probably also plays a central role in defense response of Norway spruce against H. parviporum [28]. In addition, salicylic acid (SA) signaling pathway can, likewise, have an important role in this pathosystem as the SA-mediated hypersensitive response (HR) can be facilitated by necrotrophic pathogen infection [30]. Understanding the similarities and differences in transcriptional response between variable environmental conditions can thus provide useful information over how Norway spruce -H. parviporum pathosystems adapt to climate change.
We could show previously that H. parviporum benefit from host stress (i.e. drought) [18]. In this study, we performed RNA-seq analysis on samples derived from the same trees as in   [18]. We concentrated on Norway spruce, to find out which genetic pathways are involved in the specific responses to unfavorable environmental stress (drought) with and without the pathogen, H. parviporum. In addition, we also built a transcriptional study based on previous studies [28,[31][32][33][34] to investigate the interplay between H. parviporum infection and defense-related gene expression mediated by drought stress.

Plant and Fungal Material
The setting used in this experiment is from Terhonen et al. 2019. Summarily, plant material consisted of 18 three-year-old, apparently healthy and vital, Norway spruce saplings purchased from the nursery Schlegel and Co Gartenprodukte (provenience of the saplings: No. 840 11, Thüringer Wald and Frankenwald, montane zone 600 m). Seedlings were potted into 3-liter plastic pots filled with fertilized peat (Flora gard, TKS®2 Instant Plus, Hermann Meyer KG, Rellingen, Germany). The potted saplings were then acclimatized to the greenhouse conditions for 16-days prior to the water experiment, during which they received tap water, as required, to maintain moist soil. No additional fertilization was given during the experiment.
H. parviporum strain was received from the North-West German Forest Research Institute, collected by G. Langer and colleagues. H. parviporum (strain NW-FVA 0459) was isolated from P. abies in 2010 (Germany, Lower Saxony, forest department Oerrel, Bobenwald). DNA of H. parviporum was extracted from 150 mg of the homogenized mycelium sample using the method developed by Keriö et al. (2020) and the fungal strain used in this study was tested for the specificity with species-specific primers developed for H. parviporum (KJ-F and KJ-R; [19]). In brief, DNA template (100ng), buffer (KCl extra buffer, 1X), 1.5 mM MgCl2, primers KJ-F and KJ-R (each concentration of 0.5 µM), a dNTP-mix (each deoxynucleotide in a concentration of 200 µM) and 20 U/mL of DNA-polymerase (VWR) was adjusted to 25 µL reaction with autoclaved Milli-Q H2O. The cycling conditions used were: Initial denaturation 10 min at 95 °C, followed by a 3-step cycling: Denaturation 30 s at 95 °C; annealing 35 s at 67 °C; extension 1 min at 72 °C, for 40 cycles. A final extension step for 7 min at 72 °C was applied. Taq DNA polymerase (Qiagen) was used for PCR amplification of ITS regions with the primer pair, ITS1-F and ITS4. Briefly, the PCR protocol was as follows: 1X CoralLoad PCR Buffer, 200 µM dNTP, 0.5 µM primer 1, 0.5 µM primer 2, 100ng template DNA, 0.2 U/µL DNA polymerase; the reaction was adjusted to 25 µL with autoclaved Milli-Q H2O. The PCR conditions used were 94 °C for 3 min; 30 cycles of 94 °C for 30 s, 55 °C for 1 min, 72 °C for 1 min; a final step at 72 °C for 10 min was applied before storing the PCR products at 4 °C. Possible contaminations were determined with a negative control using sterile water as a template in both PCR protocols. StainIN™ RED Nucleic Acid Stain was used to confirm DNA amplicons on a 2% agarose gel and the visual detection was made by ultraviolet transillumination. The expected band (fragment size around 350 bp) with KJ-F and KJ-R primers was visualized and the ITS region PCR products were purified and sequenced using the ITS4 primer at Microsynth SEQLAB (Göttingen, Germany). The FASTA files thus obtained were checked with BioEdit to confirm that the pathogen was not contaminated with other fungi strains. ITS region sequencing indicated no contamination in the culture of H. parviporum (see [18]).

Experimental Design
The following experiment is described in detail in [18]. In summary, the study was conducted at the Forest Botany and Tree Physiology greenhouses, in Göttingen, Germany (51°33'28.4" N 9°57'30.5" E) from early April until early August 2018. From this experiment, we chose 18 seedlings (three per water treatment group) from randomly block-assigned to waterproof tables with either high water (high water group = +group) or low water (low water group = −group) availability treatments, that were experimentally inoculated with H. parviporum, mock-inoculated controls or left entirely untreated. The water treatment experiment was running for 35 days before the inoculations were performed. H. parviporum isolate (NW-FVA 0459) was inoculated using a puncher (5mm) through the bark to reach the sapwood surface. Equal sized plugs from pure culture of H. parviporum or control (2% MEA) were placed onto the exposed surface and sealed with Parafilm ® [18]. The inoculation experiment was run for 70 days [18]. At the end of the experiment, we collected 30 samples aiming to sequence 18 samples. The stems of Norway spruce saplings (5cm both directions: above and below the infection point was collected, around 11 cm) were sampled into liquid nitrogen and stored individually into 15-ml Falcon tubes at -80 °C for long-term storage.

RNA extraction
A total of 30 Norway spruce stems were selected, comprising of five sample replicates for each of the six treatment groups for small-scale totRNA extractions; a total of 24 (with four sample replicates for each of the six treatment groups) were selected for upscale totRNA extractions. Only 18 (three per treatment group, 3*6=18, were sent for sequencing). For the first round of extractions, the surface of inoculation point (bark, phloem, sapwood) was scraped with a scalpel and the sample was ground into fine powder (Mixer Mill MM 400 from Retsch GmbH with a set program of 25.0 Hz for 20 s). The samples were handled with liquid nitrogen throughout the entire milling process. The ground product was then stored in 15-ml Falcon tubes at -80 °C for later use in the extractions. The extraction protocols used were modified from Zeng et al. (2018) [35] and unfortunately, the first try out with 0.1-0.5 g of sample transferred to 2-ml Eppendorf tubes was unsuccessful. We subsequently re-did the scraping from the infection point and increased the area to collect more material (the sample included bark, phloem and sapwood). The new samples were similarly ground into fine powder and extracted for totRNA, of which they were successful. From this, we establish the following extraction protocol, which we had modified from Zeng et al. (2018) [35] for our study.
Day 1: 2-3 g of sample was transferred to a 50-ml Falcon tube, where 15 ml of Extraction buffer (2% CTAB, 2% PVP, 100 mM Tris-HCl, 25 mM EDTA, 2 M NaCl, DEPC-treated H2O, pre-heated to 65 °C) was added. 100 µl of β-mercaptoethanol was added to each tube and vigorously vortexed before incubating the tubes at 65 °C for 5 min. An equal volume of chloroform:isoamyl alcohol (24:1) was added to each tube and vortexed vigorously before centrifuging the tubes at 10,000 g for 15 min. The upper aqueous phase was transferred to a new sterile 50-ml Falcon tube. The addition of chloroform:isoamyl alcohol, then centrifuging the tubes and subsequently transferring the upper aqueous phase was repeated. 1/4 volume of 10 M LiCl was added to each tube and inverted 8 times before precipitating the RNA overnight at 4 °C Day 2: The tubes were centrifuged using the benchtop Eppendorf 5810R Centrifuge (Rotor A-4-81 with bucket for 7 x 50 ml conical tubes) with a speed of approximately 6480 rcf at 4 °C for 30 min. The supernatant was collected in a new 50-ml Falcon tube and set aside for gDNA extraction. The remaining pellet (containing totRNA) was briefly dried in a laminar flow hood. 1 ml of SSTE (1 M NaCl, 0.5% SDS, 10 mM Tris-HCl, 1 mM EDTA, DEPC-treated H2O, preheated to 65 °C) was then added to dissolve the pellet. Thereafter, an equal volume of chloroform:isoamyl alcohol was added and vortexed vigorously before centrifuging at 10,000 g for 15 min. The upper aqueous phase (maximum 500 µl) was transferred to a new 2-ml Eppendorf tube and three volumes of cold 100% absolute EtOH were added to each tube. The resultant mixture is incubated at -20 °C for 1-2 hours before transferring to -80 °C overnight for RNA precipitation.
Day 3: totRNA samples are centrifuged at maximum centrifugal speed at 4 °C for 30 min. The supernatant is pipetted out and the remaining totRNA pellet is dried in a sterile hood for approximately 3 min. The pellet is then resuspended in 100 µl of DEPC-treated H2O, where it can be subsequently used in molecular analyses or stored at -80 °C for longterm storage.

Quantitation and purity analyses of extracted totRNA
A NanoPhotometer (NP80) from Implen GmbH was used to quantitate totRNA. A Qubit 3.0 Fluorometer from Life Technologies (Catalogue number: Q33216) was used to quantitate totRNA. The analysis was performed according to the manufacturer's instructions. Three biological replicates from each treatment (2 µg of total RNA) were sent for RNA sequencing (RNAseq) to Novogene. RNA integrity and quantitation were assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Purified RNA was used for library construction using rRNA depletion method (GlobinZero kit) for Illumina® (NEB, USA), followed by Illumina paired-end (2 × 150 bp) RNAseq in the facilities of the Novogene (HK) company limited. The raw reads have been submitted to NCBI SRA under Bioproject PRJNA761217.

RNA-seq analysis
Quality control was made by Novogene through in-house scripts. Briefly, clean data (clean reads) were obtained by removing reads containing adapter and poly-N sequences and reads with low quality from raw data. The filtering process was as follows: (1) Discard reads with adapter contamination; (2) Discard reads when uncertain nucleotides constitute more than 10% of either read (N > 10%); (3) Discard reads when low quality nucleotides (base quality less than 20) constitute more than 50% of the read. At the same time, Q20, Q30 and GC content of the clean data were calculated. All downstream analyses were based on the clean data with high quality. The quality was assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reference genome and gene model annotation files were downloaded from genome website browser (the processed reads were mapped against the genome assembly of Norway spruce (v 1.0) [36] downloaded from ftp://plantgenie.org/Data/ConGenIE/Picea_abies/) directly. Paired-end clean reads were mapped to the reference genome using HISAT2 software [37]. The raw read counts table was loaded into R studio version 3.5.1 [38] and differential expression analysis between two conditions (three biological replicates per condition) was performed using DESeq2 R package [39]. Conditions were: normal water versus low water in each treatment (treatments: non-treated, mock inoculated or H. parviporum infection). To track down only drought related gene pathways, we compared all treatments to non-treated (+) water treatment. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the False Discovery Rate (FDR). Genes with an adjusted P value < 0.05 and the cutoff value for FDR < 0.05 were assigned as differentially expressed. We searched for statistically expressed genes between different conditions/treatments with Venny 2.1 (Oliveros 2007-2015). Results were compared on genes expressed in the study by Chaudhary et al. (2021). When we found interesting genes, we visualized the expression in each treatment. Similarly, we studied how H. parviporum resistance candidate gene laccase PaLAC5 (MA_97119g0010 and MA_97119g0020) were expressed in each treatment [24].

GO enrichment test
The GO functions and enrichment analysis was performed using ConGenIE.org (http://congenie.org/). Annotate GeneList workflow was created on the ConGenIE platform selecting P. abies as the species of interest. The gene ID's were pasted in the search field and annotated. The enrichment analysis was carried out using the functional enrichment tools. The GO slims was run and based on the annotations the transcripts were categorized in the biological process and molecular function groups with their respective descriptions. The significantly enriched GO terms were then categorized into Biological Process, Molecular Function and Cellular Component.

Transcriptional responses validated with Quantitative real-time polymerase chain reaction (qRT-PCR)
For cDNA preparation we followed the protocol described by Raffaello and Asiegbu (2017) [40]. Summarily, 1 µg sample of total RNA was treated with DNase I (Thermo Scientific). The mixture was incubated at 37 °C for 30 min and then 65 °C for 10 min to inactivate the DNase. Before the second incubation, 1 µl of 50 mM EDTA (Thermo Scientific) was added to terminate the reaction. After that, 1 µl random hexamers primers (0.2 µg, Thermo Scientific) and 0.5 µl nuclease-free water were added to the RNA sample. The above mixture was heated at 65 °C for 5 min and cooled down on ice immediately. Then, a 7.5 µl mixture containing 4 µl 5× RevertAid reverse transcriptase buffer (Thermo Scientific), 2 µl dNTPs 27 (10 mM), 0.5 µl RiboLock RNase inhibitor (Thermo Scientific) and 1 µl RevertAid reverse transcriptase (200U/µl, Thermo Scientific), was added and the reaction was incubated at 25 °C for 10 min, then at 42 °C for 60 min. Finally, the reaction was stopped by incubating at 70 °C for 5 min. Synthesized cDNA was purified using GenElute PCR Clean-Up Kit (Sigma-Aldrich, Finland) referring to instructions provided by the manufacturer.
In this experiment, we have prepared three biological replicates from each of the six treatment groups to make six pooled cDNA mixes for the qRT-PCR. We have also prepared three technical replicates per treatment per gene tested. Four genes were tested: LOX [27], ERF1 [28], CHIIV [34] and p/DIR32 [33]. A housekeeping gene, Elongation factor 1α (ELF1α) [31], was used as a reference against the tested genes. To prepare the loading samples for qRT-PCR, gene mastermixes comprising of 0.5 µM each of forward and reverse primers for each selected gene, qPCR Supermix (SsoFast™ EvaGreen® Supermix (BIO-RAD)), and DEPC-H2O as a diluent were made. A total of 20 µl per qRT-PCR reaction (including the pooled cDNA) were then loaded on to a 96-well V-bottom plate and sealed tight with plate foil. The qRT-PCR (iQ5™ Multicolor Real-Time PCR Detection System (BIO-RAD)) cycling conditions ran were set for an initial 2 min at 95 °C, then 40 cycles of 20 s at 95 °C, 30 s at 60 °C, and 30 s at 72 °. According to a primer efficiency test prior to running the qRT-PCR, we decided on a cDNA concentration of 100 ng per pooled cDNA. The relative expressions were subsequently calculated with the method 2-ΔΔCT [41].

RNA-seq
Our analysis of RNA-seq data identified 114 transcripts with significantly different expression between high (+) and low (-) water groups under pathogen attack (Supplementary file 1).
Comparing H. parviporum infection under drought stress, the biological processes (in upregulated) significantly enriched were metabolic process, transport, protein folding, and translation. Molecular function GO terms in pathogen-infected seedlings were significant in transferase activity, hydrolase activity, and binding. Cellular component GO term was significantly related to plant cell organelles (apoplast, chloroplast, Golgi apparatus, mitochondrion, cell wall, plasma membrane) (Supplementary file 1). The roles played by these genes are expressed as GO terms.
The drought increased substantially the gene expression under H. parviporum inoculation (Table 1). Exclusively 377 genes were downregulated and 129 genes upregulated only in H. parviporum-infected (-) group (Supplementary file 1). The upregulated Gene Ontology (GO) terms assigned to the "Biological process" in the (-) water group are: protein folding, response to stress, oxidation-reduction process, rRNA processing, primary metabolic process, phosphorylation, pentose-phosphate shunt and Golgi organization. The terms assigned to the "Molecular function" are: ATP binding, oxidoreductase activity, nucleotide binding, unfolded protein binding, transferase activity and protein kinase activity. The downregulated genes expressed resulted in different gene functions and the most abundant within this group are metabolic process, oxidation-reduction process, response to wounding, response to fungus, phosphorylation, defense response to fungus, response to water deprivation, and response to oxidative stress. Only the Heterobasisioninfected (-) groups were found to be functionally enriched (Supplementary file 1). Table 1. Statistically differentially expressed genes of each treatment when compared to non-treated (+) water group We found four genes that were expressed significantly only under drought stress MA_189802g0010, MA_10265000g0010, MA_10435878g0010 and MA_178006g0010 (Figure 1, Table 2). These genes were not observed to be significantly expressed in normal water treatments. Only gene DEG MA_10435878g0010 had a role in the metabolic process (GO:0008152) ( Table 2).

qRT-PCR
Normalization between technical triplicate samples was achieved by determining the mean fold change of each gene (i.e. LOX, ERF1, CHIIV, p/DIR32) against the housekeeping gene, ELF1, for each treatment scheme ( Table 3). The difference in gene expression between technical triplicates were mostly with a mean fold change value (2 ∆∆ ) of 1 irrespective of water treatment groups and inoculation status, indicating little non-biological variation. This is further illustrated with Figure 4, which generally demonstrated consistency between triplicates except for two treatment groups in p/DIR32 (SD: H. parviporum, normal = 0.788; mock control, low = 0.822) and one group in ERF1 (SD: mock-control, low = 0.565).  To obtain meaningful interpretation of our results, we also compared the fold change in gene expression mediated by H. parviporum and physical wounding as singular and combinatorial factors to examine how Norway spruce respond differently under drought conditions (Table 4).  On the other hand, with the exception of CHIIV, the tested genes were downregulated when the plants received normal water treatment (LOX = 0.27, ERF1 = 0.87, p/DIR32 = 0.08). No change was observed for ERF1 even at low water treatment condition (ERF1 = 1.06). These observations are further represented in Figure 5.
In general, LOX and p/DIR32 showed an increase in gene expression levels under drought stress when combinatorically influenced by H. parviporum infection and physical wounding in comparison to their respective singular factors that are independent of each other. Under normal water conditions, ERF1 are dominantly expressed in H. parviporuminfected Norway spruce; this is the same for CHIIV in Norway spruce that are physically wounded. This observation is flipped under low water conditions ( Figure 5).   [42] identified several candidate genes related to Heterobasidion sp. resistance from which we found only four (MA_110169g0010, MA_14707g0010, MA_15852g0010, MA_10427673g0020) in H. parviporum infected plants to be statistically expressed. The accumulation of these genes was induced in drought treatment in almost all cases. The same trend could be found for Palac5 genes (MA_97119g0010 and MA_97119g0020) [24] as they accumulated more in H. parviporum infected plants under drought, and in nontreated (-) group plants these genes were at the minimum level. These results highlight the need for combined factor studies: abiotic and biotic stress in plant gene expression studies to resolve this combination effect as some pathways can be induced by the environment.

Changing environment: drought as a new factor in forest pathosystems
Climatic extremes such as drought affect trees directly, by inducing water deficiency, and indirectly, by making trees more susceptible to fungal pathogens [10,18]. This leads to unknown behavior of host-fungi interactions, resulting in increased disease outbreaks caused by native fungal pathogens [43]. Prolonged water and nutrient deficiencies are causing severe damages in Norway spruce and the drought disturbances are expected to increase in following years [2]. Norway spruce provenances react differentially to drought [44]. Trujillo-Moya et al. (2018) [44] found drought response indicators (24 genes, 29 SNPs) of Norway spruce genotypes. One of these 24 genes, MA_588952g0010 (no similarity to any known gene), was found exclusively to be significantly expressed in H. parviporum (-) treatment. Our analysis of RNA-seq data identified only four genes (MA_189802g0010, MA_10265000g0010, MA_10435878g0010 and MA_178006g0010) significantly expressed in drought stressed trees. The gene (MA_178006g0010) is a member of the BURP domain family (RD22-like). Members of this family have been described to be associated with several plant species subjected to stress conditions [45]. In Arabidopsis thaliana, the AtRD22 and AtUSPL1 (both members of the Arabidopsis thaliana BURP) are upregulated as part of the abscisic acid (ABA)-mediated moisture stress response [45]. ABA is a phytohormone that adjusts the drought stress response in plants. Additional ABA application enhances tolerance of wheat [46] and tomato [47] seedlings to drought and ABA is accumulating under drought in maize [48].
Gene MA_10265000g0010 is highly homologous to intracellular pathogenesis-related PR10 proteins. These genes produce proteins that are toxic to invading fungal pathogens and act in defense signaling. PR10 gene down-regulation was observed in maize early in drought-induced plants [49]. PR-10 protein was also noted to accumulate in needles of drought-stressed Pinus pinaster seedlings [50]. The sequence of gene MA_10265000g0010 has a 100% query coverage and 98.97% identity to Picea glauca putative intracellular pathogenesis-related protein (picg4) [51]. Ekramoddoullah (2004) [51] showed that the expression of PR-10 (including picg4) genes is induced by both biotic (white pine blister rust, Cronartium ribicola) and abiotic (cold) stress. Therefore, we suggest that the gene MA_10265000g0010 could be a possible marker of abiotic and biotic factors (i.e. drought stress and pathogen attack) for transcriptional response studies in conifers [51] as the expression seems to be accumulating in the combination of these stressors.
Gene MA_10435878g0010 plays a role in metabolism and it belongs to the 2OG-Fe(II) oxygenase superfamily. It has been found to be significantly expressed in biosynthesis of secondary metabolites and in metabolic pathways under different light spectra in Norway spruce [52]. It appears that this gene reacts to environmental change and is accumulating more when we combine biotic stress (H. parviporum) with abiotic stress (drought). MA_189802g0010 belongs to the thaumatin family and it responds to fungal organisms and salinity; the association of this gene with P. abies was previously unknown.
In addition to the RNA-seq analysis, another aim of this study was to uncover preliminary insights into the defense responses of H. parviporum-infected Norway spruce under drought stress. LOX [31], ERF1 [28], CHIIV [34], and p/DIR32 [33] have previously been demonstrated to be activated and expressed when Norway spruce is infected with H. parviporum; thus far, our observations from LOX and ERF1 through qRT-PCR concur with these studies. On the other hand, although we had similar observations [33] from p/DIR32 showing greater expression in wounded Norway spruce than H. parviporum-infected ones, they were only induced under drought stress conditions. This similar observation was also the same in CHIIV [34]. Plausible reasons for these findings are unfortunately unknown to date. Nevertheless, we discover that this defense-related gene response is more complex when additional abiotic stress factors are present. Summarily, we deduce that water availability plays a critical role in the gene switch between up-/downregulation strategies employed by H. parviporum-infected Norway spruce. This is evident when p/DIR32 regulation is heavily influenced by water supply. Although CHIIV is more responsive to the presence of physical wounding under normal water conditions, this response can be easily overtaken by H. parviporum infection when subjected to drought stress. This suggests a preferential defense response of Norway spruce to H. parviporum infection over physical wounding, depending on water availability conditions. Likewise, this preferential response can also be influenced by the type of damage and infection present in the plant. For example, when affected by physical wounding, CHIIV is expressed under normal water conditions while ERF1 is expressed under drought stress conditions; this is the other way around when Norway spruce is infected with H. parviporum. Through these observations, the complexity of defense-related gene response in H. parviporum-infected Norway spruce is demonstrated, and piques new research questions: within the Norway spruce -H. parviporum forest pathosystem, how does the host make its choice in activating a specific defense-related gene, and how does it relate back to the various abiotic signaling pathways involved? Currently, due to the explorative nature of our study, the RNA-seq and qRT-PCR dataset we have gathered thus far is, unfortunately, not sufficiently comprehensive to make sensible inferences to the above-mentioned questions. Therefore, opportunities for further in-depth research in this direction may be prompted.

Norway spruce resistance in a changing climate
Norway spruce is not only highly economic but also possesses ecological value in Europe, particularly in Northern and Central European countries. The current breeding activities concentrate on developing more productive and climate resilient trees in Scandinavia. Similarly, several markers are suggested for Heterobasidion annosum s.l. resistance [24,42,53,54]. These traits should all be studied for suitability in breeding programs as due to the long tree generations, Norway spruce might not keep up with the changes in the local environment, making them more susceptible for infection by the Heterobasidion species. In addition, with the new discovery of preferential defense gene response in Norway spruce mediated by water availability, understanding how specific defense-related genes are regulated in the midst of a changing climate will be pivotal to support successful breeding of these trees. The expression of these gene regions would be good to be tested in multifactorial analysis. Necrosis caused by fungal pathogens are increasing under host stress [3,10,18], and revealing the molecular determinants of resistance/virulence in a changing environment of forest pathosystems would be needed to mitigate this detrimental response. The major threats to the sustainable supply of forest tree products are adverse climate that benefit the pests and diseases at the cost of the host. Changes in the environment claim for the inclusion of new traits in breeding programs. Increasing the availability of Norway spruce and H. parviporum transcriptome data thus creates new opportunities to discover genes that underlie this tree-pathogen interaction on the projected climate change continuum.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1 Author Contributions: For research articles with several authors, a short paragraph specifying their individual contributions must be provided. The following statements should be used "Conceptualization, X.X. and Y.Y.; methodology, X.X.; software, X.X.; validation, X.X., Y.Y. and Z.Z.; formal analysis, X.X.; investigation, X.X.; resources, X.X.; data curation, X.X.; writing-original draft preparation, X.X.; writing-review and editing, X.X.; visualization, X.X.; supervision, X.X.; project administration, X.X.; funding acquisition, Y.Y. All authors have read and agreed to the published version of the manuscript." Please turn to the CRediT taxonomy for the term explanation. Authorship must be limited to those who have contributed substantially to the work reported.
Funding: This research was funded by the Faculty of Forest Sciences and Forest Ecology, Georg-August-Universität Göttingen, Germany.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data is contained within the article, the supplementary material and at the NCBI SRA under Bioproject PRJNA761217.