Genome-Wide Exon-Capture Approach Identifies Genetic Variants of Norway Spruce Genes Associated With Susceptibility to Heterobasidion parviporum Infection

Root and butt rot caused by members of the Heterobasidion annosum species complex is the most economically important disease of conifer trees in boreal forests. Wood decay in the infected trees dramatically decreases their value and causes considerable losses to forest owners. Trees vary in their susceptibility to Heterobasidion infection, but the genetic determinants underlying the variation in the susceptibility are not well-understood. We performed the identification of Norway spruce genes associated with the resistance to Heterobasidion parviporum infection using genome-wide exon-capture approach. Sixty-four clonal Norway spruce lines were phenotyped, and their responses to H. parviporum inoculation were determined by lesion length measurements. Afterwards, the spruce lines were genotyped by targeted resequencing and identification of genetic variants (SNPs). Genome-wide association analysis identified 10 SNPs located within 8 genes as significantly associated with the larger necrotic lesions in response to H. parviporum inoculation. The genetic variants identified in our analysis are potential marker candidates for future screening programs aiming at the differentiation of disease-susceptible and resistant trees.

Introduction

Norway spruce [Picea abies (L.) H. Karst.] is one of the dominant species in European boreal forests. As a major source of timber, it is one of the most commercially important species in Nordic countries. However, a combination of biotic and abiotic factors pose a threat to sustainable timber supply from Norway spruce plantations (Lindroth and St. Clair, 2013; Oliva et al., 2013; Robert et al., 2013; Lind et al., 2014; Wang et al., 2015). Among microbial pathogens of Norway spruce, necrotrophic fungi of the Heterobasidion annosum species complex occupy an outstanding position due to their devastating effect on commercial plantations (Keriö et al., 2014). Members of this species complex are causative agents of root and butt rot of conifer trees throughout the boreal and temperate zones of Northern hemisphere (Abu et al., 2004; Asiegbu et al., 2005a; Garbelotto and Gonthier, 2013; Gunulf et al., 2013; Skrøppa et al., 2015). They affect growth rate, mortality, and timber quality of conifer trees as well as steadily increase the risk of windfall (Gori et al., 2013; Keriö et al., 2014). The annual losses due to Heterobasidion infections are estimated at €50 million for Finnish forest industry (Gori et al., 2013; Piri and Hamberg, 2015), and at €790 million for European forest owners (Asiegbu et al., 2005a). Based on morphology, ecology, and molecular characters, H. annosum species complex is classified into five species, three of which occur in Europe: H. annosum sensu stricto, H. parviporum, and H. abietinum (Asiegbu et al., 2005a; Garbelotto and Gonthier, 2013; Hansson et al., 2014). Their preferred hosts are Scots pine, Norway spruce, and Silver fir, respectively (Asiegbu et al., 2005a). H. parviporum occurs primarily in forest areas where its major host, Norway spruce, grows naturally (Asiegbu et al., 2005a,b; Lind et al., 2014). The tree infection might occur either through wounds and damaged area on stem or, alternatively, pathogen might move from infected to healthy trees via root contacts (Asiegbu et al., 2005a).

Norway spruce trees possess a repertoire of defenses against microbial pathogens and pests. These defense reactions are commonly classified into physical and chemical defenses. An example of physical defense mechanisms is provided by thick outer bark, which acts as an efficient barrier that most pathogens cannot penetrate. It protects both dead tissues (heartwood) and living tissues (phloem, cambium and sapwood) (Kovalchuk et al., 2013). Chemical defense of spruce trees is achieved via production of diverse secondary metabolites such as flavonoids, lignans, phenolics, stilbenes, and terpenes. Many of them possess antimicrobial activities against potential pathogens, and they could delay or prevent the establishment of pathogens within the tree (Danielsson et al., 2011). Upon recognition of potential pathogens or tissue damage caused by microorganisms or insects, trees activate induced defenses that include cell wall reinforcement, production of phenolics, terpenes and pathogenesis-related (PR) proteins, and formation of traumatic resin ducts (TD), and polyphenolic parenchyma (PP) cells (Arnerup et al., 2013).

Growth rate and wood properties (Cunningham et al., 2006; Zubizarreta-Gerendiain et al., 2009; Luoranen and Viiri, 2016; Levkoev et al., 2017) and impact of climate change (Kolár et al., 2017) have been the major economic traits selected for in Norway spruce breeding programs. Recent studies have also aimed at identifying genes implicated in Norway spruce resistance against H. parviporum infection (Lind et al., 2014). Their results led to the identification of the gene PaLAR3, encoding leucoanthocyanidin reductase, an enzyme of the flavonoid biosynthetic pathway. A low frequency allele of this gene is associated with higher resistance to H. parviporum and increased content of (+)-catechin, a compound showing a fungistatic effect on H. parviporum (Nemesio Gorriz et al., 2016). However, Norway spruce resistance to the infection is a quantitative character, and it is likely to be controlled by multiple loci (Lind et al., 2014). Thus, further studies are required to identify additional genetic determinants of resistance against H. parviporum infection in Norway spruce.

The genome-based discoveries open the way to better understanding the effects of host genes associated with resistance to microbial pathogens. Availability of the draft sequence of Norway spruce genome (Nystedt et al., 2013) provides opportunities for high-throughput identification of genetic variants associated with the resistance to H. parviporum infection. Sequence capture or exon capture is a cost-efficient approach to determine genetic variation in species with a large genome size or without a reference genome having a wide range of scale of the capture from several targeted loci to a million target regions (Grover et al., 2012; Müller et al., 2015). This approach eliminates the need to set up and optimize thousands of PCR reactions for candidate genes, allowing instead for the parallel enrichment and sequencing of thousands of targeted regions. Exon capture has been successfully applied for such species as poplar, eucalyptus, polyploid wheat, polyploid switchgrass, and loblolly pine (Pavy et al., 2016). In the present study, we used the exon capture approach for the identification of Norway spruce genetic loci associated with the resistance or susceptibility to H. parviporum infection. The objectives of this study were:

