Combining QTL Mapping and Transcriptomics to Decipher the Genetic Architecture of Phenolic Compounds Metabolism in the Conifer White Spruce

Abstract

Conifer forests worldwide are becoming increasingly vulnerable to the effects of climate change. Although the production of phenolic compounds (PCs) has been shown to be modulated by biotic and abiotic stresses, the genetic basis underlying the variation in their constitutive production level remains poorly documented in conifers. We used QTL mapping and RNA-Seq to explore the complex polygenic network underlying the constitutive production of PCs in a white spruce (Picea glauca) full-sib family for 2 years. QTL detection was performed for nine PCs and differentially expressed genes (DEGs) were identified between individuals with high and low PC contents for five PCs exhibiting stable QTLs across time. A total of 17 QTLs were detected for eight metabolites, including one major QTL explaining up to 91.3% of the neolignan-2 variance. The RNA-Seq analysis highlighted 50 DEGs associated with phenylpropanoid biosynthesis, several key transcription factors, and a subset of 137 genes showing opposite expression patterns in individuals with high levels of the flavonoids gallocatechin and taxifolin glucoside. A total of 19 DEGs co-localized with QTLs. Our findings represent a significant step toward resolving the genomic architecture of PC production in spruce and facilitate the functional characterization of genes and transcriptional networks responsible for differences in constitutive production of PCs in conifers.

Introduction

Empirical evidence indicates that global warming will increase pressure on conifer forests via the intensification of drought events and the introduction of new insects and pathogens (Dale et al., 2001; Sturrock et al., 2011). As long-lived species, trees largely depend on their mechanical and chemical defenses to survive and reproduce (Davis and Shaw, 2001). Trees can protect themselves by producing secondary metabolites such as terpenes and phenolic compounds including flavonoids, monolignols, phenolic acids, stilbenes, and coumarins (Berini et al., 2018). The constitutive production of phenolic compounds (PCs) is usually thought to confer broad resistance to pathogens and abiotic factors in trees (Francheschi et al., 2005; Kessler, 2015). In spruces (Picea spp.), PCs are frequently associated with defense responses against pathogenic fungi (Brignolas et al., 1995; Hammerbacher et al., 2011, 2013, 2014) or pest insects (Brignolas et al., 1998; Faccoli and Schlyter, 2007; Delvas et al., 2011; Schiebe et al., 2012; Warren et al., 2015), and a recent study also highlighted their involvement in drought resistance in Picea abies (Schiop et al., 2017). Hence, PCs likely play a major role in climate adaptation in conifers. However, the genetic basis underlying their quantitative variation remains poorly understood (Ralph et al., 2006a; Warren et al., 2015), mainly due to the occurrence of large and complex gene families in conifers, combined with a limited knowledge of their functions (Hamberger and Bohlmann, 2006; Prunier et al., 2016).

The regulatory mechanism of phenylpropanoid biosynthesis is complex and involves several families of transcription factors (TFs) that regulate the expression of downstream genes (Yang et al., 2012). Recent advances showed that the PCs pathway in plants is also controlled at different branches by R2R3-MYB repressors (Cavallini et al., 2015; Ma and Constabel, 2019). The regulation of late biosynthetic anthocyanin and proanthocyanin genes is orchestrated by the ternary MBW complex, involving transcription factors from the R2R3-MYB, basic helix–loop–helix (bHLH) and WD40 classes (Xu et al., 2015; Ma et al., 2018). Although the MBW network is evolutionary conserved in plants (Ramsay and Glover, 2005), the considerable expansion of TF families in plant lineages (Feller et al., 2011) could lead to functionally divergent genes between angiosperms and gymnosperms. The complex regulatory system of PCs metabolism, i.e., genes and their networks of interactions, remains largely undeciphered in conifer species. Recent work evidenced protein-protein interactions among different MYB, bHLH and WDR (Nemesio-Gorriz et al., 2017) and NAC gene family (Dalman et al., 2017) in Norway spruce, while other studies identified a variety of gene families (Hamberger and Bohlmann, 2006; Ralph et al., 2006a; Warren et al., 2015), genes (Celedon et al., 2017) and transcription factors (Bomal et al., 2008; Bedon et al., 2010; Deng and Lu, 2017) playing a role in the regulation of PC pathway.

The majority of PC traits are known to be predominantly polygenic (Külheim et al., 2011; Francisco et al., 2016; Ganthaler et al., 2017), hence the need to deploy genome-wide approaches to study their synthesis. Quantitative trait loci (QTL) mapping was used successfully to uncover the genetic architecture of PCs production in crops (Li et al., 2016; Czyczyło-Mysza et al., 2019) and tree species (Verdu et al., 2014; Caseys et al., 2015). QTL analyses have also been successfully conducted in the Pinaceae, using traits relative to growth and phenology (Pelgas et al., 2011; Prunier et al., 2013) and biotic stress resistance (Porth et al., 2012; Lind et al., 2014). However, these analyses suffer from low resolution in conifers because of the low density of genetic linkage maps currently available and because of the considerable time and resources needed to grow, phenotype and genotype progenies. Recently, several studies in angiosperms showed that combining QTL mapping and transcriptome profiling is a robust approach to identify candidate genes underlying complex traits (Xu et al., 2015; Ye et al., 2017; Zhang et al., 2017; Jian et al., 2019). Transcriptome profiling can be used in combination with genetic mapping to narrow down the number of candidate genes identified by QTL and pointing out key genes involved in the mechanisms of interest. Hence, we used a similar approach to identify key genes involved in the PC pathway in white spruce [Picea glauca (Moench) Voss], a conifer species with high genome complexity and for which considerable genomic resources such as a draft genome sequences (Birol et al., 2013; Warren et al., 2015), an annotated expressed gene catalog (Rigault et al., 2011), high-throughput genotyping arrays (Pavy et al., 2013), and a high-density genetic linkage map (Pavy et al., 2017) are currently available.

In this study, we first conducted QTL mapping to assess the genomic architecture of the constitutive production of nine candidate PCs in a white spruce full-sib family. Three flavonoids, one neolignan and one stilbenoid, for which we identified QTLs stable across 2 years of measurement, were further analyzed with RNA-Seq in order to identify key genes involved in the regulatory network of PC metabolism. Differentially expressed genes (DEGs) among individuals displaying contrasting phenotypes (high vs. low metabolite content) were identified and compared with genes located in QTLs.

Materials and Methods

Workflow

In this study, QTL and RNA-Seq approaches were used to identify genes involved in the metabolism of nine PCs in a random subset of 1,976 siblings previously used to produce a high-resolution genetic map for white spruce (Pavy et al., 2017). The nine PCs, all related to various biotic and abiotic stress responses in spruce species as well as in other plants (Table 1), were quantified for 164 and 202 siblings in two different years (Figure 1). QTL analyses were conducted with phenotypes obtained for the same progeny used for genetic mapping. From there, five metabolites showed a significant underlying genetic component for both years of assessment and were deemed of prime interest to be further investigated at the transcriptomic level using an RNA-Seq approach. For each of these metabolites, at least 20 siblings presenting contrasting phenotypes (i.e., highest and lowest metabolite content) were selected for individual transcriptome profiling and analyzed to highlight sets of differentially expressed genes.

Table 1

Metabolite	Molecule skeleton	Class	MW (g mol−1)a	Formula	Biotic stressb	Abiotic stressc
Astringin		Stilbenoid	406.1	C20H22O9	Fungus Hammerbacher et al. (2011)	Not studied
Catechin		Flavonoid	290.3	C15H14O6	Fungus Bahnweg et al. (2000)	Chobot et al. (2009)
Gallocatechin		Flavonoid	306.3	C15H14O7	Fungus Hammerbacher et al. (2018)	Wang et al. (2016)
Isorhapontin		Stilbenoid	420.4	C21H24O9	Fungus Hammerbacher et al. (2011)	Not studied
Neolignan-2		Lignan	400.5	C23H28O6	Insect Schiebe et al. (2012)	Moura et al. (2010)
Piceid		Stilbenoid	390.4	C20H22O8	Insect and fungus Brignolas et al. (1998)	Villangó et al. (2016)
Procyanidin B1		Flavonoid	578.5	C30H26O12	Insect Rohde et al. (1996)	Varela et al. (2016)
Taxifolin		Flavonoid	304.2	C15H12O7	Fungus Evensen et al. (2000)	Not studied
Taxifolin glucoside		Flavonoid	466.4	C21H22O12	Insect and fungus Brignolas et al. (1995)	Not studied

