As trait complexity grows larger in organisms, gene families encoding regulatory proteins, such as transcription factors (TFs) are generally expanded and this expansion contributes to the evolutionary diversification of traits, e.g., developmental fate or specialized metabolism (SM) (Lang et al., 2010). The concerted action of multiple classes of TFs in regulatory protein complexes may allow for further trait diversification (Singh, 1998). A well-known example of a specific regulatory protein complex driving the diversification of multiple traits in plants is the R2R3-MYB/basic helix-loop-helix (bHLH)/WD-Repeat (MBW) complex (Pesch et al., 2015; Xu et al., 2015). The MBW members cooperatively participate to form a transcriptional activating complex in which the bHLH member plays a central role and interacts with both, the MYB transcription factor and the WD40 protein. Collectively, these three members are required for the regulation of specialized metabolic pathways and of specific developmental programs in higher plants (Baudry et al., 2004; Feller et al., 2011; Tominaga-Wada and Wada, 2014; Xu et al., 2015). The MBW complex is essentially modular in its organization in diverse plant systems (Xu et al., 2015) and is well studied in the model plant Arabidopsis thaliana. For instance, Arabidopsis WD40 protein is represented by the single-copy gene Transparent Testa Glabra 1 (TTG1) that controls all traits associated with MBW complexes (Tominaga-Wada et al., 2011). Three bHLH proteins, TT8, EGL3 and GL3, act in a partially redundant manner in the Arabidopsis' MBW complex (Nesi et al., 2001; Heim et al., 2003; Toledo-Ortiz et al., 2003; Baudry et al., 2004; Pires and Dolan, 2010). Finally, the MYB TFs regulate the expression of specific genes involved in the biosynthesis of specialized metabolites (SM) (Dias et al., 2003; Baudry et al., 2004; Gonzalez, 2009). Angiosperms and gymnosperms diverged around 300 million years ago (Savard et al., 1994; Lu et al., 2014). The most widespread among the extant gymnosperms are the conifers, which dominate large forest areas, especially in boreal and montane forests in the Northern Hemisphere. Despite their evolutionary history and their ecological importance, our knowledge of the MBW complex in conifers is episodic and focused on R2R3-MYB TFs, which have been relatively well studied regarding their regulatory roles in the biosynthesis of different classes of SMs. Xue et al. (2003) identified 10 black spruce (Picea mariana) R2R3-MYBs. PmMBF1, a member of R2R3-MYBs family, is capable of inducing pigment accumulation during transient overexpression in maize cell lines and caused transactivation of the anthocyanidin-related Bz2 promoter in spruce and larch cell lines. 18 R2R3-MYBs have been discovered in white spruce (Picea glauca) and loblolly pine (Pinus taeda) (Bedon et al., 2007). Tissue specific expression patterns suggest potential roles of particular members of R2R3-MYBs in lignin biosynthesis. Indeed, the roles of at least two of the pine R2R3-MYB TFs in lignification was later confirmed as white spruce lines overexpressing these two pine R2R3-MYB TFs showed increased lignin accumulation in cell walls and induction of genes in the shikimate and monolignol biosynthetic pathways (Bomal et al., 2008). In addition, white spruce lines overexpressing the white spruce TF PgMYB14 showed increased terpene and flavonoid accumulation (Bedon et al., 2010). Similar to Arabidopsis, there appears to be a high degree of complexity associated with the R2R3-MYB component of the MBW complex in conifers, resulting in specific transcriptional responses. It has been recently shown that a similar set of genes exhibit opposing transcriptional outputs in the transgenic lines overexpressing diverse members of the R2R3-MYB TF family (Bomal et al., 2008, 2014), suggesting distinct actions of closely related R2R3-MYB TFs in conifers. In two independent studies, Bedon et al. (2010) and Xue et al. (2003) identified bHLH-binding motifs in the upstream regulatory regions of R2R3-MYB TFs, suggesting bHLH-MYB regulatory interactions that carry out a downstream gene activation cascade. However, there are very few reports on potential coniferous bHLH or WD40 members of the MBW complex. In addition, there are no reports involving conifer WD40 genes, while a stress-induced bHLH gene in Norway spruce (Picea abies) encoding a predicted protein with similarity to Arabidopsis TT8 was superficially reported by Lundén et al. (2015).