1) to assess the genetic polymorphism among selected clonal lines of Norway spruce and

2) to identify genetic variants associated with the control of lesion lengths caused by H. parviporum inoculation.

Materials and Methods

Plant and Fungal Material

The study material comprised 533 Norway spruce saplings. These represented 64 clonal lines which are a sample of the spruce breeding material for Central Finland (7 to 10 ramets per clonal line) (Table 1). The ortet trees were selected in two 10-year-old Norway spruce progeny trials in Pieksämäki (N 62°23′, E 27°17′) for their favorable phenotype (mainly for growth and late flushing). The ortets were propagated by rooted cuttings, which were then grown three seasons at the breeding nursery of the Natural Research Institute Finland. All the clones were untested and previously unexploited in research. Prior to the inoculation experiment, they were cultivated in a greenhouse (12 h daylight, day 22°C; night 18°C) at the Viikki campus of the University of Helsinki for 2 months.

TABLE 1

Table 1. Detailed information on the studies clones.

A H. parviporum isolate 04009 (a courtesy of Kari Korhonen, LUKE) was grown on 2% malt extract agar (MEA) plates supplemented with Norway spruce sawdust. The plates were inoculated with H. parviporum and incubated at 20°C for 4 weeks. For the mock inoculation, the sterile 2% MEA plates with 2% sawdust were used as a control.

Phenotyping for Susceptibility to H. parviporum Infection

After the 2-month acclimatization, 4 to 7 saplings of each clone were inoculated with H. parviporum, and three saplings of each clone were mock-inoculated to perform as a control. The inoculation method used (Figure 1) was similar to the previously described one (Swedjemark and Stenlid, 1995, 1996; Sun et al., 2011; Keriö et al., 2014). The stem surface was sterilized with 70% ethanol before a hole was made through the bark with a 70% ethanol-sterilized puncher (3 mm diameter) to reach the xylem surface. The distance of the holes from the stem base was ~5 cm. Plugs of H. parviporum-colonized sawdust-containing malt extract agar of the same diameter as the holes were placed in the holes and then covered with Parafilm (Pechiney Plastic Packaging, Chicago, IL, USA). For control saplings, holes were inoculated with uncolonized sawdust-containing malt extract agar plugs. Saplings were harvested 2 months after inoculation. Harvested material was stored at −20°C until used. The necrotic lesions in phloem and xylem tissues were measured after removing the periderm tissues using a sterilized knife. The saplings dimensions (stem diameter, height, and volume) were also documented.

FIGURE 1

Figure 1. (A) illustration of inoculation method used in Norway spruce seedlings experiment. (B) Inoculated seedling stem with sawdust-agar pre-colonized by pathogen H. parviporum (nr. 04009). (C) Covering the inoculum with Parafilm. (D) Measurement horizontal and vertical necrosis on phloem (E) and xylem.

Statistical Analysis of Phenotyping Data

Data were analyzed using the SPSS version 25.0 (IBM Corporation, New York, USA) and Microsoft Excel 2016. A one-way analysis of variance (ANOVA) was performed to detect differences among the necrosis sizes of each clone and among the inoculation and mock treatments. The normality of the error variances was checked with Levene's test and by skewness and kurtosis values for the data distribution, the values for lesion size were converted to the logarithmic scale (Supplementary Files S1, S2). The effect of the clone on the lesion size was estimated with the following general linear model:

$\begin{array}{ll}{Y}_{ij}=\mu +c_i+e_{ij}& \left(1\right)\end{array}$

where Y_ij is the value of ijth observation, μ is the overall mean, c_i is the fixed effect of the clone i, and e_ij is the random residual term associated with the ith clone and jth ramet. Differences between clones were tested with the least significant difference (LSD) test between the means.

The growth rates of H. parviporum were assessed by t-tests for two-sample assuming unequal variances. Pearson's correlation analysis was used to determine correlations of diameter, height, and volume of the tested clonal lines with the sizes of necrotic lesions caused by H. parviporum. Differences were considered as statistically significant if p-value was below the threshold of 0.05.

Estimation of Heritability

Broad-sense heritability was estimated both for the growth traits (height, diameter, and volume of the saplings) and the lengths of various lesions using the following linear model:

$\begin{array}{ll}{\text{y}}_{\text{ij}}=\mu +c_{\text{i}}+{\text{w}}_{\text{ij}}& \left(2\right)\end{array}$

where y_ij is the ijth phenotypic observation of the trait y, μ is the general mean, c_i is the random effect ith clone and w_ij is the residual effect due to the jth sapling of the ith clone. The variance components of the model effects ${\sigma }_{\text{c}}^2$ and ${\sigma }_{\text{w}}^2$ were computed using procedure Mixed of the SAS package. Heritability was first calculated separately for both treatments (fungal vs. mock inoculation) as a ratio of the clonal variance to the phenotypic variance:

$\begin{array}{ll}{\text{h}}^2={\text{r }{\sigma }^2}_{\text{c}}{({\sigma }_{\text{c}}^2+{\sigma }_{\text{w}}^2)}^{-1}& \left(3\right)\end{array}$

and its standard error as sqrt(Var (${\sigma }_{\text{c}}^2$))/(${\sigma }_{\text{c}}^2$ + ${\sigma }_{\text{w}}^2$)⁻¹.

Secondly, to analyze the magnitude of clone by treatment interaction in the lesion traits we used the mixed effect model:

$\begin{array}{ll}{\text{y}}_{\text{ijk}}=\mu +{\text{c}}_{\text{i}}+{\text{T}}_{\text{j}}+{\text{c}}_{\text{i}}{\text{T}}_{\text{j}}+{\text{w}}_{\text{ijk}}& \left(4\right)\end{array}$