Bark phenolic compounds analyzed in a white spruce full-sib family.

a

Molecular weight of the phenolic compound (g mol⁻¹) and chemical formula are indicated.

b

Studies demonstrating the role of the phenolics tested in tree response to biotic stresses (insect and pathogens) in Picea sp. are reported.

c

Studies demonstrating the role of the phenolics tested in tree response to abiotic stresses (drought) in plant species are reported.

Figure 1

Plant Material and Phenotypic Data Measurements

Twigs were collected on 164 and 202 full-sib siblings in August 2014 (age 15) and August 2017 (age 18), respectively, from a QTL mapping population previously used in several studies (cross C94-1-2516 between ♀77111 × ♂2388; Figure 1; see Pelgas et al., 2011; Pavy et al., 2012, 2017) and raised in a common garden located in Saint-Antonin, Québec (47°45N; 69°28 W). Seventy-nine siblings were common to both sampling years. Twigs (3-4 copies each) corresponding to the two parents originally crossed were also sampled in August 2015, 2016, and 2017 at the Cap-Tourmente Arboretum of Natural Resources Canada, Québec (47°05 N; −70°47 W). All sampled trees were originally produced from grafted material with a normal growth. For each tree, current-year twigs (18 cm long) were collected at breast height on the same side. Each sample was cut longitudinally with a razor blade to separate the woody inner stem tissues, the pith, from the bark tissues. Each sample was then cut into two parts, immediately frozen in liquid nitrogen, and stored at−80°C until further use: one part for metabolite quantification and the other one for transcriptomic analyses.

Metabolite Quantification: Extraction of Phenolic Compounds

Extraction of PCs, namely astringin, catechin, gallocatechin, isorhapontin, 4-[1,3-dihydroxy-2-[2-hydroxy-4-(3-hydroxypropyl)phenoxy]propyl]-β-D-xylopyranoside (neolignan-2; Supplementary Figure 1), piceid, procyanidin B1, taxifolin and taxifolin glucoside, was performed at the Max Planck Institute in Germany. Briefly, PCs were extracted from ground bark samples and quantified by LC-tandem mass spectrometry on an Agilent 1200 HPLC system (Agilent, Santa Clara, CA, United States) coupled to an API 3200 mass analyzer (Sciex, Darmstadt, Germany). Analyst v1.5 software (Applied Biosystems) was used for data acquisition and processing. Linearity of metabolite detection for quantification was verified by external calibration curves for catechin, taxifolin, astringin and procyanidin B1. Other metabolite concentrations were determined relative to the calibration curve of the metabolite most closely resembling it. Detailed procedures are given in Supplementary Methods 1 in Supplementary Material. Descriptive statistics and phenotypic distributions of metabolites are reported in Supplementary Table 1 and Supplementary Figure 2. Pairwise correlations between metabolite concentrations measured in 2014 and 2017 were calculated with the software R 3.5 (R Core Team, 2013) using the rcor function in the Hmisc R package (Harrell, 2014).

QTL Analyses

QTL analyses were conducted in order to identify genomic regions accounting for variation in metabolite concentrations among siblings and underlying candidate genes putatively carrying causal variants. QTLs detection was performed for the nine PCs investigated (Table 1). Data were analyzed on a yearly basis using the two parental linkage maps consisting in 2,774 (female) and 2,308 (male) high-quality SNPs mapping in as many separated genetic bins (see details in Pavy et al., 2017), resulting in four distinct QTL analyses per metabolite. All individuals included herein were part of the progeny originally used to generate the consensus linkage map consisting of 8,793 distinct gene loci (Pavy et al., 2017). QTL analyses were conducted with the software MapQTL v6.0 (van Ooijen and Kyazma, 2009), using interval mapping, which is a robust method against deviations from normality (van Ooijen and Kyazma, 2009). A QTL with a LOD score ≥ 3.1 was considered significant (P < 0.05 based on 1,000 genome-wide permutations of the markers). QTLs detected on each parental map were then positioned on the consensus map. Genes within 1 LOD of either side of the QTL peak were considered as candidate genes. Given the extremely high LOD score associated with neolignan-2 QTL, the window was expanded to 15 cM, which corresponded to the average window size of other significant QTLs detected in this study.

RNA-Seq Experiment

QTL analyses revealed that five metabolites, namely gallocatechin, neolignan-2, piceid, procyanidin B1 and taxifolin glucoside, had significant QTLs across both years surveyed (Figure 1). The expression of genes involved in the biosynthesis of those five metabolites were further investigated using RNA-Seq to identify sets of differentially expressed genes (DEGs) between two groups of individuals that showed contrasting concentrations of metabolites. Out of the 202 siblings collected in 2017, we selected 80 individuals showing high or low concentrations in gallocatechin, piceid, procyanidin B1 and taxifolin glucoside metabolites (Figure 1), some of them showing extreme phenotypes for more than one metabolite. The transcriptomic profiles of 20 phenotypically divergent individuals per PC (10 high-metabolite content vs. 10 low-metabolite content individuals) were then compared. Since the neolignan-2 displayed a bimodal distribution of concentrations (Figure 2), each of the two contrasted groups of individuals included 35 siblings instead of 10.

Figure 2

RNA Extraction, Libraries Preparation and Sequencing

For the 80 individuals selected for RNA-Seq, total RNA was isolated from frozen bark tissues following the method of Chang et al. (1993), with the protocol modifications described in Pavy et al. (2008), as detailed in Supplementary Methods 2 in Supplementary Material. For each sibling, a cDNA library was generated from 1 μg of total RNA using the KAPA stranded mRNA-Seq Kit (KAPA Biosystems, Roche Sequencing solutions, Canada). This kit contains all buffers and enzymes required for the poly(A) mRNA capture and the construction of stranded mRNA-Seq libraries of 100 ng−4 μg of intact total RNA. A PCR of 13 cycles was then performed for each cDNA library having a specific adapter (Illumina TruSeq HT). The quality of total RNA was assessed using the Agilent 2100 Bioanalyzer with RNA 6000 Nano LabChips (Agilent Technologies Inc, Santa Clara, CA, USA) and RNA concentration was determined with a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, U.S.A). For all samples, RNA concentration exceeded 100 ng μl⁻¹ and RNA Integrity Number (RIN) exceeded 7.8. Libraries were first tagged individually, and then merged into a single sequencing pool at equimolar concentrations. The pool was sequenced at McGill University and Genome Quebec Innovation Center (Montreal, Quebec, Canada) using two lanes of Illumina NovaSeq 6000 S2 PE100.

Pre-processing and Differential Gene Expression Analyses

For each of the five metabolites tested, gene expression was compared between groups of individuals showing high and low PC contents. Quality of raw sequence data was first checked using FASTQC v0.11.2 (Andrews et al., 2014) for each sibling sequenced. Adapter sequences were trimmed using Cutadapt (Martin, 2011) and checked again with FASTQC v0.11.2 to ensure proper trimming. For each metabolite, reads of individuals making up the high and low PC content groups were pooled. High-quality reads were then aligned against the most recent version of the white spruce reference transcriptome (Rigault et al., 2011) using Salmon v0.11.0 with the following options: –Gcbias, –seqBias, –validatingMappings and minScoreFraction (0.8) (Patro et al., 2017).