Given the previous reports of spruce transcripts with similarity to Arabidopsis TT8 and TT2 genes where two different genes have been identified as potential spruce AtTT2 homologs (Xue et al., 2003; Arnerup et al., 2011; Lundén et al., 2015), the aim of this study was to survey and identify potential members of the MBW complex in Norway spruce (Picea abies). We hypothesized that the Norway spruce genome includes a complete MBW complex repertoire, including WD40 proteins homologous to TTG1. We also hypothesized that Norway spruce would harbor several potential AtTT2 homologs in the R2R3-MYBs of subgroup S5 as seen in other species (Stracke et al., 2001; Yoshida et al., 2008; Schaart et al., 2013; Soler et al., 2015). Consistent with our hypotheses, one AtTTG1 homolog, two AtTT8 homologs and six putative AtTT2 homologs were identified in the Norway spruce genome, including the previously described TT2 (Arnerup et al., 2011). As we could identify two TT8 and six TT2 paralogs, we hypothesized that these sets of paralogs would display differentiation in protein function and/or expression patterns consistent with the theory of sub-functionalization (Innan and Kondrashov, 2010). Therefore, we quantified the expression patterns of the paralogs in different tissues and in response to various hormone analogs simulating abiotic and biotic stress and assessed their biophysical interaction with each other and the Norway spruce AtTTG1 homolog in a Yeast-two-Hybrid (Y2H) assay. In angiosperms, members of the MBW complex are known to regulate late biosynthetic genes in the flavonoid biosynthetic pathway when overexpressed in combination or as individual genes (Baudry et al., 2004; Mellway et al., 2009; Matus et al., 2010; Schaart et al., 2013). Thus, in order to elucidate if the previously identified spruce AtTT2 homologs (Xue et al., 2003; Arnerup et al., 2011) control secondary metabolism, we created transgenic Norway spruce cell lines overexpressing these two R2R3-MYB TFs and two of the newly identified R2R3-MYB TFs that showed transcriptional activity. We measured the expression levels of genes in the flavonoid biosynthetic pathway as well as the flavonoid, neolignan and stilbene contents of the transgenic lines.

Plants depend on (TFs) to respond and adapt to environmental changes. For this reason, plants possibly possess a higher TF diversity than other living organisms (Shiu et al., 2005). The amount of TFs varies between plant species, but higher plants generally have expanded TF gene families (Lang et al., 2010). Even though previous studies suggested that conifers have a lower number of TFs than angiosperms (Rigault et al., 2011; Canales et al., 2014), the publication of their genomes has shown that the TF number is similar but that TF families evolved differently in conifers, showing expansion or contraction of specific subgroups within TF families (Bedon et al., 2010; Nystedt et al., 2013). In this study we explored the specific diversity of a MYB-bHLH-WDR (MBW) complex in Norway spruce and compare our results with the present knowledge in conifers.

The WD40 protein is portrayed as the ubiquitous member of modular MBW complexes (Ramsay and Glover, 2005; Xu et al., 2015). In Arabidopsis, the WD40 member is represented by the stably expressed single-copy gene AtTTG1, which influences all traits associated with the MBW complex (Tominaga-Wada et al., 2011). Consistent with these reports, we find a single potential homolog of AtTTG1 in Norway spruce, PaWD40-1, which is constitutively expressed in all tissues and in response to various types of abiotic stressors as expected of an AtTTG1 Homolog. The predicted PaWD40-1 protein sequence shows substantial identity and similarity to AtTTG1, particularly in the C-terminal region, which has been shown to be important for the interaction between AtTTG1 with the bHLH TF member of the complex, AtTT8 (Matsui and Ohme-Takagi, 2010). Furthermore, consistent with the conservation of the predicted interaction domain, PaWD40-1 interacts with the tested Norway spruce bHLH proteins included in the yeast two-hybrid assay. We could not detect direct physical interactions between PaWD40-1 and any of the PaMYB proteins and we do not preclude that the proteins may still interact in vivo. This also has been reported in Arabidopsis (Baudry et al., 2004), where direct interactions between TTG1 and TT2 homologs were not observed, and in strawberry (Schaart et al., 2013). Directionality appears to be relatively common in yeast two-hybrid assays involving members of MBW complexes and it has been attributed to partial misfolding possibly producing false negatives (Baudry et al., 2004; Schaart et al., 2013). Taken together, our data suggest that PaWD40-1 is the TTG1 homolog in Norway spruce.

R/B-like subgroup IIIf bHLH proteins participate in MBW complexes and some of these proteins in subgroup IIIf are involved in the control of pigmentation and flavonoid biosynthesis in angiosperms (Park et al., 2007; Park, 2012; Schwinn et al., 2014). Our study gives a first insight into this subgroup of bHLH TFs in conifers. The phylogenetic analyses of our bHLH candidates suggested that the three paralogous bHLH proteins are the closest homologs to AtTT8, as the three bHLH candidates showed a shorter phylogenetic distance to AtTT8 than to the other Arabidopsis bHLH subgroup IIIf members. Further analysis of the predicted amino acid sequences suggested that PabHLH-1 and PabHLH-2 encode functional bHLH proteins and that PabHLH-3 is, most likely, a pseudogene. PabHLH-3 lacks large parts of the ID, AD and bHLH domains, which are essential for TT8 function (Pattanaik et al., 2008; Feller et al., 2011). Thus, we conclude that the Norway spruce genome contains two functional AtTT8 homologs, PabHLH-1 and PabHLH-2. However, the level of sequence divergence between the paralogs and the identification of orthologous sequences to both genes in the white spruce and loblolly pine genomes indicate that it is not a recent duplication but that it predates the divergence between Picea and Pinus, approximately 90–100 Mya (Lu et al., 2014). Even though PabHLH-1 and PabHLH-2 show similar expression patterns in most tissues and in response to abiotic stress, their protein interaction with subgroup 5 R2R3-MYB is largely dissimilar. PabHLH-2 interacts with the MYB proteins of the subgroup S5, PaMYB29, PaMYB31, and with PaMYB32 and PaMYB33, while PabHLH-1 only interacts with PaMYB33 (Figure 5). Thus, the two paralogs appear to have diverged functionally since their separation.

In Arabidopsis, the R2R3-MYB TF subgroup 5 is represented by the single member AtTT2 (Nesi et al., 2001), but in the legume model plant Lotus japonicus there are three members of this subgroup capable of restoring function in tt2 mutants (Yoshida et al., 2008). When co-expressed with LjTT8 and LjTTG1, the three LjTT2s show different activation of late biosynthetic gene transcription in the flavonoid biosynthesis pathway (Yoshida et al., 2010) and this variation in activation strength is associated with substitutions in the amino acid sequences of the LjTT2s. Previously, two different genes with similarity to AtTT2 have been reported from the genus Picea (Xue et al., 2003; Arnerup, 2011), leading us to hypothesize that the conifer R2R3-MYB TF subgroup 5 also would contain multiple members with different functions. We isolated the cDNA sequence from seven Norway spruce R2R3-MYB TFs. Based on phylogenetic analysis, we suggest that six of these genes PaMYB29, PaMYB30, PaMYB31, PaMYB32, PaMYB33, and PaMYB35 are members of, or highly similar to members of, subgroup 5 of the R2R3-MYB TF family. The phylogeny also indicated that the R2R3-MYB TFs have different evolutionary paths; PaMYB29, PaMYB30, and PaMYB31 cluster together with AtTT2 (Nesi et al., 2001) and its homologs in angiosperms (Pazares et al., 1987; Mellway et al., 2009; Terrier et al., 2009) forming the subgroup S5, while PaMYB35 shows similarity to VvMBPA1 and WPS-II SG MYBs of woody angiosperms (Bogs et al., 2007; Soler et al., 2015). However, PaMYB32 and PaMYB33 have no known homologs outside the gymnosperms, indicating that they represent a gymnosperm-specific subgroup of the R2R3-MYB TF family, which we named CS-I. Examples of R2R3-MYB family subgroup expansions in conifers compared to Arabidopsis have already been reported as it is the case for the expanded R2R3-MYB family subgroup 4, which was shown to be involved in the regulation of early flavonoid biosynthesis genes in conifers (Bedon et al., 2010; Bomal et al., 2014). Soler et al. (2015) reported that subgroup 5 is one of three R2R3–MYB TF subgroups, all controlling proanthocyanidin and anthocyanin pathways, which appear to be expanded in woody perennial angiosperms compared to herbaceous angiosperms. Thus, the formation of woody tissues may depend on a more refined regulation of these pathways, and that some of these subgroups may predate the separation of angiosperm gymnosperm lineages while others appear lineage specific.

Grotewold (2005) proposed a model to explain the functional divergence of recently duplicated regulatory genes and how that relates to metabolic diversity in plants. In that model, “the duplication of a regulatory gene for a metabolic pathway is followed by mutations that result in a TF with a partial loss of function” modulating the expression of some but not all genes in a pathway and allowing intermediates to accumulate. This sub-functionalization of the recently duplicated regulatory genes is then followed by its apparent neo-functionalization of metabolic pathways. Thus, under this model, the paralogous R2R3-MYB TFs identified in Norway spruce would regulate genes in the flavonoid pathway differently. To test this, we generated transgenic lines overexpressing PaMYB29, PaMYB32, PaMYB33, and PaMYB35. Overexpression of either of these TFs lead to higher expression levels of PaLAR3 compared to wild type lines. PaLAR3 catalyzes the last step in the flavonoid biosynthesis pathway to produce catechin and it is a key factor in resistance of Norway spruce to fungal pathogens (Hammerbacher et al., 2014; Nemesio-Gorriz et al., 2016). Interestingly, PaMYB32, the Norway spruce ortholog of PmMBF1 (Xue et al., 2003), appeared to also specifically regulate PaLAR3 as well as PaLAR4, while overexpression of its close paralog, PaMYB33, activated the expression of PaPAL1 and PaANR2, PaANR3, PaANR5, PaLAR1, PaLAR3, and PaLAR4 encoding enzymes in the late flavonoid biosynthesis. Together with their contrasting transcriptional responses to hormones and abiotic stress, our observations suggest that the gene duplication that gave rise to PaMYB32 and PaMYB33 was followed by a sub-functionalization of the paralogs, as predicted by the model on functional divergence of recently duplicated regulatory genes (Grotewold, 2005). Similar evolutionary changes have been observed for R2R3-MYB family TFs in several angiosperms (Dias et al., 2003; Chai et al., 2014; Zhao and Bartley, 2014).

The contrasting expression pattern of the gene pair PaMYB32 and PaMYB33 is the most obvious example, in this study, of distinct separation of tissue- or stress-dependent expression patterns between the R2R3–MYB TFs, but all studied genes show individual expression patterns. Such variation in expression patterns is consistent with the proposed modularity of the MBW complex where an exchange, or the tissue specificity, of the R2R3-MYB TFs would specify the regulation of individual traits (Dias et al., 2003; Baudry et al., 2004; Gonzalez, 2009; Xu et al., 2015) in combination with particular bHLH proteins (Ramsay and Glover, 2005). The gene expression profiles of the Norway spruce R2R3-MYB TFs and the expression patterns of late genes in flavonoid biosynthesis in the overexpression lines could indicate that the previously isolated genes PaMYB29 and PaMYB32 (Xue et al., 2003; Arnerup, 2011) control specific functions irrespective of tissue. The two genes were expressed at similar levels in most tested tissues, responded to abiotic and biotic stress, and PaMYB29 and PaMYB32 overexpressing lines showed an upregulation of specific genes in flavonoid biosynthesis. Interestingly enough, the gene expression patterns in the PaMYB29 overexpression lines suggest that this particular TF may contribute to the induction of the genes PaANR3 and PaLAR3 in response to biotic stress. The more generalized late-biosynthesis-gene upregulation seen in PaMYB33 and PaMYB35 overexpression lines and the expression pattern of PaMYB30, PaMYB31, PaMYB33, and PaMYB35, which showed a higher degree of tissue specificity and more restricted responses to abiotic stress compared to PaMYB29 and PaMYB32, suggests that these R2R3-MYB family TFs may perform their regulatory role in particular organs or distinct cell types. Such organ specificity is reported for Arabidopsis' subgroup 7 R2R3-MYB members (Stracke et al., 2007), for paralogous R2R3-MYB genes, PdMYB2 and PdMYB20, in poplar (Chai et al., 2014) and to some extent among conifer subgroup 4 R2R3-MYB family TFs (Bedon et al., 2007, 2010).

Previously, transcriptional profiling studies have shown that overexpression of conifer R2R3-MYBs acting on the shikimate and monolignol pathways can regulate specific sets of genes depending on the MYB that is overexpressed (Bedon et al., 2014). This regulation could potentially channel carbon flux into different directions leading to the accumulation of different metabolites, however such changes in metabolite patterns was not confirmed by metabolite profiling. In this study we attempted to correlate activation of genes in the flavonoid pathway with the soluble flavonoid, stilbene and lignan levels in the overexpression lines. The R2R3-MYB overexpressing lines showed higher expression levels of members of the PaLAR gene family compared to WT controls and overall this was reflected in their chemical profiles. All R2R3-MYB overexpressing lines, except lines overexpressing PaMYB35, had significantly higher catechin levels. This apparent incongruence between gene expression and metabolite accumulation in lines overexpressing PaMYB35 could be an effect of substrate unavailability for particular enzymatic processes, due to the general activation of the late flavonoid biosynthesis pathway, leading to lower metabolite accumulation despite high levels of gene expression (Treutter, 2005). Furthermore, the possibility that the products made by the activated flavonoid biosynthesis pathway were laid down as insoluble lignin or proanthocyanidins affecting the profiles of soluble phenolics cannot be excluded. In fact, Bomal et al. (2008) observed reduced levels of soluble phenolics and an increased deposition of lignin and in PtMYB1 and PtMYB8 overexpressing white spruce lines. Despite the difficulties to find simple correlations of specific specialized metabolites with the overexpression of specific Norway spruce subgroup 5 R2R3-MYB TFs, the general impression is that several of the Norway spruce subgroup 5 R2R3-MYB TFs regulate the activity of the genes in the late flavonoid biosynthesis pathway and the subsequent accumulation flavonoids.

Comparative analyses of gymnosperm and angiosperm TF gene families provided evidence of monophyletic gene family expansions that occurred in gymnosperms after the angiosperm–gymnosperm split (Guillet-Claude et al., 2004; Liu and Ekramoddoullah, 2004; Bedon et al., 2010) and illustrates the different evolutionary trajectories of the gene families in seed plants (Bedon et al., 2010). Our study on selected members of the MBW complex in Norway spruce further adds to this picture as we identified a bHLH gene pair orthologous to AtTT8 and several members of the R2R3-MYB TF family subgroup 5 in the conifer, Norway spruce. Taken together, previous studies on conifer MBW members and this study indicate that conifer Subgroups 4 and 5 control early and late gene expression in the flavonoid pathway as well as stilbene and neolignan biosynthetic pathways. To fulfill these functions, it is likely that they form complexes with bHLH and WD40 proteins as suggested by Zhao et al. (2013). Testing interactions of PabHLH-1 and PabHLH-2 with members of R2R3-MYB subgroups 4 and 5 would shed light on the involvement of the bHLH paralogs in the control of the different biosynthetic modules. Furthermore, overexpression of members of R2R3-MYB under inducible promoters or silencing expression by RNA-interference might lead to a clearer elucidation of their functions.

MN and ME conceived the study. MN isolated the candidate genes, performed the phylogenetic analyses, prepared the vectors for yeast and spruce transformation, extracted the RNA and performed the expression studies on plant tissues and cell lines. PB performed the protein interaction experiment and analyzed the data together with SM. JA performed the expression study to test biotic and abiotic stress on spruce plants. KD conducted the experiment to test abiotic stress on spruce cell lines. AH performed the chemical analysis of the cell lines. MN drafted the manuscript and coordinated the writing. PB, KD, AH, JS, SM, and ME contributed with advice on data analyses and advice during the process of writing the manuscript. MN wrote the final manuscript and all co-authors read and approved the final version of the manuscript.