where T_j is the fixed effect of the jth treatment (fungal inoculation or mock inoculation), and c_iT_j is the random interaction effect between ith clone and jth treatment. This enabled us to separate the clone and clone by environment (treatment) interaction components of variance in the estimation of broad-sense heritability:

$\begin{array}{ll}{\text{h}}^2={\sigma }_{\text{c}}^2{({\sigma }_{\text{c}}^2+{\sigma }_{\text{cT}}^2+{\sigma }_{\text{w}}^2)}^{-1}& \left(5\right)\end{array}$

The relative magnitude of the interaction variance vs. the total genetic variance (R_cT) was calculated as a simple ratio:

$\begin{array}{l}{\text{R}}_{\text{cT}}={\sigma }_{\text{cT}}^2{({\sigma }_{\text{c}}^2+{\sigma }_{\text{cT}}^2)}^{-1}\end{array}$

Library Preparation and Target Enrichment Sequencing

Total genomic DNA was extracted from the needles using a standard cetyl–trimethyl ammonium bromide (CTAB) method with modifications described in Terhonen et al. (2011). DNA concentrations were estimated with PicoGreen dsDNA quantification assay (ThermoFisher Scientific, USA) and DNA integrity was analyzed by visualizing the DNA on a 0.8% w/v agarose electrophoresis gel. A total of 500 ng was utilized for construction of libraries compatible with Illumina sequencing machines. In summary, the DNA was mechanically sheared to a mean fragment size of 300 bp, followed by repair of the ends of the molecules, phosphorylation and adenylation. Common adapters suited for Illumina sequencing were ligated on each side of the molecules containing 8 bp indexes (i5 and i7). The ligated libraries were amplified with 10 cycles of PCR with common primers to enrich the libraries for properly ligated molecules and the resulting libraries were quantified with PicoGreen. A set of 80,000 probes was utilized to capture genic regions of the Norway spruce genome (Vidalis et al., in review). In summary, these probes target 34761 genes with a median of two probes per gene, and were selected based on their uniqueness on the genome. Probe sequences were designed based on the P. abies v1.0 genome (Nystedt et al., 2013) and RNA-based, 120 nucleotide long probes were synthesized by Agilent Technologies. Sixteen libraries were equimolarly combined to a total of 750 ng for target enrichment, which was performed following SureSelect Target Enrichment System (Agilent Technologies). Enriched libraries were sequenced using Illumina HiSeq3000 instrument on a 2 × 100 bp sequencing mode. The library construction and enrichment sequencing were performed by RAPiD Genomics (USA). The generated sequence reads were submitted to the GenBank Short Reads Archive (SRA) and are available under the BioProject accession number PRJNA450911.

Variant Identification and Genome-Wide Association Analysis

The generated sequencing data was processed to identify the SNPs segregating in the population. Short reads were filtered for quality and aligned to the scaffolds that contain probes of the P. abies v1.0 genome using Mosaik version 2.1 as described before (Neves et al., 2013). The resulting alignment BAM files were used to identify SNPs with Freebayes v1.0.2 (Garrison and Marth, 2012) on a populational level (parameters: –max-complex-gap 1 –theta 0.01 –no-complex –no-mnps –no-indels). The resulting SNPs were filtered using VCF.filter v0.1.13 (Müller et al., 2017) to keep markers with a minimum quality of 10, an average sequencing depth between 15 and 150, a minimum allele frequency of 0.01, a missing data lower than 0.4 and assigning genotypes with sequencing depth lower than 3 as missing. The VCF file containing the information on the identified SNPs was deposited to the Dryad repository (Mukrimin et al., 2018).