Read counts were further normalized using the DESeq2 Bioconductor package in R (Anders and Huber, 2010; Love et al., 2014). Differential gene expression analyses were then carried out using the default pipeline implemented in the function “DESeq.” The DEGs were then determined using a threshold of < 0.05 (adjusted p-value after applying the Benjamini-Hochberg's multiple-testing correction; Benjamini and Hochberg, 1995). Expression level of DEGs was further classified as high (log₂fc ≥ 1) or low (log₂fc ≤ −1) according to the log₂ fold change values (log₂fc) obtained when comparing expression in high vs. low metabolite content individuals. DEGs with a total read count lower than 10 were discarded. Venn diagrams were plotted to highlight unique and shared DEGs among metabolites using a bio-analytic web tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Functional Characterization and Enrichment Tests of Differentially Expressed Genes

Using the BLAST2GO PRO suite (Gotz et al., 2008), functional annotations of the DEGs were retrieved by performing a BLASTX search against the RefSeq database (RefSeq release 91, using a cut-off e-value of ≤10⁻¹⁰). Homologous protein domains from translated sequences were then identified by searching against the Interpro database. GO annotations (molecular functions and biological processes) were obtained for each DEG and slimmed down into the more general GO plant slim terms. For each metabolite, enrichment tests were performed by comparing GO associated with DEGs with those associated with the full set of expressed genes, using Fisher's exact tests (P < 0.05) in the BLAST2GO PRO suite. The transcription factor database PlantRegMap/PlantTFDB v5.0 was used to identify potential transcription factors among DEGs (TFs; http://planttfdb.cbi.pku.edu.cn/; Jin et al., 2017; Tian et al., 2020). Finally, DEG functions were visualized in MapMan v3.5.1 (Thimm et al., 2004) using DEG mapping files generated by Mercator (v3.6; Lohse et al., 2014). MapMan classifies genes into various relevant functional categories, also referred as BIN classes herein.

Results

Variation in Phenolic Compounds and Correlations

The nine PCs examined in this study belonged to the flavonoid, the stilbenoid, and the lignan biosynthesis pathways (Table 1). Astringin and isorhapontin were the most abundant metabolites (62.35 and 22.20 mg g⁻¹ DW, respectively), followed by taxifolin (12.55 mg g⁻¹ DW). The other PCs were relatively less abundant, with concentrations <3.28 mg g⁻¹ DW (Supplementary Table 1). All metabolites followed a unimodal and continuous distribution (Supplementary Figure 2), except for neolignan-2, which displayed a clear bimodal distribution in parents (Figure 2A) and progeny (Figure 2B), suggesting the possibility of a major gene effect.

Most metabolite concentrations within individuals were positively correlated (Figure 3). An exception was found for gallocatechin and taxifolin glucoside, which appeared negatively correlated in 2014 (Figure 3). As a general trend, PCs synthesized from the same biosynthetic pathway were most highly correlated, as observed for the two stilbenoids astringin and piceid (r = 0.60 in 2014 and r = 0.70 in 2017), and for the flavonoids catechin and procyanidin B1 (r > 0.73). However, one other exception was the strong correlation between piceid (a stilbenoid) and catechin (a flavonoid) (r > 0.78 in 2017). The concentration of neolignan-2 only moderately correlated with other metabolites (0.14 < r < 0.36), but variations in concentration among individuals were stable across time (r = 0.94, P < 0.05; based on the 79 individuals measured in both years; data not shown).

Figure 3

QTL Analyses

A total of 17 significant QTLs (LOD score ≥ 3.1) were detected for eight of the nine metabolites tested (Table 2). Individual QTLs included between 23 and 160 genes, resulting in a list of 871 unique candidate genes (Supplementary Table 2). QTLs for astringin, catechin and taxifolin were only detected from 2017 data (Table 2). Significant QTLs for gallocatechin, neolignan-2, piceid, and procyanidin B1 were observed for both years on at least one co-localizing region of a unique parent. In addition, two distinct QTLs were detected for taxifolin glucoside in 2014 and 2017, on different parents and genomic regions. The phenotypic variance explained (PVE) by single QTLs ranged between 8.0% (which is approximately corresponding to the genome-wide significance threshold for this experiment) and 91.3%. Interestingly, very high LOD scores (86.3 and 75.9) and PVE values (91.3 and 82.4%) were obtained for the only QTL detected for neolignan-2 in both years (Table 2; Figure 2C). Those observations highlight the occurrence of a genetic effect with strong impact on this phenotype on linkage group 4 (Figure 2C). Overall, QTLs were located on seven out of the 12 linkage groups (LG), corresponding to the 12 chromosomes of white spruce. Co-localization of QTLs for both years of measurement was observed for four PC (gallocatechin on LG8; neolignan-2 and piceid on LG4; and procyanidin B1 on LG1). This indicates that genotype-phenotype linkages were robust to inter-annual climate variations, which strengthen the involvement of these regions in the control of the studied phenotypes. These compounds were thus selected for transcriptomic analysis along with taxifolin glucoside, for which QTLs were detected for both years of measurement, but on different parents (Table 2; Figure 1).

The functional annotation of genes located within QTL regions relied on BLAST2GO annotations, GO terms classification, and protein family signatures. While most QTL genes had no known direct role in the phenylpropanoid pathway, we found two genes previously reported as involved in the PC biosynthesis in Picea abies, namely a leucoanthocyanidin reductase-like gene (LAR; GQ03701_M12) for procyanidin B1, and a flavonoid 3′,5′-hydroxylase 2-like gene (F3′5′H; GQ03712_G11) for taxifolin glucoside (Supplementary Table 2).

Relative Functions of Differentially Expressed Genes Between Groups of High and Low PC Content Individuals

The five selected metabolites (gallocatechin, neolignan-2, piceid, procyanidin B1, and taxifolin glucoside) were further investigated using RNA-Seq to compare the transcriptomic profiles of individuals showing contrasting phenotypes (high vs. low metabolite concentrations for each PC). A summary of RNA-Seq statistics is presented in Supplementary Table 3. A total of 603 unique DEGs were identified at a rate of 3 to 372 DEGs per metabolite (Figure 4A; Supplementary Table 4). DEGs overlap among metabolites was very limited to non-existent, except for gallocatechin and taxifolin glucoside, which shared 137 DEGs (Figure 4B). Except for one gene, these DEGs showed higher expression in high taxifolin glucoside content individuals and lower expression in high gallocatechin content individuals (Figure 4C). Among the 603 DEGs, 436 were successfully annotated with a GO term associated to molecular functions and biological processes (Supplementary Figure 3). The GO enrichment analyses revealed that DEGs were particularly overrepresented in monosaccharide metabolic (2.5%), secondary metabolic (1.9%) and phenylpropanoid metabolic (1.9%) processes (Supplementary Figure 3). Overrepresented molecular functions included stress-related functions such as catalytic (59.2%), oxidoreductase (16.8%), lyase (7.6%), and transferase (5.7%) activities, and functions related to the phenylpropanoid pathway such as UDP-glycosyltransferase and O-methyltransferase activities (Supplementary Figure 3).

Figure 4

MapMan analyses were conducted to classify DEGs into relevant functional bins (or functional categories; Supplementary Figure 4; Supplementary Table 5). DEGs were assigned to 29 functional bins. While an important proportion of the 603 DEGs were not assigned (37%), a significant proportion of DEGs was related to secondary metabolism (8.9%), stress (5.9%) and signaling (5%) (Supplementary Figure 4). The stress-related genes included 20 unique DEGs involved in abiotic stress response (Figure 5A). The most relevant functions associated to the 144 unique DEGs involved in biotic stress response included pathogenesis-related (PR) proteins, heat shock proteins, transcription factors as well as various genes involved in hormone signaling and defense (Figure 5A). The 59 defense-related DEGs identified by MapMan (Figure 5A) included 11 putative PR-proteins, and 48 genes involved in secondary metabolism, among which 36 genes were involved in the shikimate/phenylpropanoid pathway (Figure 5B).

Figure 5

Identification of Candidate Genes Specifically Associated to the Phenylpropanoid Pathway

For all metabolites considered, the combination of MapMan classification and Blast2GO functional annotations allowed for the identification of 50 DEGs involved in the phenylpropanoid pathway associated with monophenols, phenylpropanoids, flavonoids, as well as lignins and lignans metabolism (Table 3; Figure 5B). From this subset, 31 and 30 DEGs were respectively associated to gallocatechin and taxifolin glucoside, including 15 genes shared between the two DEG sets. In addition, four other genes were associated with neolignan-2 and procyanidin B1. For all metabolites, the most represented gene families were the putative dirigent proteins (DIR) (9 genes), laccases (7 genes) and plant peroxidases (6 genes) (Table 3). Among the 50 DEGs identified, 36 DEGs were differentially expressed between the groups of high and low PC content individuals. A large proportion of putative DIR showed higher expression, with three genes associated with taxifolin glucoside (GQ03808_J11, GQ03806_D05, GQ03307_E08), two with gallocatechin (GQ036606_J07, GQ03815_M16), one with neolignan-2 (WS00740_J05) and one with procyanidin B1 (GQ01301_K10). Two DIR showing lower expression were detected for gallocatechin (GQ03806_D05, GQ03803_O03). Most laccases associated to the phenylpropanoid pathway were also differentially expressed, with four highly expressed and three lower expressed genes in individuals with high taxifolin glucoside and high gallocatechin content, respectively. In addition, three putative peroxidases (GQ0202_L09, GQ03322_C02, GQ02016_E21), two O-methyltransferases (WS00740_E09, GQ03507_F11), one cytochrome P450 (WS00736_D10) and one NmrA-like (GQ03009_B07) putative protein showed higher expression in high taxifolin glucoside content individuals.

In this study, genes known to be involved in the phenylpropanoid pathway in white spruce or other plant species, and genes that were found differentially expressed were designated as key genes. Among the 50 DEGs associated to the phenylpropanoid pathway (Table 3), we identified five key DEGs, namely the 4-coumarate-CoA ligase (4CL), the phenylalanine ammonia-lyase (PAL), cinnamoyl-CoA reductase (CCR), p-coumaroyl shikimate/quinate 3′–hydroxylase (C3H) and caffeic O-methyltransferase (COMT) (Figure 6). Six other key enzymes were associated specifically to the flavonoid pathway, and included the flavonol synthase (FLS), the UDP-dependent glucosyl transferase (UGT), the flavonoid 3′5′-hydroxylase (F3′5′H - CYP75A), the bifunctional dihydroflavonol 4-reductase/flavanone 4-reductase (DFR), the leucoanthocyanidin reductase (LAR) and the anthocyanidin reductase (ANR) (Figure 6). It should be noted that FLS, UGT and F3H displayed opposite gene expression patterns in individuals producing high levels of taxifolin glucoside, compared to those producing high levels of gallocatechin (see Supplementary Table 4 for details).

Figure 6

Differentially Expressed Transcription Factors as Potential Regulators of the Phenylpropanoid Pathway

A total of 41 transcription factors (TFs) were differentially expressed among the high and low PCs content groups (Supplementary Table 6). Among those, MYB (10 genes), WRKY (5 genes), NAC and AP2-EREBP (3 genes) were the most represented families. No significant TFs were found for piceid, one for procyanidin B1, while only three TFs were detected for neolignan-2. The remaining TFs were associated to taxifolin glucoside and gallocatechin. One bHLH (GQ01301_E17), one WD40 (GQ03304_I14), one WRKY (GQ03304_E19) and eight MYBs genes showed high and low expression in individuals producing high levels of taxifolin glucoside and gallocatechin, respectively.

Combining QTL and Transcriptomic Analyses to Identify Candidate Genes for Phenolic Compound Variation

A total of 19 DEGs that were associated to gallocatechin, neolignan-2, and taxifolin glucoside, were also located within QTL regions (Supplementary Table 7). Those genes can be safely designated as candidate genes for PC variation. They included five hydrolases, four uncharacterized putative proteins as well as two genes encoding methylesterases and transferases. None of them were formally associated to the phenylpropanoid pathway.

Discussion

Neolignan: A Likely Rare Case of Monogenic Trait in Conifers

Two types of phenotypic distributions were observed in the present study (Figure 2; Supplementary Figure 2). Stilbenoids and flavonoids displayed a continuous and unimodal distribution of concentrations typical for quantitative traits controlled by multiple genes with small effects (Supplementary Figure 2), as generally observed for Pinaceae taxa (Routaboul et al., 2012; Wahyuni et al., 2014; Ganthaler et al., 2017). The distribution of their 15 QTLs across 7 linkage groups (Table 2), the moderate proportion of variance explained by single QTLs (10.9% ≤ PVE ≤ 13%; Table 2), and the distribution of 197 DEGs across the 12 spruce linkage groups (Supplementary Table 2) further support the scenario of a polygenic control of stilbenoid and flavonoid production. In contrast, neolignan-2 displayed a bimodal distribution of metabolite content (Figure 2). This distribution suggests a monogenic control of metabolite concentration, consistent with the identification of a single major QTL on LG4 (Figure 2C) explaining up to 91.3% of the phenotypic variance (Table 2). Until now, studies reporting the occurrence of neolignans in the Picea genus, and more generally in conifers, are still scarce [but see Hong et al. (2014); for review see Tanase et al. (2019)]. The neolignan analyzed in the current study (neolignan-2) was recently identified in Norway spruce (Nemesio-Gorriz et al., 2017), and found under the regulation of several MYB genes. Neolignans have diverse physiological roles in angiosperm plants, including defense against fungus and insects (Choi et al., 2009; Saguez et al., 2013), but little is known regarding their specific role in conifers. Monogenic traits usually represent prime candidates for marker-assisted selection in breeding programs. In this sense, further research should focus on the characterization and physiological role of neolignan-2 in white spruce, and densifying the QTL involved with more gene markers, given that less than one third of the transcriptome has been mapped in white spruce (Pavy et al., 2017).

Metabolite Data and Expression Profiles Reveal New Insights Into the Transcriptional Control of Flavonoids

Transcriptomic data highlighted a set of 137 genes showing opposite expression patterns between high gallocatechin and high taxifolin glucoside individuals (r = −0.88, P < 0.001; Figure 4). While further investigations are needed to determine which of these genes play a role in the constitutive production of these two flavonoids, such a strong negative relationship suggests that these 137 genes are likely co-regulated, as commonly observed for genes related to the biosynthesis of flavonoids (Honda et al., 2002; Tian et al., 2017). This opposite expression pattern also points to a possible feedback mechanism influencing the production of these two metabolites. However, the lack of significant negative correlation between gallocatechin and taxifolin glucoside concentrations in contrasted phenotypes (Supplementary Figure 5) indicates that gene expression and metabolite levels are decoupled. This does not rule out the possibility of a feedback mechanism, given that metabolite concentrations can be affected by metabolic network connectivity, when metabolites are involved in several metabolic pathways (Wegner and Kummer, 2005; Zelezniak et al., 2014). Alternatively, the activation of these genes, or a subset of them, may simply enhance the production of taxifolin glucoside, without affecting the production of gallocatechin, while their repression would yield an opposite pattern. Since both pathways require the same substrate, this scenario would imply that substrate availability is not a limiting factor for the synthesis of these two metabolites.

Trancriptomics Highlighted Key Genes Involved in the Phenylpropanoid Pathway and Environmental Stress Response

Transcriptomics revealed that members of gene families generally associated with the phenylpropanoid pathway in plants (Reinprecht et al., 2017) such as DIR, PER-LAC, cytochrome P450, NmrA-like and O-methyltransferases were differentially expressed in white spruce (Table 3; Figure 6). Some key genes involved in the flavonoid pathway (FLS, UGT and F3H; Figure 6) displayed opposite gene expression patterns in individuals producing high levels of taxifolin glucoside and gallocatechin. This observation indicates that these genes are important actors of the transcriptional machinery governing the production of flavonoids. Two genes showing higher expression in trees producing high taxifolin glucoside and gallocatechin contents (GQ03519_N09 and GQ0082_B18; Figure 6) encoded the flavanone 3-hydroxylase (F3H), an enzyme mediating the production of taxifolin in spruce (Hammerbacher et al., 2019). A previous study reported a negative feedback regulation of F3H expression by catechins in Camellia sinensis (Singh et al., 2008). Similar mechanism may occur in white spruce, since catechin and taxifolin metabolite levels were strongly positively correlated in trees producing high or low levels of gallocatechin and taxifolin glucoside (Supplementary Figure 5). Expression level of UGT85A8 was higher in individuals producing high levels of taxifolin glucoside, in line with the commonly accepted hypothetical biosynthesis pathway in spruce where UGT facilitated the glucosylation of taxifolin to taxifolin glucoside (Hammerbacher et al., 2019). One gene encoding a flavonoid 3′5′-hydroxylase (F3′5′H; GQ04107_C21) had a higher expression level in high gallocatechin content individuals (Table 3; Figure 6). F3′5′H is an important branch point enzyme in flavonoid biosynthesis that catalyzes the conversion of flavonols into 3′,4′,5′-hydroxylated derivatives and allows the formation of gallocatechin (Deng and Lu, 2017; Figure 6). In Norway spruce, the formation of gallocatechin preferentially leads through catechin rather than dihydromyricetin and leucodelphinidin, and was induced following bark beetle-fungus infection (Hammerbacher et al., 2018, 2019). Finally, one gene located in the procyanidin B1 QTL (GQ03701_M12) encoded the leucoanthocyanidin reductase (LAR) (Figure 6), a bifunctional enzyme catalyzing the reduction of leucocyanidin and controlling the degree of polymerization of proanthocyanidins (Jun et al., 2018; Yu et al., 2019). In Norway spruce, the homologous gene of GQ03701_M12 (Hammerbacher et al., 2014) and other LAR genes (Oliva et al., 2015; Nemesio-Gorriz et al., 2016) were found to be involved in response to fungus infection. Interestingly, one dirigent protein (DIR, WS00740_J05) was highly expressed in individuals with high neolignan-2 content (Table 3; Figure 6). Dirigent proteins are involved in the stereoselective reaction forming the lignan pinoresinol from coniferyl alcohol (Davin and Lewis, 2005), and could also act in the formation of lignin (Burlat et al., 2001). These proteins are usually induced by wounding as well as weevil and budworm herbivory attacks (Ralph et al., 2006b; Lippert et al., 2007).

Phenylpropanoid biosynthesis is a common plant response to biotic and abiotic stress. In line with this idea, several DEGs identified herein were previously found involved in climate adaptation (Supplementary Table 8; Hornoy et al., 2015), cold hardening and cold acclimation in various spruce species (Supplementary Table 8; Holliday et al., 2008; Kayal et al., 2011; Pelgas et al., 2011). We also identified DEGs related to plant defense pathways other than the phenylpropanoid pathway, such as the terpenoid pathway (Figure 5B; Supplementary Tables 4, 8), possibly indicating signaling crosstalk between secondary metabolites production and stress response (Jacobo-Velázquez et al., 2015; Isah, 2019).

Identification of Transcription Factors Potentially Involved in the Transcriptional Control of Flavonoids

We identified 40 differentially expressed TFs from families (i.e., MYB, bHLH, WD40, WRKY and AP2/EREBP; Supplementary Table 6) known as involved in spruce PCs pathway (e.g., Bomal et al., 2008; Bedon et al., 2010; Nemesio-Gorriz et al., 2017). This is consistent with the fact that in higher plants, several R2R3-MYB proteins are known to activate the early steps of flavonoid biosynthesis, whereas late biosynthetic genes are rather controlled by the MYB-bHLH-WD40 (MBW) complex (Xu et al., 2015; Ma et al., 2018). In addition, recent studies reported that the action of the MBW complex can be modified by WRKY TFs (Lloyd et al., 2017). In this study, we found two R2R3-MYBs (GQ03719_G10, white spruce homolog of PaMYB33; and GQ04002_F03, white spruce homolog of PaMYB31) showing higher expression in high taxifolin glucoside content individuals (Supplementary Table 6), and for which protein-protein interactions with bHLH members were shown to play a role in the regulation of the flavonoid pathway in Norway spruce (Nemesio-Gorriz et al., 2017). Interestingly, the overexpression of PaMYB33 in transgenic cell lines also activated the expression of genes encoding PAL, ANR and LAR enzymes (Nemesio-Gorriz et al., 2017), suggesting that similar mechanisms may govern the late flavonoid biosynthesis in white spruce. The 12 TFs displaying opposite expression patterns between gallocatechin and taxifolin glucoside in our study (Supplementary Table 6) might be directly involved in the formation of the MBW complex or might be downstream target genes of the MBW complex. Given that the regulatory mechanism involving the MBW complex appears highly conserved in higher plants (Xu et al., 2015), our findings provide a relevant framework for studying the complex transcriptional network governing the biosynthesis of flavonoids, in white spruce and other conifer species.

Integrating QTL Mapping and Transcriptomic Approaches: Potential, Limits and Perspectives

While RNA-Seq provided valuable insights into genes involved in white spruce phenylpropanoid pathway and their regulation, the combination of both QTL mapping and RNA-Seq approaches yielded mixed results. We identified 19 candidate DEGs colocalizing within QTLs, but none of them were formally associated with the PC pathway (Supplementary Table 7). However, considering that the landscape of genes essential for the regulation of phenylpropanoids in plants is still not exhaustive (for a review see Biala and Jasiński, 2018), we should not exclude the possibility that these candidate genes could be indirectly linked to the phenylpropanoid biosynthesis pathway. Further, the fact that one third (203) of the 603 DEGs identified herein have not been positioned yet on the current white spruce genetic map (Pavy et al., 2017) suggests that some additional DEGs could have colocalized with the detected QTLs. Association studies focusing on these DEGs would also help identify more candidate genes. Transcriptomics alone proved powerful to investigate the combinatorial gene regulation of flavonoids in white spruce. It allowed the identification of 137 genes likely co-regulated, along with several candidate regulators, and 50 genes encoding key enzymes of the white spruce phenylpropanoid pathway. These results may be further explored by focusing on the discovery of variable genomic regions responsible for variations in gene expression (i.e., eQTL studies), given that PCs regulatory networks appear relatively conserved in plants (e.g., MBW complex; see Xu et al., 2015).

References

• 1

  AndersS.HuberW. (2010). Differential expression analysis for sequence count data. Genome Biol.11:106. 10.1186/gb-2010-11-10-r106

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  AndrewsK. R.HohenloheP. A.MillerM. R.HandB. K.SeebJ. E.LuikartG. (2014). Trade-offs and utility of alternative RADseq methods: reply to Puritz et al. Mol. Ecol.23, 5943–5946. 10.1111/mec.12964

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  BahnwegG.SchubertR.KehrR. D.Müller-StarckG.HellerW.LangebartelsC.et al. (2000). Controlled inoculation of Norway spruce (Picea abies) with Sirococcus conigenus: PCR-based quantification of the pathogen in host tissue and infection-related increase of phenolic metabolites. Trees14, 435–441. 10.1007/s004680000058

  • CrossRef
  • Google Scholar

• 4

  BedonF.BomalC.CaronS.LevasseurC.BoyleB.MansfieldS. D.et al. (2010). Subgroup 4 R2R3-MYBs in conifer trees: gene family expansion and contribution to the isoprenoid- and flavonoid-oriented responses. J. Exp. Bot.61, 3847–3864. 10.1093/jxb/erq196

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  BenjaminiY.HochbergY. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B. Stat. Methodol. 57, 289–300. 10.1111/j.2517-6161.1995.tb02031.x

  • CrossRef
  • Google Scholar

• 6

  BeriniJ. L.BrockmanS. A.HegemanA. D.ReichP. B.MuthukrishnanR.MontgomeryR. A.et al. (2018). Combinations of abiotic factors differentially alter production of plant secondary metabolites in five woody plant species in the boreal-temperate transition zone. Front. Plant Sci.9:1257. 10.3389/fpls.2018.01257

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  BialaW.JasińskiM. (2018). The phenylpropanoid case – it is transport that matters. Front. Plant Sci.9:1610. 10.3389/fpls.2018.01610

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  BirolI.RaymondA.JackmanS. D.PleasanceS.CoopeR.TaylorG. A.et al. (2013). Assembling the 20 Gb white spruce (Picea glauca) genome from whole-genome shotgun sequencing data. Bioinformatics29, 1492–1497. 10.1093/bioinformatics/btt178

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  BomalC.BedonF.CaronS.MansfieldS. D.LevasseurC.CookeJ. E. K.et al. (2008). Involvement of Pinus taeda MYB1 and MYB8 in phenylpropanoid metabolism and secondary cell wall biogenesis: a comparative in planta analysis. J. Exp. Bot.59, 3925–3939. 10.1093/jxb/ern234

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  BrignolasF.LacroixB.LieutierF.SauvardD.DrouetA.ClaudotA.-C.et al. (1995). Induced responses in phenolic metabolism in two norway spruce clones after wounding and inoculations with Ophiostoma polonicum, a bark beetle-associated fungus. Plant Physiol.109, 821–827. 10.1104/pp.109.3.821

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  BrignolasF.LieutierF.SauvardD.ChristiansenE.BerrymanA. A. (1998). Phenolic predictors for Norway spruce resistance to the bark beetle Ips typographus (Coleoptera: Scolytidae) and an associated fungus, Ceratocystis polonica. Can. J. For. Res.28, 720–728. 10.1139/x98-037

  • CrossRef
  • Google Scholar

• 12

  BurlatV.KwonM.DavinL. B.LewisN. G. (2001). Dirigent proteins and dirigent sites in lignifying tissues. Phytochemistry57, 883–897. 10.1016/S0031-9422(01)00117-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  CaseysC.StrittC.GlauserG.BlanchardT.LexerC. (2015). Effects of hybridization and evolutionary constraints on secondary metabolites: the genetic architecture of phenylpropanoids in European Populus species. PLoS ONE10:e0128200. 10.1371/journal.pone.0128200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  CavalliniE.MatusJ. T.FinezzoL.ZenoniS.LoyolaR.GuzzoF.et al. (2015). The phenylpropanoid pathway is controlled at different branches by a set of R2R3-MYB C2 repressors in grapevine. Plant Physiol.167, 1448–1470. 10.1104/pp.114.256172

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  CeledonJ. M.YuenM. M. S.ChiangA.HendersonH.ReidK. E.BohlmannJ. (2017). Cell-type- and tissue-specific transcriptomes of the white spruce (Picea glauca) bark unmask fine-scale spatial patterns of constitutive and induced conifer defense. Plant J.92, 710–726. 10.1111/tpj.13673

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  ChangS.PuryearJ.CairneyJ. (1993). A simple and efficient method for isolating RNA from pine trees. Plant Mol. Biol. Rep.11, 113–116. 10.1007/BF02670468

  • CrossRef
  • Google Scholar

• 17

  ChobotV.HuberC.TrettenhahnG.HadacekF. (2009). (±)-Catechin: chemical weapon, antioxidant, or stress regulator?J. Chem. Ecol.35, 980–996. 10.1007/s10886-009-9681-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  ChoiN. H.ChoiG. J.MinB. S.JangK. S.ChoiY. H.KangM. S.et al. (2009). Effects of neolignans from the stem bark of Magnolia obovata on plant pathogenic fungi. J. Appl. Microbiol.106, 2057–2063. 10.1111/j.1365-2672.2009.04175.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  Czyczyło-MyszaI. M.CyganekK.DziurkaK.QuarrieS.SkrzypekE.MarcińskaI.et al. (2019). Genetic parameters and QTLs for total phenolic content and yield of wheat mapping population of CSDH lines under drought stress. Int. J. Mol. Sci.20:6064. 10.3390/ijms20236064

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  DaleV. H.JoyceL. A.McnultyS.NeilsonR. P.AyresM. P.FlanniganM. D.et al. (2001). Climate change and forest disturbances: climate change can affect forests by altering the frequency, intensity, duration, and timing of fire, drought, introduced species, insect and pathogen outbreaks, hurricanes, windstorms, ice storms, or landslides. Bioscience51, 723–734. 10.1641/0006-3568(2001)051[0723:CCAFD]2.0.CO;2

  • CrossRef
  • Google Scholar

• 21

  DalmanK.WindJ. J.Nemesio-GorrizM.HammerbacherA.LundénK.EzcurraI.et al. (2017). Overexpression of PaNAC03, a stress induced NAC gene family transcription factor in Norway spruce leads to reduced flavonol biosynthesis and aberrant embryo development. BMC Plant Biol.17:6. 10.1186/s12870-016-0952-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  DavinL. B.LewisN. G. (2005). Lignin primary structures and dirigent sites. Curr. Opin. Biotechnol.16, 407–415. 10.1016/j.copbio.2005.06.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  DavisM. B.ShawR. G. (2001). Range shifts and adaptive responses to quaternary climate change. Science292, 673–679. 10.1126/science.292.5517.673

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  DelvasN.BauceÉ.LabbéC.OllevierT.BélangerR. (2011). Phenolic compounds that confer resistance to spruce budworm. Entomol. Exp. Appl.141, 35–44. 10.1111/j.1570-7458.2011.01161.x

  • CrossRef
  • Google Scholar

• 25

  DengY.LuS. (2017). Biosynthesis and regulation of phenylpropanoids in plants. CRC Crit. Rev. Plant Sci.36, 257–290. 10.1080/07352689.2017.1402852

  • CrossRef
  • Google Scholar

• 26

  EvensenP. C.SolheimH.HøilandK.StenersenJ. (2000). Induced resistance of Norway spruce, variation of phenolic compounds and their effects on fungal pathogens. For. Pathol.30, 97–108. 10.1046/j.1439-0329.2000.00189.x

  • CrossRef
  • Google Scholar

• 27

  FaccoliM.SchlyterF. (2007). Conifer phenolic resistance markers are bark beetle antifeedant semiochemicals. Agric. For. Entomol.9, 237–245. 10.1111/j.1461-9563.2007.00339.x

  • CrossRef
  • Google Scholar

• 28

  FellerA.MachemerK.BraunE. L.GrotewoldE. (2011). Evolutionary and comparative analysis of MYB and bHLH plant transcription factors. Plant J.66, 94–116. 10.1111/j.1365-313X.2010.04459.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  FrancheschiV. R.KrokeneP.ChristiansenE.KreklingT. (2005). Anatomical and chemical defenses of conifer bark against bark beetles and other pests. New Phytol.167, 353–376. 10.1111/j.1469-8137.2005.01436.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  FranciscoM.AliM.FerreresF.MorenoD. A.VelascoP.SoengasP. (2016). Organ-specific quantitative genetics and candidate genes of phenylpropanoid metabolism in Brassica oleracea. Front. Plant Sci.6:1240. 10.3389/fpls.2015.01240

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  GanthalerA.StögglW.KrannerI.MayrS. (2017). Foliar phenolic compounds in Norway spruce with varying susceptibility to Chrysomyxa rhododendri: analyses of seasonal and infection-induced accumulation patterns. Front. Plant Sci.8:1173. 10.3389/fpls.2017.01173

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  GotzS.Garcia-GomezJ. M.TerolJ.WilliamsT. D.NagarajS. H.NuedaM. J.et al. (2008). High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res.36, 3420–3435. 10.1093/nar/gkn176

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  HambergerB.BohlmannJ. (2006). Cytochrome P450 mono-oxygenases in conifer genomes: discovery of members of the terpenoid oxygenase superfamily in spruce and pine. Bioch. Soc. Trans.34, 1209–1214. 10.1042/BST0341209

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  HammerbacherA.KandasamyD.UllahC.SchmidtA.WrightL. P.GershenzonJ. (2019). Flavanone-3-Hydroxylase plays an important role in the biosynthesis of spruce phenolic defenses against bark beetles and their fungal associates. Front. Plant Sci.10:208. 10.3389/fpls.2019.00208

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  HammerbacherA.PaetzC.WrightL. P.FischerT. C.BohlmannJ.DavisA. J.et al. (2014). Flavan-3-ols in Norway spruce: biosynthesis, accumulation, and function in response to attack by the bark beetle-associated fungus Ceratocystis polonica. Plant Physiol.164, 2107–2122. 10.1104/pp.113.232389

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  HammerbacherA.RaguschkeB.WrightL. P.GershenzonJ. (2018). Gallocatechin biosynthesis via a flavonoid 3′, 5′-hydroxylase is a defense response in Norway spruce against infection by the bark for beetle-associated sap-staining fungus Endoconidiophora polonica. Phytochemistry148, 78–86. 10.1016/j.phytochem.2018.01.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  HammerbacherA.RalphS. G.BohlmannJ.FenningT. M.GershenzonJ.SchmidtA. (2011). Biosynthesis of the major tetrahydroxystilbenes in spruce, astringin and isorhapontin, proceeds via resveratrol and is enhanced by fungal infection. Plant Physiol.157, 876–890. 10.1104/pp.111.181420

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  HammerbacherA.SchmidtA.WadkeN.WrightL. P.SchneiderB.BohlmannJ.et al. (2013). A common fungal associate of the spruce bark beetle metabolizes the stilbene defenses of Norway spruce. Plant Physiol.162, 1324–1336. 10.1104/pp.113.218610

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  HarrellF. E.Jr. (2014). Hmisc Package v3.14–4. Available online at: https://cran.r-project.org/web/packages/Hmisc/index.html (accessed March 15, 2018).

  • Google Scholar

• 40

  HollidayJ. A.RalphS. G.WhiteR.BohlmannJ.AitkenS. N. (2008). Global monitoring of autumn gene expression within and among phenotypically divergent populations of Sitka spruce (Picea sitchensis). New Phytol.178, 103–122. 10.1111/j.1469-8137.2007.02346.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  HondaC.KotodaN.WadaM.KondoS.KobayashiS.SoejimaJ.et al. (2002). Anthocyanin biosynthetic genes are coordinately expressed during red coloration in apple skin. Plant Physiol. Biochem.40, 955–962. 10.1016/S0981-9428(02)01454-7

  • CrossRef
  • Google Scholar

• 42

  HongS. S.JeongW.KimJ. K.KwonJ. G.LeeJ. Y.AhnE.-K.et al. (2014). Neolignan inhibitors of antigen-induced degranulation in RBL-2H3 cells from the needles of Pinus thunbergii. Fitoterapia99, 347–351. 10.1016/j.fitote.2014.10.015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  HornoyB.PavyN.GérardiS.BeaulieuJ.BousquetJ. (2015). Genetic adaptation to climate in white spruce involves small to moderate allele frequency shifts in functionally diverse genes. Genome Biol. Evol.7, 3269–3285. 10.1093/gbe/evv218

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  IsahT. (2019). Stress and defense responses in plant secondary metabolites production. Biol. Res.52:39. 10.1186/s40659-019-0246-3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  Jacobo-VelázquezD. A.González-AgüeroM.Cisneros-ZevallosL. (2015). Cross-talk between signaling pathways: the link between plant secondary metabolite production and wounding stress response. Sci. Rep.5:8608. 10.1038/srep08608

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  JianH.ZhangA.MaJ.WangT.YangB.ShuangL. S.et al. (2019). Joint QTL mapping and transcriptome sequencing analysis reveal candidate flowering time genes in Brassica napus L. BMC Genomics20:21. 10.1186/s12864-018-5356-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  JinJ.TianF.YangD.-C.MengY.-Q.KongL.LuoJ.et al. (2017). PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants. Nucleic Acids Res.45, D1040–D1045. 10.1093/nar/gkw982

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  JunJ. H.XiaoX.RaoX.DixonR. A. (2018). Proanthocyanidin subunit composition determined by functionally diverged dioxygenases. Nat. Plants4, 1034–1043. 10.1038/s41477-018-0292-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  KayalW. E.AllenC. C.JuC. J.-T.AdamsE. R. I.King-JonesS.ZahariaL. I.et al. (2011). Molecular events of apical bud formation in white spruce, Picea glauca. Plant Cell. Environ.34, 480–500. 10.1111/j.1365-3040.2010.02257.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  KesslerA. (2015). The information landscape of plant constitutive and induced secondary metabolite production. Curr. Opin. Insect Sci.8, 47–53. 10.1016/j.cois.2015.02.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  KülheimC.YeohS. H.WallisI. R.LaffanS.MoranG. F.FoleyW. J. (2011). The molecular basis of quantitative variation in foliar secondary metabolites in Eucalyptus globulus. New Phytol.191, 1041–1053. 10.1111/j.1469-8137.2011.03769.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  LiM.-W.MuñozN. B.WongC.-F.WongF.-L.WongK.-S.WongJ. W.-H.et al. (2016). QTLs regulating the contents of antioxidants, phenolics, and flavonoids in soybean seeds share a common genomic region. Front. Plant Sci.7:854. 10.3389/fpls.2016.00854

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  LindM.KällmanT.ChenJ.MaX.-F.BousquetJ.MorganteM.et al. (2014). A Picea abies linkage map based on SNP markers identifies QTLs for four aspects of resistance to heterobasidion parviporum infection. PLoS ONE9:e101049. 10.1371/journal.pone.0101049

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  LippertD.ChowriraS.RalphS. G.ZhuangJ.AeschlimanD.RitlandC.et al. (2007). Conifer defense against insects: proteome analysis of Sitka spruce (Picea sitchensis) bark induced by mechanical wounding or feeding by white pine weevils (Pissodes strobi). Proteomics7, 248–270. 10.1002/pmic.200600525

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  LloydA.BrockmanA.AguirreL.CampbellA.BeanA.CanteroA.et al. (2017). Advances in the MYB–bHLH–WD repeat (MBW) pigment regulatory model: addition of a WRKY factor and co-option of an anthocyanin MYB for betalain regulation. Plant Cell Physiol.58, 1431–1441. 10.1093/pcp/pcx075

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  LohseM.NagelA.HerterT.MayP.SchrodaM.ZrennerR.et al. (2014). Mercator: a fast and simple web server for genome scale functional annotation of plant sequence data. Plant Cell Environ.37, 1250–1258. 10.1111/pce.12231

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15:550. 10.1186/s13059-014-0550-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  MaD.ConstabelC. P. (2019). MYB repressors as regulators of phenylpropanoid metabolism in plants. Trends Plant Sci.24, 275–289. 10.1016/j.tplants.2018.12.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  MaD.ReicheltM.YoshidaK.GershenzonJ.ConstabelC. P. (2018). Two R2R3-MYB proteins are broad repressors of flavonoid and phenylpropanoid metabolism in poplar. Plant J.96, 949–965. 10.1111/tpj.14081

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 60

  MartinM. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J.17, 10–12. 10.14806/ej.17.1.200

  • CrossRef
  • Google Scholar

• 61

  MouraJ. C. M. S.BonineC. A. V.de Oliveira Fernandes VianaJ.DornelasM. C.MazzaferaP. (2010). Abiotic and biotic stresses and changes in the lignin content and composition in plants. J. Integr. Plant Biol.52, 360–376. 10.1111/j.1744-7909.2010.00892.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  Nemesio-GorrizM.BlairP. B.DalmanK.HammerbacherA.ArnerupJ.StenlidJ.et al. (2017). Identification of Norway spruce MYB-bHLH-WDR transcription factor complex members linked to regulation of the flavonoid pathway. Front. Plant Sci.8:305. 10.3389/fpls.2017.00305

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  Nemesio-GorrizM.HammerbacherA.IhrmarkK.KällmanT.OlsonÅ.LascouxM.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. 10.1104/pp.16.00685

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 64

  OlivaJ.RommelS.FossdalC. G.HietalaA. M.Nemesio-GorrizM.SolheimH.et al. (2015). Transcriptional responses of Norway spruce (Picea abies) inner sapwood against Heterobasidion parviporum. Tree Physiol.35, 1007–1015. 10.1093/treephys/tpv063

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 65

  PatroR.DuggalG.LoveM. I.IrizarryR. A.KingsfordC. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods14, 417–419. 10.1038/nmeth.4197

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 66

  PavyN.GagnonF.RigaultP.BlaisS.DeschênesA.BoyleB.et al. (2013). Development of high-density SNP genotyping arrays for white spruce (Picea glauca) and transferability to subtropical and Nordic congeners. Mol. Ecol. Resour.13, 324–336. 10.1111/1755-0998.12062

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  PavyN.LamotheM.PelgasB.GagnonF.BirolI.BohlmannJ.et al. (2017). A high-resolution reference genetic map positioning 8.8 K genes for the conifer white spruce: structural genomics implications and correspondence with physical distance. Plant J.90, 189–203. 10.1111/tpj.13478

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 68

  PavyN.PelgasB.BeauseigleS.BlaisS.GagnonF.GosselinI.et al. (2008). Enhancing genetic mapping of complex genomes through the design of highly-multiplexed SNP arrays: application to the large and unsequenced genomes of white spruce and black spruce. BMC Genomics9:21. 10.1186/1471-2164-9-21

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  PavyN.PelgasB.LarocheJ.RigaultP.IsabelN.BousquetJ. (2012). A spruce gene map infers ancient plant genome reshuffling and subsequent slow evolution in the gymnosperm lineage leading to extant conifers. BMC Biol.10:84. 10.1186/1741-7007-10-84

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 70

  PelgasB.BousquetJ.MeirmansP. G.RitlandK.IsabelN. (2011). QTL mapping in white spruce: gene maps and genomic regions underlying adaptive traits across pedigrees, years and environments. BMC Genomics12:145. 10.1186/1471-2164-12-145

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  PorthI.WhiteR.JaquishB.AlfaroR.RitlandC.RitlandK. (2012). Genetical genomics identifies the genetic architecture for growth and weevil resistance in spruce. PLoS ONE7:e44397. 10.1371/journal.pone.0044397

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 72

  PrunierJ.PelgasB.GagnonF.DespontsM.IsabelN.BeaulieuJ.et al. (2013). The genomic architecture and association genetics of adaptive characters using a candidate SNP approach in boreal black spruce. BMC Genomics14:368. 10.1186/1471-2164-14-368

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  PrunierJ.VertaJ.-P.MacKayJ. J. (2016). Conifer genomics and adaptation: at the crossroads of genetic diversity and genome function. New Phytol.209, 44–62. 10.1111/nph.13565

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 74

  R Core Team (2013). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for statistical computing. Available online at: http://www.R-project.org

  • Google Scholar

• 75

  RalphS.ParkJ.-Y.BohlmannJ.MansfieldS. D. (2006b). Dirigent proteins in conifer defense: gene discovery, phylogeny, and differential wound- and insect-induced expression of a family of DIR and DIR-like genes in spruce (Picea spp.). Plant Mol. Biol.60, 21–40. 10.1007/s11103-005-2226-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 76

  RalphS. G.YuehH.FriedmannM.AeschlimanD.ZeznikJ. A.NelsonC. C.et al. (2006a). Conifer defence against insects: microarray gene expression profiling of Sitka spruce (Picea sitchensis) induced by mechanical wounding or feeding by spruce budworms (Choristoneura occidentalis) or white pine weevils (Pissodes strobi) reveals large-scale changes of the host transcriptome. Plant Cell Environ.29, 1545–1570. 10.1111/j.1365-3040.2006.01532.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  RamsayN. A.GloverB. J. (2005). MYB–bHLH–WD40 protein complex and the evolution of cellular diversity. Trends Plant Sci.10, 63–70. 10.1016/j.tplants.2004.12.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 78

  ReinprechtY.PerryG. E.PaulsK. P. (2017). A comparison of phenylpropanoid pathway gene families in common bean. Focus on P450 and C4H genes in The Common Bean Genome, eds Pérez de la VegaM.SantallaM.MarsolaisF. (Cham: Springer International Publishing), 219–261. 10.1007/978-3-319-63526-2_11

  • CrossRef
  • Google Scholar

• 79

  RigaultP.BoyleB.LepageP.CookeJ. E. K.BousquetJ.MacKayJ. J. (2011). A white spruce gene catalog for conifer genome analyses. Plant Physiol.157, 14–28. 10.1104/pp.111.179663

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 80

  RohdeM.WaldmannR.LunderstädtJ. (1996). Induced defence reaction in the phloem of spruce (Picea abies). and larch (Larix decidua) after attack by Ips typographus and Ips cembrae. Forest Ecol. Manag.86, 51–59. 10.1016/S0378-1127(96)03802-9

  • CrossRef
  • Google Scholar

• 81

  RoutaboulJ.-M.DubosC.BeckG.MarquisC.BidzinskiP.LoudetO.et al. (2012). Metabolite profiling and quantitative genetics of natural variation for flavonoids in Arabidopsis. J. Exp. Bot.63, 3749–3764. 10.1093/jxb/ers067

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 82

  SaguezJ.AttoumbréJ.GiordanengoP.Baltora-RossetS. (2013). Biological activities of lignans and neolignans on the aphid Myzus persicae (Sulzer). Arthropod Plant Interact.7, 225–233. 10.1007/s11829-012-9236-x

  • CrossRef
  • Google Scholar

• 83

  SchiebeC.HammerbacherA.BirgerssonG.WitzellJ.BrodeliusP. E.GershenzonJ.et al. (2012). Inducibility of chemical defenses in Norway spruce bark is correlated with unsuccessful mass attacks by the spruce bark beetle. Oecologia170, 183–198. 10.1007/s00442-012-2298-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 84

  SchiopS. T.Al HassanM.SestrasA. F.BoscaiuM.SestrasR. E.VicenteO. (2017). Biochemical responses to drought, at the seedling stage, of several Romanian Carpathian populations of Norway spruce (Picea abies L. Karst). Trees31, 1479–1490. 10.1007/s00468-017-1563-1

  • CrossRef
  • Google Scholar

• 85

  SinghK.RaniA.KumarS.SoodP.MahajanM.YadavS. K.et al. (2008). An early gene of the flavonoid pathway, flavanone 3-hydroxylase, exhibits a positive relationship with the concentration of catechins in tea (Camellia sinensis). Tree Physiol.28, 1349–1356. 10.1093/treephys/28.9.1349

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 86

  SturrockR. N.FrankelS. J.BrownA. V.HennonP. E.KliejunasJ. T.LewisK. J.et al. (2011). Climate change and forest diseases. Plant Pathol.60, 133–149. 10.1111/j.1365-3059.2010.02406.x

  • CrossRef
  • Google Scholar

• 87

  TanaseC.CoşarcăSMunteanD.-L. (2019). A critical review of phenolic compounds extracted from the bark of woody vascular plants and their potential biological activity. Molecules24:1182. 10.3390/molecules,24061182

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 88

  ThimmO.BläsingO.GibonY.NagelA.MeyerS.KrügerP.et al. (2004). MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J.37, 914–939. 10.1111/j.1365-313X.2004.02016.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 89

  TianF.YangD.-C.MengY.-Q.JinJ.GaoG. (2020). PlantRegMap: charting functional regulatory maps in plants. Nucleic Acids Res.48, D1104–D1113. 10.1093/nar/gkz1020

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 90

  TianJ.ZhangJ.HanZ.-Y.SongT.-T.LiJ.-Y.WangY.-R.et al. (2017). McMYB12 transcription factors co-regulate proanthocyanidin and anthocyanin biosynthesis in Malus crabapple. Sci. Rep.7:43715. 10.1038/srep43715

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 91

  van OoijenJ. W.KyazmaB. V. (2009). MapQTL® 6, Software for the Mapping of Quantitative Trait Loci in Experimental Populations of Diploid Species. Wageningen: Kyazma B.V.

  • Google Scholar

• 92

  VarelaM. C.ArslanI.ReginatoM. A.CenzanoA. M.LunaM. V. (2016). Phenolic compounds as indicators of drought resistance in shrubs from Patagonian shrublands (Argentina). Plant Physiol. Biochem.104, 81–91. 10.1016/j.plaphy.2016.03.014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 93

  VerduC. F.GuyotS.ChildebrandN.BahutM.CeltonJ.-M.GaillardS.et al. (2014). QTL analysis and candidate gene mapping for the polyphenol content in cider apple. PLoS ONE9:e107103. 10.1371/journal.pone.0107103

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 94

  VillangóS. Z.SzekeresA.BencsikO.LáposiR.PálfiZ.ZsófiZ. S. (2016). The effect of postveraison water deficit on the phenolic composition and concentration of the Kékfrankos (Vitis vinifera L.) berry. Sci. Hortic.209, 113–116. 10.1016/j.scienta.2016.06.010

  • CrossRef
  • Google Scholar

• 95

  WahyuniY.Stahl-HermesV.BallesterA.-R.de VosR. C. H.VoorripsR. E.MaharijayaA.et al. (2014). Genetic mapping of semi-polar metabolites in pepper fruits (Capsicum sp.).: towards unravelling the molecular regulation of flavonoid quantitative trait loci. Mol. Breed.33, 503–518. 10.1007/s11032-013-9967-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 96

  WangW.XinH.WangM.MaQ.WangL.KaleriN. A.et al. (2016). Transcriptomic analysis reveals the molecular mechanisms of drought-stress-induced decreases in Camellia sinensis leaf quality. Front. Plant Sci.7:385. 10.3389/fpls.2016.00385

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 97

  WarrenR. L.KeelingC. I.YuenM. M. S.RaymondA.TaylorG. A.VandervalkB. P.et al. (2015). Improved white spruce (Picea glauca) genome assemblies and annotation of large gene families of conifer terpenoid and phenolic defense metabolism. Plant J.83, 189–212. 10.1111/tpj.12886

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 98

  WegnerK.KummerU. (2005). A new dynamical layout algorithm for complex biochemical reaction networks. BMC Bioinformatics6:212. 10.1186/1471-2105-6-212

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 99

  XuE.VaahteraL.HõrakH.HinchaD. K.HeyerA. G.BroschéM. (2015). Quantitative trait loci mapping and transcriptome analysis reveal candidate genes regulating the response to ozone in Arabidopsis thaliana. Plant Cell Environ.38, 1418–1433. 10.1111/pce.12499

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 100

  YangC.-Q.FangX.WuX.-M.MaoY.-B.WangL.-J.ChenX.-Y. (2012). Transcriptional regulation of plant secondary metabolism. J. Integr. Plant Biol.54, 703–712. 10.1111/j.1744-7909.2012.01161.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 101

  YeJ.YangY.ChenB.ShiJ.LuoM.ZhanJ.et al. (2017). An integrated analysis of QTL mapping and RNA sequencing provides further insights and promising candidates for pod number variation in rapeseed (Brassica napus L.). BMC Genomics18:71. 10.1186/s12864-016-3402-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 102

  YuK.JunJ. H.DuanC.DixonR. A. (2019). VvLAR1 and VvLAR2 are bifunctional enzymes for proanthocyanidin biosynthesis in grapevine. Plant physiol.180, 1362–1374. 10.1104/pp.19.00447

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 103

  ZelezniakA.SheridanS.PatilK. R. (2014). Contribution of network connectivity in determining the relationship between gene expression and metabolite concentration changes. PLoS Comput. Biol.10:e1003572. 10.1371/journal.pcbi.1003572

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 104

  ZhangD.ZhangH.ChuS.LiH.ChiY.Triebwasser-FreeseD.et al. (2017). Integrating QTL mapping and transcriptomics identifies candidate genes underlying QTLs associated with soybean tolerance to low-phosphorus stress. Plant Mol. Biol.93, 137–150. 10.1007/s11103-016-0552-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar