Introduction
============

Outbreaks of insect herbivores can cause widespread mortality of forest trees and drastically alter landscapes. Examples of forest pests with large outbreak dynamics and severe impacts on their environments include the spruce budworm (SBW) (Gray and MacKinnon, [@b15]), the mountain pine beetle (Kurz *et al*., [@b25]; Raffa *et al*., [@b37]) and the emerald ash borer (Poland and McCullough, [@b36]). The SBW \[*Choristoneura fumiferana* (Clem)\] is a lepidopteron native to eastern North America that defoliates conifers, primarily mature spruce (*Picea*) and fir (*Abies*) trees (Blais, [@b3]) (Figure[1](#fig01){ref-type="fig"}a). The SBW outbreak of 1950--1993 covered an area of 850 000 km^2^ in Canada, killed 45--58% of tree hosts in highly affected areas and decreased wood yields by 300--6800 m^3^ km^−2^ (Gray and MacKinnon, [@b15]). Efforts to predict future occurrences and model the potential impacts of SBW outbreaks indicate that they will be exacerbated by changing environmental conditions (Gray and MacKinnon, [@b15]), similar to the recent spread of the mountain pine beetle epidemic (Kurz *et al*., [@b25]; Raffa *et al*., [@b37]). An understanding of mechanisms of resistance to the SBW, including its molecular underpinnings, may enhance our ability to predict and potentially mitigate the impacts of SBW and evaluate the ability of host trees to acclimate and adapt to changing conditions.

![Experimental system and expression of the *Pgβglu-1* gene.(a) Distribution of white spruce (*Picea glauca*) and spruce budworm (*Choristoneura fumiferana*; SBW). \[Distribution map accessed at <https://cfs.nrcan.gc.ca/projects/107> (30 March 2014) and reproduced with permission from the Natural Resources Canada/Canadian Forest Service.\] The black circle and square show the test trees site and the population survey area, respectively. In the SBW box: larval stage 6 and adult moth. In the white spruce box: resistant (R) and recovering non-resistant (N-R) phenotypes shown in 2010 about 3 years after a local spruce budworm outbreak.(b) Genes with largest differences in transcript abundance in the foliage of R (*n* = 7) and N-R trees (*n* = 7), identified by using an oligonucleotide microarray.(c) The relative expression of glycosyl hydrolase family 1 transcripts was compared across eight tissues: vegetative buds (buds), foliage, xylem mature (m), xylem juvenile (j), phelloderm, roots, megagametophyte (MGG) and embryonic culture (EC). Putative glycosyl hydrolase transcripts are annotated with the closest Arabidopsis match found in the TAIR database (<http://www.arabidopsis.org/>). Expression in foliage was measured in 4-year-old trees not previously exposed to SBW.](tpj0081-0068-f1){#fig01}

Host resistance against SBW has been linked to growth phenology and the chemical composition of foliage (Clancy, [@b8]; Daoust *et al*., [@b10]; Delvas *et al*., [@b11]), but the underlying mechanisms and population patterns of variation are largely unknown. Constitutive and inducible defence mechanisms against insects have been identified in herbaceous plants like *Nicotiana attenuata* and *Arabidopsis thaliana* (Kessler and Baldwin, [@b22]; Todesco *et al*., [@b48]; Züst *et al*., [@b54]) and in forest trees such as spruce (Hall *et al*., [@b16]) and poplar (Irmisch *et al*., [@b21]). In forest trees, an understanding of the molecular mechanisms underlying variations in insect resistance has begun to develop (Robert *et al*., [@b41]; Hall *et al*., [@b16]) and is expected to be key to determining their ability to survive and adapt, when considering their high levels of variability across populations, wide-ranging distributions and changing environmental conditions.

As forward and reverse genetic techniques can be difficult to apply in non-model species and are poorly developed for mature forest trees, we used an alternative approach based on comparative transcriptome profiling of resistant (R) and non-resistant (N-R) trees in a forest setting to discover genes involved in SBW resistance in white spruce (*Picea glauca* (Voss.) Moench.). Mature white spruce trees with contrasting resistance levels were previously identified during a localized SBW outbreak (1998--2007) (Figure[1](#fig01){ref-type="fig"}a) (Daoust *et al*., [@b10]). Average annual defoliation intensities served to identify trees as R (0--20% defoliation) or N-R (30--70% defoliation). At the end of the outbreak, R trees had remained healthy whereas N-R trees were dead or had recovered from defoliation but had sparse foliage (Figure[1](#fig01){ref-type="fig"}a). The foliage of R trees accumulated high levels of the acetophenone metabolites piceol and pungenol as well as the corresponding glucosides, picein and pungenin (Figure[2](#fig02){ref-type="fig"}). In contrast, the foliage of N-R trees contained predominantly the acetophenone glucosides and very low levels of the aglycons (Delvas *et al*., [@b11]). Compared with controls, laboratory-reared SBW fed with piceol and pungenol at concentrations found in the needles of R trees showed survival rates reduced by up to 50%, slower development and lower pupal mass (Delvas *et al*., [@b11]).

![The structures of biologically active acetophenone aglycons, pungenol and piceol, and their respective glucosides, pungenin and picein.](tpj0081-0068-f2){#fig02}

Comparative transcriptome screening of mature R and N-R white spruce trees, reported here, revealed a glucosyl hydrolase gene, *Pgβglu-1*, highly expressed in R trees. We identified and characterized the corresponding enzyme function and variation in gene expression, which explained differences in the acetophenone aglycon levels in R and N-R trees. We monitored *Pgβglu-1* transcript profiles over time and between years, its inheritance and its geographic origins within a white spruce population. In summary, we report a β-glucosidase gene of plant acetophenone metabolism, whose expression and function explain differences in acetophenone levels associated with SBW resistance in white spruce.

Results
=======

Transcriptome profiling identifies *Pgβglu-1* as highly expressed in SBW-resistant trees
----------------------------------------------------------------------------------------

We screened the expression of nearly 24 000 different white spruce genes in the foliage of R and N-R trees. Transcripts of 236 genes were more abundant in R trees, while transcripts of 250 genes were more abundant in N-R trees (*P* \< 0.05) (Table S1 in Supporting Information). Transcript levels for most of the genes varied less than 10-fold between R and N-R trees. In contrast, transcripts corresponding to a gene identified here as *Pgβglu-1* (GenBank BT114253) were 770 times more abundant in current-year foliage of R trees than in N-R trees (Figure[1](#fig01){ref-type="fig"}b). Transcripts of only one other gene of unknown function varied more than 20-fold (Figure[1](#fig01){ref-type="fig"}b). Validation by quantitative reverse transcriptase polymerase chain reaction assay (RT-qPCR) showed *Pgβglu-1* transcript abundance about 1000-fold higher in R compared with N-R trees, with the same samples as used with the microarray (Mann--Whitney--Wilcoxon, *P* = 0.0015).

Full-length cDNA and genomic sequences of *Pgβglu-1*
----------------------------------------------------

We cloned and sequenced the full-length (FL) *Pgβglu-1* complementary DNA (cDNA) (GenBank KJ780719) and the corresponding genomic DNA. The translated *Pgβglu-1* FL cDNA sequence encodes a PgβGLU-1 protein 506 amino acids long (Figure S2) which contains a glycosyl hydrolase family 1 domain characteristic for enzymes that hydrolyse phenolic glucosides, disaccharides or other glycosylated substrates (Hill and Reilly, [@b17]) and shares 52% of amino acids with the *A. thaliana* enzyme β-glucosidase 40.

*Pgβglu-1* is preferentially expressed in foliage
-------------------------------------------------

The white spruce gene catalogue, which represents 27 720 unique sequences (Rigault *et al*., [@b39]), has 242 sequences containing a glycosyl hydrolase domain (33 different PFAM domains); 22 of these contained a glycosyl hydrolase family 1 domain most similar to plant β-glucosidases. Of these sequences, *Pgβglu-1* had the highest expression among those that appeared to be selectively expressed in foliage (Figure[1](#fig01){ref-type="fig"}c; data from the PiceaGeneExpress database; Raherison *et al*., [@b38]).

Assessment of *Pgβglu-1* sequence variation in R and N-R trees
--------------------------------------------------------------

We identified intact *Pgβglu-1* genes in both R and N-R trees. Sequence variations were screened over 3598 bp including exons (13) and introns (12) as well as upstream and downstream regions. Comparing genomic *Pgβglu-1* sequences from seven R and seven N-R trees, we identified a total of 19 polymorphic sites in the translated protein sequences. However, none of the variations were specific to either R or N-R trees (Table S2), and none of the variations affected the conserved functional sites (Figure S2). Thus, these variations would not easily explain the observed chemical differences in R and N-R trees. None of the trees were homozygous nonsense *Pgβglu-1* mutants.

Differential *Pgβglu-1* expression correlates with levels of biologically active acetophenones
----------------------------------------------------------------------------------------------

We monitored levels of *Pgβglu-1* transcripts, the acetophenones piceol and pungenol and the acetophenone glucoside picein over two growth seasons, 3 and 6 years after SBW outbreak. Levels of *Pgβglu-1* transcripts, piceol and pungenol were higher in R than N-R trees in June or August in both years (Figure[3](#fig03){ref-type="fig"}a--c). *Pgβglu-1* expression was positively correlated with piceol (Pearson\'s correlation, *r* = 0.66, *P* \< 0.001) and pungenol (Pearson\'s correlation, *r* = 0.68, *P* \< 0.001) and negatively correlated with picein (Pearson\'s correlation, *r* = −0.39, *P* = 0.02). Piceol and pungenol levels were positively correlated (Pearson\'s correlation, *r* = 0.84, *P* \< 0.001). These correlations are as expected for putative substrate and products of the predicted β-glucosidase enzyme activity of the *Pgβglu-1* gene product.

![Differential *Pgβglu-1* expression and acetophenones levels were maintained over time. Foliar phenotypes of test trees during 2010 (*n* = 8) and 2013 (*n* = 10) of resistant (R) and non-resistant trees (N-R). We report results of mixed model analysis with year, month, resistance phenotype as fixed effects and with their interactions.(a) *Pgβglu-1* transcript levels were higher in R than N-R trees (*F*~1,26~ = 92.17, *P *\< 0.001) and in June than August (*F*~1,26~ = 10.58, *P *= 0.003). (B--D) Acetophenone concentrations from R and N-R trees. Piceol and pungenol levels are higher in R than in N-R trees (piceol *F*~1,28~ = 37.77, *P *\< 0.001, pungenol *F*~1,28~ = 50.42, *P *\< 0.001). Picein levels did not differ between resistance phenotypes (*F*~1,28~ = 0.01, ns) but varied between months (*F*~1,28~ = 5.04, *P* = 0.03). All panels show mean ± SEM.](tpj0081-0068-f3){#fig03}

Predicted properties and structural features of PgβGLU-1
--------------------------------------------------------

The predicted PgβGLU-1 protein has a theoretical isoelectric point (pI) of 5.57, an N-terminal signal peptide (score of 0.915 out of 1) and a cleavage site between amino acids 25 and 26 (Figure S2) (SignalP 4.1 server; Petersen *et al*., [@b35]), and appears to be targeted to the secretory pathway ([target]{.smallcaps}P v1.1 software; secretory pathway score 0.966 out of a maximum of 1; scores for mitochondrion, chloroplast and others were below 0.124; secretory pathway reliability class = 1 out of 5, where 1 is highest; Emanuelsson *et al*., [@b12]). [wolfpsort]{.smallcaps} and [multiloc]{.smallcaps} software (Höglund *et al*., [@b18]; Horton *et al*., [@b20]) indicated that PgβGLU-1 is a vacuolar protein (vacuolar score of 0.92 out of 1; scores for extracellular matrix and endoplasmic reticulum were below 0.05). These predictions suggest that PgβGLU-1 is targeted to the vacuole and may be released to its lumen after cleavage of a hydrophobic signal peptide.

We used protein modelling for initial assessment of possible interactions of PgβGLU-1 with the glucosylated acetophenones picein and pungenin. PgβGLU-1 was predicted to have a (β/α)~8~-barrel fold similar to other GH1 enzymes (Figure[4](#fig04){ref-type="fig"}a) with two catalytic glutamate residues (Sansenya *et al*., [@b46]). Amino acids of the modelled active site (Figure[4](#fig04){ref-type="fig"}b) interacting with the sugar moiety are highly conserved (Marana, [@b29]), while amino acids interacting with the aglycone component are less conserved. The modelled Michaelis complex with picein suggested H-bonding with the catalytic glutamate residues (Figure[4](#fig04){ref-type="fig"}b). The methyl and oxygen functionalities of the phenolic group of picein are positioned in the hydrophobic and positively charged regions, respectively, of the predicted GH1 aglycone subsite (Figure[4](#fig04){ref-type="fig"}c).

![Structural model of the *Pgβglu-*1 gene product.(a) Cartoon representation of the overall structure of the PgβGLU-1 protein.(b) Predicted PgβGLU-1 Michaelis complex with picein. Catalytic site glutamic acid residues (boxed, E199, E413): nucleophile E413 Oe1 was hydrogen bonded (dashed lines) with O5 (2.92 Å) and the distance between Oe2 of the catalytic acid/base residue E199 and O1 of picein is 3.35 Å. Three amino acids (underlined, T202, Y206, R293) predicted to form hydrogen bonds with the phenolic moiety of picein.(c) Electrostatic potential surface representation, negatively and positively charged areas in red and blue, respectively.](tpj0081-0068-f4){#fig04}

PgβGLU-1 protein catalyses the formation of biologically active acetophenones
-----------------------------------------------------------------------------

To substantiate putative functions with biochemical evidence, recombinant PgβGLU-1 protein was produced in, and purified from, *Nicotiana benthamiana* leaves, which can serve as a reliable expression system for plant proteins that are difficult to produce in microbial hosts. PgβGLU-1 enzyme was active with both picein and pungenin, forming the respective aglycons (Figure[5](#fig05){ref-type="fig"}a). The enzyme kinetic parameters of PgβGLU-1 were determined in assays with 0.01--2.0 m[m]{.smallcaps} picein. *K*~m~, *k*~cat~ and *k*~cat~/*K*~m~ were 1.28 ± 0.53 m[m]{.smallcaps}, 0.41 sec^−1^ and 0.32 m[m]{.smallcaps}^−1^ sec^−1^, respectively. The kinetic parameters are within the range of previously characterized *A. thaliana* β-glucosidases BGLU45 and BGLU46 (Escamilla-Treviño *et al*., [@b13]). PgβGLU-1 showed high level of specificity for picein compared with the following glucosides with similar structures: 4-nitrophenyl-glucoside, vanillin 4-*O*-β-glucoside, salicin, arbutin and benzyl-β-glucoside (Figure[5](#fig05){ref-type="fig"}b). Together, these results established the glucosides picein and pungenin as substrates and the aglycons piceol and pungenol as products of PgβGLU-1 enzyme activity.

![Enzyme activity of recombinant PgβGLU-1 protein.(a) Liquid chromatography mass spectrometry (LCMS) profiles for PgβGLU-1 enzyme assays. Recombinant PgβGLU-1 was assayed by incubating purified protein with an extract from white spruce needles (top) or 0.1 m[m]{.smallcaps} picein (bottom). Note that pungenin was not commercially available for testing. The LCMS profiles are shown as extracted ion chromatograms in positive mode for the parent ion of pungenol (peak 3) and pungenin (peak 2) (152 + 1; dashed line) or piceol (peak 4) and picein (peak 1) (136 + 1; continuous line). Controls were performed with no protein or recombinant GFP protein.(b) The relative activity of PgβGLU-1 on other phenylpropanoid glucosides was compared to picein. Little or no activity was observed with compounds other than picein.](tpj0081-0068-f5){#fig05}

Variation of resistance traits in the white spruce population and inheritance of expression level
-------------------------------------------------------------------------------------------------

Variability of *Pgβglu-1* expression, piceol, pungenol and picein was analysed with an independent set of genetically unrelated trees from 23 locations covering a large geographic area in the eastern range of white spruce (De Lafontaine *et al.,* [@b26]; Figure[1](#fig01){ref-type="fig"}a, Table S3) which were grown in a common garden with no observed defoliation by SBW. The levels of *Pgβglu-1* transcripts and acetophenones varied widely (Figure[6](#fig06){ref-type="fig"}a--d). The highest observed levels surpassed those of test R trees by 1.3-fold and two-fold for *Pgβglu-1* expression and both of the aglycons, respectively, and the lowest levels of the same traits were similar to the N-R trees. Several trees had intermediate levels of transcripts or acetophenones, consistent with continuous distributions of the traits in the population. Abundance of *Pgβglu-1* transcripts was positively correlated with piceol (Pearson\'s correlation, *r* = 0.43, *P* = 0.011) and pungenol levels (Pearson\'s correlation, *r* = 0.45, *P* = 0.006) as observed in test R and N-R trees but no correlation was observed between levels of *Pgβglu-1* transcripts and picein (Pearson\'s correlation, *r* = −0.23, not significant). Piceol and pungenol levels were also positively correlated (Pearson\'s correlation, *r* = 0.51, *P* = 0.002).

![Resistance trait variations in white spruce population at the end of the growth season. Distribution of *Pgβglu-1* expression (a), piceol (b), pungenol (c) and picein (d) into abundance classes determined as six equal fractions of the maximal level for each trait. The analysis used 15-year-old trees from a common garden not previously exposed to spruce budworm.](tpj0081-0068-f6){#fig06}

Inheritance of *Pgβglu-1* gene expression levels was evaluated in 1-year-old progeny of R and N-R test trees grown in controlled conditions with no exposure to SBW. Transcript accumulation was on average 20-fold higher in R progeny than N-R progeny (Figure[7](#fig07){ref-type="fig"}). Expression levels varied more widely in the R tree progeny, as expected for plants grown from wind-pollinated seeds where the paternal parent is unknown. Expression levels in young trees were also much lower than observed in mature parent trees, which was not surprising as SBW defoliation occurs mostly in mature trees.

![Inheritance of *Pgβglu-1* expression in young progeny of resistant (R) and non-resistant (N-R) parent trees.Differential expression in the foliage of 1-year-old trees in wind-pollinated families derived from R (*n* = 80 trees, four families) and N-R (*n* = 80 trees, four families) mother trees (Mann--Whitney--Wilcoxon test, *P* = 4.13 × 10^−5^). Trees were grown in the greenhouse and not exposed to spruce budworm.](tpj0081-0068-f7){#fig07}

Discussion
==========

Using a genomics approach including transcriptome and gene expression analyses and biochemical functional characterization, we discovered a white spruce *Pgβglu-1* gene and identified its function as a control point in accumulation of acetophenone in SBW resistance. Prior to their implication in SBW resistance (Delvas *et al*., [@b11]), biological activities of acetophenones in tree defence were not known, although previous work linked the ratio of piceol to picein to general defences and abiotic stress in spruce (Hoque, [@b19]; Lokke, [@b28]) and piceol to pathogen resistance in some herbaceous plant species (Curir *et al*., [@b9]). One of the acetophenone glucosides, pungenin, was tested on SBW and found to be a modest feeding deterrent for sixth-instar larvae, but no effect was found on development or mortality (Strunz *et al*., [@b47]). In contrast, the aglycons retarded development and increased mortality (Delvas *et al*., [@b11]). The upstream biosynthetic route for acetophenones was previously proposed (Negrel and Javelle, [@b32]) (Figure[8](#fig08){ref-type="fig"}, step 1), but the enzymes or genes of this pathway had not been described. Variation in the level of accumulation of biologically active acetophenones in SBW R and N-R white spruce trees appears to be determined by glucosyl hydrolysis of acetophenone glucosides, which is controlled by *Pgβglu-1* transcript levels and the encoded PgβGLU-1 enzyme activity (Figure[8](#fig08){ref-type="fig"}, step 3). The PgβGLU-1 enzyme catalyses the release of both piceol and pungenol, and their levels of accumulation were correlated both to each other and to the levels of *Pgβglu-1* transcript in the natural population, indicating that they are likely to be products of the same enzyme activity. Availability of the glucoside precursor (Figure[8](#fig08){ref-type="fig"}, step 2) does not appear to be rate-limiting for a release of aglycons, since in N-R trees levels of piceol are low despite high levels of picein (Figure[3](#fig03){ref-type="fig"}). Full assessment of the stoichiometry of acetophenone aglycon and glucoside pools requires future discovery and characterization of the corresponding glucosyltransferase(s).

![Schematic of a proposed biosynthetic pathway for the hydroxylacetophenones piceol and pungenol and their respective glucosides picein and pugenin.The proposed biosynthesis of piceol and pungenol involves, respectively, coumaroyl-CoA or caffeolyl CoA as precursors (1) (Negrel and Javelle, [@b32]). This is followed by glucosylation (2) to produce the glucosides picein and pungenin which are found in resistant and non-resistant trees. Deglucosylation by PgβGLU-1 releases the aglycons in resistant trees. The schematic does not account for different compartments, membranes or different cell types potentially involved in the formation, transport and accumulation of the glucosides and release of aglycons.](tpj0081-0068-f8){#fig08}

The tissue and cellular localization of piceol, pungenol and their glucosides, as well as the localization of PgβGLU-1 and relevant glucosyltransferase enzymes are not known. Many plant defence molecules are glycosylated and stored in the vacuole, while the β-glucosidases that may activate and release them are localized in a different part of the cell or tissue (Morant *et al*., [@b31]). Unlike, for example, cyanogenic glucosides or glucosinolates, which come into contact with their bio-activating β-glucosidases upon tissue disruption, acetophenone aglycons were found in spruce tissues that had not been exposed to mechanical damage prior to sampling. This suggests that the acetophenone glucosides come into contact with β-glucosidase enzymes without tissue disruption. Whether the glucosides and PgβGLU-1 are both contained in the same compartment, potentially the vacuole, or transported to a specialized cell type or some other location is unknown and warrants further investigation. High levels of acetophenone aglycon may be toxic to the tissues producing them (Hoque, [@b19]). It is therefore possible that acetophenones are first produced as glucosides allowing them to be transported to the vacuole and are then released by PgβGLU-1. This type of scheme where the bulk of acetophenones are initially glucosylated and then released may allow for tighter control over the concentration and containment of potentially harmful aglycons in the cell.

Here we showed that genetically heritable variation of expression of a single *Pgβglu-1* gene with preferential expression in foliage underpins the resistance of white spruce to defoliation by SBW. The accumulation of acetophenones in newly formed foliage of mature trees was found several years after previous exposure to the SBW, and was also observed in trees of broad geographic origins for which there was no obvious prior exposure to SBW. Therefore, accumulation of the acetophenone aglycons as a defence appears to be constitutive in white spruce in contrast to other glucosyl hydrolase-mediated defence systems, which invoke glycoside hydrolysis upon tissue damage for the release of active defence compounds from precursors such as cyanogenic glucosides or glucosinolates (Morant *et al*., [@b31]). Insect defence mechanisms characterized in several other conifer species are both constitutive and inducible. Resistance of lodgepole pine against the mountain pine beetle involves quantitative multigenic control (Yanchuk *et al*., [@b53]), and may be based on a suite of chemical and physical defences (Franceschi *et al*., [@b14]; Bohlmann, [@b4]). In Sitka spruce, resistance to the white pine weevil has been associated with variations of stone cells and resin ducts (King *et al*., [@b23]) and constitutive and induced terpenoid levels and their underlying terpenoid synthases (Hall *et al*., [@b16]; Roach *et al*., [@b40]). The present study showed that levels of both *Pgβglu-1* transcripts and acetophenone aglycons active against the SBW are elevated in early summer when the most destructive SBW larval stages are active (Miller, [@b30]). The constitutive production of acetophenone aglycons provides an immediate defence against ravenous SBW larvae. We did not investigate whether *Pgβglu-1* transcript or protein expression is inducible and if induction could afford phenotypic plasticity of SBW resistance.

The differences of expression of *Pgβglu-1* in R and N-R trees suggest the presence of variation of this gene in *cis* regulatory sequences, or a variation of *trans*-acting factors. Variation in the coding sequences of *Pgβglu-1* could not explain the observed phenotypic differences. Phenotypic variation within and among species result from coding sequence and regulatory mutations but some types of variations are more likely to result from regulatory mutations because of the inherent properties of transcription (Wray, [@b52]). Transcription is modular and allows fine-tuning. Regulatory mutations tend to be dominant or additive. Hence, expression variation is usually less pleiotropic and more directly visible to selection than coding variation (Wray, [@b52]; Galego Romero *et al*., [@b42]). Recent studies have shown that variations in gene expression contribute to evolutionary changes in systems ranging from crop plants (Chen *et al*., [@b7]; Koenig *et al*., [@b24]) to animals and human diseases (Romero *et al*., [@b42]). Through transcriptional variation and high outbreeding, tree populations may accumulate diverse alleles that enable adaptation to changing conditions over long periods of time.

Outbreaks of SBW have increased in frequency and scale in the last two centuries (Blais, [@b3]) and a potentially large outbreak is currently emerging in eastern Canada (NRCAN, <http://www.nrcan.gc.ca/forests/insects-diseases/13383>). Combating SBW with insecticides may be a short-term solution to reduce damage but is costly to deploy and has met with environmental concerns. While hundreds of millions of white spruce trees are planted in North America every year for forest renewal and to boost wood production, an outbreak of SBW could cancel out potential yield improvements within its range. The molecular underpinning of resistance to SBW with differential expression of a *Pgβglu-1* gene and its encoded substrate-specific β-glucosidase activity described here could enable breeding for resistance to SBW. Incorporation of knowledge about resistance mechanisms and their molecular underpinning may also improve existing models of SBW spread and risk assessment.

Experimental procedures
=======================

Plant material and sampling
---------------------------

### Test trees

*Picea glauca* trees were identified as phenotypically R or N-R to SBW defoliation in a plantation established in 1963 and located at Saint-Cyrille-de-Wendover, Québec, Canada (45°93′N, 72°52′W) (Daoust *et al*., [@b10]). Plantation trees were grown from open-pollinated seeds from local white spruce populations. Resistance phenotypes were based on levels of SBW defoliation (R, 0--20%; N-R, 30--70%) during a local outbreak (1998--2007) and field rearing of SBW in sleeve cages on the trees (Daoust *et al*., [@b10]). A total of seven R and eight N-R trees were selected from the study of Delvas *et al*. ([@b11]). All trees were 47 years old at the start of the study in 2010. Current-year foliage was sampled in 2010 on 17 June and 6 August. For each tree, three different foliage samples were taken from different parts of the live crown and used for separate analyses. One sample was from the top third of the live crown and two samples from different branches in the central third of the crown. In 2013, five R and five N-R trees were sampled on 18 June and 13 August and current-year foliage was sampled from the central third of the crown.

### Progeny of test trees

Half-sib progeny of four R and four N-R trees were obtained by growing plants from open-pollinated seeds that were collected in September 2011 (120--200 seeds per tree), cleaned and stored at −20°C. Seeds were germinated and seedlings were grown in high-density polyethylene 110 ml pots, a potting mix of peat, perlite and vermiculite (60, 20, 20%, respectively) and fertilized weekly with 20/20/20. Plants were grown in a plastic greenhouse from June 2012 to April 2013, which represents an extended growth period, and then placed in cooled growth rooms for a simulated dormancy period of 13 weeks. Plants were transferred to 1.7 L pots in July 2013 and placed in an open-end plastic house for growth. Two weeks after the transfer, foliage samples comprising newly formed foliage and 1-year-old needles were taken from the eight open-pollinated families (20 trees per family).

### Population trees

White spruce trees originating from 23 different locations were sampled in a common garden established in 1999 in Valcartier, Québec, Canada (46°56′N, 71°29′W) (Table S3). A total of 360 genotypes were available in the common garden and had been selected in genetic progeny tests of white spruce established with seeds from locations throughout the Province of Québec and are part of a tree breeding programme. In the population survey, we selected one tree randomly for each of the locations. In 2013, current-year leaves were collected from the central third of the crown of 79 trees on 3 July (early during growth season) and 1 October (late during growth season).

All foliage samples were flash frozen in liquid N~2~ immediately after removal from the trees and stored at −80°C. Foliage was ground to a fine powder using a MixerMill 300 (Retsch, <http://www.retsch.com/>) and steel grinding balls cooled in nitrogen. Powdered tissue was stored at −80°C until extraction for metabolite and RNA analyses.

Transcriptome analyses
----------------------

Total RNA was extracted as described (Chang *et al*., [@b6]) with modifications (Pavy *et al*., [@b34]) and stored at −80°C. The total RNA concentration was determined using a NanoDrop 1000 (Thermo Scientific, <http://www.thermoscientific.com/>) and assessed for quality with an Agilent 2100 Bioanalyzer and Agilent RNA 6000 Nano Kit LabChips (Agilent Technologies Inc., <http://www.agilent.com/>). Transcriptome profiling was carried out with RNA extracted from seven R and seven N-R trees sampled in 2012 (17 June). All three samples per tree (total *n* = 42) were analysed with a custom microarray comprising oligonucleotide probes for 23 853 unique *P. glauca* gene sequences (Raherison *et al*., [@b38]). Hybridizations were performed using HS 400® Pro Hybridization Stations (TECAN, <http://www.tecan.com/>). Before use in hybridization, microarray slides were washed with 0.5× saline-sodium citrate (SSC), 0.1% sodium dodecyl sulphate (SDS) at 37°C and with 2× SSC, 0.5% SDS at 50°C for 20 sec each. They were pre-hybridized for 1 h with medium agitation at 65°C in a pre-hybridization buffer (5× SSC, 0.1% SDS, 0.2 mg ml^−1^ BSA and 0.1 mg ml^−1^ herring sperm DNA) to prevent non-specific binding. Slides were washed twice with 2× SSC, 0.5% SDS. The RNA probes were prepared for hybridization by using Amino Allyl Message Amp® II aRNA Amplification Kit following the manufacturer\'s instructions (Applied Biosystems/Ambion, <http://www.lifetechnologies.com/>) and labelled with Alexa Fluor® 555 (Molecular Probes Inc., <http://www.lifetechnologies.com/>). The number of dye molecules incorporated per 1000 nucleotides was calculated by measuring absorbance on a NanoDrop 1000 (Thermo Scientific) (<http://www.nanodrop.com/>). Only samples with more than 30 dye molecules per 1000 nucleotides proceeded to the next steps. Labelled RNA samples were fragmented using RNA Fragmentation Reagents following the manufacturer\'s instructions (Applied Biosystems/Ambion), denatured for 2 min at 65°C, cooled on ice for 1 min and centrifuged at 13 000 ***g*** for 30 sec. A 120 μl volume of hybridization buffer (50% formamide, 5× SSC, 0.1% SDS and 0.1 mg ml^−1^ herring sperm DNA) pre-heated to 46°C was added to each sample for hybridization and samples were kept at 46°C in a heating block until they were added to the slides. Hybridization solutions were injected onto each slide and hybridized at 45°C for 16 h with medium agitation. Slides were washed twice with 2× SSC, 0.5% SDS at 45°C, three times with 0.5× SSC, 0.1% SDS, respectively at 45, 37 and 23°C, twice with 0.1× SSC at 23°C and once with Nanopure water at 23°C. Slides were then dried using pressurized nitrogen at 39 p.s.i. for 2 min 30 sec. Slides were scanned the same day in a PowerScanner® (TECAN) at 5 μm resolution. Laser power was set to 50% for all of the slides but the photomultiplier (PMT) was adjusted for each slide using an autogain function. Hybridization signals were converted to raw intensities with ArrayPro® Analyser version 6.3 software (MediaCybernetics, <http://www.mediacy.com/>). Median intensity and local background were obtained for each spot on the array. Bad spots resulting from dust, high local background or saturation were flagged. Net intensity was calculated by subtracting the trimmed background from median intensity and was used for further analysis.

Data analysis was carried out using the [flexarray]{.smallcaps} 1.6 software package (<http://genomequebec.mcgill.ca/FlexArray>). A-quantile normalization was used to normalize intensity distributions between arrays. Statistically significant variation was determined by univariate analysis of variance, corrected for multiple testing (Benjamini and Hochberg, [@b2]). The list of differentially expressed genes is given in Table S1. Gene identification and annotations were as previously described (Rigault *et al*., [@b39]). Tissue expression of white spruce sequences similar to known plant β-glucosidase genes was extracted from the PiceaGeneExpress database (Raherison *et al*., [@b38]).

Full length cloning of the *Pgβglu-1* cDNA
------------------------------------------

A β-glucosidase gene sequence was identified by microarray transcriptome profiling (*Pgβglu-1*); however, the longest available cDNA clone (GenBank BT114253) contained only a partial protein-coding sequence based on similarity searches to public databases. Analysis of RNA-Seq data (Rigault *et al*., [@b39]) produced a putative full-length (FL) RNA transcript of 2118 nucleotides (nt) with a 5′ untranslated region (UTR) of 177 nt and a 3′ UTR of 423 nt. Primers were designed to match the ends of the putative FL RNA sequence (Table S4); a cDNA fragment was amplified by PCR from one of the R trees and cloned using a TOPO TA Cloning® Kit for sequencing (Invitrogen, <http://www.invitrogen.com/>). Briefly, Platinum® Taq DNA Polymerase High Fidelity (Invitrogen) was used following the manufacturer\'s instructions. The amplified product was purified using the PCR Purification Kit (Invitrogen). Then, 1.25 units of Platinum® Taq DNA Polymerase (Invitrogen), 1× Platinum® Taq PCR buffer minus Mg, 1.5 m[m]{.smallcaps} MgCl~2~ and 1 m[m]{.smallcaps} of dATP were added, and the reaction was incubated for 10 min at 72°C for post-amplification addition of 3′ A- overhangs. The sample was placed on ice and the TOPO® cloning reaction was started immediately following the manufacturer\'s instructions (Invitrogen) with 15 min incubation at 24°C. Transformation was done with OneShot® Chemically Competent *Escherichia coli* (Invitrogen). The DNA from positive clones was purified using QIAprep® (Qiagen, <http://www.qiagen.com/>) and a colony-PCR was performed to verify the length of the insert before sequence analysis. Alignment with previously obtained sequences used the BioEdit Sequence Alignment Editor and sequence similarity searches were carried out with BLAST. The resulting cDNA was 1761 nt and the predicted protein coding sequence was 1518 nt or 506 amino acids, a 5′ UTR of 58 nt and a 3′ UTR of 183 nt (GenBank KJ780719). It contained a coding sequence that was identical to previous sequences for the gene.

Reverse transcriptase-qPCR analysis
-----------------------------------

Reverse transcriptase-qPCR with gene-specific primers (Table S4) was used to validate microarray results in the mature R (*n* = 5) and N-R (*n* = 8) test trees. Primers were designed using Primer3Plus software (<http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi>), self-complementarity of the newly designed primers was verified using Oligo Calc: Oligonucleotide Properties Calculator software (<http://www.basic.northwestern.edu/biotools/oligocalc.html>) and specificity was verified against the *P. glauca* gene catalogue (Rigault *et al*., [@b39]). Complementary DNA synthesis used 500 ng of total RNA and the Superscript® First-Strand cDNA synthesis system for RT-PCR (Invitrogen). Resulting cDNAs were diluted 1:4 in RNase-free water before qPCR quantification. The PCR mixtures were composed with a QuantiFast® SYBR® Green PCR kit (Qiagen) as follows: 1× master mix, 300 n[m]{.smallcaps} of 5′ and 3′ primers and 5 μl of cDNA in a final volume of 15 μl. Amplifications were carried out in a LightCycler® 480 (Roche, <http://www.roche.com/>) as described in Boyle *et al*. ([@b5]). The LRE method (Rutledge and Stewart, [@b43]) adapted for Excel (Boyle *et al*., [@b5]) was used to calculate the number of transcript molecules, which was normalized to a ratio calculated by the geometric mean of three reference genes: elongation factor 1a (EF1-α) (BT102965), cell division cycle 2 (CDC2) (BT106071) and ribosomal protein L3A (BT115036) as described elsewhere (Beaulieu *et al*., [@b1]).

Amplification specificity was assed based on melting curve analyses and by amplicon sequencing. For RT-qPCR in N-R trees, preliminary analyses showed that amplicons were produced which contained intronic sequences. Intron-containing amplicons were specific to N-R trees and the *Pgβglu-1* gene, and were greatly decreased by DNase I treatment of the RNA samples. To prevent potential confounding effects associated with residual genomic DNA, all of the primers used to assess transcript levels were positioned such as to span the junction of two exons.

To verify that RT-qPCR assay was robust when analysing genetically diverse individuals, three different primer pairs were designed and tested in test trees and in our survey of natural population trees. The results from all three primers pairs were highly correlated for the population trees, indicating high repeatability of results (Figure S1).

Identification of the genomic sequence of *Pgβglu-1*
----------------------------------------------------

Genomic DNA was extracted using DNeasy Plant Mini Kit (Qiagen). The DNA concentrations were determined using a NanoDrop 1000 (Thermo Scientific) and DNA integrity was assessed with agarose gel electrophoresis. Primers for the *Pgβglu-1* gene (BT114253) were designed as described for RT-PCR. Using genomic DNA sampled in 2010 for all the test trees, PCR amplification was conducted with Platinum Taq DNA polymerase High Fidelity (LifeTechnologies, <https://www.lifetechnologies.com/>). The PCR programme was: 5 min at 94°C followed by 35 cycles of 15 sec at 94°C and 1 min at 62°C and a final elongation step of 1 min at 68°C. Amplicons were purified with ExoSAP-IT® for PCR Product Cleanup (Affymetrix, <http://www.affymetrix.com/>) and Sanger sequenced. Sequences were aligned and analysed with BioEdit software (<http://www.mbio.ncsu.edu/bioedit/bioedit.html>). The sequence assembly was performed using the [chromaspro]{.smallcaps} 1.5 software (<http://www.technelysium.com.au/ChromasPro.html>). A unique consensus sequence was obtained by retaining the most frequent nucleotide at each site. Sites within sequences were identified as polymorphic with at least one variant out of 14 trees.

Protein modelling and substrate docking
---------------------------------------

A PgβGLU-1 protein model was built based on the crystal structure of Os4BGLU12 (Sansenya *et al*., [@b46]) using the modelling software nest from the Jackal package (<http://wiki.c2b2.columbia.edu/honiglab_public/index.php/Software:Jackal>). Model quality was assessed by analysis of a Ramachandran plot through PROCHECK (Laskowski *et al*., [@b27]). The one- or two-dimensional structures of picein, pungenin, piceol and pungeol were obtained from the PDB, PubChem and Zinc databases. Their three-dimensional structure was modelled with Babel software (<http://openbabel.org/wiki/Main_Page>). Substrate docking of picein and pungenin in the three-dimensional PgβGLU-1 protein model was performed with the software [autodock vina]{.smallcaps} (Trott and Olson, [@b49]) and results interpreted based on reports for substrate specificity in family 1 glycoside hydrolases and plant β-glucosidases (Verdoucq *et al*., [@b50]; Marana, [@b29]) by using a grid with dimensions of 30 Å × 30 Å × 30 Å centred on the catalytic acid/base (E199) and nucleophile (E413) residues. Images were generated using [pymol]{.smallcaps} (<http://www.pymol.org/>). The theoretical pI of PgβGLU-1 was determined (<http://web.expasy.org/compute_pi/>). Targeting of the PgβGLU-1 protein was assessed by predicting signal peptides with the software [signal]{.smallcaps}P (Peterson *et al*., 2001), whether it is secreted with [target]{.smallcaps}P (Emanuelsson *et al*., [@b12]) and its localization with [wo]{.smallcaps}LF PSORT (Horton *et al*., [@b20]) and [multiloc]{.smallcaps} (Höglund *et al*., [@b18]).

Protein expression and purification
-----------------------------------

*Pgβglu-1* without the signal peptide was cloned into the GoldenGate ready-expression vector pEAQ-GG, and fused to secreted SlCYS8-tag with a Factor Xa cleavable linker (Sainsbury *et al*., [@b44], [@b45]). This construct, as well as pEAQ-GG-GFP, was transformed into *Agrobacterium tumefaciens* stain AGL1 for *N. benthamiana* expression. For tissue infiltration, bacteria grown at 28°C under agitation in kanamycin-containing Luria--Bertani (LB) culture medium were harvested by centrifugation at 2400 ***g*** for 10 min, and resuspended in filtration medium \[10 m[m]{.smallcaps} 2-(N-morpholine)-ethanesulphonic acid (MES) buffer, pH 5.8, 10 m[m]{.smallcaps} MgCl~2~\] obtaining an OD~600~ of 0.5. The bacterial suspensions were blended with an equal volume of *A. tumefaciens* cells harbouring protein p19, a gene-silencing suppressor from the tomato bushy stunt virus (Voinnet *et al*., [@b51]). The bacterial mixtures were put in 1-ml needle-free syringes and gently forced into the abaxial side of fully expanded *N. benthamiana* leaves. After 5 days, leaves were harvested, flash frozen and stored at −80°C. Leaves were weighed and ground in liquid nitrogen. Five millilitres of extraction buffer (50 m[m]{.smallcaps} MES-KOH, pH 6; 50 m[m]{.smallcaps} NaCl) per gram fresh weight was added to the tissue and tissue was incubated at 4°C with shaking for 1 h. Supernatant was collected by centrifuging for 15 min at 3000 ***g*** and freezing overnight at −80°C to precipitate Rubisco protein. The protein extract was defrosted in a water bath and centrifuged for 15 min at 3000 ***g***. Imidazole (20 m[m]{.smallcaps}) and DTT (1 m[m]{.smallcaps}) were added to supernatant. Prepared Ni-NTA resin (Qiagen) was added and extracts were incubated at 4°C with shaking for 1 h. Protein was purified using a gravity column. The column was washed with buffer (50 m[m]{.smallcaps} MES-KOH, pH 6; 150 m[m]{.smallcaps} NaCl; 20 m[m]{.smallcaps} imidazole; 1 m[m]{.smallcaps} DTT) and eluted by incubating the resin for 10 min with elution buffer (50 m[m]{.smallcaps} MES-KOH, pH 6; 150 m[m]{.smallcaps} NaCl; 150 m[m]{.smallcaps} imidazole; 1 m[m]{.smallcaps} DTT) before collecting protein. Protein was desalted using a PD-25 desalting columns (GE Lifesciences, <http://www.gelifesciences.com/>) and desalting buffer (50 m[m]{.smallcaps} MES-KOH, pH 6; 10% glycerol). Protein was quantified using Bradford reagent (Bio-Rad, <http://www.bio-rad.com/>). The secretion SlCYS8 tag was cleaved by incubating the protein overnight at 16°C with Factor Xa (1 μl per 250 μl of protein).

Βeta-glucosidase enzyme assays
------------------------------

Enzyme assays were performed using reaction buffer (20 m[m]{.smallcaps} MES-KOH, pH 6), 0.01--2.0 m[m]{.smallcaps} picein (Toronto Research Chemicals Inc., <http://www.trc-canada.com/>) or foliar methanolic extract, 10 μl purified protein and 100 μ[m]{.smallcaps} benzoic acid as an internal standard in a total reaction volume of 200 μl. Assays were incubated for 1 h at 37°C then stopped with ethyl acetate and vortexed well. After centrifugation for 10 min at 1000 ***g***, 300 μl of the top layer was place in a new vial and dried under nitrogen gas. Assays were resuspended in 100 μl of 100% methanol and analysed using LC-MS. Relative activity assays were performed as above except that 0.025 m[m]{.smallcaps} of the various substrates was used. A LC-MS analysis was performed using a LC-MSD-Trap-XCT_plus with a 15-cm SB-C18 column (Agilent Technologies Inc.). Solvent A was water with 0.2% (v/v) formic acid; solvent B was 100% (v/v) acetonitrile with 0.2% (v/v) formic acid. The following gradient was used: increase to 5% solvent B from 0 to 0.5 min; increase to 22% solvent B from 0.5 to 5.0 min; increase to 35% solvent B from 5.0 to 10.0 min; increase to 50% solvent B from 10.0 to 13.0 min; increase to 95% solvent B from 13.0 to 16.0 min; holding 95% solvent B from 16.0 to 17.0 min; decrease to 5% solvent B from 17.0 to 17.1 min. Column flow rate was 0.8 ml min^−1^. Picein, piceol, pungenin and pungenol were identified using the parent mass for the aglycons piceol and pungenol, 137(+) and 153(+) respectively. Also, authentic standards of picein, piceol and pungenol were used to verify the identities (Toronto Research Chemicals Inc.; Sigma-Aldrich; TCI America, Portland, <http://www.tcichemicals.com/>). Pungenin is not commercially available.

Extraction and HPLC analysis of phenolic compounds
--------------------------------------------------

Phenolic compounds were extracted using 50--100 mg of fine needle powder with 70% aqueous methanol HPLC grade as the solvent. Benzoic acid at 1 mg ml^−1^ was used as an internal standard. Aqueous methanol (600 μl) was added to the needle powder and incubated at 4°C on an agitation plate for 48 h. Aqueous methanol was removed (after centrifugation at 13 000 ***g*** for 10 min) and preserved (−80°C) after 6, 24 and 48 h of incubation. After 6 and 24 h of incubation, a fresh 600 μl of aqueous methanol was added to each sample to continue the extraction. Extracts obtained after 6, 24 and 48 h were pooled as a single extract. The Varian Prostar HPLC used for quantification was equipped with a UV detector 325, autosampler 420 and a solvent delivery module 240. Acetophenones were separated through a pre-column Polaris MetaGuard 4.6 mm and a column Polaris 250 mm × 4.6 mm C18-A (Agilent Technologies Inc.) both heated at 62°C. The solvent and solvent gradient were identical to those used in the enzyme assay. The column flow rate was 1.5 ml min^−1^. Ten microlitres of extract diluted 1:2 in solvent A was injected. Quantification was done using calibration curves for picein, piceol and pungenol.

The authors thank P. Huron, G. Tétrault, D. St-Pierre, M. Lamara, P. Lenz and L. Madilao for technical assistance; C. Goulet and F. Sainsbury for pEAQ-GG constructs and the *Agrobacterium* strain; B. Schneider, C. Paetz, M. Fischer and S. Withers for their assistance in chemical structure confirmation; and D. Baulcombe for the pEAQ-p19 construct. The work was supported by the IFOR Research Consortium; the Natural Sciences and Engineering Research Council of Canada (NSERC; Strategic Project Grant to EB, JB, JJM); and funds from Genome Canada, Genome BC and Genome Québec for the SMarTForests Project (to JB, JJM). JB is a UBC Distinguished University Scholar.

Conflict of interest
====================

The authors confirm that they have no conflict of interest to be declared in accordance with the journal policy.

Accession numbers
=================

Sequence data for *Pgβglu-1* are available in GenBank (GenBank KJ780719) and transcriptome profiling data are available in the Gene Expression Omnibus (GEO, acc. GSE57301).

Supporting Information
======================

Additional Supporting Information may be found in the online version of this article.

###### 

Validation of RT-quantitative PCR assay.

###### 

Predicted amino acid sequence of PgβGLU-1 protein.

###### 

Differentially expressed genes between resistant (R) and non-resistant (N-R) trees.

###### 

Variations in *Pgβglu-1* amino acid sequence of resistant (R) and non-resistant (N-R) trees.

###### 

White spruce provenances sampled for the population survey.

###### 

Primer sequences used for quantitative PCR analysis, molecular cloning, and sequencing assays.

[^1]: These authors contributed equally to the work as co-first authors.