The filtered markers for the 64 clonal lines were subsequently used for genome-wide association analysis (GWAS). The average xylem lesion size for each clone was used as a quantitative phenotype for the GWAS, which was carried out using EMMAX (Kang et al., 2010) including the BN kinship matrix in the model to control for relatedness in the population. To control for multiple testing, the Bonferroni and the Benjamini and Hochberg (1995) procedures were used (α = 0.05). The Manhattan and Quantile-Quantile plot (Q-Q plot) were calculated using the qqman package on R (Turner SD. Qqman: an R Package For Visualizing Gwas Results Using Q-Q And Manhattan Plots. 2014; https://doi.org/10.1101/005165). Since the spruce genome is not ordered and fragmented, the p-values plotted on the Manhattan plot do not represent any physical or genetic distance, but rather each SNP was simply considered a coordinate on the x-axis for visualization purposes (Supplementary File S3).

Results

Variation in Lesion Length in Response to Fungal Inoculation

In our experiment, the lesions formed in response to inoculation with H. parviporum were significantly larger than the lesions developing after mock treatment. The differences between mock and fungal inoculation treatments were observed both in xylem and in phloem, and they were statistically significant (p < 0.001) in both cases (Table 2).

TABLE 2

Table 2. Size of necrotic lesions in spruce saplings (mean ± SE) for inoculation and mock treatments.

The lesions developed faster in vertical than in horizontal direction, as indicated by significant differences between vertical and horizontal lesion lengths in all treatments [e.g., F_{(1, 680)} = 96.58, F_{(1, 680)} = 547.82, F_{(1, 380)} = 130.52, and F_{(1, 381)} = 402, respectively at p < 0.001]. Horizontal lesion length in phloem was significantly higher than in xylem (p < 0.001), but there was no significant differences between the tissues in the vertical lesion length.

The analyzed spruce clones displayed significant differences in size of necrotic lesions both in phloem and in xylem [F_{(63, 277)} = 2.84 and F_{(63, 277)} = 2.77, respectively at p < 0.001]. The distribution of lesion length among the analyzed clones was similar to normal (Figure 2). The data for the individual clones are summarized in Figure 3. Clonal lines V31094 and V31406 had the largest lesion sizes and were classified as the most susceptible to the infection.

FIGURE 2

Figure 2. LogNormal distributions of vertical necrotic stem lesion size in phloem (A, C) and xylem (B, D) based on number of clones and individual saplings, respectively.

FIGURE 3

Figure 3. Clone means of lesion size in (A) phloem vertical necrosis; (B) phloem horizontal necrosis; (C) xylem vertical necrosis; and (D) xylem horizontal necrosis.

Based on Pearson correlation coefficients (r), there were significant negative correlations between the growth traits (stem diameter, seedling height, and volume) and the observed lesion sizes. The strongest negative correlation was found in inoculated xylem that differed significantly at p < 0.01 (Table 3).

TABLE 3

Table 3. Pearson correlation coefficients (r) according to both phloem and xylem lesion size to growth traits.

The heritability estimates, indicating the proportion of clonal variation of the total phenotypic variation in each treatment, were moderate to large, ranging from 0.18 to 0.80. The levels of heritability were consistently larger for the clones subjected to mock treatment than the respective estimates for inoculated trees. The genetic coefficients of variation displayed the similar tendency (Table 4). The clonal variance used in calculating these heritability values was in fact composed of two components of variation: the true genetic (clonal) and genotype by treatment interaction variances. When these two sources of variation were distinguished by Equation (4), much smaller heritability estimates were obtained both for the lesion traits (0.15–0.23) and the growth traits (0.07–0.13) (Table 5). The relative magnitude of the genotype by treatment interaction variance (Equation 5) is generally low, but markedly high for phloem horizontal necrosis.

TABLE 4

Table 4. Estimates of among-clone and within-clone (residual) variance components and their standard errors.

TABLE 5

Table 5. Clone, clone by treatment and residual variance components, broad-sense heritability (h2) and the relative magnitude of the clone by treatment variance (RcT).

Sequence Capture and Sequencing of Gene Space in Norway Spruce

We genotyped 64 clonal lines of Norway spruce by resequencing of predicted exons corresponding to 34 761 gene models. On average, 6.3 million reads were generated per clonal line. Reads were aligned to Norway spruce reference genome, with the median experiment-wide sequencing depth of 40.35×. The efficiency of the sequence enrichment strategy was high, capturing 79 170 out of 80 000 probes. Calling of genetic variants identified 373 384 high-quality biallelic SNPs.

Identification of Sequence Variants Associated With the Necrotic Lesion Size

Genome-wide association analysis identified 34 genes (36 SNPs) that were significantly associated (at 5% FDR significance level) with the size of necrotic lesions formed in response to fungal inoculation (Supplementary File S4). As we had a limited number of clonal lines available for GWAS analysis, we additionally selected only those variants, which occurred in at least 4 clonal lines. The resulting set included 10 SNPs located within 8 gene models (Table 6). All of the identified variants were associated with larger necrotic lesions and belonged to low frequency alleles (frequency among studied clonal lines <5%). Three of the identified genes were of particular interest because of their potential involvement in the regulation of plant defense responses: the gene MA_940838g0010 with a similarity to Arabidospis thaliana ILITYHIA, which is required for the establishment of systemic acquired resistance (Monaghan and Li, 2010), the gene MA_4047g0010 encoding a predicted subtilase (subtilisin-like serine peptidase), and the gene MA_3905g0010 with a similarity to the plant-specific class of HD2-type histone deacetylases.

TABLE 6

Table 6. SNPs identified by GWAS as significantly associated with the size of necrotic lesions formed in response to fungal inoculation.

Discussion

In this study, we characterized a set of Norway spruce clonal lines to explore the genetic control of susceptibility to H. parviporum infection. As an indicator of susceptibility, we used size of necrotic lesions developing in response to fungal inoculation. The artificial fungal inoculation of seedlings and saplings has been widely practiced under controlled greenhouse conditions to investigate tree susceptibility to Heterobasidion infection (Swedjemark and Stenlid, 1996; Swedjemark et al., 2001; Sun et al., 2011; Krokene et al., 2012). This approach is much faster than the field experiments (Swedjemark and Stenlid, 1996), but yet much more complex than a potential selection via genetic markers. Inoculation with Heterobasidion results in the pathogen spread in host tissues and in the development of necrotic lesions around the inoculation point.

In this study, the necrotic lesions significantly differed among the tested clones, which concurs with the results of previous studies (Swedjemark and Stenlid, 1996; Swedjemark et al., 2001; Keriö et al., 2014; Skrøppa et al., 2015). Our findings show that saplings with lower values of growth parameters (stem diameter, height, and volume) appeared considerably more sensitive to fungal inoculation and displayed larger necrotic lesions in response to H. parviporum infection. This is in accordance with previous results that reported a strong negative correlation between growth traits and H. parviporum resistance (Karlsson et al., 2008). The investigated clones represented a practically random sample of genotypes from a single breeding population aimed at a specific region (Central Finland). We could not detect any common variables characterizing these clones that would explain the differences observed in their susceptibility.

Heritability is often used in quantitative genetics, behavior genetics, selective breeding, among other fields to quantify the ratio of genetic to phenotypic variation in a specific population and a specific environment. In this study, all the traits showed moderate to high heritability, suggesting a relatively small role of environmental differences in the trait expression when the analyses were performed separately for the two treatments. Interestingly, heritability estimates for traits measured after the mock treatment were consistently higher than the respective estimates on clones subjected to inoculation, suggesting these traits could be indeed considered to be genetically different. In the second analysis for the combined data (including a fixed treatment effect and a random clone by treatment variance) all the traits studied showed moderate heritability. The relative magnitude of the interaction was high for a single trait, phloem horizontal necrosis, but small to moderate for the rest of the lesion traits and very small for all of the three growth traits. This kind of discrepancy is interesting, yet difficult to explain. The finding needs to be verified in further studies to make sure that it was not just a random occurrence.

GWAS became an important tool in the identification of genetic loci controlling economically important traits in various plant species (Kump et al., 2011; Tian et al., 2011; Li et al., 2013; Wei et al., 2015). Our analysis identified 10 SNPs associated with the control of the size of necrotic lesions developed in response to H. parviporum inoculation. These SNPs correspond to 8 genes. None of the identified genes was previously reported as being implicated in control of resistance of Norway spruce to Heterobasidion infection. At the same time, none of the SNP markers associated with the Norway spruce resistance to H. parviporum in a previous screening program (Lind et al., 2014) was recovered in our analysis. Possible explanations for the lack of overlap of our results with the previous ones could be differences in the employed marker identification strategies; the populations used, whereas in the case of our work a multi-pedigree breeding population was used vs. a single bi-parental QTL cross used in the referenced study; and the relatively small number of samples employed in both studies.

Our analysis identified several genes, which might be involved in the regulation of plant defense responses in Norway spruce. The identification of a gene showing similarity to A. thaliana ILITYHIA (ILA) in our study is particularly noteworthy, as the corresponding protein was demonstrated to play a key role in the induced defense reactions in Arabidopsis (Monaghan and Li, 2010). Furthermore, by controlling chloroplast development, it might be implicated in the regulation of the ROS accumulation and programmed cell death in response to a pathogen attack (Faus et al., 2018). However, nothing is known about the function of this protein in other plant species, including conifer trees. The role of this gene in defense reactions of Norway spruce against H. parviporum infection deserves further studies.

The gene MA_4047g0010 encodes a predicted subtilase (subtilisin-like serine peptidase). This protein family is largely expanded in plant genomes, and their physiological functions include the control of growth and development, regulation of programmed cell death and senescence and plant responses to biotic and abiotic stressors (Schaller et al., 2018). Several plant subtilases are important modulators of plant immune responses (Tornero et al., 1996; Ramírez et al., 2013; Serrano et al., 2016). A specific group of subtilases, called phytaspases, are involved in the regulation of programmed cell death (hypersensitive response) triggered by the pathogen attack (Chichkova et al., 2004; Schaller et al., 2018). Our transcriptomic profiling of asymptomatic and symptomatic Norway spruce trees naturally infected by Heterobasidion sp. showed that the gene MA_4047g0010 had higher expression level in symptomatic trees, with the p value (p = 0.06) close to the significance threshold (Supplementary File S5). However, taking into account a large number of the predicted subtilase genes in spruce genome, it is difficult to predict the specific role of the identified gene, and its characterization would require further experimental work.

The MA_3905g0010, identified in our analysis, shows similarity to the plant-specific class of HD2 histone deacetylases. Members of this class play a role, among others, in ABA and abiotic stress responses (Sridha and Wu, 2006; Luo et al., 2012; Han et al., 2016). They might potentially contribute to the ABA-mediated regulation of programmed cell death, which, in turn, might be one of the factors determining sensitivity or resistance to necrotrophic pathogens (Coll et al., 2011).

The connection of the remaining genes with the responses of spruce trees to fungal infection are less clear. No obvious link between their deduced function and the response of spruce trees against fungal infection could be established. However, it is possible that identified markers are physically linked with neighboring genes, located in the adjacent chromosomal loci. Unfortunately, the draft genome assembly of Norway spruce is unordered and highly fragmented (N50 = 4,869 bp; Nystedt et al., 2013), which makes it not always possible to find genes located adjacent to the identified SNP markers.

The identified genetic variants represent low frequency alleles (occurrence frequency among analyzed clonal lines <5%). This observation indicates that more extensive genotyping programs involving several hundreds of individuals might be required to identify the most relevant alleles controlling responses of Norway spruce to H. parviporum infection. In parallel, genomic selection approaches can be used to advance the selection of tolerant trees in the breeding program, as these methods do not rely on the dissection of quantitative traits to rank the most favorable trees.

This study is the first attempt to identify genes involved in the control of Norway spruce resistance to H. annosum infection using sequence capture enrichment strategy. Obtained results illustrate the suitability of this approach for genotyping of Norway spruce, a species with a very large genome and a high content of repetitive elements. Some of the identified genes are promising candidates, and further analysis should investigate their role in spruce defense reactions. However, our study also illustrates that even larger efforts might be required to identify genetic variants controlling economically important traits and occurring in low frequency in natural populations.

Author Contributions

Phenotyping of Norway spruce saplings was performed by MM and EJ; total genomic DNA extraction was done by MM and EJ; genetic and genome-wide association analyses were done by LN; heritability analysis was carried out by MH; MM, AK, and LN analyzed the data; MM and AK prepared the manuscript draft; MK and FA conceived the study and contributed to the experiment design. All authors read and approved the final version of the manuscript.

Funding

This study was supported by funding from the Academy of Finland (Academy of Finland grant nr. 276862).

Conflict of Interest Statement

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

Acknowledgments

We greatly acknowledge the Ministry of Research, Technology, and Higher Education (RISTEKDIKTI) of the Republic of Indonesia for the financial support of MM's Ph.D. project.

Supplementary Material

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

References

Abu, S. M., Li, G., and Asiegbu, F. O. (2004). Identification of Heterobasidion annosum (S-type) genes expressed during initial stages of conidiospore germination and under varying culture conditions. FEMS Microbiol. Lett. 233, 205–213. doi: 10.1016/j.femsle.2004.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnerup, J., Nemesio-Gorriz, M., Lundén, K., Asiegbu, F. O., Stenlid, J., and Elfstrand, M. (2013). The primary module in Norway spruce defence signalling against H. annosum s.l. seems to be jasmonate-mediated signalling without antagonism of salicylate-mediated signalling. Planta 237, 1037–1045. doi: 10.1007/s00425-012-1822-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Asiegbu, F. O., Adomas, A., and Stenlid, J. (2005a). Conifer root and butt rot caused by Heterobasidion annosum (Fr.) Bref. s.l. Mol. Plant Pathol. 6, 395–409. doi: 10.1111/j.1364-3703.2005.00295.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Asiegbu, F. O., Nahalkova, J., and Li, G. (2005b). Pathogen-inducible cDNAs from the interaction of the root rot fungus Heterobasidion annosum with Scots pine (Pinus sylvestris L.). Plant Sci. 168, 365–372. doi: 10.1016/j.plantsci.2004.08.010

CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate - a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300.

Google Scholar

Chichkova, N. V., Kim, S. H., Titova, E. S., Kalkum, M., Morozov, V. S., Rubtsov, Y. P., et al. (2004). A plant caspase-like protease activated during the hypersensitive response. Plant Cell 16, 157–171. doi: 10.1105/tpc.017889

PubMed Abstract | CrossRef Full Text

Coll, N. S., Epple, P., and Dangl, J. L. (2011). Programmed cell death in the plant immune system. Cell Death Differ. 18:1247. doi: 10.1038/cdd.2011.37

PubMed Abstract | CrossRef Full Text | Google Scholar

Cunningham, C., Zimmermann, N. E., Stoeckli, V., and Bugmann, H. (2006). Growth response of Norway spruce saplings in two forest gaps in the Swiss Alps to artificial browsing, infection with black snow mold, and competition by ground vegetation. Can. J. For. Res. 36, 2782–2793. doi: 10.1139/x06-156

CrossRef Full Text | Google Scholar

Danielsson, M., Lundén, K., Elfstrand, M., Hu, J., Zhao, T., Arnerup, J., et al. (2011). Chemical and transcriptional responses of Norway spruce genotypes with different susceptibility to Heterobasidion spp. infection. BMC Plant Biol. 11:154. doi: 10.1186/1471-2229-11-154

PubMed Abstract | CrossRef Full Text | Google Scholar

Faus, I., Niñoles, R., Kesari, V., Llabata, P., Tam, E., Nebauer, S. G., et al. (2018). Arabidopsis ILITHYIA protein is necessary for proper chloroplast biogenesis and root development independent of eIF2α phosphorylation. J. Plant Physiol. 224–225, 173–182. doi: 10.1016/j.jplph.2018.04.003

CrossRef Full Text | Google Scholar

Garbelotto, M., and Gonthier, P. (2013). Biology, epidemiology, and control of Heterobasidion species worldwide. Annu. Rev. Phytopathol. 51, 39–59. doi: 10.1146/annurev-phyto-082712-102225

PubMed Abstract | CrossRef Full Text | Google Scholar

Garrison, E., and Marth, G. (2012). Haplotype-based variant detection from short-read sequencing. arXiv preprint arXiv:1207.3907.

Google Scholar

Gori, Y., Cherubini, P., Camin, F., and La Porta, N. (2013). Fungal root pathogen (Heterobasidion parviporum) increases drought stress in Norway spruce stand at low elevation in the Alps. Eur. J. For. Res. 132, 607–619. doi: 10.1007/s10342-013-0698-x

CrossRef Full Text | Google Scholar

Grover, C. E., Salmon, A., and Wendel, J. F. (2012). Targeted sequence capture as a powerful tool for evolutionary analysis. Am. J. Bot. 99, 312–319. doi: 10.3732/ajb.1100323

PubMed Abstract | CrossRef Full Text | Google Scholar

Gunulf, A., Wang, L., Englund, J. E., and Rönnberg, J. (2013). Secondary spread of Heterobasidion parviporum from small Norway spruce stumps to adjacent trees. For. Ecol. Manage. 287, 1–8. doi: 10.1016/j.foreco.2012.09.011

CrossRef Full Text | Google Scholar

Han, Z., Yu, H., Zhao, Z., Hunter, D., Luo, X., Duan, J., et al. (2016). AtHD2D gene plays a role in plant growth, development, and response to abiotic stresses in Arabidopsis thaliana. Front. Plant Sci. 7:310. doi: 10.3389/fpls.2016.00310

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansson, D., Wubshet, S., Olson, Å., Karlsson, M., Staerk, D., and Broberg, A. (2014). Secondary metabolite comparison of the species within the Heterobasidion annosum s.l. complex. Phytochemistry 108, 243–251. doi: 10.1016/j.phytochem.2014.08.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, H. M., Sul, J. H., Service, S. K., Zaitlen, N. A., Kong, S. Y., Freimer, N. B., et al. (2010). Variance component model to account for sample structure in genome-wide association studies. Nat. Genet. 42:348. doi: 10.1038/ng.548

PubMed Abstract | CrossRef Full Text | Google Scholar

Karlsson, B., Tsopelas, P., Zamponi, L., Capretti, P., Soulioti, N., and Swedjemark, G. (2008). Susceptibility to Heterobasidion parviporum in Picea abies clones grown in different environments. For. Pathol. 38, 83–89. doi: 10.1111/j.1439-0329.2008.00543.x

CrossRef Full Text | Google Scholar

Keriö, S., Niemi, S. M., Haapanen, M., Daniel, G., and Asiegbu, F. O. (2014). Infection of Picea abies clones with a homokaryotic isolate of Heterobasidion parviporum under field conditions. Can. J. For. Res. 45, 227–235. doi: 10.1139/cjfr-2014-0247

CrossRef Full Text | Google Scholar

Kolár, T., Cermák, P., Trnka, M., Žid, T., and Rybníček, M. (2017). Temporal changes in the climate sensitivity of Norway spruce and European beech along an elevation gradient in Central Europe. Agric. For. Meteorol. 239, 24–33. doi: 10.1016/j.agrformet.2017.02.028

CrossRef Full Text | Google Scholar

Kovalchuk, A., Keriö, S., Oghenekaro, A. O., Jaber, E., Raffaello, T., and Asiegbu, F. O. (2013). Antimicrobial defenses and resistance in forest trees: challenges and perspectives in a genomic era. Annu. Rev. Phytopathol. 51, 221–244. doi: 10.1146/annurev-phyto-082712-102307

PubMed Abstract | CrossRef Full Text | Google Scholar

Krokene, P., Lahr, E., Dalen, L. S., Skrøppa, T., and Solheim, H. (2012). Effect of phenology on susceptibility of Norway spruce (Picea abies) to fungal pathogens. Plant Pathol. 61, 57–62. doi: 10.1111/j.1365-3059.2011.02487.x

CrossRef Full Text | Google Scholar

Kump, K. L., Bradbury, P. J., Wisser, R. J., Buckler, E. S., Belcher, A. R., Oropeza-Rosas, M. A., et al. (2011). Genome-wide association study of quantitative resistance to southern leaf blight in the maize nested association mapping population. Nat. Genet. 43, 163–168. doi: 10.1038/ng.747

PubMed Abstract | CrossRef Full Text | Google Scholar

Levkoev, E., Kilpeläinen, A., Luostarinen, K., Pulkkinen, P., Mehtätalo, L., Ikonen, V., et al. (2017). Differences in growth and wood density in clones and provenance hybrid clones of Norway spruce. Can. J. For. Res. 47, 389–399. doi: 10.1139/cjfr-2016-0285

CrossRef Full Text | Google Scholar

Li, H., Peng, Z., Yang, X., Wang, W., Fu, J., Wang, J., et al. (2013). Genome-wide association study dissects the genetic architecture of oil biosynthesis in maize kernels. Nat. Genet. 45, 43–50. doi: 10.1038/ng.2484

PubMed Abstract | CrossRef Full Text | Google Scholar

Lind, M., Källman, T., Chen, J., Ma, X. F., Bousquet, J., Morgante, M., et al. (2014). A Picea abies linkage map based on SNP markers identifies QTLS for four aspects of resistance to Heterobasidion parviporum infection. PLoS ONE 9:e101049. doi: 10.1371/journal.pone.0101049

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindroth, R. L., and St. Clair, S. B. (2013). Adaptations of quaking aspen (Populus tremuloides Michx.) for defense against herbivores. For. Ecol. Manage. 299, 14–21. doi: 10.1016/j.foreco.2012.11.018

CrossRef Full Text | Google Scholar

Luo, M., Wang, Y. Y., Liu, X., Yang, S., Lu, Q., Cui, Y., et al. (2012). HD2C interacts with HDA6 and is involved in ABA and salt stress response in Arabidopsis. J. Exp. Bot. 63, 3297–3306. doi: 10.1093/jxb/ers059

PubMed Abstract | CrossRef Full Text | Google Scholar

Luoranen, J., and Viiri, H. (2016). Deep planting decreases risk of drought damage and increases growth of Norway spruce container seedlings. New For. 47, 701–714. doi: 10.1007/s11056-016-9539-3

CrossRef Full Text | Google Scholar

Monaghan, J., and Li, X. (2010). The HEAT repeat protein ILITYHIA is required for plant immunity. Plant Cell Physiol. 51, 742–753. doi: 10.1093/pcp/pcq038

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukrimin, M., Kovalchuk, A., Neves, L. G., Jaber, E., Haapanen, M., Kirst, M., et al. (2018). Data from: genome-wide exon-capture approach identifies genetic variants of Norway spruce genes associated with susceptibility to Heterobasidion parviporum infection. Dryad Digit. Repos. doi: 10.3389/fpls.2018.00793. [Epub ahead of print].

CrossRef Full Text

Müller, H., Jimenez-Heredia, R., Krolo, A., Hirschmugl, T., Dmytrus, J., Boztug, K., et al. (2017). VCF. Filter: interactive prioritization of disease-linked genetic variants from sequencing data. Nucleic Acids Res. 45, W567–W572. doi: 10.1093/nar/gkx425

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, T., Freund, F., Wildhagen, H., and Schmid, K. J. (2015). Targeted re-sequencing of five Douglas-fir provenances reveals population structure and putative target genes of positive selection. Tree Genet. Genomes 11:816. doi: 10.1007/s11295-014-0816-z

CrossRef Full Text | Google Scholar

Nemesio Gorriz, M., Hammerbacher, A., Ihrmark, K., Källman, T., Olson, Å., Lascoux, M., et al. (2016). Different alleles of a gene encoding leucoanthocyanidin reductase (PaLAR3) influence resistance against the fungus Heterobasidion parviporum in Picea abies. Plant Physiol. 171, 2671–2681. doi: 10.1104/pp.16.00685

PubMed Abstract | CrossRef Full Text | Google Scholar

Neves, L. G., Davis, J. M., Barbazuk, W. B., and Kirst, M. (2013). Whole-exome targeted sequencing of the uncharacterized pine genome. Plant J. 75, 146–156. doi: 10.1111/tpj.12193

PubMed Abstract | CrossRef Full Text | Google Scholar

Nystedt, B., Street, N. R., Wetterbom, A., Zuccolo, A., Lin, Y. C., Scofield, D. G., et al. (2013). The Norway spruce genome sequence and conifer genome evolution. Nature 497, 579–584. doi: 10.1038/nature12211

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliva, J., Bernat, M., and Stenlid, J. (2013). Heartwood stump colonisation by Heterobasidion parviporum and H. annosum s.s. in Norway spruce (Picea abies) stands. For. Ecol. Manage. 295, 1–10. doi: 10.1016/j.foreco.2013.01.005

CrossRef Full Text | Google Scholar

Pavy, N., Gagnon, F., Deschênes, A., Boyle, B., Beaulieu, J., and Bousquet, J. (2016). Development of highly reliable in silico SNP resource and genotyping assay from exome capture and sequencing: an example from black spruce (Picea mariana). Mol. Ecol. Resour. 16, 588–598. doi: 10.1111/1755-0998.12468

PubMed Abstract | CrossRef Full Text | Google Scholar

Piri, T., and Hamberg, L. (2015). Persistence and infectivity of Heterobasidion parviporum in Norway spruce root residuals following stump harvesting. For. Ecol. Manage. 353, 49–58. doi: 10.1016/j.foreco.2015.05.012

CrossRef Full Text | Google Scholar

Ramírez, V., López, A., Mauch-Mani, B., Gil, M. J., and Vera, P. (2013). An extracellular subtilase switch for immune priming in Arabidopsis. PLoS Pathog. 9:e1003445. doi: 10.1371/journal.ppat.1003445

PubMed Abstract | CrossRef Full Text | Google Scholar

Robert, J. A., Pitt, C., Bonnett, T. R., Yuen, M. M., Keeling, C. I., Bohlmann, J., et al. (2013). Disentangling detoxification: Gene expression analysis of feeding mountain pine beetle illuminates molecular-level host chemical defense detoxification mechanisms. PLoS ONE 8:e77777. doi: 10.1371/journal.pone.0077777

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaller, A., Stintzi, A., Rivas, S., Serrano, I., Chichkova, N. V., Vartapetian, A. B., et al. (2018). From structure to function – a family portrait of plant subtilases. New Phytol. 218, 901–915. doi: 10.1111/nph.14582

PubMed Abstract | CrossRef Full Text | Google Scholar

Serrano, I., Buscaill, P., Audran, C., Pouzet, C., Jauneau, A., and Rivas, S. (2016). A non canonical subtilase attenuates the transcriptional activation of defence responses in Arabidopsis thaliana. Elife 5:e19755. doi: 10.7554/eLife.19755.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Skrøppa, T., Solheim, H., and Hietala, A. (2015). Variation in phloem resistance of Norway spruce clones and families to Heterobasidion parviporum and Ceratocystis polonica and its relationship to phenology and growth traits. Scand. J. For. Res. 30, 103–111. doi: 10.1080/02827581.2014.963144

CrossRef Full Text | Google Scholar

Sridha, S., and Wu, K. (2006). Identification of AtHD2C as a novel regulator of abscisic acid responses in Arabidopsis. Plant J. 46, 124–133. doi: 10.1111/j.1365-313X.2006.02678.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, H., Paulin, L., Alatalo, E., and Asiegbu, F. O. (2011). Response of living tissues of Pinus sylvestris to the saprotrophic biocontrol fungus Phlebiopsis gigantea. Tree Physiol. 31, 438–451. doi: 10.1093/treephys/tpr027

PubMed Abstract | CrossRef Full Text | Google Scholar

Swedjemark, G., and Stenlid, J. (1995). Susceptibility of conifer and broadleaf seedlings to Swedish S and P strains of Heterobasidion annosum. Plant Pathol. 44, 73–79. doi: 10.1111/j.1365-3059.1995.tb02717.x

CrossRef Full Text | Google Scholar

Swedjemark, G., and Stenlid, J. (1996). Variation in spread of Heterobasidion annosum in clones of Picea abies grown at different vegetation phases under greenhouse conditions. Scand. J. For. Res. 11, 137–144. doi: 10.1080/02827589609382921

CrossRef Full Text | Google Scholar

Swedjemark, G., Stenlid, J., and Karlsson, B. (2001). Variation in growth of Heterobasidion annosum among clones of Picea abies incubated for different periods of time. For. Pathol. 31, 163–175. doi: 10.1046/j.1439-0329.2001.00238.x

CrossRef Full Text | Google Scholar

Terhonen, E., Marco, T., Sun, H., Jalkanen, R., Kasanen, R., Vuorinen, M., et al. (2011). The effect of latitude, season and needle-age on the mycota of scots pine (Pinus sylvestris) in Finland. Silva Fennica 45, 301–317. doi: 10.14214/sf.104

CrossRef Full Text | Google Scholar

Tian, F., Bradbury, P. J., Brown, P. J., Hung, H., Sun, Q., Flint-Garcia, S., et al. (2011). Genome-wide association study of leaf architecture in the maize nested association mapping population. Nat. Genet. 43, 159–162. doi: 10.1038/ng.746

PubMed Abstract | CrossRef Full Text | Google Scholar

Tornero, P., Mayda, E., Gómez, M. D., Cañas, L., Conejero, V., and Vera, P. (1996). Characterization of LRP, a leucine-rich repeat (LRR) protein from tomato plants that is processed during pathogenesis. Plant J. 10, 315–330. doi: 10.1046/j.1365-313X.1996.10020315.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Gunulf, A., Pukkala, T., and Rönnberg, J. (2015). Simulated Heterobasidion disease development in Picea abies stands following precommercial thinning and the economic justification for control measures. Scand. J. For. Res. 30, 174–185. doi: 10.1080/02827581.2014.978887

CrossRef Full Text | Google Scholar

Wei, X., Liu, K., Zhang, Y., Feng, Q., Wang, L., Zhao, Y., et al. (2015). Genetic discovery for oil production and quality in sesame. Nat. Commun. 6:8609. doi: 10.1038/ncomms9609

PubMed Abstract | CrossRef Full Text | Google Scholar

Zubizarreta-Gerendiain, A., Peltola, H., Pulkkinen, P., and Kellomaki, S. (2009). Effects of genetic entry and competition by neighbouring trees on growth and wood properties of cloned Norway spruce (Picea abies). Ann. For. Sci. 66:806. doi: 10.1051/forest/2009075

CrossRef Full Text | Google Scholar