# PLANT SECONDARY COMPOUNDS IN FOREST ECOSYSTEMS UNDER GLOBAL CHANGE: FROM DEFENSE TO CARBON SEQUESTRATION

EDITED BY : Bartosz Adamczyk and Judy Simon PUBLISHED IN : Frontiers in Plant Science

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88963-010-3 DOI 10.3389/978-2-88963-010-3

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# PLANT SECONDARY COMPOUNDS IN FOREST ECOSYSTEMS UNDER GLOBAL CHANGE: FROM DEFENSE TO CARBON SEQUESTRATION

Topic Editors: Bartosz Adamczyk, Natural Resources Institute Finland (LUKE), Finland Judy Simon, University of Konstanz, Germany

Citation: Adamczyk, B., Simon, J., eds. (2019). Plant Secondary Compounds in Forest Ecosystems Under Global Change: From Defense to Carbon Sequestration. Lausanne: Frontiers Media. doi: 10.3389/978-2-88963-010-3

# Table of Contents


Jarmo K. Holopainen, Virpi Virjamo, Rajendra P. Ghimire, James D. Blande, Riitta Julkunen-Tiitto and Minna Kivimäenpää

*16 Interannual and Seasonal Dynamics of Volatile Organic Compound Fluxes From the Boreal Forest Floor*

Mari Mäki, Juho Aalto, Heidi Hellén, Mari Pihlatie and Jaana Bäck

*30 Combinations of Abiotic Factors Differentially Alter Production of Plant Secondary Metabolites in Five Woody Plant Species in the Boreal-Temperate Transition Zone*

John L. Berini, Stephen A. Brockman, Adrian D. Hegeman, Peter B. Reich, Ranjan Muthukrishnan, Rebecca A. Montgomery and James D. Forester

*47 Fertilization Changes Chemical Defense in Needles of Mature Norway Spruce (*Picea abies*)*

Line Nybakken, Marit H. Lie, Riitta Julkunen-Tiitto, Johan Asplund and Mikael Ohlson

*56 Flavanone-3-Hydroxylase Plays an Important Role in the Biosynthesis of Spruce Phenolic Defenses Against Bark Beetles and Their Fungal Associates*

Almuth Hammerbacher, Dineshkumar Kandasamy, Chhana Ullah, Axel Schmidt, Louwrance P. Wright and Jonathan Gershenzon


Mette Sørensen, Federico Cozzi, Berin A. Boughton, Allison Maree Heskes and Elizabeth Heather Jakobsen Neilson

# Editorial: Plant Secondary Compounds in Forest Ecosystems Under Global Change: From Defense to Carbon Sequestration

Judy Simon<sup>1</sup> \* and Bartosz Adamczyk <sup>2</sup>

*<sup>1</sup> Plant Interactions Ecophysiology Group, Department of Biology, University of Konstanz, Konstanz, Germany, <sup>2</sup> Natural Resources Institute Finland (LUKE), Helsinki, Finland*

Keywords: climate change, boreal, temperate, eucalypts, plant secondary metabolites, herbivory, phenolics, volatiles

**Editorial on the Research Topic**

#### **Plant Secondary Compounds in Forest Ecosystems Under Global Change: From Defense to Carbon Sequestration**

Forests are crucial for the future of human society, not only by their provision of fundamental ecosystems services, such as timber, food, and recreation, but also for mitigating climate change via the capture and long-term storage of carbon (MEA, 2005; Sabatini et al., 2019). At the ecosystem level, plant secondary compounds (PSCs) drive key processes, such as litter degradation, and thus nutrient cycling (Chomel et al., 2016; Adamczyk et al., 2018). At the plant level, PSCs are employed to counter herbivore attacks, cope with changes in the environment, and ultimately ensure plant survival. PSCs are chemically diverse and their abundance shifts depending on biotic and abiotic stressors, but also across species and even organs within a species. The research highlighted in this collection combines short- and long-term studies in the greenhouse and in the field with a focus on phenolics and terpenes, both at the tree and ecosystem level.

Climate change factors, such as elevated atmospheric CO<sup>2</sup> levels and/or enhanced temperature might lead to counteracting responses of trees and forest ecosystems. Holopainen et al. reviewed how forest ecosystems change depending on the influence of abiotic and biotic factors on PSC synthesis. Different aspects of climate change show counteracting effects on PSCs: for example, in response to rising CO<sup>2</sup> total phenolics increased and terpenes decreased in the leaves, whereas warming showed the opposite effect (Holopainen et al.). According to Mäki et al. the temporal dynamics of PSCs such as monoterpenes in forest ecosystems also depend on temperature and rainfall and are linked to the degradation of litter. Berini et al. used an untargeted metabolomics approach for five tree species and concluded that effects of the environment are species-specific and differ depending on the abiotic factors in question. However, indirect consequences of measures to mitigate global change are likely to cause further problems. For example, fertilization with nitrogen to increase carbon sequestration might result in increased insect herbivory in the future. The phenolic levels of mature Picea abies trees were reduced after 13 years of nitrogen deposition which—in combination with higher nitrogen content in the needles—might severely decrease tree resistance to herbivory (Nybakken et al.). In contrast, in Picea abies, levels of the PSCs dihydroflavonol and flavan-3-ol increased after inoculation with a fungus related to bark beetle attack and negatively affected both bark beetle and fungus (Hammerbacher et al.).

The allocation of resources to PSCs is however not just regulated by environmental factors, but also by inherent aspects. For example, the genetic variation of Pinus pinaster in their

#### Edited and reviewed by:

*Henry D. Adams, Oklahoma State University, United States*

\*Correspondence: *Judy Simon judy.simon@uni-konstanz.de*

#### Specialty section:

*This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science*

Received: *23 May 2019* Accepted: *07 June 2019* Published: *26 June 2019*

#### Citation:

*Simon J and Adamczyk B (2019) Editorial: Plant Secondary Compounds in Forest Ecosystems Under Global Change: From Defense to Carbon Sequestration. Front. Plant Sci. 10:831. doi: 10.3389/fpls.2019.00831*

**4**

constitutive and induced resistance against herbivory was explained by PSC profiles rather than single compounds (López-Goldar et al.). Stark and Martz investigated whether dioecious Juniperus communis shows a gender specific response to a reduction in biomass because females might allocate more resources to defense (i.e., PSCs), whereas males invest more into growth. However, genders did not differ significantly in their PSC levels, most likely due to a high demand of PSCs to adapt to abiotic stresses (Stark and Martz). Moreover, Marques dos Santos et al. used a new method combining UHPLC, DAD, and MS/MS analyses to show that formylated phloroglucinols were found in glands embedded subdermally in leaves, flowers, and flower buds in eight eucalypts and one Corymbia species.

In conclusion, PSCs are important drivers of ecosystem processes and directly affect the interactions between trees and other organisms as well as the environment. Considering that most studies investigate tree seedlings or only the short-term

#### REFERENCES


effects on PSCs (Holopainen et al.), the new knowledge gained from the studies presented in this research topic is a further step in understanding the responses of trees to biotic and abiotic stressors and the underlying mechanisms at the tree and forest level, and ultimately in developing and evaluating sustainable management strategies to adapt forests to global change.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

# ACKNOWLEDGMENTS

We would like to thank the authors, reviewers, and the Frontiers in Plant Science Editorial Office for their support in creating this research topic.

biodiversity in European temperate forests. Global Change Biol. 25, 536–548. doi: 10.1111/gcb.14503

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

Copyright © 2019 Simon and Adamczyk. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Climate Change Effects on Secondary Compounds of Forest Trees in the Northern Hemisphere

Jarmo K. Holopainen<sup>1</sup> \*, Virpi Virjamo<sup>2</sup> , Rajendra P. Ghimire<sup>1</sup> , James D. Blande<sup>1</sup> , Riitta Julkunen-Tiitto<sup>2</sup> and Minna Kivimäenpää<sup>1</sup>

<sup>1</sup> Department of Environmental and Biological Sciences, Kuopio Campus, University of Eastern Finland, Kuopio, Finland, <sup>2</sup> Department of Environmental and Biological Sciences, Joensuu Campus, University of Eastern Finland, Joensuu, Finland

Plant secondary compounds (PSCs), also called secondary metabolites, have high chemical and structural diversity and appear as non-volatile or volatile compounds. These compounds may have evolved to have specific physiological and ecological functions in the adaptation of plants to their growth environment. PSCs are produced by several metabolic pathways and many PSCs are specific for a few plant genera or families. In forest ecosystems, full-grown trees constitute the majority of plant biomass and are thus capable of producing significant amounts of PSCs. We summarize older literature and review recent progress in understanding the effects of abiotic and biotic factors on PSC production of forest trees and PSC behavior in forest ecosystems. The roles of different PSCs under stress and their important role in protecting plants against abiotic and biotic factors are also discussed. There was strong evidence that major climate change factors, CO<sup>2</sup> and warming, have contradictory effects on the main PSC groups. CO<sup>2</sup> increases phenolic compounds in foliage, but limits terpenoids in foliage and emissions. Warming decreases phenolic compounds in foliage but increases terpenoids in foliage and emissions. Other abiotic stresses have more variable effects. PSCs may help trees to adapt to a changing climate and to pressure from current and invasive pests and pathogens. Indirect adaptation comes via the effects of PSCs on soil chemistry and nutrient cycling, the formation of cloud condensation nuclei from tree volatiles and by CO<sup>2</sup> sequestration into PSCs in the wood of living and dead forest trees.

#### Keywords: CO2, drought, ozone, phenolics, temperature, terpenes, UV-B, VOCs

#### INTRODUCTION

Plant secondary compounds (PSCs) have high chemical diversity with an estimated 200,000 compounds (Dixon and Strack, 2003). In higher plants, terpenoids (approximately 30,000 known compounds, Lämke and Unsicker, 2018), alkaloids (21,000) (Wink, 2010), and phenolic compounds (8,000) (Munné-Bosch, 2012) are the most diverse PSC groups. The majority of PSCs are related to plant chemical defense and herbivore pressure, which are considered the main drivers of PSC diversification (Lämke and Unsicker, 2018). Terpenoids originate directly from glycolysis of glucose via the mevalonate pathway or the methylerythritol phosphate (MEP) pathway (Bartram et al., 2006). Phenolic compounds originate from the shikimic acid pathway, which is related to the metabolism of carbohydrates and aromatic amino acids (Seigler, 1998; Lindroth, 2012).

#### Edited by:

Bartosz Adamczyk, University of Helsinki, Finland

#### Reviewed by:

Carsten Kulheim, Michigan Technological University, United States Ryo Nakabayashi, RIKEN, Japan

> \*Correspondence: Jarmo K. Holopainen jarmo.holopainen@uef.fi

#### Specialty section:

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

Received: 09 May 2018 Accepted: 12 September 2018 Published: 02 October 2018

#### Citation:

Holopainen JK, Virjamo V, Ghimire RP, Blande JD, Julkunen-Tiitto R and Kivimäenpää M (2018) Climate Change Effects on Secondary Compounds of Forest Trees in the Northern Hemisphere. Front. Plant Sci. 9:1445. doi: 10.3389/fpls.2018.01445

**6**

Most of the alkaloids are derived from amino acid precursors (Wink, 2010), but coniferous alkaloids are synthesized via the polyketide pathway (Seigler, 1998). Terpenes and phenolics are the most extensively studied PSCs in forest trees (Lindroth, 2012; Lämke and Unsicker, 2018), while alkaloids have been considered in relatively few studies (Virjamo et al., 2014). Alkaloids and phenolics may compete more with protein synthesis than terpenoids (Koricheva, 2002). Up to 10% of recently fixed carbon can be allocated to volatile organic compounds (VOCs) in stressed plants (Peñuelas and Staudt, 2010), while stored terpenoids typically require secretory structures such as resin canals and glandular trichomes to avoid autotoxicity (Goodger et al., 2013).

At the global level, climate change (CC) is strongly linked to increasing CO<sup>2</sup> concentration in the atmosphere, which affects temperature globally (IPCC, 2014). Depletion of the stratospheric ozone (O3) layer is partly affected by elevated CO<sup>2</sup> and leads to increased UV-B radiation (IPCC, 2014). In the lower atmosphere (troposphere) O<sup>3</sup> acts as a greenhouse gas and is phytotoxic to plants (Vapaavuori et al., 2009; Lindroth, 2010). At the regional scale, global warming affects arctic and boreal areas most strongly and will include drastic changes, e.g., in precipitation, drought, and cloudiness (shading), affecting temperature-dependent ecosystem processes in forests such as nutrient availability (IPCC, 2014).

Climate change directly affects the abiotic conditions of forest trees and their growth, physiology and defense including induced PSC production (Metlen et al., 2009; **Figure 1**). The responses of other organisms, such as microbial plant pathogens, herbivorous insects and arrival of exotic invasive herbivores and pathogens will alter the biotic stress burden of trees (Donnelly et al., 2012). The richness of chemical compounds originating from PSCs in forest ecosystems is affected by PSC reactivity (Kundu et al., 2012). Many of the volatile reaction products form secondary organic aerosols (SOA), which may affect cloud formation (Zhao et al., 2017; Rose et al., 2018; Scott et al., 2018) and provide important ecosystem – atmosphere feedback, thus helping vegetation to adapt to CC (e.g., Niinemets, 2018). SOA particles may become wet or dry deposited on forest vegetation and accumulate in soil (Holopainen et al., 2017) together with litter PSCs (Smolander et al., 2012).

We summarize the major results of literature reviews and meta-analyses on the quantitative responses of tree PSCs to specific CC-related abiotic factors, with restriction to literature up to 2015 in **Figure 2**. Recent results of PSC (after 2015) and plant VOC (after 2012) research are then reviewed. We then discuss the roles of multiple CC-related stresses on PSCs of forest trees and the potential of non-volatile and volatile PSCs to modulate plant adaptation to CC and their involvement in climate feedback.

#### SUMMARY OF OLDER REVIEWS

Terpenoids and phenolic compounds in leaves and terpenoids in emissions are the most frequently studied PSCs of forest trees in response to CC-related factors (**Figure 2**). Major CC factors such as elevated CO<sup>2</sup> and warming have the most distinctive, and often contrasting effects on PSCs. Greater availability of CO<sup>2</sup> increases phenolics in foliage and reduces terpenoids both in foliage and in emissions, while warming reduces phenolics in foliage and increases terpenoids in foliage and emissions (Zvereva and Kozlov, 2006; Peñuelas and Staudt, 2010). CO2+warming increased phenolics in foliage, but reduced phenolics in woody tissues of woody plants (Zvereva and Kozlov, 2006). O<sup>3</sup> effects on foliar concentrations of terpenes and volatile emissions have been variable, but increased phenolic compounds have been found. One reason for variable responses of terpenes could be withinspecies variation in O<sup>3</sup> sensitivity (Valkama et al., 2007) and variable O<sup>3</sup> reactivity of terpenes (Blande et al., 2014). Major terpenoid groups such as monoterpenes (MT) and sesquiterpenes (SQT) in foliage have been less studied separately under most of the other CC related stresses.

#### RESPONSES TO ABIOTIC FACTORS – PSC CONCENTRATIONS IN TREES

Certain plant PSC groups such as flavonoids (Lavola et al., 2013) and MTs (Semiz et al., 2007) are under strong genotypic control. This may increase the potential for forest tree chemical defenses to evolutionarily adapt to environmental changes (Lavola et al., 2013). However, high within-species variation in PSCs is problematic for short-term experimental set-ups and may mask the effects of environmental factors on PSC concentrations and emissions.

## Elevated CO<sup>2</sup> and Warming

Total leaf phenolics, particularly condensed tannins (Julkunen-Tiitto et al., 2015), follow the traditional carbon-nutrient balance hypothesis (Bryant et al., 1983), which predicts that photosynthesised carbon will be diverted to carbon-rich PSCs, if photosynthesis or CO<sup>2</sup> availability are at a high level and nutrient availability limits carbon allocation to plant growth (Koricheva et al., 1998). However, many flavonoids and other phenolics, and particularly terpenoids, do not support that hypothesis (Lindroth, 2012). Recent studies have confirmed a trend with foliar phenolics; doubled ambient CO<sup>2</sup> increased salicylates and phenolic acids (Sobuj et al., 2018), anthocyanins and flavonoids (Vanzo et al., 2015; Nissinen et al., 2016), while less than doubled CO<sup>2</sup> did not affect phenolics (McKiernan et al., 2012; Randriamanana et al., 2018). Furthermore, some phenolic groups, such as phenolic glycosides, are not as responsive (Jamieson et al., 2017). In Eucalyptus spp., there was no observed effect of elevated CO<sup>2</sup> (McKiernan et al., 2012) or elevated CO<sup>2</sup> decreased (Bustos-Segura et al., 2017) stored foliar terpenes.

Reduction in foliar phenolics at elevated temperatures (Julkunen-Tiitto et al., 2015) has been explained by a dilution effect whereby more carbon was allocated to growth supporting structures with turnover of phenolics slow or even absent (Kosonen et al., 2012). Recently, it has been reported that warming reduces flavonoids (Nissinen et al., 2016; Zhang et al., 2018), salicylates and phenolic acids (Sivadasan et al., 2018) in deciduous trees and total phenolics in a conifer (Zhang et al., 2018). Flavonoids, alkaloids, and saponins (triterpenes) were

increased in Robinia pseudoacacia grown in warmer open-top chambers compared to seedlings grown in an open field (Zhao et al., 2016). However, enclosure may significantly affect foliar PSCs (Peltonen et al., 2005).

Warming induced increase in foliar terpenes (Valkama et al., 2007; Julkunen-Tiitto et al., 2015) is possibly explained by improved plant growth and larger storage space for stored terpenes such as MTs and resin acids (Sallas et al., 2003). Recent results showing higher MT concentrations in xylem resin (Ferrenberg et al., 2017) and increased number of needle resin canals (Kivimäenpää et al., 2016) in Pinus spp. in warmer environments supports this view. Warming (+2 ◦C) increased concentrations of Picea abies needle alkaloids (Virjamo et al., 2014) and parental temperature range correlated with conifer alkaloid concentrations (Gerson et al., 2009; Virjamo and Julkunen-Tiitto, 2016).

# O<sup>3</sup>

O<sup>3</sup> is phytotoxic and reduces plant photosynthesis and growth (Bueker et al., 2015; Li et al., 2017) and induces substantial changes in protein activity, gene expression, signaling pathways, and metabolism even before any tissue damage can be detected (Bueker et al., 2015; Vainonen and Kangasjärvi, 2015). Activated PSC metabolism (increased biosynthesis and reduced turnover) explains the increased concentrations of phenolic and terpenoid PSCs in tree foliage found in a meta-analysis (Valkama et al., 2007). In response to O<sup>3</sup> treatment, elevated levels of total phenolics (Gao et al., 2016), condensed tannins (Couture et al., 2017) and flavonoids (Cotrozzi et al., 2018) were found in foliage of deciduous trees, but this was not the case in needles of coniferous trees (Riikonen et al., 2012).

# Drought

Drought in warm periods is expected to become the most frequent widespread climatic extreme that negatively affects terrestrial ecosystems (Schwalm et al., 2017). Forests in southern Europe will face stronger summer drought than those in the north (Ruosteenoja et al., 2018). The observed trends in PSCs are not strong, although reduction of some terpenes and phenolics have been reported (Koricheva et al., 1998; Julkunen-Tiitto et al., 2015). In Quercus ilex leaves, moderate drought periods resulted in higher concentrations of total polyphenolics (Nogues et al., 2014; Rivas-Ubach et al., 2014). Drought stress did not affect alkaloid concentrations in conifers (Gerson and Kelsey, 2004). Increasing severity of drought increased resin canal density and resin acid content in the stem wood of Pinus sylvestris (Heijari et al., 2010). Contents of most MTs and sesquiterpenoids in needles increased with increasing severity of water deficit, but


FIGURE 2 | Summary of the results in comprehensive reviews and meta-analyses of plant PSCs responses in their concentration or emissions in forest trees under climate change related stresses. This tabulation covers nearly 400 original scientific articles. Type of stresses: O<sup>3</sup> = elevated ozone, CO<sup>2</sup> = elevated carbon dioxide, UV-B = elevated UV-B radiation. The tabulation gives the direction of PSC responses to single environmental factors, except for the combined O3+CO<sup>2</sup> and CO2+warming effects. Arrow direction (increase or decrease) shows the effect size (significant or highly significant effect in meta-analysis) or the amount of evidence in reviews. Horizontal bidirectional arrow indicates the studies with non-significant results and up and down arrows in the same cell indicate significant results in both directions dominate. Question mark indicates that the reviews and meta-analyses did not provide enough data. Gray arrows in CO2, warming and CO2+warming columns indicate woody tissue responses as specified in reference 2. Small numbers in the upper-right corner of reach cell indicate the source reference: 1 = Koricheva et al., 1998; 2 = Zvereva and Kozlov, 2006; 3 = Valkama et al., 2007; 4 = Peñuelas and Staudt, 2010; 5 = Lindroth, 2012; 6 = Julkunen-Tiitto et al., 2015; 7 = Robinson et al., 2012; 8 = Li et al., 2017. References 1–3 and 7–8 are meta-analyses and 4–6 literature reviews. Meta-analyses give the size of the effect, and the reviews give the frequency of observations in each response categories (–, 0, +).

resin acids showed the opposite trend and polyphenols did not respond (Sancho-Knapik et al., 2017). Drought slightly increased the needle MT contents of the coastal, but not of the interior Pseudotsuga menziesii provenance (Kleiber et al., 2017).

#### Shading

Shading limits photosynthesis and it is predicted that less carbon is allocated to PSCs in shaded conditions (Koricheva et al., 1998). Giertych et al. (2015) tested the hypothesis with light-demanding and shade-tolerant deciduous tree species. All species had higher total phenolic content of leaves in sunny conditions, but shade tolerant species invested more carbon into growth when more light was available.

#### UV-B

Short-wave (280–315 nm) UV-B radiation induces tree foliage to produce more phenolic acids and flavonoids as protective pigments (Julkunen-Tiitto et al., 2015; Kaling et al., 2015). In Populus tremula seedlings, two flavonoid compounds increased under enhanced UV-B radiation (Nissinen et al., 2017). Filtering out solar UV-B radiation from reaching P. tremula in mountains lead to lower flavonoid concentration, but increased salicin in leaves, and salicylates and total phenolics in stem (Stromme et al., 2018). However, a three-year exposure of Salix myrsinifolia to 32% elevated UV-B in a field site did not result in significant changes in flavonoid concentrations, and in male trees total salicylates were even reduced (Ruuhola et al., 2018). Enhanced levels of UV-B radiation did not affect alkaloid concentrations in conifers (Virjamo et al., 2014).

#### RESPONSES TO ABIOTIC FACTORS – EMISSION OF VOLATILE PSCs

#### Warming

Warming is one of the most intensively studied factors in volatile PSC studies (Peñuelas and Staudt, 2010), due to VOC emissions being strongly light and temperature dependent and some, e.g., isoprene not being stored in leaves (Sharkey and Yeh, 2001). At the global scale, isoprene represents nearly 50% of the estimated biogenic VOC emission, and emissions of MTs and SQTs together comprise about 15 and 3%, respectively (Guenther et al., 2012). Warming increased emissions of volatile terpenoids most consistently (Peñuelas and Staudt, 2010) by stimulating biosynthesis (Loreto and Schnitzler, 2010; Fini et al., 2017) and temperature-dependent release from leaf storage structures (Grote and Niinemets, 2008; Peñuelas and Staudt, 2010). VOC emissions from storage structures are high in conifers, e.g., over

40% in Pinus, and very low in deciduous trees (Ghirardo et al., 2010). Differences in emission responses of MTs and SQTs to temperature might be caused by temporal storage of less volatile SQTs in leaf surface waxes (Joensuu et al., 2016), and their temperature-dependent adherence and release (Himanen et al., 2015).

Warming of +3 ◦C increased (∼70%) emission rates of isoprene-emitting and MT-emitting Quercus spp. (Staudt et al., 2017). Isoprene-emitting Populus spp. have an emission maximum at +45◦C even if grown in elevated CO<sup>2</sup> (Potosnak et al., 2014; Niinemets and Sun, 2015). Long-term warming of 1 ◦C+ambient increased MT (fivefold), SQT (fourfold) and green leaf volatile (40%) emissions of Betula pendula (Hartikainen et al., 2012). Ontogeny of needle development in terpene-storing conifers affects variability in MT emission rate responses to warming (Esposito et al., 2016; Kivimäenpää et al., 2016; Ghimire et al., 2017; Tiiva et al., 2018). Some specific MTs of forest trees, particularly β-ocimene are known to be heat-stress indicators (Jardine et al., 2017).

# CO<sup>2</sup>

There is clear evidence that exposure to elevated CO<sup>2</sup> reduces terpenoid-based plant VOC emissions (Peñuelas and Staudt, 2010; Potosnak et al., 2014; Niinemets and Sun, 2015; Mochizuki et al., 2017) and uncouples isoprene emissions from photosynthesis. Inhibition of isoprene synthesis becomes stronger with increasing CO<sup>2</sup> concentration (Wilkinson et al., 2009; Rasulov et al., 2018) and isoprene is more responsive than MTs (Peñuelas and Staudt, 2010). There is still no consensus of the mechanisms that lead to the contrasting responses of photosynthesis and isoprene biosynthesis (Fini et al., 2017) and emission at elevated CO<sup>2</sup> (Rasulov et al., 2018). In MTstoring gymnosperms the reduced MT needle concentration at elevated CO<sup>2</sup> (Sallas et al., 2003) may partly explain reduced MT emissions.

# O<sup>3</sup>

Plant VOCs easily react with O<sup>3</sup> and samples may be destroyed if ozone scrubbers are not used in sampling (Blande et al., 2014). O<sup>3</sup> exposure experiments have shown that plant VOC responses are variable and depend on the O<sup>3</sup> concentration (Peñuelas and Staudt, 2010). Acute ozone exposure can induce similar VOCs to herbivore damage (Blande et al., 2014). Moderately elevated O<sup>3</sup> concentrations increase emissions of MTs in coniferous (Kivimäenpää et al., 2013, 2016; Ghimire et al., 2017) and deciduous trees (Carriero et al., 2016), but decrease isoprene emission in poplar (Yuan et al., 2016) and MT emissions in Larix were not affected (Mochizuki et al., 2017).

#### Drought

Warm and dry periods lead to drought stress; mild drought may increase VOC emissions (Peñuelas and Staudt, 2010; Bourtsoukidis et al., 2014; Yuan et al., 2016), not affect (Nogues et al., 2014) or have increasing-decreasing trend (Simpraga et al., 2011), while severe drought results in a reduction of emissions (Peñuelas and Staudt, 2010; Rodriguez-Calcerrada et al., 2013; Luepke et al., 2016; Marino et al., 2017; Saunier et al., 2017; Staudt et al., 2017). Drought stimulated herbivore-induced GLV and MT emissions in Alnus glutinosa (Copolovici et al., 2014).

#### Shading

Shading by branch position can reduce isoprene and MT emissions from forest trees (Esposito et al., 2016; Juran et al., 2017; van Meeningen et al., 2017). This suggests that cloudiness that reduces penetration of solar radiation to the canopy (Niinemets, 2018) may reduce carbon allocation to volatile PSCs. However, in the canopy of Fagus sylvatica, MT emissions were highest in the semi-shaded leaves (Simpraga et al., 2013).

## COMBINED ABIOTIC FACTORS VERSUS NON-VOLATILE AND VOLATILE PSCs

The effects of eutrophication [nitrogen (N) deposition and forest fertilization] on tree VOCs are poorly known, although it appears to slightly increase isoprene emissions (**Figure 2**; Peñuelas and Staudt, 2010). Generally, exposure to additional abiotic factors may mitigate, intensify or not affect foliar chemical defense (Sallas et al., 2001; Zvereva and Kozlov, 2006; Vapaavuori et al., 2009) or VOC (Kivimäenpää et al., 2016; Ghimire et al., 2017; Tiiva et al., 2017) responses to the first factor. For example elevated O<sup>3</sup> and N availability had additive effects on terpenoid emissions from a deciduous tree (Carriero et al., 2016) and conifer (Kivimäenpää et al., 2016), whereas N in addition to warming suppressed the emissions of non-oxygenated MTs from a conifer (Ghimire et al., 2017). In Populus sp. N addition (Yuan et al., 2017) or drought (Yuan et al., 2016) did not mitigate the negative effects of O<sup>3</sup> on isoprene emission.

An elevated CO2-related decrease of terpenoids in conifer needles was mitigated by elevated temperature and CO2-induced increase in phenolics was intensified by elevated temperature (Zvereva and Kozlov, 2006). Furthermore, elevated CO<sup>2</sup> may mitigate a typical O<sup>3</sup> response such as ozone-induced increase of some phenolics in B. pendula (Vapaavuori et al., 2009), but in Larix CO<sup>2</sup> and O<sup>3</sup> did not have interactive effects on volatiles (Mochizuki, et al., 2017). Along altitudinal gradients, the total phenolics in Quercus robur foliage were found to decrease at lower elevations where temperature was warmer (Abdala-Roberts et al., 2016). Our analysis shows that in addition to a warminginduced reduction in total phenolics, elevated UV-B radiation may increase flavonoids and phenolic acids in foliage (Julkunen-Tiitto et al., 2015). However, recent combined warming and UV-B exposure outdoors has confirmed that warming has a stronger increasing effect on leaf phenolics (Nissinen et al., 2017) and terpenoid emissions (Maja et al., 2016) than elevated UV-B radiation.

# PLANT PSCs AND BIOTIC STRESS

Abiotic stresses affect growth, reproduction and behavior of herbivorous insects and severity of microbial plant pathogens and their biotic stress on trees (**Figure 1**). Furthermore, arrival of exotic invasive herbivores and pathogens with warming

related to CC may add another threat, because local tree species have not adapted to use PSCs to defend against these invasive species (Donnelly et al., 2012). Increased accumulation of PSCs in plant tissues under elevated atmospheric CO<sup>2</sup> resulted in compensation feeding and more severe defoliation (Docherty et al., 1996). Warming may reduce constitutive PSCs in foliage (Stark et al., 2015) and improve (Grainger and Gilbert, 2017) or reduce (Ghimire et al., 2017) insect performance and lead to induced production of non-volatile (Rubert-Nason et al., 2015) or volatile PSCs (Blande et al., 2010; Holopainen and Gershenzon, 2010).

#### PSC TRANSMITTED ECOSYSTEM FEEDBACKS TO CLIMATE CHANGE

Forests are a globally important sink of atmospheric CO<sup>2</sup> (Pan et al., 2011; Pukkala, 2018). Carbon is sequestered in soil carbon (44%), living biomass (42%), deadwood (8%), and litter (5%) (Pan et al., 2011). Reducing forest harvesting and adding dead wood as large trees is the fastest way to increase CO<sup>2</sup> sequestration in boreal forests (Pukkala, 2018). Reduced decay rate of dead wood decreases release of CO<sup>2</sup> back to the atmosphere. Concentrations of both resin acids (Karppanen et al., 2007) and stilbenes (Ioannidis et al., 2017) in heartwood of Pinus sp. can exceed 10% of dry weight and significantly reduce decay rate of dead wood by fungal rot (Karppanen et al., 2007) and by woodborers (Nerg et al., 2004).

The PSCs are released into the soil (Kainulainen et al., 2003) and atmosphere (Aaltonen et al., 2011) faster from the needle litter than dead wood. In conifer forests, MTs contribute over 90% of litter emissions peaking in early summer and during the needle fall in autumn (Aaltonen et al., 2011). After 19-months of decomposition, pine needles have lost 96, 79, and 86 percent of initial concentrations of MTs, resin acids and total phenolics, respectively (Kainulainen et al., 2003). Warm periods in the boreal conifer forests increase reactive MT emissions to the forest atmosphere and stimulate formation of SOA, particularly in evenings (Rose et al., 2018) and potentially stimulate cloudiness (Joutsensaari et al., 2015; Zhao et al., 2017), thereby reducing penetration of solar radiation to vegetation and mitigating warming and drought stress (Niinemets, 2018).

Volatile PSCs have several functions in the species interactions affecting herbivores and their natural enemies and activating defenses in healthy plants (Blande et al., 2014). CC factors also affect herbivore-induced volatiles, but studies on forest trees and woody plants remain limited (Holopainen and Gershenzon,

#### REFERENCES


2010). Non-volatile PSCs released from plant leaf and needle litter accumulate in soil and affect C or net N mineralization in forest soils (e.g., Smolander et al., 2012). Growing three seasons at elevated CO<sup>2</sup> and O<sup>3</sup> did not affect concentrations of total phenolics and MTs in P. sylvestris needles, but at elevated O<sup>3</sup> resin acids were increased. However, during decomposition this difference in needle litter was lost (Kainulainen et al., 2003). In leaf litter of B. pendula growing at elevated CO<sup>2</sup> and O<sup>3</sup> concentrations, several phenolic compounds were elevated, and the growth of juvenile earthworms was reduced (Kasurinen et al., 2007) suggesting that CC reduced quality of food for soil animals. In Alnus incana litter, UV-A and UV-B exclusion affected concentrations of phenolic groups variably, whereas in birch litter there were no significant differences in phenolic compounds (Kotilainen et al., 2009).

#### CONCLUDING REMARKS

There is strong evidence that CO<sup>2</sup> increases phenolic compounds in foliage, but limits terpenoids in foliage and emissions, while warming decreases phenolic compounds in foliage, but increases terpenoids in foliage and emissions. CO<sup>2</sup> in combination with other abiotic stresses has more variable effects. As current PSC information is mostly from tree seedlings, long-term experiments covering the whole ontogeny of trees in forest ecosystems are needed to better understand CC effects on PSCs in trees. Assessment of the role of volatile PSCs in sun-screening and nonvolatile PSCs in control of decay rate of carbon storage in dead wood will be necessary to improve our knowledge of the capacity of PSCs to mitigate CC. Transgenic trees with silenced or overexpressed genes for key volatile PSCs (Kaling et al., 2015; Vanzo et al., 2015) will be important tools in this task.

#### AUTHOR CONTRIBUTIONS

JH drafted the manuscript. All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

This work was funded by University of Eastern Finland BORFOR Top Research Area and Academy of Finland, projects nos. 278424, 267360, 305343, and 309425.

oak species. Am. J. Bot. 103, 2070–2078. doi: 10.3732/ajb.16 00310






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

Copyright © 2018 Holopainen, Virjamo, Ghimire, Blande, Julkunen-Tiitto and Kivimäenpää. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Interannual and Seasonal Dynamics of Volatile Organic Compound Fluxes From the Boreal Forest Floor

Mari Mäki1,2 \*, Juho Aalto1,3, Heidi Hellén<sup>4</sup> , Mari Pihlatie1,5,6 and Jaana Bäck1,2,6

1 Institute for Atmospheric and Earth System Research/Forest Sciences, Helsinki, Finland, <sup>2</sup> Department of Forest Sciences, Faculty of Agriculture and Forestry, University of Helsinki, Helsinki, Finland, <sup>3</sup> Department of Physics, Faculty of Science, University of Helsinki, Helsinki, Finland, <sup>4</sup> Finnish Meteorological Institute, Helsinki, Finland, <sup>5</sup> Department of Agricultural Sciences, Faculty of Agriculture and Forestry, University of Helsinki, Helsinki, Finland, <sup>6</sup> Viikki Plant Science Centre, University of Helsinki, Helsinki, Finland

#### Edited by:

Judy Simon, Universität Konstanz, Germany

#### Reviewed by:

Jürgen Kreuzwieser, University of Freiburg, Germany Ivika Ostonen, University of Tartu, Estonia

> \*Correspondence: Mari Mäki mari.maki@helsinki.fi

#### Specialty section:

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

Received: 26 October 2018 Accepted: 05 February 2019 Published: 22 February 2019

#### Citation:

Mäki M, Aalto J, Hellén H, Pihlatie M and Bäck J (2019) Interannual and Seasonal Dynamics of Volatile Organic Compound Fluxes From the Boreal Forest Floor. Front. Plant Sci. 10:191. doi: 10.3389/fpls.2019.00191 In the northern hemisphere, boreal forests are a major source of biogenic volatile organic compounds (BVOCs), which drive atmospheric processes and lead to cloud formation and changes in the Earth's radiation budget. Although forest vegetation is known to be a significant source of BVOCs, the role of soil and the forest floor, and especially interannual variations in fluxes, remains largely unknown due to a lack of long-term measurements. Our aim was to determine the interannual, seasonal and diurnal dynamics of boreal forest floor volatile organic compound (VOC) fluxes and to estimate how much they contribute to ecosystem VOC fluxes. We present here an 8 year data set of forest floor VOC fluxes, measured with three automated chambers connected to the quadrupole proton transfer reaction mass spectrometer (quadrupole PTR-MS). The exceptionally long data set shows that forest floor fluxes were dominated by monoterpenes and methanol, with relatively comparable emission rates between the years. Weekly mean monoterpene fluxes from the forest floor were highest in spring and in autumn (maximum 59 and 86 µg m−<sup>2</sup> h −1 , respectively), whereas the oxygenated VOC fluxes such as methanol had highest weekly mean fluxes in spring and summer (maximum 24 and 79 µg m−<sup>2</sup> h −1 , respectively). Although the chamber locations differed from each other in emission rates, the inter-annual dynamics were very similar and systematic. Accounting for this chamber location dependent variability, temperature and relative humidity, a mixed effects linear model was able to explain 79–88% of monoterpene, methanol, acetone, and acetaldehyde fluxes from the boreal forest floor. The boreal forest floor was a significant contributor in the forest stand fluxes, but its importance varies between seasons, being most important in autumn. The forest floor emitted 2–93% of monoterpene fluxes in spring and autumn and 1–72% of methanol fluxes in spring and early summer. The forest floor covered only a few percent of the forest stand fluxes in summer.

Keywords: biogenic volatile organic compound, flux, forest floor, temperature, seasonality, vegetation, decomposition

# INTRODUCTION

fpls-10-00191 February 20, 2019 Time: 18:19 # 2

Global forest ecosystems are the largest existing source of biogenic volatile organic compounds (BVOCs) (Guenther et al., 1995). BVOC emissions are dominated by isoprene and monoterpenes (Lathière et al., 2005; Guenther et al., 2012), and they have a crucial role in the atmosphere as their oxidation products drive secondary organic aerosol (SOA) formation (Kulmala et al., 1998, 2013). SOA formation from volatile organic compounds (VOCs) through oxidation processes and SOA particles can change the Earth's radiation budget as a cooling feedback mechanism (Virtanen et al., 2010; Kourtchev et al., 2016). The cooling feedback mechanism could become stronger if the temperature were to increase in high latitudes by an estimated 6–8◦C by 2100 (IPCC Fifth Assessment Report (AR5), 2014), boosting SOA load and cloud formation by increasing VOC emissions from the boreal biosphere (Kulmala et al., 2014), which is a well-documented VOC source (Rinne et al., 2007). VOCs affect hydroxyl (OH) radical and ozone formation and oxidation processes (Jacob et al., 2005; Hüve et al., 2007). Ozone formation in photochemical reactions requires NOx and reactive VOCs (Crutzen, 1979; Logan, 1985), and for this reason, it is important to quantify seasonal VOC emission fluxes from different sources. These sources should be quantified more accurately, because several studies have shown that measured and modeled ozone deposition fluxes and OH radical reactivities include significant differences (Mogensen et al., 2011; Wolfe et al., 2011; Rannik et al., 2012; Zhou et al., 2017). A gap in atmospheric oxidant sinks between measurements and models can be decreased by including forest floor emissions in these models.

The total VOC load in the atmosphere remains largely unknown because VOCs are difficult to measure and the boreal biosphere emits at least 25 different VOCs (Schallhart et al., 2018). VOCs are released by tree shoots (Aalto et al., 2014), stems (Vanhatalo et al., 2015) and forest floor (Aaltonen et al., 2013), and more specifically by decomposition processes, vegetative litter, root stores and plant VOC synthesis (Hayward et al., 2001; Leff and Fierer, 2008; Faubert et al., 2010). VOC synthesis of plants and microbes together with VOC volatilization are strongly affected by temperature (Guenther et al., 1993; Kesselmeier and Staudt, 1999). Temperature dependent VOC fluxes from the boreal forest floor (Aaltonen et al., 2013; Mäki et al., 2017) can increase when soil temperatures increase in a warming climate (Fan et al., 2014). Warming can also change vegetation cover and directly influence the abundance of the different VOC sources hence affect the amount and blend of emitted VOCs. In this study, we tested whether we could find evidence of such changes in an 8-year continuous data set of forest floor VOC fluxes. Previously reported measurements have been rather shortterm, and only few cover more than one growing season (Hellén et al., 2006; Asensio et al., 2007, 2008; Faubert et al., 2010).

Volatile organic compound fluxes between the forest floor and the atmosphere are dominated by isoprenoids (e.g., isoprene, monoterpenes, sesquiterpenes) and oxygenated VOCs (alcohols, aldehydes, and ketones) (Hellén et al., 2006; Aaltonen et al., 2011, 2013; Mäki et al., 2017). VOC fluxes from the boreal forest floor are highest during spring and autumn (Hellén et al., 2006; Aaltonen et al., 2011, 2013; Mäki et al., 2017). High spring fluxes are likely connected to snowmelt, which increases soil moisture and accelerates microbial decomposition and leaf growth by exposing soil and leaves to radiation. In the northern ecosystems, VOCs are produced by evergreen and deciduous dwarf shrubs (Rinnan et al., 2013), mosses (Hanson et al., 1999), roots (Hayward et al., 2001), and decomposing litter (Isidorov and Jdanova, 2002; Isidorov et al., 2016). A warmer climate will lead to an earlier snowmelt and growing season start, which can cause increased VOC flux in spring and possibly affect the emission blend. A warmer climate can also change precipitation patterns and affect VOC emissions, because soil water content impacts plant metabolism and microbial decomposition (Staudt et al., 2002; Davidson and Janssens, 2006). High monoterpene load from the soil surface into the atmosphere in autumn is likely released by decomposing and decaying pine litter (Mäki et al., 2017), which contains monoterpene storages (Kainulainen and Holopainen, 2002). VOCs can also be released by freeze– thaw and dry–wet cycles (Asensio et al., 2007, 2008; Insam and Seewald, 2010) in spring and autumn.

In this study, our aim was to use the long-term field measurements with automated chambers to determine the dynamics of forest floor VOC fluxes and to estimate the contribution of forest floor fluxes to the whole forest ecosystem fluxes. We analyzed the effects of environmental and sitespecific variation (i.e., air and soil temperature, growing season length, soil moisture, and vegetation cover) on forest floor VOC exchange in order to develop a statistical model that could explain the interannual, seasonal, and daily patterns in emission rates.

# MATERIALS AND METHODS

#### Continuous Chamber Flux Measurements and Supporting Data

We measured the forest floor VOC fluxes in the southern boreal forest (established in 1962) at the SMEAR II (Station for Measuring Ecosystem-Atmosphere Relations) station (61◦ 51<sup>0</sup> N, 24◦ 17<sup>0</sup> E, 180 m above sea level). Canopy basal area is covered by Pinus sylvestris (75%), Picea abies (15%), and broadleaf species (10%) such as Betula pendula and Sorbus aucuparia. The forest floor is covered by ericoid shrubs (35.4%) such as Vaccinium vitisidaea, Vaccinium myrtillus, and Calluna vulgaris, mosses (67.8%) such as Pleurozium schreberi, Dicranum polysetum, Dicranum scoparium, and Hylocomium splendens, tree seedlings (0.2%) such as Sorbus aucuparia, Pinus sylvestris, Betula pendula, and Picea abies, and grasses (8.4%) such as Deschampsia flexuosa and Melampyrum sylvaticum (Mäki et al., 2017). The soil type is Haplic podzol, and organic soil contains the largest proportion of carbon (356 mg g−<sup>1</sup> ) and nitrogen (13 mg g−<sup>1</sup> ) compared with mineral soil (5–32 and ∼1 mg g−<sup>1</sup> , respectively).

We measured the VOC fluxes using three automated dynamic (flow-through) chambers (80 cm × 40 cm × 25 cm) installed on permanent stainless-steel soil collars (80 cm × 40 cm × 10 cm) every spring (Aaltonen et al., 2013). The shared volume of chamber and soil collar was 112 L. The automated chambers

were connected to a quadrupole proton transfer reaction mass spectrometer (quadrupole PTR-MS, Ionicon Analytik, Innsbruck, Austria, de Gouw and Warneke, 2007) and VOC fluxes were measured from April 20, 2010 to November 7, 2017. Chamber frames were made of aluminum and the sides and top of the chamber were covered by a transparent fluorinated ethylene propylene (FEP) film (0.05 mm) (Aaltonen et al., 2013). Mixing the air with the chamber closed was done with two fans (wind speed 2.2 m s−<sup>1</sup> and air volume of flowthrough 167 L min−<sup>1</sup> ) and the chamber opening direction was toward the south. We drew air samples from the chambers using a flow of 1.1–4.0 L min−<sup>1</sup> through a 64 m heated FEP tube (inner diameter 4 mm) and determined the VOC concentrations from a side flow (0.1 L min−<sup>1</sup> ) through a polytetrafluoroethylene (PTFE) tube (Aaltonen et al., 2013) (**Figure 1**). Chamber closure time was 15 min. To avoid underpressure inside the chambers, we substituted the sample flow with well-mixed ambient air using a 30–100% higher flow rate compared with the sampling flow. The VOC concentrations of the substitute air were measured immediately before every chamber enclosure. We measured the masses (amu) of methanol (33), acetaldehyde (45), acetone (59), isoprene (69), benzene (79), monoterpene fragment (81), methyl butenol (87), toluene (93), hexenal (99), hexanal (101), monoterpenes (137), and methyl salicylate (153). We measured all the masses every year, except benzene, methyl salicylate, hexenal, and hexanal between 2010 and 2012 and methyl butenol and toluene between 2013 and 2017. The calibration gas (Apel-Riemer Environmental, Inc., Denver, CO, United States) included methanol, acetaldehyde, acetone, isoprene, benzene, toluene, hexanal (except in summer 2010), and α-pinene. The chamber system was automated and every chamber was measured eight times a day (Aaltonen et al., 2013). The instrument was calibrated 1–4 times per month based on the scheme presented by Taipale et al. (2008). The quadrupole PTR-MS has a high sensitivity and a

short response time (Taipale et al., 2008), suitable to quantify oxidized VOCs (Aaltonen et al., 2013), but is unable to separate compounds with the same molecular mass such as different monoterpenes.

We measured the forest floor carbon dioxide (CO2) and H2O fluxes from the same chambers using infrared light absorption analyzers (URAS4, Hartman and Braun, Frankfurt, Germany), until spring 2013, when the instrument was replaced with a LICOR LI-840A (LI-COR, Lincoln, NE, United States). We estimated the H2O and CO<sup>2</sup> fluxes using the same mass balance equation (Kolari et al., 2012) that was used for the VOC flux calculations. Precipitation was monitored with an FD12P weather sensor (Vaisala Oyj, Helsinki, Finland) at a height of 18.0 m. Soil volumetric water content and soil temperature were measured from four to five pits and at each pit on every soil horizon (H horizon, eluvial E horizon, illuvial B horizon, and parent material C horizon). The mean thickness of the organic soil is 6.0 cm, the eluvial E horizon 2.0 cm and the illuvial B horizon 16 cm at the SMEAR II stand (Mäki et al., 2017). The soil moisture and soil temperature sensors were placed in the middle of each soil horizon. Soil temperatures were measured at 15-min intervals using silicon temperature sensors (Philips KTY81–110, Philips Semiconductors, Eindhoven, Netherlands). Volumetric water contents were recorded at 60-min intervals using the timedomain reflectometry method (TDR-100, Campbell Scientific, Ltd., United Kingdom). Photosynthetically active radiation (PAR) was measured using an LI-190 quantum sensor (Li-Cor Biosciences, Lincoln, NE, United States) above the canopy. The monthly total litterfall and the fraction of needles was determined with 21 L collectors (Mäki et al., 2017). We determined plant coverage (%) of ericoid shrubs such as Vaccinium vitis-idaea and Vaccinium myrtillus, mosses such as Pleurozium schreberi, Dicranum polysetum, Dicranum scoparium, and Hylocomium splendens, and other plant species such as Linnaea borealis and Deschampsia flexuosa for each measurement chamber based

on visual assessment in 2010 and 2017 (**Table 1**). Total plant coverage of soil collars was 55–94% in 2017.

#### Continuous Ecosystem Flux Measurements

fpls-10-00191 February 20, 2019 Time: 18:19 # 4

In order to assess the importance of forest floor fluxes to the total VOC budget of the stand, we compared the measured VOC fluxes from the forest floor with the total ecosystem fluxes. The data set of ecosystem fluxes was published in Rantala et al. (2014, 2015). Ecosystem fluxes of 27 different VOCs were measured with the profile method using the quadrupole PTR-MS (Ionicon Analytik GmbH, Innsbruck, Austria) with a 2-s sampling time from two heights in the canopy (4.2 and 8.4 m) and four heights above the canopy (16.8, 33.6, 50.4, and 67.2 m) (Rantala et al., 2014, 2015). Measurements were performed eight times a day using continuous flow (33 L min−<sup>1</sup> ) through a 100 m long sample tube (PTFE, inner diameter 8 mm). The air sample entered the PTR-MS as a side flow through the 4 m long sampling tube (PTFE, inner diameter 1.6 mm) (Rantala et al., 2014). The quadrupole PTR-MS completed a whole measurement cycle within 1 h, in which the quadrupole PTR-MS ran through each measurement height nine times, with one individual repetition at each height taking 1 min (Rantala et al., 2014). Losses of monoterpenes and oxygenated VOCs in the tube walls are typically small (Kolari et al., 2012), because flow rate is relatively high and the ratio of the tube wall area and tube diameter is relatively small. We measured VOC fluxes for the whole ecosystem from May 28, 2010 to September 8, 2014. The calculations of 30-min average volume mixing ratios were described by Taipale et al. (2008). The VOC fluxes were estimated from the profile measurements using the surface layer profile method based on the Monin-Obukhov theory, and calculations were described in detail by Rantala et al. (2014). For the flux calculations, the ambient temperature was measured from each height using Pt-100 sensors, and turbulence parameters were determined using the three-dimensional acoustic anemometer (Solent 1012R2, Gill Instruments, Ltd., United Kingdom) at a height of 23 m.

#### Chamber Flux Calculations

Flux calculation was based on the mass balance equation (Kolari et al., 2012). The rate of change of concentrations during the first 400 s as well as the known concentrations of substitute air were taken into account. We quantified the VOC fluxes according to the VOC concentration change (C) during the chamber closure, which is derived from the mass balance equation (Hari et al., 1999), Eq. 1:

$$V\frac{dc}{dt} = E + F(C\_i - C) \tag{1}$$

where V is the chamber volume, E is the VOC emission rate (positive flux) or the VOC uptake rate (negative flux), F is the flow rate of ingoing air, and C<sup>i</sup> is the VOC concentration of the air used to replace the sampled air volume. Calculating Eq. 1 for VOC concentration C as a function of time t after chamber closure leads to solution (Kolari et al., 2012):

$$C(t) = C\_0 + (\frac{C\_i - C\_0}{V} + \frac{E}{F})(1 - e^{-\frac{E}{V}}) \tag{2}$$

where C<sup>0</sup> is the VOC concentration (µg m−<sup>3</sup> ) at the time when the chamber was closed. The emission rate E (µg m−<sup>2</sup> h −1 ) corresponds to the rate of change in concentration during the first 400 s after the chamber was closed (Eq. 2).

There was a negative exponential correlation between relative humidity (%) in the chamber headspace and water soluble VOC fluxes, and for this reason we only used methanol, acetone, and acetaldehyde fluxes measured under 75% relative humidity. This threshold led to removal of 74– 81% of the methanol, acetone, and acetaldehyde fluxes for the chambers. All the data analyses for methanol, acetone, and acetaldehyde were performed using filtered data. Relative humidity inside the chamber headspace was measured right before the closure.

We tested the normality of different VOC fluxes measured from the automated chambers (n = 3) with the Kolmogorov–Smirnov test (df = 2). As the data were nonnormally distributed, the statistical differences between years (**Table 2**) and between chambers (**Supplementary Table S1**) were determined using the non-parametric Kruskal–Wallis test (df = 2) at a significance level of < 0.05 (**Supplementary Table S1**).

Weekly mean monoterpene fluxes were estimated using the whole data set (**Figure 2**) and weekly mean methanol fluxes using daytime measurements from 9 a.m. to 8 p.m. (n = 3) (**Figure 3**), because nighttime measurements were mainly filtered away with a 75% relative humidity threshold. We determined

TABLE 1 | Total plant coverage (%) and coverage (%) of Vaccinium vitis-idaea, Vaccinium myrtillus, Linnaea borealis, and mosses of three measurement chambers in 2010 and 2017.


Plant coverage was determined visually from photographs. <sup>∗</sup>Aaltonen et al. (2013).

TABLE 2 | Start and end dates of volatile organic compound (VOC) flux measurements and total annual litterfall (gDW m−<sup>2</sup> ); environmental conditions calculated from March to October: annual mean temperature (◦C), temperature sum (sum of daily mean temperatures > 5 ◦C), mean soil moisture in the O horizon (m<sup>3</sup> m−<sup>3</sup> ), and precipitation sum (mm); timing of snowmelt and first snow in each year; annual mean monoterpene, methanol, acetone, and acetaldehyde fluxes (µg m−<sup>2</sup> h −1 ) with standard deviation of three measurement chambers between 2010 and 2017 at the SMEAR II station.


Snowmelt date is when snow cover is mainly melt.

weekly mean methanol fluxes using nighttime measurements from 9 p.m. to 8 a.m. (**Supplementary Figure S1**). The mean monoterpene, methanol, and acetone fluxes were estimated at different times of day from the boreal forest floor in spring, summer, and autumn (**Figure 4**) based on the assumption that summer starts when daily mean temperature is over 10◦C, and autumn begins when daily mean temperature is below 10◦C. We determined the effect of chamber temperature and relative

FIGURE 4 | Mean (A) monoterpene, (B) methanol, and (C) acetone fluxes (n = 3 chambers), and standard deviation (µg m−<sup>2</sup> h −1 ) during different times of the day from the boreal forest floor in spring, summer, and autumn, from 2010 to 2017. Summer starts when daily mean temperature is over 10◦C, and autumn begins when daily mean temperature is below 10◦C.

humidity on VOC fluxes from the forest floor (**Figure 5** and **Supplementary Figure S2**).

The cumulative sum (µg m−<sup>2</sup> h −1 ) of methanol (33 amu), acetaldehyde (45 amu), acetone (59 amu), and monoterpenes (137 amu) for each 2-week period was calculated for each chamber and for the whole ecosystem using an unfiltered data set. The mean of the cumulative sums of three chambers was compared with the total ecosystem fluxes by calculating the proportion of soil fluxes relative to the whole ecosystem fluxes for each 2-week period (**Figure 6**).

#### The Mixed Effects Linear Model

We developed a mixed effects linear model to determine which parameters can be used to explain the fluxes of the different VOCs. Model parameters are presented in **Supplementary Table S2**. We used the effect of chamber temperature, relative humidity, and above-canopy PAR to model monoterpene, methanol, acetone, and acetaldehyde fluxes. VOC fluxes (V) were modeled by the mixed effects linear model (3):

$$V = B\_0 + B\_t + B\_R + B\_{tR} + \in,\tag{3}$$

where B<sup>0</sup> represents a fixed intercept parameter, B<sup>t</sup> represents fixed unknown parameters associated with chamber temperature, B<sup>R</sup> represents fixed unknown parameters associated with relative humidity of the chamber, and BtR represents fixed parameters for interaction of the measurement chamber with relative humidity and temperature. In the model (3), the error term ∈ is presumed to have the form:

$$
\epsilon \in \infty\_{CR} + \infty\_{CT} + \mu,\tag{4}
$$

where ∝ CR represents random parameters related to interaction of the measurement chamber (1, 2, and 3) and relative humidity of the chamber, ∝ CT represents random parameters related to interaction of the measurement chamber and chamber temperature, and u is an unobservable random error term.

#### Accuracy of the Chamber Measurement Method

We tested the chamber wall effects in laboratory conditions by flushing gas containing the study compounds into the chamber headspace, and the concentrations of ingoing and outgoing air were determined by sampling air into the Tenax TA-Carboback-B adsorbent tubes. Concentration of the calibration gas in the ingoing air was from 1.7 to 3.1 µg m−<sup>3</sup> , varying with the compounds. The chamber was placed on a flat surface covered with a transparent FEP film (0.05 mm); air (flow rate 3 dm<sup>3</sup> min−<sup>1</sup> ) was moving into and out of the chamber through FEP tubes (inner diameter 6 mm), and we drew air samples from a side flow using a flow of 0.1–0.15 dm<sup>3</sup> min−<sup>1</sup> . We determined the concentrations of individual monoterpenes (α-pinene, camphene, β-pinene, 13-carene, p-cymene, 1,8-cineol, limonene, terpinolene, linalool,

and myrcene) and sesquiterpenes (longicyclene, isolongifolene, β-caryophyllene, aromadendrene, and α-humulene) from the adsorbent tubes by using a thermodesorption instrument (PerkinElmer TurboMatrix 650, Waltham, MA, United States) coupled with a gas chromatograph (PerkinElmer Clarus 600, Waltham, MA, United States) with a mass selective detector (PerkinElmer Clarus 600T, Waltham, MA, United States) (Mäki et al., 2017). Relative standard deviation (RSD, %) was calculated as a standard deviation between the four parallel samples taken during the chamber enclosure in the laboratory using constant concentrations of VOCs of ingoing air. RSD shows that the error of sampling and analytical method (thermal desorption– gas chromatography–mass spectrometry) was relatively low for both monoterpenes (5–14%) and sesquiterpenes (6–9%) (**Supplementary Table S3**).

#### RESULTS

## Interannual, Seasonal, and Diurnal Dynamics in Soil VOC Exchange

#### Interannual Dynamics

Forest floor VOC exchange was dominated by monoterpenes and oxygenated VOCs such as methanol, acetone, and acetaldehyde, while isoprene fluxes were mainly close to zero. Annual mean monoterpene fluxes ranged from 7 to 19 µg m−<sup>2</sup> h −1 , methanol fluxes from −1 to 3 µg m−<sup>2</sup> h −1 , acetone fluxes from 0 to 2 µg m−<sup>2</sup> h −1 , and acetaldehyde fluxes from 1 to 2 µg m−<sup>2</sup> h −1 (**Table 2**). The VOC fluxes showed comparable emission rates between years.

#### Seasonal Dynamics

Interestingly, the seasonal dynamics of monoterpene fluxes from the forest floor differed from those of oxygenated VOCs. The highest weekly mean monoterpene fluxes were observed from the forest floor in spring (May–June) and in autumn (September– October) (maximum 59 and 86 µg m−<sup>2</sup> h −1 , respectively) (**Figure 2**). Weekly mean methanol fluxes were highest in spring and summer (May–August) (**Figure 3**: maximum 24 and 79 µg m−<sup>2</sup> h −1 , respectively) based on daytime measurements from 9 a.m. to 8 p.m. We found that acetone and acetaldehyde fluxes followed a similar seasonal pattern to that of methanol. Acetone and acetaldehyde fluxes were clearly lower than monoterpene fluxes but were at the same level as methanol fluxes (**Table 2**).

#### Diurnal Dynamics

We determined the diurnal dynamics of the forest floor VOC fluxes in spring, summer, and autumn. The diurnal maximum fluxes for monoterpenes, methanol, and acetone were observed in spring (19, 4, and 3 µg m−<sup>2</sup> h −1 ), summer (33, 27, and 9 µg m−<sup>2</sup> h −1 ), and autumn (28, 10, and 1 µg m−<sup>2</sup> h −1 ), respectively (**Figure 4**). The monoterpene and methanol fluxes were strongest during the daytime in all seasons, while acetone and acetaldehyde fluxes showed clear diurnal dynamics only in spring and summer. Relatively high monoterpene, methanol, and acetone fluxes were also observed from 6 p.m. to 8 p.m. in summer.

#### Temperature and Relative Humidity Effect on Forest Floor VOC Exchange

Our results showed that monoterpene, methanol, and acetone fluxes correlate with chamber temperature (**Supplementary Figure S2**). Chamber temperature explained 14–61% of methanol and 25–57% of acetone fluxes in spring and in summer (**Supplementary Figure S2**), but not in autumn, except in chamber 1 (55 and 51%). Monoterpene emissions showed correlation (7–50%) with chamber temperature from spring to summer and weak correlation in autumn (2–10%) (**Supplementary Figure S2**). Relative humidity had a significant effect on monoterpene and methanol fluxes (**Figure 5**).

#### Mixed Effects Linear VOC Emission Model

The model analysis showed that the VOC flux rates from the forest floor are indeed strongly affected by temperature and relative humidity. The model explained 79–88% of monoterpene, methanol, acetone, and acetaldehyde fluxes from the boreal forest floor (**Figures 7**, **8** and **Supplementary Figure S3**). This fairly simple statistical modeling approach is able to capture the individual behavior of each measurement chamber.

#### Importance of Forest Floor to Ecosystem VOC Exchange

We compared forest floor VOC emissions to the whole ecosystem emissions measured by the flux gradient method between 2010 and 2014 (Rantala et al., 2014, 2015). The forest floor accounted for 2–93% of monoterpene fluxes relative to forest stand fluxes in spring and autumn and 1–72% of methanol fluxes in spring and early summer (**Figure 6**). The role of forest floor in forest stand fluxes was only a few percent in summer. Fluxes were compared using the whole data set without filtering with 75% relative humidity. VOC deposition dominated in the forest floor or forest stand during some 2-week periods (negative values in **Figure 6**).

#### Uncertainties in VOC Exchange Measurements From the Forest Floor

The chambers differed from each other in emission rates but showed very similar and systematic temporal dynamics. Monoterpene and methanol fluxes were highest from chamber 1, where temperature was the highest (**Supplementary Table S1**). The acetaldehyde and acetone fluxes were mainly highest from chamber 1, but the differences between chambers were small. Chamber temperature was strongly correlated with ambient air.

# DISCUSSION

#### Interannual Dynamics

Our 8-year long data series of forest floor VOC exchange showed relatively small variation between the years. Annual mean monoterpene fluxes were very similar to earlier studies performed on the boreal forest floor (Aaltonen et al., 2011: α-pinene 0–14 µg m−<sup>2</sup> h −1 ; Mäki et al., 2017: total monoterpenes 23 µg m−<sup>2</sup> h −1 ; Wang et al., 2018: total monoterpenes 3–10 µg m−<sup>2</sup> h −1 ). Also, mean fluxes of oxygenated VOCs, methanol, acetone, and acetaldehyde, were similar to earlier measurements at the same site in 2010 (Aaltonen et al., 2013: methanol −0.6 to 7.2 µg m−<sup>2</sup> h −1 ; acetone −0.8 to 2.2 µg m−<sup>2</sup> h −1 ; acetaldehyde 0.8–2.2 µg m−<sup>2</sup> h −1 ).

Most of the interannual variability in monoterpene exchange resulted from variations in the temperature sum and precipitation between the years, while a similar trend was not observed for methanol, acetone, and acetaldehyde. Periods of high precipitation decreased monoterpene fluxes, which were seen also in low annual mean fluxes in 2011, 2012, and 2017. The effect of precipitation was assumed to be due to increased relative humidity and formation of water films on soil and leaf surfaces, which decreases VOC evaporation from surfaces. Annual mean monoterpene flux was also high in 2016 due to high monoterpene fluxes in October, which were likely released by fresh decomposing litter (Mäki et al., 2017; Wang et al., 2018).

# Seasonal Dynamics

The strong seasonal variation on monoterpene and VOC fluxes most probably results from seasonality in the production of VOCs via plant VOC synthesis (Aaltonen et al., 2011; Faubert et al., 2012) and litter decomposition (Hayward et al., 2001; Greenberg et al., 2012; Mäki et al., 2017). Temperature regulates plant VOC synthesis (Guenther et al., 1993), and hence this could explain the higher annual mean monoterpene fluxes during the warm summers (2010 and 2013) compared to the cold summers (2012 and 2017) (**Table 2**).

At our site, Vanhatalo et al. (2015) observed bursts of monoterpene from Scots pine stems during spring recovery in April. Similarly, a release of monoterpenes from belowground storages such as roots (Hayward et al., 2001) could explain the high forest floor monoterpene emissions during spring. Based on our study, it remains unclear whether monoterpenes were mostly emitted by litter or by roots, while the production in plants was likely small due to the fact the monoterpene emissions from all three chambers were similar in magnitude despite the differences in plant species coverage. Vegetation in chamber 1 was dominated by mosses that emit mostly isoprene and oxygenated VOCs and only minimal amounts of monoterpenes (Hanson et al., 1999; Hellén et al., 2006), while vegetation in chambers 2 and 3 were mostly dominated by Vaccinium spp. that emits monoterpenes (Aaltonen et al., 2011; Faubert et al., 2012). There were hardly any VOC flux measurements directly after snowmelt, while Hellén et al. (2006) found that seasonally highest monoterpene fluxes from the boreal forest floor were measured right after snowmelt at our measurement

site. High monoterpene emissions in autumn are linked to decomposing litter (Greenberg et al., 2012; Mäki et al., 2017) and the organic soil, which contains easily available carbon for microbial decomposition.

Seasonality of oVOC fluxes in the forest floor was comparable to those of the forest stand fluxes. Gross primary production of ground vegetation is generally highest from mid-June to mid-August at our site (Kolari et al., 2006). The maximum forest floor methanol fluxes coincide with the maximum biomass production, which suggests that new surface vegetation growth is a major methanol source at our measurement site. The seasonality of forest floor methanol emissions follows the seasonality of canopy methanol exchange (Aalto et al., 2014) indicating that the processes and sources of methanol are the same both in the canopy and in the forest floor. Other evidence is that oxygenated VOC fluxes were lowest from chamber 3 with the lowest vegetation cover from June to August (**Supplementary Table S1**). Methanol can also be released by microbes that synthesize VOCs (Bäck et al., 2010; Mancuso et al., 2015)especially in summer, when decomposer activity is high due to temperature driven enzyme activity (Davidson and Janssens, 2006). Methanol fluxes in early autumn are probably related to plant senescence – leaf litter is a strong methanol source (Warneke et al., 1999) and the senescence of Vaccinium myrtillus leaves is in early autumn (Kulmala et al., 2008).

#### Diurnal Dynamics

The diurnal dynamics in forest floor VOC exchange were clearly affected by temperature and radiation changes between the daytime and nighttime, which corresponds with the earlier observations at the site (Aaltonen et al., 2013). Monoterpene synthesis in plants is typically light and temperature dependent (Kesselmeier and Staudt, 1999), and similar light and temperature dependent VOC production has been demonstrated with ericoid shrubs, grasses, and mosses (Hellén et al., 2006; Faubert et al., 2010; Aaltonen et al., 2013; Lindwall et al., 2015; Kramshøj et al., 2016). Monoterpene and methanol fluxes are affected by temperature due to the volatility of these compounds (Guenther et al., 1993) and by radiation that heats leaf surfaces. Methanol fluxes are also strongly dependent on stomatal conductance (Nemecek-Marshall et al., 1995), which is high during daytime when plants maintain their photosynthesis by taking up CO<sup>2</sup> through stomata and maintain water transport through transpiration.

Nighttime monoterpene emissions observed in spring and autumn (**Figure 4**) were likely affected by soil processes independent from light availability. Oxygenated VOC fluxes were close to zero in the nighttime, but we also observed deposition on leaf surfaces and chamber walls. Night-time deposition and daytime emissions of methanol were also observed from the bare temperate cropland in spring (Bachya et al., 2018). In our study, the nighttime deposition of methanol, acetone, and acetaldehyde is linked to lower temperatures and higher humidity in the chamber headspace. VOCs were likely released from the aqueous layer on plant and soil surfaces in the morning, when air humidity drops.

# Modeling Soil VOC Exchange

Chamber temperature explained monoterpene, methanol, and acetone fluxes from the forest floor (**Supplementary Figure S2**), similar to methanol and acetone fluxes from Pinus sylvestris shoots, which correlated with temperature (Aalto et al., 2014). Temperature dependence of monoterpenes is not as clear, because monoterpenes are emitted immediately from synthesis or with a delay from storage structures in resin ducts, glandular trichomes, or other storage structures (Laothawornkitkul et al., 2009). Litter also contains monoterpene storage structures and accelerates microbial decomposition (Hayward et al., 2001; Isidorov and Jdanova, 2002; Leff and Fierer, 2008; Gray et al., 2010; Isidorov et al., 2016). Monoterpene release from litter is likely directly linked to temperature due to missing monoterpene synthesis. However, low correlation between chamber temperature and monoterpene fluxes gives an indication that monoterpenes are released from and consumed simultaneously in the forest floor. Ramirez et al. (2010) showed that soil can be a sink for VOCs released by litter, like methanol that is used as a carbon source by bacteria via the ribulose monophosphate cycle and the serine cycle (Quayle and Ferenci, 1978). Different abiotic and biotic stresses also affect plant VOC emissions (Baldwin et al., 2006; Loreto and Schnitzler, 2010; Niinemets et al., 2013), which can make temperature responses of forest floor fluxes more complex.

The mixed effects linear model accounting for location dependent variability, temperature, and relative humidity was able to explain 79–88% of monoterpene, methanol, acetone, and acetaldehyde fluxes from the boreal forest floor (**Figures 7, 8** and **Supplementary Figure S3**). It seems that with such a multiannual data set, a fairly simple statistical modeling approach is able to capture individual behavior of each measurement chamber, while the challenge remains how to capture the spatial variation of soil VOC exchange. This model is not directly suitable to estimate forest floor VOC exchange in other ecosystems based on ambient temperature and relative humidity. The model clearly underestimates high methanol, acetone, and acetaldehyde fluxes, indicating that environmental factors other than temperature and relative humidity are the main drivers of these fluxes in the chambers. One likely reason is deposition of these water soluble molecules on moist leaf and chamber surfaces. Deposition should be modeled with a different mixed effects linear model, because the processes behind VOC deposition are different than the biological and physicochemical processes, which regulate production and evaporation of VOCs.

Volatile organic compound fluxes seemed to be strongly stimulated by temperature when the relative humidity effect is included (**Figures 7**, **8**). VOC fluxes were observed to decrease with increasing relative humidity (**Figure 5**). Therefore our results give an indication that VOC fluxes from the boreal forest floor will likely increase in a warming climate. The global mean surface temperature is expected to increase by a minimum of 0.3–1.7◦C under Representative Concentration Pathway (RCP) 2.6 and a maximum of 2.6–4.8◦C under RCP8.59 by 2100 (IPCC Fifth Assessment Report (AR5), 2014).

A temperature increase of 2◦C can lead to a twofold increase in monoterpene fluxes and a five-fold increase in sesquiterpene fluxes in subarctic ecosystems (Valolahti et al., 2015). In our study, monoterpene fluxes correlate with temperature in spring and summer, but the flux rate rise was smaller with increasing temperature, possibly because vegetation is less sensitive to temperature fluctuations in boreal ecosystems compared with arctic ecosystems.

In boreal ecosystems, a warming climate is expected to increase ericoid shrub cover and decrease moss cover, which can affect leaf litter quality and decomposition rates (De Long et al., 2016). Decomposition activity is expected to increase due to temperature dependent enzyme activity of microbes, and this can affect belowground VOC production. Warming can also affect plant VOC synthesis, release of evaporated VOCs through stomata, and the abundance of different plant species. The boreal forest floor is prevailingly covered by monoterpeneemitting Vaccinium spp. (Faubert et al., 2012; Aaltonen et al., 2013) at our site. Vaccinium myrtillus has shown a 36% increase in aboveground biomass by warming soil (Anadon-Rosell et al., 2014), which could make understory vegetation less diverse (Dawes et al., 2011) and impact monoterpene production in a warming climate.

#### Forest Floor Affects Ecosystem Fluxes

The comparison indicated that the forest floor is a significant contributor to the whole ecosystem fluxes, but its importance varies between seasons, being most important in autumn. The forest floor covered only a few percent of the forest stand fluxes in summer. Our results are in line with those of Aaltonen et al. (2013), who showed that the forest floor covered from several percent to 10s of percent of the total ecosystem fluxes depending on the compound and the season. Our method of comparing the observed forest floor level VOC emissions with ecosystem level VOC fluxes does not include oxidation of VOCs taking place during the transport from the forest floor to the abovecanopy atmosphere. Therefore our results presented in **Figure 7** should be considered as order of magnitude estimates rather than definitive values determining the proportion of forest floor oriented VOCs.

The seasonal dynamics of oxygenated VOCs (methanol, acetaldehyde, and acetone) were similar between the forest floor and the whole ecosystem (data not shown), indicating that vegetation is likely a significant source of oxygenated VOCs from the forest floor. In trees, methanol production is connected to growth and is produced mainly in leaves and needles in spring, and in stem and root expansion in summer. Soil VOC production can also influence the canopy fluxes, because methanol produced in the stem and roots during their growth, can be transported to the canopy via transpiration stream (Rissanen et al., 2018), and hence be emitted through stomata (Folkers et al., 2008). Peaking acetone and acetaldehyde fluxes above a hardwood forest in autumn were speculated to result from leaf senescence and decaying biomass (Karl et al., 2003). Acetone and acetaldehyde are also produced in the air from the oxidation of other VOCs.

Our results show that monoterpene fluxes from the forest floor contributed the highest load (9–93%) to the whole ecosystem fluxes in autumn due to low canopy fluxes. Monoterpenes are likely released from decomposing litter and plant senescence, because VOC synthesis of ground vegetation likely decreased similarly to shoot VOC emissions. This can provide a significant contribution to the OH sink in boreal forest air. VOCs affect formation and oxidation processes of OH radicals and ozone (Jacob et al., 2005; Hüve et al., 2007). Monoterpenes, isoprene, and other organic compounds were calculated to cover about 24% of the total OH reactivity in August based on a onedimensional vertical chemistry-transport model (Mogensen et al., 2011). Later, it was estimated that monoterpenes cover ∼14% and oxidized VOCs ∼44% of OH radical sinks (Mogensen et al., 2015). Ozone reacts mainly with inorganic compounds and only ∼3% with monoterpenes and 6% with sesquiterpenes and is deposited on the soil and vegetation surfaces (Mogensen et al., 2015; Zhou et al., 2017).

## Uncertainties

Temporal dynamics of forest floor VOC exchange can rather well be captured with three chambers, while several factors may affect and increase the uncertainties in chamber measurements. Our results indicate that vegetation cover affects BVOC emission rates and together with varying temperature creates spatial variation in the fluxes. The contribution of light to soil VOC fluxes is relatively small due to poor light availability under a closed canopy. Monoterpenes can also stick on lipophilic plant leaves (Mäki et al., 2017), which could explain why the highest fluxes were observed from chamber 1 with the lowest Vaccinium spp. coverage.

The quadrupole-PTR-MS is an instrument with high sensitivity and short response time and it has been shown to be able to perform online-measurements of reactive trace gasses from grasslands (Pape et al., 2009). Flux measurements of water soluble VOCs such as methanol, acetone, and acetaldehyde are sensitive to biases, because they tend to be adsorbed on moist surfaces (vegetation, soil, and chamber walls). For this reason, soil chamber data of water soluble VOCs was filtered with 75% relative humidity. A relative humidity of 70% in the chamber has been recommended by Kolari et al. (2012). Relative humidity of the headspace is often over 75% when measurements are performed during the nighttime or rainy days, meaning that water soluble VOCs are dissolved in moist surfaces on chamber walls and leaves. Seasonal fluxes of oxygenated VOCs were slightly overestimated, because the flux measurements performed during the nighttime or rainy days are underrepresented in the data. There were also gaps in the data due to technical analyzer problems, which affected the data coverage between years. Data were missing from mid-June to October in 2012, from September to October in 2014, and from April to mid-August in 2016.

We found hardly any correlation between the VOC fluxes and soil temperature. This may be caused by the spatial variation in soil temperature, and the lack of soil temperature sensors installed right next to the soil chambers, and in the top-most litter layer, where most of the soil-emitted VOCs are released (Hayward et al., 2001; Greenberg et al., 2012). Moreover, the effect of soil temperature on VOC fluxes may be hidden by the complexity of different VOC sources and their different responses to temperature.

## CONCLUSION

fpls-10-00191 February 20, 2019 Time: 18:19 # 12

We used the 8-year data set to assess the interannual, seasonal, and diurnal dynamics of forest floor VOC fluxes and compared them with the simultaneously measured total ecosystem fluxes. The forest floor affects ecosystem VOC exchange, emitting 2– 93% of monoterpene fluxes in spring and autumn and 1–72% of methanol fluxes in spring and early summer. Oxygenated VOC fluxes showed similar seasonal dynamics between the forest stand and the forest floor.

The highest weekly mean monoterpene fluxes were observed from the forest floor in spring and autumn (maximum 59 and 86 µg m−<sup>2</sup> h −1 , respectively) and the highest weekly mean methanol fluxes in spring and summer (maximum 24 and 79 µg m−<sup>2</sup> h −1 , respectively). The seasonal dynamics indicate that litter and ground vegetation were the dominant sources of monoterpenes and oxygenated VOC from the boreal forest floor. Forest floor VOC exchange was dominated by monoterpenes and methanol and the flux dynamics was relatively similar between the 8 years, whereas the emission rates differed between the chamber locations. Accounting for location dependent variability and temperature and relative humidity, a mixed effects linear model was able to explain 79– 88% of monoterpene, methanol, acetone, and acetaldehyde fluxes from the boreal forest floor. Forest floor VOC exchange was measured in the prevailing climate, but based on temperature responses of monoterpenes, and especially oxygenated VOCs, it seems that forest floor VOC fluxes will likely increase in warming climate.

Forest floor VOC exchange should be measured using continuous long-term measurements in the different ecosystems

#### REFERENCES


to define the contribution of soils to ecosystem VOC exchange and hence their effect on atmospheric processes globally.

#### DATA AVAILABILITY

The data of forest floor VOC fluxes for this study can be found in the **Supplementary Material**.

## AUTHOR CONTRIBUTIONS

MM was responsible for preparing the manuscript. All authors participated in planning the experiment, analyzing data, preparing the manuscript, and reviewed the manuscript.

# FUNDING

This research was funded by the Academy of Finland Centre of Excellence Programme (Grant No. 307331), Academy Research Fellow project (Grant No. 2884941), the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 757695), and by the Jenny and Antti Wihuri Foundation.

#### ACKNOWLEDGMENTS

The SMEAR II station staff are acknowledged for their help in installing and maintaining the continuous field measurements.

#### SUPPLEMENTARY MATERIAL

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


spectrometry (PTR-TOF-MS). Appl. Soil Ecol. 86, 182–191. doi: 10.1016/j. apsoil.2014.10.018


above a boreal forest. Atmos. Chem. Phys. 18, 815–832. doi: 10.5194/acp-18- 815-2018


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

Copyright © 2019 Mäki, Aalto, Hellén, Pihlatie and Bäck. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Combinations of Abiotic Factors Differentially Alter Production of Plant Secondary Metabolites in Five Woody Plant Species in the Boreal-Temperate Transition Zone

John L. Berini1,2 \*, Stephen A. Brockman<sup>3</sup> , Adrian D. Hegeman<sup>3</sup> , Peter B. Reich2,4,5 , Ranjan Muthukrishnan<sup>6</sup> , Rebecca A. Montgomery2,4 and James D. Forester2,6

<sup>1</sup> Conservation Biology Graduate Program, University of Minnesota, St. Paul, MN, United States, <sup>2</sup> Institute on the Environment, University of Minnesota, St. Paul, MN, United States, <sup>3</sup> Department of Horticultural Science, The Microbial and Plant Genomics Institute, University of Minnesota, St. Paul, MN, United States, <sup>4</sup> Department of Forest Resources, University of Minnesota, St. Paul, MN, United States, <sup>5</sup> Hawkesbury Institute for the Environment, Western Sydney University, Penrith, NSW, Australia, <sup>6</sup> Department of Fisheries, Wildlife, and Conservation Biology, University of Minnesota, St. Paul, MN, United States

#### Edited by:

Bartosz Adamczyk, University of Helsinki, Finland

#### Reviewed by:

Raquel Esteban, University of the Basque Country (UPV/EHU), Spain Adam Matkowski, Wroclaw Medical University, Poland

> \*Correspondence: John L. Berini beri0015@umn.edu

#### Specialty section:

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

Received: 21 May 2018 Accepted: 09 August 2018 Published: 05 September 2018

#### Citation:

Berini JL, Brockman SA, Hegeman AD, Reich PB, Muthukrishnan R, Montgomery RA and Forester JD (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. doi: 10.3389/fpls.2018.01257 Plant secondary metabolites (PSMs) are a key mechanism by which plants defend themselves against potential threats, and changes in the abiotic environment can alter the diversity and abundance of PSMs. While the number of studies investigating the effects of abiotic factors on PSM production is growing, we currently have a limited understanding of how combinations of factors may influence PSM production. The objective of this study was to determine how warming influences PSM production and how the addition of other factors may modulate this effect. We used untargeted metabolomics to evaluate how PSM production in five different woody plant species in northern Minnesota, United States are influenced by varying combinations of temperature, moisture, and light in both experimental and natural conditions. We also analyzed changes to the abundances of two compounds from two different species – two resin acids in Abies balsamea and catechin and a terpene acid in Betula papyrifera. We used permutational MANOVA to compare PSM profiles and phytochemical turnover across treatments and non-metric multidimensional scaling to visualize treatment-specific changes in PSM profiles. We used linear mixed-effects models to examine changes in phytochemical richness and changes in the abundances of our example compounds. Under closed-canopy, experimental warming led to distinct PSM profiles and induced phytochemical turnover in B. papyrifera. In open-canopy sites, warming had no influence on PSM production. In samples collected across northeastern Minnesota, regional temperature differences had no influence on PSM profiles or phytochemical richness but did induce phytochemical turnover in B. papyrifera and Populus tremuloides. However, warmer temperatures combined with open canopy resulted in distinct PSM profiles for all species and induced phytochemical turnover in all but Corylus cornuta. Although neither example compound in A. balsamea was influenced by any of the abiotic conditions, both compounds in B. papyrifera exhibited

**30**

significant changes in response to warming and canopy. Our results demonstrate that the metabolic response of woody plants to combinations of abiotic factors cannot be extrapolated from that of a single factor and will differ by species. This heterogeneous phytochemical response directly affects interactions between plants and other organisms and may yield unexpected results as plant communities adapt to novel environmental conditions.

Keywords: phytochemical turnover, PSM diversity, untargeted metabolomics, balsam fir, beaked hazel, paper birch, red maple, trembling aspen

## INTRODUCTION

Plant secondary metabolites (PSMs) are one of the primary ways in which plants respond to environmental variability, and regulation of PSM production is strongly influenced by the local environment (Wink, 1988; Bennett and Wallsgrove, 1994; Bray et al., 2000; Hirt and Shinozaki, 2003). Many interactions between plants and other organisms are mediated by PSMs (Farmer, 2001; Karban et al., 2006; Karban, 2008), and thus, the biochemical mechanisms that influence these interactions are modulated, at least in part, by the presence, absence, or magnitude of various environmental factors (DeLucia et al., 2012; Jamieson et al., 2012). For instance, changes in the amount and seasonality of precipitation have been shown to influence concentrations of cyanogenic glycosides (Gleadow and Woodrow, 2002; Vandegeer et al., 2013), and elevated concentrations of atmospheric CO<sup>2</sup> often result in increased concentrations of condensed tannins (Lindroth, 2012). Evidence is mounting that recent warming may also influence the production of PSMs (Kuokkanen et al., 2001).

Studies investigating the influence of warming on PSM production suggest that temperature-induced changes to PSMs may be species, compound, or even context dependent. For example, warming has been shown to have no effect on levels of phenolics in red maple (Acer rubrum, Williams et al., 2003), Norway spruce (Picea abies, Sallas et al., 2003), and Scots pine (Pinus sylvestris, Sallas et al., 2003) but resulted in decreased levels of phenolics in dark-leaved willow (Salix myrsinifolia, Veteli et al., 2006) and silver birch (Betula pendula, Kuokkanen et al., 2001). Additionally, warming has been shown to increase levels of terpene-based compounds in Norway spruce (Sallas et al., 2003), Ponderosa pine (Pinus ponderosa, Constable et al., 1999), and Scots pine (Sallas et al., 2003) but has been shown to both increase (Constable et al., 1999) and decrease (Snow et al., 2003) levels of monoterpenes in Douglas fir (Pinus menziesii). While evidence of warming-induced changes to phytochemistry is important to our understanding of how plants will respond to future climates, in natural settings, elevated temperature often combines with other abiotic conditions to influence PSM production and potentially modulate any changes to phytochemistry that might otherwise be induced by warming alone.

As temperatures continue to rise, global precipitation patterns are expected to shift (Hurrell, 1995; Alexander et al., 2006; IPCC, 2014) and light availability to understory plants will likely be altered due to changes in the frequency and intensity of forest disturbance patterns (Canham et al., 1990; Dale et al., 2001). While variability in each of these environmental factors has been shown to influence production of PSMs on their own (Bryant et al., 1983; Dudt and Shure, 1994; Pavarini et al., 2012), combinations of factors can have a distinct effect (Rizhsky et al., 2002, 2004; Mittler, 2006; Zandalinas et al., 2018). Moreover, plant responses to combinations of abiotic factors can be either synergistic or antagonistic (Bonham-Smith et al., 1987; Mittler, 2006; Zandalinas et al., 2018). For example, drought has been shown to enhance cold tolerance (Cloutier and Andrews, 1984), but also exacerbate a plant's intolerance of high temperatures (Rizhsky et al., 2002). Further, different combinations of salinity and high temperatures have been shown to have both positive and negative influences on the metabolism of reactive oxygen species and stomatal response (Zandalinas et al., 2018). Regardless, the vast majority of current research remains focused on the influences of individual conditions rather than considering potential interactions among them.

Until recently, the majority of studies investigating the potential influence of different abiotic factors largely considered the effects of these factors on individual compounds or small groups of compounds. However, individual metabolites rarely, if ever, function in isolation (Gershenzon et al., 2012). Rather, the influence of any one compound is dependent on conditions within the local environment, as well as the relative abundance of numerous other metabolites within a plant's array of chemical constituents (Dyer et al., 2003; Richards et al., 2010; Gershenzon et al., 2012; Jamieson et al., 2015). Thus, understanding how changes in the abiotic environment will influence a plant's metabolic profile is important for interpreting how these changes will influence the abundance and biological role of individual compounds as well.

Phytochemical diversity influences how effective plants are when defending against a range of threats (Gershenzon et al., 2012; Frye et al., 2013; Richards et al., 2015). Compounds may act synergistically, thereby forming mixtures that can provide enhanced protection against potential hazards (Gershenzon, 1984; Harborne, 1987; Gershenzon et al., 2012). Indeed, recent evidence suggests that the number of individual compounds comprising a plant's phytochemical profile can even influence local biological diversity via the influence of changes in toxicity on rates of herbivory (Richards et al., 2015). Increased diversity of secondary metabolites may also allow for more precise communication between plants, thereby allowing for more robust protection against a range of conditions (Iason et al., 2005; Poelman et al., 2008; Gershenzon et al., 2012; Frye et al., 2013). Two metrics that are useful for assessing changes in phytochemical diversity are "phytochemical richness" (i.e., the

absolute number of compounds produced) and "phytochemical turnover" (i.e., the degree of overlap among the compounds produced), as both measures provide different insights into the metabolic response of plants to a range of abiotic conditions.

Variability in phytochemistry, even within the same species, may influence ecosystem structure and function through an array of chemically driven ecological effects (Bukovinszky et al., 2008; Gillespie et al., 2012; Sedio et al., 2017). The growthdifferentiation balance hypothesis (GDBH) suggests that as the local environment becomes increasingly stressful, growth processes will become limited and the production of PSMs will increase until the point that PSM production also becomes limited by resource acquisition/availability (Lerdau et al., 1994). While phytochemical diversity has not been explicitly tested in light of the GDBH, studies have shown that herbivore-induced secondary chemistry can be completely suppressed in some woody species under a range of abiotic conditions (Lewinsohn et al., 1993), rendering them vulnerable to further invasion by pests and pathogens. While the number of studies investigating the effects of warming and other abiotic conditions on PSM production is rapidly growing, we currently have a limited understanding of how different abiotic factors may interact to influence phytochemical diversity (Bidart-Bouzat and Imeh-Nathaniel, 2008; Jamieson et al., 2012, 2015). The objective of this study was to determine how elevated temperatures may influence the production of PSMs and to evaluate how the addition of other abiotic factors may modulate this effect.

While a targeted approach uses standard model compounds to identify and observe changes in specific compounds selected a priori, an untargeted (i.e., global) approach makes no assumptions regarding specific metabolites, and therefore, allows one to observe global changes across the entire metabolic profile. Thus, the strength of an untargeted approach lies in the potential to discover unanticipated changes in metabolic profiles as a result of environmental perturbations (Crews et al., 2009). Although untargeted metabolomics have been used in medicine for clinical diagnosis of various diseases, including numerous forms of cancer (Sreekumar et al., 2009; Jain et al., 2015), this study is among the first to apply this method to an ecological setting.

We used an untargeted metabolomics approach to evaluate how the phytochemical profiles of five different woody plant species are influenced by temperature, soil moisture, and light. Specifically, we tested the hypothesis that elevated temperatures alter the production of PSMs by leading to phytochemical profiles that are distinct from those found at ambient temperature (H1) and that warming will change phytochemical diversity via reductions in phytochemical richness or a high degree of turnover (H2). We also tested the hypothesis that the addition of other abiotic factors, specifically high light and drought, will either magnify or nullify temperature-induced changes in phytochemical profiles and PSM diversity (H3). Finally, because individual compounds may vary greatly in response to heterogeneity in the abiotic environment, we identified two 'example compounds' from balsam fir (Abies balsamea – two unspecified di-terpene resin acids) and paper birch (Betula papyrifera – catechin and another unspecified diterpene resin acid) and analyzed the effects of different sets of

abiotic factors (high-temperature, light, and drought) on their relative abundance. Specifically, we tested the hypothesis that individual compounds will respond to different conditions and combinations of conditions by either increasing or decreasing in relative abundance, potentially in a non-uniform and unpredictable manner (H4).

#### MATERIALS AND METHODS

#### Experimental Design

The Boreal Forest Warming at an Ecotone in Danger (B4WarmED) project is an ecosystem experiment that simulates both above- and below-ground warming in a boreal forest community. The experiment was conducted at Cloquet Forestry Center (CFC; Cloquet, MN, United States) and was initiated in 2008. The experimental design consists of a 2 (overstory – open and closed) × 3 (warming – ambient, ambient +1.7◦C, and ambient +3.4◦C) × 2 (precipitation – ambient and ambient −40%) factorial design with six replicates (two per block) per treatment combination, for a total of 72 – 7.1 m<sup>2</sup> plots (Rich et al., 2015). Within each plot, 121 seedlings of 11 tree species were planted into the remaining herbaceous vegetation in a gridded design (Rich et al., 2015). Above-ground biomass was warmed using a Temperature Free-Air-Controlled Enhancement System (T-FACE) and below-ground biomass was warmed via buried resistance-type heating cables (Rich et al., 2015). Above- and below-ground temperatures have been monitored and logged at 15-min intervals since spring 2008. In 2012, event-based rain exclosures were installed on nine plots in the open overstory replicates of the warming experiment, which allowed for safe and reliable removal of rainfall. Mean annual rainfall exclusion from June to September ranges from 42 to 45%.

We collected plant samples from the B4WarmED project during two different time periods. On July 14, 2013, we collected samples of balsam fir and paper birch that were grown under closed overstory and three warming treatments, and on July 15, 2014, we collected samples of balsam fir, paper birch, trembling aspen (Populus tremuloides), and red maple (Acer rubrum) grown under open overstory in the three warming treatments and two precipitation treatments. Where possible, we collected recentgrowth foliar tissue from two plants per species within each replicate plot. However, some replicates contained either one or no plants with enough leaf tissue to sample. Samples sizes were particularly small during 2014, so we were forced to group individual warming treatments (ambient, +1.7◦C, +3.4◦C) into a binary response (ambient temperature vs. elevated temperature). All plant samples were collected within a 2-h time period. Upon collection, samples were flash frozen with dry ice, and subsequently stored in a −80◦C freezer to minimize chemical degradation. We broadly refer to samples collected from the B4WArmED project as our "experimental" samples.

To investigate how temperature and light conditions may interact to influence phytochemical production in a natural forest environment, we collected samples of balsam fir, paper birch, trembling aspen, and beaked hazel (Corylus cornuta) from open and closed canopy environments across two regions

in northeastern Minnesota (**Figure 1**). These regions exhibit differences in mean-maximum summer temperature (maximum daily temperature averaged across June, July, and August) of approximately 5.5◦C (**Supplementary Table S1**). On July 14, 2015, we collected a minimum of 3 biological replicates from each species within each set of abiotic conditions. The sampling design consists of a 2 (overstory – open and closed) × 2 (temperature – warm and cool) design with three plot replicates per treatment combination, for a total of 12 – 400 m<sup>2</sup> plots. Open-canopy plots allowed us to evaluate high-light conditions on production of PSMs and were located in areas that were clear-cut in 2006 (i.e., open overstory), while closed-canopy plots were located in areas that experienced no known overstory disturbance since at least 1985 (i.e., closed overstory). Thus, light conditions for all plots were based on whether the overstory was open (i.e., high light) or closed (i.e., low light). Temperature logger data collected for a parallel study from similar plot types suggest that average high temperatures from May 1, 2015 to July 14, 2015 ranged from 30.4◦C in low-light plots in the cool region to 36.6◦C in high-light plots in the warm region. All field samples were collected on the same day, within an 8-h period. Upon collection, samples were flash frozen with dry ice, and subsequently stored in a −80◦C freezer. For brevity, we occasionally refer to samples collected throughout northeast Minnesota as "observational" samples.

#### Study Organisms

Balsam fir is a mid- to large-sized species of conifer, growing to 26 m in height, with shallow roots (Smith, 2008). It is highly vulnerable to drought, fire, and spruce budworm (Choristoneuro fumiferana) infestations (Engelmark, 1999), and modest climate warming has been shown to decrease net photosynthesis and growth by as much as 25% (Reich et al., 2015). Paper birch can grow to 28 m in height (Smith, 2008) and is drought and shade intolerant (Iverson and Prasad, 1998; Iverson et al., 2008). While it can grow rapidly and live to be 250 years of age, seedlings need significant light to prosper (Kneeshaw et al., 2006). Elevated temperatures have been shown to influence foliar nitrogen, lignin, and condensed tannins in both paper birch and trembling aspen with the specific response varying as a function of species and climate (Jamieson et al., 2015). Trembling aspen is one of the most widespread tree species in North America and occurs on a wide-range of soil types and in various climatic conditions (Smith, 2008). It is sensitive to both drought and shade (Iverson and Prasad, 1998; Iverson et al.,

2008) and may become increasingly vulnerable to other potential stressors under conditions of drought and high temperatures (Worrall et al., 2008). Red maple is a moderately large tree, growing to 29 m in height and is known to be tolerant to a wide-range of precipitation conditions, from drought to seasonal flooding (Smith, 2008). While this species is expected to prosper under future climate scenarios (Iverson and Prasad, 1998; Iverson et al., 2008) and performed well under experimental warming (Reich et al., 2015), both prolonged flooding and severe drought have been shown to result in senescence and decreased growth, respectively (Nash and Graves, 1993). Beaked hazel, a shadetolerant shrub that can grow to 4 m tall, is a common understory species in both conifer and deciduous forests and occurs almost exclusively in fire prone habitats (Smith, 2008). Beaked hazel is highly sensitive to fire and previous work suggests that growth may be limited by soil moisture (Johnston and Woodard, 1985).

#### Metabolite Analysis

Tissue samples were lyophilized for 72 h and then homogenized and extracted using 25 mg (+/−2.5 mg) of each sample. Homogenization and extraction were performed for 5 min at a frequency of 1500 Hz with 1 ml of 70% isopropyl alcohol at –20◦C using a bead mill and 2.8 mm tungsten carbide beads (Sped Sample Prep GenoGrinder 2010, Metuchen, NJ, United States). Samples were then subjected to centrifugation at 16,000 × g for 5 min. The supernatant was then removed and subjected to an additional centrifugation step, 16,000 × g for an additional 5 min, and the supernatant was collected for subsequent analysis. Finally, 20 µL of each supernatant sample was removed and pooled to use as a control. All samples were then stored at –80◦C.

We analyzed samples with liquid chromatography mass spectrometry (LC-MS) using an Ultimate 3000 UHPLC (ultrahigh-performance liquid chromatography) system coupled to a Q Enactive hybrid quadrupole-Orbitrap mass spectrometer with a heated electrospray ionization (HESI) source (Thermo Fisher Scientific, Bremen, Germany). We injected 1 µL of each sample per analysis onto an ACQUITY UPLC HSS T3 column, 100 Å, 1.8 µm, 2.1 mm × 100 mm (Waters, Milford, MA, United States) using a gradient composed of solvents A: 0.1% formic acid and B: acetonitrile. Specifically, 0–2 min, 2% B; 6 min, 24% B; 9 min, 33% B; 12 min, 65% B; 16 min, 80% B; 20 min 93% B; 21 min 98% B; 22 min 98% B; 23 min 2% B; 23–25 min 2% B. Samples were analyzed in a randomized order to minimize systematic bias from instrument variability and carryover. Fullscan analysis was performed using positive/negative ion polarity switching, a 115–1500 m/z scan range, a resolution of 70,000 (at m/z 200), maximum fill times of 100 ms, and target automatic gain control (AGC) of 1 × 10<sup>6</sup> charges. Ion fragmentation was performed using a higher-energy collision dissociation (HCD) cell and resulting MS/MS data were collected using a resolution of 17,500, maximum fill times of 100 ms, and an AGC target of 2 × 10<sup>5</sup> charges. Normalized collision energies (NCE) ranged from 10 to 45 in increments of 5. All data were collected using Xcalibur version 2.2 (Thermo Fisher Scientific, Bremen, Germany).

### Example Compounds

To determine which chemical features varied consistently and significantly among each treatment and species group, we aligned, smoothed, background subtracted, and analyzed all chromatographic data using analysis of variance (α = 0.001) via Genedata 7.1 (Genedata, Basel, Switzerland). We assigned putative metabolite identities only to the features found to be significantly abundant (ANOVA, α = 0.001) with an exact mass and higher-energy collisional dissociation (HCD) MS/MS fragmentation spectra. We determined molecular formulae by using exact mass to calculate the most probable elemental composition for each feature (**Supplementary Table S2**). We then manually interpreted HCD spectra collected at numerous collision energies (**Supplementary Figures S1–S3**), and compared these to the MassBank database using MetFusion (Gerlich and Neumann, 2013). Where possible, we confirmed the identity of individual compounds via comparison to an authenticated standard (Sigma-Aldrich) and assigned other putative identities by matching molecular formulae to those of previously observed metabolites in Betula (Julkunen-Tiitto et al., 1996) and Abies (Otto and Wilde, 2001). Specifically, we analyzed changes in the relative abundance of catechin and an unspecified terpene acid in paper birch and two unspecified diterpene resin acids in balsam fir. The identification of catechin was confirmed by comparison of accurate mass, LC-retention and MS/MS fragmentation properties of commercially available standard compounds for both catechin and its frequently associated isomer epicatechin which were distinguishable by chromatographic separation. There has been a great deal of work investigating the biological and ecological activity of catechin and terpenoid-based metabolites (Tahvanainen et al., 1985; Gershenzon and Croteau, 1992; Berg, 2003; Stolter et al., 2005); and as a result, we expect our results regarding these compounds to be broadly relevant.

#### Data Processing and Statistical Analysis

Data processing and statistical analyses were conducted using R 3.5.0 (R Core Team, 2017). To initiate data processing, we used the xcmsRaw function in the xcms package (Smith et al., 2006; Tautenhahn et al., 2008; Benton et al., 2010) to read our raw mzML files into R. After separating our data by polarity using the split function in the base package, we used the findPeaks.centwave function for peak detection, which we parameterized as follows: ppm = 2, peakwidth = c(5,20), prefilter = c(1,15000000), mzCenterFun = "apex," integrate = 1, mzdiff = −0.001, fitgauss = F, snthresh = 10. Once peak detection was complete, we trimmed the resulting polarity-specific data frames based on retention time and retained only those peaks detected between 1 and 21 min.

A major shortfall of employing LC-MS to perform "untargeted profile analysis," as we did here, is the production of two independent but partially overlapping datasets resulting from ion polarity switching. While polarity switching is useful for detection of features that can only be detected via either positive or negative ionization, some features are detectable under both ionization modes, therefore resulting in two independent data sets containing a small subset of common features. Moreover,

interpretation of statistical results is challenging due to the presence of parallel sets of analyses with common features contributing to both. To alleviate these issues, we combined positive and negative polarities using the find.matches function in the Hmisc package (Harrell and Dupont, 2018). The find.matches function allows one to identify which rows in a data matrix align with those in a separate, identically formatted matrix by allowing the user to define a tolerance level for the numerical columns in each matrix. Thus, to determine our common features in the positive and negative ionization datasets that result from LC-MS, we created two matrices for positive and negative polarity, containing three separate columns – the mass of each detected peak, an assigned name for each peak, and retention time. To ensure that corresponding features from each ionization mode were capable of alignment, we subtracted 2.1046, roughly the mass of two protons, from all masses in the positive polarity dataset. For those features identified as common among both ionization modes, we retained peak data from the polarity exhibiting greatest mean intensity across all samples. We then assigned new peak names to identify which peaks were present in either positive or negative polarity vs. those that were found in both. All output created using the find.matches function was manually checked to ensure that all peaks identified as having a match in one polarity, had their mate identified as a match in the other.

We used permutational MANOVA (perMANOVA, Anderson, 2001) to compare PSM profiles between abiotic conditions. When analyzing PSM profiles, differences were estimated using Canberra dissimilarity matrices (Dixon et al., 2009). Analysis was performed with the adonis function (from the vegan package, Oksanen et al., 2015), which allowed us to account for our blocked sampling design via the strata argument. Both differences in the centroids among conditions or differences in multivariate dispersion can lead to statistically significant results when using perMANOVA. To determine whether differences among centroids were contributing to perMANOVA results, we created mean dissimilarity matrices using the meandist function and we used the betadisper function to assess multivariate homogeneity of variance (i.e., dispersion, Oksanen et al., 2015). We used non-metric multidimensional scaling (NMDS, Kruskal, 1964) to visualize differences in PSM profiles among conditions, which we performed using the metaMDS function in the vegan package (Oksanen et al., 2015). We set our dimensionality parameter (k) to 2 and projected condition-specific effects onto NMDS plots using the ordiellipse function to plot 95% confidence ellipses based on standard error (Oksanen et al., 2015).

To evaluate treatment-induced changes to PSM diversity, we calculated phytochemical richness based on the presence and absence of individual compounds, then tested the main effect of treatment on richness with block (experimental samples) or site ID (observational samples) as our random effect using linear mixed-effects models (lme function within the nlme package, Pinheiro et al., 2015). To analyze phytochemical turnover (i.e., the degree of overlap between the phytochemical profiles of individual plants across and between conditions), we created dissimilarity matrices based on binary datasets representing the presence or absence of each feature using Jaccard's Index. We evaluated condition-specific differences in phytochemical turnover using perMANOVA via the adonis function, and evaluated the influence of multivariate centroids and homogeneity of variance on perMANOVA results as detailed above (Anderson, 2001; Oksanen et al., 2015).

Weather data from CFC shows that ambient air temperature, cumulative precipitation from 1 January to collection date, and leaf surface temperature were not statistically different between 2012 and 2013 or between specific sample sets (2013 – closed overstory, 2014 – open overstory). However, soil moisture and soil temperature vary strongly between years and sample sets, and differences between experimental and observational samples are likely to be even greater. Thus, samples collected during different years were analyzed independently of one another as individual data sets.

For analytical and visualization purposes, the condition or set of conditions assumed to impart the least amount of metabolic change during each year was labeled as our reference group, to which all other conditions were compared for that sample year. For Year 1 (2013), we designated "ambient" as our reference category, while samples grown under ambient temperature and ambient precipitation were designated as our reference category for Year 2 (2014). We designated samples collected from cold region, low-light conditions as our reference category for Year 3 (2015). To help visualize how different abiotic conditions may influence PSM production in different species, we calculated the number of chemical features that increased and decreased by ≥ 75%, relative to our reference category and created scaled Venn Diagrams representing these relationships.

Finally, we used linear mixed-effects models to test the main effect of abiotic condition on the relative abundance of our example compounds, with sample block as our random effect for experimental samples and plot ID as our random effect for observational samples (lme function within the nlme package, Pinheiro et al., 2015). These models tested whether combinations of abiotic factors influence the abundance of our known example compounds.

# RESULTS

#### Temperature

The influence of temperature was both species and context dependent. In closed overstory (Year 1), when compared to ambient, warming-induced changes to the phytochemical profile of balsam fir were not statistically significant, whereas paper birch exhibited warming-induced shifts to phytochemical profiles, thereby leading to distinct PSM profiles for the warming treatment. Analysis of multivariate dispersion and mean-dissimilarity matrices both suggest that differences in paper birch were due to temperature-induced changes in the centroid rather than dispersion (**Table 1**). NMDS plots reveal minor overlap between temperature conditions in paper birch, and balsam fir grown under moderate and high-temperatures show strong overlap with plants grown in ambient temperatures but minor overlap with each other (**Figure 2**). Warming had no effect on phytochemical richness in either species but did


 na indicates not applicable.

were compared.

influence phytochemical turnover in paper birch (**Table 1**). In open overstory (Year 2), warming had no influence on PSM profiles or PSM diversity (i.e., richness or turnover), regardless of species (**Table 1**). NMDS plots support these findings in that there is no discernible relationship between temperature and PSM profiles, regardless of species (**Figure 3**). In observational samples collected throughout northeast Minnesota (Year 3), temperature on its own had no influence on plant PSM profiles or phytochemical richness values. However, phytochemical turnover was significantly different in plants from different temperature regions in paper birch (perMANOVA, F = 5.912, r <sup>2</sup> = 0.179, P = 0.0003) and trembling aspen (perMANOVA, F = 3.322, r <sup>2</sup> = 0.156, P = 0.0012). NMDS plots suggest that each species responds differently to combinations of temperature and light (i.e., canopy; **Figure 4**). Balsam fir produces distinct PSM profiles as a function of ambient light conditions (i.e., open vs. closed canopy), but only within the cool region, while paper birch and trembling aspen appear to have distinct PSM profiles for each combination of conditions. Conversely, beaked hazel exhibits no discernible pattern across different conditions.

Venn diagrams created to help visualize the influence of different abiotic conditions for Year 1 samples suggest that the high-temperature (+3.4◦C) treatment induced a greater response from both balsam fir and paper birch than the moderate-temperature (+1.7◦C) treatment. Specifically, the high-temperature treatment led to more features that either increased or decreased in relative abundance by 75% or more when compared to ambient or moderate-temperature treatments (**Table 2** and **Supplementary Figures S4–S6**).

#### Interactive Effects of Different Abiotic Conditions

In our Year 2 samples, the combination of drought and elevated temperature had no influence on PSM profiles or any aspect of phytochemical diversity, regardless of species (**Table 1**). These results were supported by NMDS plots (**Figure 3**). Additionally, Venn diagrams suggest large-magnitude increases or decreases in relative abundance of PSMs did not follow an obvious pattern that could be attributed to different conditions. There appears to be a high degree of overlap across conditions in those compounds that exhibit increases in relative abundance of ≥ 75%, while less overlap occurs among compounds exhibiting large declines in relative abundance. Furthermore, the influence of drought on the decline of relative abundance by ≥ 75% appears to be more distinct than that of either warming or warming and drought together (**Table 2** and **Supplementary Figures S4**–**S6**).

In observational samples from throughout northeast Minnesota (Year 3), when evaluating the effects of high temperature and light combined, balsam fir appears to create unique PSM profiles in response to different light conditions (i.e., open vs. closed canopy), but only within the cool region, while paper birch and trembling aspen appear to have distinct PSM profiles for each condition. Beaked hazel exhibits no discernible pattern (**Figure 4**). Phytochemical richness did not vary as a function of light conditions or temperature region. However, phytochemical turnover in balsam fir was significantly influenced by conditions of high light (i.e., open canopy; **Table 3**). When analyzing the interactive effects of light conditions and temperature region, all species exhibited significant differences in their PSM profile (**Table 3**), with only trembling aspen exhibiting significant differences in multivariate dispersion as a function of the combination of light condition and temperature region (**Table 3**). Although phytochemical richness was not influenced by the combined effects of temperature region and light conditions, phytochemical turnover was influenced in paper birch and trembling aspen and a marginal, non-significant trend was present in beaked hazel (**Table 3**).

Patterns in Venn diagrams detailing the influences of different conditions during Year 2 are difficult to discern, as different plant species appeared to respond to varying conditions in different ways (**Table 2** and **Supplementary Figure S5**). Drought led to more features increasing by ≥ 75% in balsam fir and paper birch, while elevated temperature led to the large-magnitude increase of more features in trembling aspen (**Table 2** and

**Supplementary Figure S5**). In red maple, the combination of drought and elevated temperature had the greatest influence on large-magnitude increases in relative abundance (**Table 2** and **Supplementary Figure S5**). The combination of drought and warming led to more large-magnitude declines in relative abundance in balsam fir and paper birch, while drought had a greater impact on red maple and trembling aspen (**Table 2** and **Supplementary Figure S5**). In observational samples (Year 3), the combination of high-light conditions and warmer temperatures led to more large-magnitude shifts in relative abundance (i.e., increasing and decreasing by 75% or more), regardless of species (**Table 2** and **Supplementary Figure S6**).

#### Example Compounds

In closed-overstory conditions (Year 1), warming resulted in significant declines in both catechin and terpene acid in paper birch but had no influence on either compound in balsam fir (**Figure 5** and **Supplementary Table S3**). In highlight conditions (Year 2), neither of the compounds in either species exhibited a significant, condition-specific change in relative abundance. However, terpene acid in paper birch was completely absent from all samples collected from high-light plots (**Figure 6** and **Supplementary Table S3**). In observational samples from throughout northeast Minnesota (Year 3), neither compound in balsam fir exhibited significant changes in relative abundance due to light conditions, temperature region, or their interaction. In paper birch, however, the interactive effects of high-light conditions and warmertemperatures resulted in a more than 250% increase in the relative abundance of catechin, while terpene acid exhibited no response, regardless of treatment (**Figure 7** and **Supplementary Table S3**).

FIGURE 4 | Non-metric multidimensional scaling (NMDS) plots detailing the influence of varying light and temperature conditions on PSM profiles of (A) balsam fir, (B) paper birch, (C) beaked hazel, and (D) trembling aspen. Ellipses represent 95% confidence intervals, based on standard error. Each species appears to respond to different abiotic conditions in a unique manner. Balsam fir appears to create unique PSM profiles in high-light conditions when compared to our reference group (closed canopy, low temperature), while paper birch and trembling aspen appear to have distinct PSM profiles for each set of conditions. Beaked hazel exhibits no discernible pattern.

# DISCUSSION

Our study is among the first to explicitly show that combinations of abiotic drivers (often potential stressors) in forest plants can lead to broad phytochemical responses that are distinct from those that result from single abiotic factors and that different species of woody plants respond to complex sets of conditions in variable ways. In our experimental samples, warming under closed canopy led to distinct PSM profiles in paper birch but not balsam fir, with paper birch exhibiting increased phytochemical turnover. Warming under open canopy had no influence on PSM profiles or any aspect of phytochemical diversity. In our observational samples collected across northeast Minnesota, warmer temperatures had no influence on PSM profiles but did lead to significant phytochemical turnover in paper birch and trembling aspen. When elevated temperature was combined with drought in Year 2 of our experimental samples, we found no influence on PSM profiles or phytochemical diversity. However, temperature variation combined with high-light conditions in our observational samples resulted in condition-specific profiles for all species and led to significant phytochemical turnover in all but beaked hazel. In general, our results indicate that the phytochemical response of plants to varying combinations of abiotic factors cannot be directly extrapolated from the response of plants to individual factors. Perhaps more importantly, our results provide evidence that heterogeneity in the abiotic environment influences secondary metabolism in woody plants via a range of complex and highly variable responses.

Few studies to date have explicitly studied the influences of heterogeneity in the abiotic environment on phytochemical diversity, and specifically, phytochemical turnover. However, it has been hypothesized that variability in which compounds


TABLE 2 | Number of chemical features that increase and decrease in relative abundance by ≥ 75% as a function the dominant stress condition.

In most scenarios, the stress condition that led to large-scale increases in relative abundance was different than that which led to large-scale decreases. "Number affected" displays the number of chemical features that either increased or decreased by ≥ 75% for the given species and stress condition.

are either present or absent may be an adaptation for variable environments, thereby decreasing vulnerability of plants to a range of potential stress conditions, including herbivory (Laitinen et al., 2000; Cheng et al., 2011). Here, we found that in some plants species, different combinations of abiotic factors can affect which compounds are either present or absent, thus leading to phytochemical turnover. For example, compounds that are absent in one set of conditions may become present within a slightly different set of conditions, or vice versa. The potential for this to occur was apparent when our example terpene acid decreased in paper birch plants subjected to experimentally elevated temperature in closed canopy but went completely undetected in plants subjected to experimental warming and drought in open canopy and exhibited no change at all in our observational samples from throughout northeast Minnesota. Suppression of individual compounds due to varying stress conditions has been observed in other studies as well. For instance, proline, which is thought to play an important role in protection from drought, is severely suppressed when plants are simultaneously subjected to drought and high temperatures (Rizhsky et al., 2004). While individual compounds can play an important role in the survival of plants subjected to a range of biotic and abiotic conditions, a plant's phytochemical profile imparts a metabolic framework that can determine the biological role and strength of individual compounds (Dyer et al., 2003; Richards et al., 2010; Gershenzon et al., 2012; Jamieson et al., 2015). Here, we show that individual compounds as well as the phytochemical context within which they operate can both be altered by variations in the abiotic environment.

Plants produce thousands of individual compounds, and variations in the relative abundance of many of these compounds can have a wide-range of effects on the biotic interactions plants have with other organisms. Catechin, which is a phenolbased precursor to proanthocyanidins (i.e., condensed tannins), is widely considered an antiherbivore defensive compound (Tahvanainen et al., 1985; Berg, 2003; Stolter et al., 2005) and can have a significant, negative impact on the development of forest pests (Roitto et al., 2009). Catechin also has antimicrobial and allelopathic effects, and plants with decreased catechin production may be at a competitive disadvantage for nutrients within the soil as it can inhibit the growth and germination of neighboring plants (Veluri et al., 2004; Inderjit et al., 2008). Terpene acids, including diterpene resin acids, are considered strong antifeedants (Ikeda et al., 1977), and the ingestion of forage with elevated concentrations of diterpenoids can result in slower development times and significantly higher mortality in herbivorous larvae (Larsson et al., 1986). Here, we show that different compounds have individualized responses based on the micro-environmental conditions that are present.

In balsam fir, warming alone led to consistent, albeit nonsignificant declines in the mean relative abundances of both resin acids. When high temperatures were combined with other abiotic factors (i.e., drought and light), resin acid 1 exhibited consistent but non-significant increases, while resin acid 2 was more variable, exhibiting no consistent trend. In paper birch, both example compounds exhibited significant changes in relative abundance as a function of different factors. While elevated temperature alone led to significant declines in catechin, the combination of elevated temperature and high light led to a more than 250% increase in relative abundance. Our example terpene acid declined with warming and was undetectable when we tried to assess the effects of drought. This particular scenario provides an example of how individual compounds may "wink in or out" due to variation in the abiotic environment.

Numerous studies have reported that high-temperature and drought interact to alter PSM production in plants (Craufurd and Peacock, 1993; Savin and Nicolas, 1996; Jiang and Huang, 2001; Rizhsky et al., 2002, 2004). Thus, we were surprised to find no interaction between drought and warming in our study. It is important to note, however, that the extremes of those treatments employed by other studies are typically greater than what we test here, sometimes increasing temperature to more than 40◦C (Rizhsky et al., 2002) and withholding water altogether for extended periods (Jiang and Huang, 2001). In our study, mean soil moisture was lower during 2014 than 2013, with mean soil temperatures being higher (unpublished


TABLE3|SummaryofresultsforobservationalsamplesusedtoassesstheinfluencesoftemperatureregionandoverstoryonPSMprofilesand

fpls-09-01257 September 4, 2018 Time: 9:5 # 12

FIGURE 5 | Relative change in abundance (%) for specific PSM compounds when compared to our reference treatment for Year 1 (ambient temperature) for (A) balsam fir and (B) paper birch in closed overstory. Neither resin acid in balsam fir was influenced by warming. In paper birch, both catechin and terpene acid declined with warming relative to ambient. Error bars represent the 95% boot-strapped confidence intervals and relative abundances significantly different than those found in the baseline treatment are identified by an asterisk (<sup>∗</sup> ).

data). Surprisingly, air temperature and leaf-tissue surface temperature during late spring/early summer (May 1 to July 15) were indistinguishable between samples years and plot types (2013 closed canopy vs. 2014 open canopy), and cumulative precipitation during the first half of each year (January 1 to July 15) was also indistinguishable (unpublished results). Combinations of abiotic factors can have one dominant factor that defines the phytochemical response of affected plants, and drought, when present, may dominate the influence of combinations of abiotic factors. Considering this, our inability to identify any treatment-specific influence on PSM profiles or phytochemical diversity may be due to low soil moisture during 2014. If plants from which samples were collected from in 2014 were experiencing some level of drought stress due to low soil moisture, this signal may have preempted any potential phytochemical response that might have occurred due to treatment.

When considering the influence of abiotic conditions on largescale shifts in relative abundance (increases or decreases ≥ 75% relative to our reference group), greater increases in temperature (+3.4◦C) appeared to have a greater influence than moderate increases (+1.7◦C). When present, drought, either alone or in combination with elevated temperature, dominated all but one of the large-scale shifts we assessed (Year 2), and in our observational samples, high-light conditions, either alone or in combination with elevated temperature, dominated all

).

of the large-scale shifts we assessed in which it was present (Year 3). As noted above, numerous studies have reported that drought has a defining impact on plants' phytochemical profiles, even when in combination with other abiotic drivers, such as elevated temperature and high light. Moreover, in our Year 1 samples, elevated temperature led to both largescale increases and large-scale decreases in relative abundance. However, the number of compounds exhibiting these shifts was substantially smaller when compared to the number of compounds influenced by the abiotic conditions evaluated in either Year 2 of our experimental samples or our observational samples (Year 3). Outside of Year 1, during which we tested only the effects of elevated temperature, it was rare when the same abiotic condition simultaneously dominated both large-scale increases and large-scale decreases in relative abundance, suggesting that different combinations of abiotic factors may influence upregulation and downregulation of different compounds.

Changes in the abundance and diversity of secondary metabolites within a plant's phytochemical profile may alter biotic interactions, potentially leading to broad-scale ecological change. For example, while some herbivores respond negatively to forage with a higher diversity of PSMs, others appear to target these plants in an effort to alleviate costs associated with external stressors via their pharmacological benefits (Forbey and Hunter, 2012). Additionally, numerous studies have reported that phytochemical diversity within a plant community is positively correlated with community diversity across multiple trophic levels (Jones and Lawton, 1991; Richards et al., 2015), influencing invertebrate predators and parasitoids, and potentially extending to vertebrate predators as well (Dicke et al., 2012).

While the consequences of different abiotic conditions on phytochemical diversity remain unpredictable, our results demonstrate that the phytochemical response of plants to combinations of abiotic factors cannot be extrapolated from that of individual factors. For instance, while warming alone may have a very specific influence on some compounds, when in combination with additional abiotic factors such as drought and light, warming may lead to highly variable and unpredictable response (Mittler, 2006), making it increasingly difficult to predict the performance of woody plants in a changing environment. Regardless, previous research suggests that changes in phytochemical production induced by variability in abiotic conditions can influence both tree resistance and pest performance traits (Jamieson et al., 2015), potentially altering the frequency and intensity of insect outbreaks (Schwartzberg et al., 2014). Elevated temperatures by themselves have been shown to reduce the competitive abilities of more southern boreal tree species when compared to co-occurring species adapted to warmer climates (Reich et al., 2015). Climate-induced changes to phytochemistry may lead to shifts in the competitive landscapes for cold-adapted trees and shrubs, potentially altering their ability to compete for resources and defend against pests and pathogens in novel climatic conditions. However, because individual compounds and the metabolic profiles of which they are a part are differentially influenced by abiotic factors and combinations of these factors, predicting how forest plants will respond to novel environmental conditions will be challenging.

The majority of biotic interactions between plants and other organisms are chemically mediated, and recent climate change has challenged our understanding of the mechanisms underlying these interactions. The primary objective of this study was to determine how warming influences plant production of secondary metabolites and how combinations of additional abiotic factors may modulate this effect.

Here, we show that heterogeneity in a range of abiotic factors broadly influence secondary chemistry in plants thereby leading to condition-specific phytochemical profiles. If our results are typical of plant responses, abiotically induced changes to secondary chemistry in woody plants could influence their rate of range expansion or contraction under novel climate scenarios. Additionally, our results contribute to current efforts to understand how continued warming will influence plants and the biotic interactions that serve as the foundation for a wide range of ecosystem processes. In the future, studies monitoring physiological changes in conjunction with global shifts in PSM profiles would provide insights into mechanisms underlying biotic interactions mediated by the local environment. As spatial and temporal patterns in the global abiotic environment continue to shift, it is imperative that we continue to learn as much as we can about the phytochemical response of plants to these changes.

#### AUTHOR CONTRIBUTIONS

JB, SB, AH, RaM, ReM, and JF formulated the study idea and developed the study methods while PR and ReM established the experimental study sites critical for the execution of this study. JB performed all the sample collection, while JB and SB performed the analytical chemistry and pre-statistical data processing. JB, RaM, and JF statistically analyzed the data. JB and SB wrote the initial draft of the manuscript. All the authors contributed to the manuscript revisions and approved the final manuscript.

#### REFERENCES


#### FUNDING

This work was supported by the Office of Science (BER), United States Department of Energy (385 Grant No. DEFG02- 07ER64456), the NSF Plant Genome Research Program (Grant Nos. IOS-0923960 and IOS-1238812), the United States Environmental Protection Agency, Science to Achieve Results Graduate Fellowship (Grant No. F13B20220), the Minnesota Environment and Natural Resources Trust Fund (Grant No. M.L. 2014, Chap. 226, Sec. 2, Subd. 05l), the Institute on the Environment at the University of Minnesota, and the Minnesota Agricultural Experiment Station Projects MIN-42-030, MIN-42-060 and MIN-41-020.

# ACKNOWLEDGMENTS

We would like to thank the Minnesota Supercomputing Institute (UMN MSI) at the University of Minnesota for software and computational support. Suggestions by L. Frelich, D. Fox, and G. DelGiudice significantly improved this study. This work would not have been possible without the logistical and technical support of K. Rice, A. Stefanski, and M. Wethington.

# SUPPLEMENTARY MATERIAL

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




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

Copyright © 2018 Berini, Brockman, Hegeman, Reich, Muthukrishnan, Montgomery and Forester. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Fertilization Changes Chemical Defense in Needles of Mature Norway Spruce (*Picea abies*)

Line Nybakken<sup>1</sup> \*, Marit H. Lie<sup>1</sup> , Riitta Julkunen-Tiitto<sup>2</sup> , Johan Asplund<sup>1</sup> and Mikael Ohlson<sup>1</sup>

<sup>1</sup> Faculty of Environmental Sciences and Natural Resource Management, Norwegian University of Life Sciences, Ås, Akershus, Norway, <sup>2</sup> Department of Environmental and Biological Sciences, University of Eastern Finland, Joensuu, Finland

*Edited by:*

Bartosz Adamczyk, University of Helsinki, Finland

#### *Reviewed by:*

Johanna Witzell, Swedish University of Agricultural Sciences, Sweden Ivika Ostonen, University of Tartu, Estonia

> *\*Correspondence:* Line Nybakken line.nybakken@nmbu.no

#### *Specialty section:*

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

*Received:* 06 April 2018 *Accepted:* 18 May 2018 *Published:* 07 June 2018

#### *Citation:*

Nybakken L, Lie MH, Julkunen-Tiitto R, Asplund J and Ohlson M (2018) Fertilization Changes Chemical Defense in Needles of Mature Norway Spruce (Picea abies). Front. Plant Sci. 9:770. doi: 10.3389/fpls.2018.00770 Nitrogen availability limits growth in most boreal forests. However, parts of the boreal zone receive significant levels of nitrogen deposition. At the same time, forests are fertilized to increase volume growth and carbon sequestration. No matter the source, increasing nitrogen in the boreal forest ecosystem will influence the resource situation for its primary producers, the plants, with possible implications for their defensive chemistry. In general, fertilization reduces phenolic compound concentrations in trees, but existing evidence mainly comes from studies on young plants. Given the role of the phenolic compounds in protection against herbivores and other forest pests, it is important to know if phenolics are reduced with fertilization also in mature trees. The evergreen Norway spruce is long-lived, and it is reasonable that defensive strategies could change from the juvenile to the reproductive and mature phases. In addition, as the needles are kept for several years, defense could also change with needle age. We sampled current and previous year needles from an N fertilization experiment in a Norway spruce forest landscape in south-central Norway to which N had been added annually for 13 years. We analyzed total nitrogen (N) and carbon (C), as well as low-molecular phenolics and condensed tannins. Needles from fertilized trees had higher N than those from controls plots, and fertilization decreased concentrations of many flavonoids, as well as condensed tannins in current year needles. In previous year needles, some stilbenes and condensed tannins were higher in fertilized trees. In control trees, the total phenolic concentration was almost five times as high in previous year needles compared with those from the current year, and there were great compositional differences. Previous year needles contained highest concentrations of acetophenone and stilbenes, while in the current year needles the flavonoids, and especially coumaroyl-astragalins dominated. Condensed tannins did not differ between current and previous year needles from control trees. In conclusion, the phenolic defense of current year needles of mature P.abies trees was strongly changed upon fertilization. This may imply that nitrogen deposition and forest fertilization leave forests less robust in a time when pests may take advantages of a changing climate.

Keywords: nitrogen, fertilization, phenolics, chemical defense, *Picea abies*, spruce, conifers

#### Nybakken et al. Fertilization Changes Spruce Defence

# INTRODUCTION

Boreal and temperate forests soils naturally have a low availability of nitrogen (N) and N is often the primary growthlimiting nutrient in these soils (Tamm, 1991; Binkley and Fisher, 2013). However, the increased nitrogen deposition rates over the last 150 years have greatly altered this in parts of these biomes (Galloway and Cowling, 2002; Meunier et al., 2016), and it is wellknown that N deposition often directly enhances forest growth (Meunier et al., 2016).Continental and northern regions of the European boreal forests, however, are relatively little affected by N deposition (e.g., Gundale et al., 2011). Yet, there is increasing interest in fertilization of these areas (e.g., Bergh et al., 2014) to increase timber volumes, and more recently, as a measure to sequester more carbon (C) (see e.g., Gundale et al., 2014).

No matter the source for the increased N availability, the consequences for the ecosystem are potentially large. N deposition has been shown to increase C uptake by trees in the range of 15–60 kg of C per kg of N deposited (Högberg et al., 2006; Pregitzer et al., 2008; Thomas et al., 2010; de Vries et al., 2014; Gundale et al., 2014). N deposition may also increase tree foliar N concentration by 22% for conifers and 13% for deciduous trees (reviewed by Throop and Lerdau, 2004). The general power of N is pointed out in the recent review by Meunier et al. (2016), who conclude that increased availability of N in the base of forest food webs results in cascade effects, which in turn substantially affect other trophical levels through alterations in food web structure and functioning.

Forest fertilization was originally introduced to increase the vitality of forests, and indeed, most studies show an increase in diameter growth after fertilization (e.g., Bergh et al., 2014). At present, however, many plant pathogens are increasing in incidence and range (Parmesan, 2006; Tylianakis et al., 2008), and it is possible that elevated N-availability is one of the drivers behind this increase. Simultaneously, N deposition has been shown to have a positive effect on insect herbivore populations (e.g., Strengbom et al., 2005), increasing the individual performance and population growth rates by up to 30% (Throop and Lerdau, 2004). Fertilization increased the palatability of bilberry (Vaccinium myrtillus) to gray sided voles (Strengbom et al., 2003).

As N availability increase the tree's C requirements for growth, it may also influence the C availability for other metabolic processes with indirect implications for ecosystem functioning. For instance, trees produce large amounts of Cbased plant secondary metabolites (PSMs). One PSM group of major ecological interest is the phenolics, with functions as diverse as herbivory defense (e.g., Bryant et al., 1983; Coley et al., 1985), pathogen resistance (Witzell and Martín, 2008), allelopathy or symbioses signaling (Inderjit, 1996; Mandal et al., 2010), frost and drought hardiness (Samanta et al., 2011), and photodamage protection (Lois, 1994; Close and McArthur, 2002). Flavonoids, stilbenes, phenolic acids, and condensed tannins are important classes of phenolics in most forest trees. According to ecological theory, the levels of C-based PSMs in plants could be expected to change with increased N. The carbon nutrient balance (CNB) hypothesis (Bryant et al., 1983), and the related growth differentiation hypothesis (GDB) Loomis, 1932; Herms and Mattson, 1992, in short postulate that the restrictions on growth caused by low N availability will lead to an accumulation of photosynthates, which can be diverted to the production of PSMs. The protein competition model (PCM) (Jones and Hartley, 1999) further claims that the underlying mechanism is that production of proteins and phenolic compounds compete for their common precursor phenylalanine. These hypotheses are corroborated in a metaanalysis by Koricheva et al. (1998), who found that plant phenolics in general decrease after fertilization. Consistent with these hypotheses, another meta-analysis showed that species adapted to resource-rich habitats grow inherently faster and invest less in defenses than species adapted to less productive habitats (Endara and Coley, 2011).

Nevertheless, surprisingly few studies have tried to test experimentally how increased N availability affects the defensive capacity of species from inherently N limited ecosystems. Of these, only a handful involve the ecologically and economically important boreal conifers. Virjamo et al. (2014) observed an increase in catechins in the needles after fertilization of 1-year old Picea abies seedlings. All other phenolics stayed unchanged. Blodgett et al. (2005) and Edenius et al. (2012) found reduced concentrations of total needle phenolics in fertilized Pinus resinosa and P.abies, respectively. On the other hand, Kytö et al. (1996) saw no effect of fertilization on the total phloem phenolics of P.abies. These previous studies have either focused on very young trees (Virjamo et al., 2014), a mixture of different aged trees and needles (Edenius et al., 2012), or have only looked at the effect on total phenolics (Kytö et al., 1996; Blodgett et al., 2005), while individual compounds or groups of compounds have not been studied. This implies that there is very little scientific knowledge on how increased N availability may influence the within-plant variation in defensive capacity. Fertilization may possibly increase growth, but leave the tree more subject to attacks from different forest herbivores or other pests. Most boreal conifers keep their needles for two or more years, and use their needles as storage compartments for carbohydrates and nutrients, with generally highest concentrations in young needles (see e.g., Ohlson, 1995 and references therein). The total levels increase and the composition of PSMs change with needle age over one season (Slimestad and Hostettmann, 1996; Slimestad, 1998; Ganthaler et al., 2017a), as well as upon infection by a rust fungus (Ganthaler et al., 2017b). However, we lack information about how such dynamics in defensive strategies e.g., over needle age, are affected by abiotic environmental changes.

Here we report on phenolic concentrations, as well as total C and N from current and previous year needles of mature spruce (P.abies), fertilized over a 13-year period. We hypothesized that increased N availability would increase N concentrations in the needles, and lead to decreased levels of PSMs, but also affect individual compounds differently, resulting in compositional changes. We further hypothesized that previous year needles would have higher defense levels and be less subject to change upon fertilization than those from the current year.

#### MATERIAL AND METHODS

#### Study Area and Fertilization Experiment

Our study area was an old and multilayered boreal Norway spruce forest landscape situated 800 m a s l near Kittilbu in SE Norway (61◦ 10′N, 09◦ 09′E). The spruce trees varied in size and age from seedlings to 20 m tall and 220 year old trees. No logging have occurred during the last 65 years, but light selective loggings have occurred in the past. The ground vegetation was dominated by bilberry V. myrtillus and bryophytes such as e.g., Pleurozium schreberi, Hylocomium splendens, and Polytrichum commune. Gauslaa et al. (2008) and Davey et al. (2017) gives information about the vegetation and the edaphic and climatic conditions in the study area.

Nitrogen has been added annually since 2003 to ten 15 × 15 m (225 m<sup>2</sup> ) experimental plots in the forest landscape at a rate of 150 kg ha−<sup>1</sup> year−<sup>1</sup> in the form of granulated pellets that contained 24.6% N, 2% P, 6% K, and trace elements (YaraMila Fullgjødsel). Ten equal sized and unfertilized plots were established to serve as control. Individual plots were located between 50 and 350 m from the next plot. Background deposition of N for the years 2007–2011 amounts to ∼6 kg ha−<sup>1</sup> year−<sup>1</sup> (http://www.environment.no/topics/air-pollution/ acid-rain/maps-deposition-of-sulphur-and-nitrogen/).

#### Chemical Extraction

We sampled spruce needles on 25 June 2016 in paper bags containing silica. From each plot, we chose three approximately equal sized and visually vital trees, from the dominating crown layer. From these c. 10 needles from both the current and the previous year's cohorts were sampled. The needle samples were taken ∼1.2 m above the ground and from the north side of the tree, to get as little variation as possible due to light exposition. The paper bags were brought to the lab the same day, and were dried in an oven at 30◦C for 48 h and thereafter stored in plastic bags in the freezer at −18◦C. The plant material was later grinded to fine powder using a Retsch MM400 ball mill (Retsch, Haag, Germany) at 30 revolutions s−<sup>1</sup> for 30–180 s. From the resulting powder, we determined total carbon (C) and nitrogen (N) with a Micro Cube (Elementar Analysen, Hanau, Germany), using 5–6 mg plant material. For phenolic analysis, further subsamples ofc. 10 mg were extracted with 600 µl methanol (MeOH) and homogenized with 3–4 zirconium oxide balls at 5,000 rpm for 20 s on a Precellys 24 homogeniser (Bertin Technologies, Montigny-le-Bretonneux, France). Samples were then cooled on ice for 15 min before being centrifuged at 15,000 rpm min−<sup>1</sup> for 3 min (Eppendorf centrifuge 5417C, Eppendorf, Hamburg, Germany). The supernatant was transferred to a 10 ml glass tube with a Pasteur pipette. The residue was again dissolved in 600 µl MeOH, homogenized, and centrifuged in the same manner as above; the supernatant was removed, and the same extraction process was conducted two more times until both the residue and the supernatant was completely colorless. The combined supernatants were evaporated in a vacuum centrifuge (Eppendorf concentrator plus; Eppendorf, Hamburg, Germany), sealed, and stored in a freezer (−18◦C) until high performance liquid chromatography (HPLC) analysis. The residues were also stored in a freezer for further analysis on MeOH—insoluble condensed tannins.

#### HPLC Analyses

The dried extracts were dissolved in 200 µl MeOH with the help of an ultra sonic cleaner (mod. no. USC200TH; VWR International LLC, Randor, USA) and diluted with 200 µl ultra-clean water (USF ELGA Maxima HPLC; Veolia Water Technologies, Saint-Maurice, France). Samples were poured into 2 ml Eppendorf tubes and centrifuged at 15000 rounds min−<sup>1</sup> for 3 min before being forced through a syringe filter (GHP Acrodisc 13 mm Syringe Filter with a 0.45µm GHP membrane; PALL Corporation, Washington, USA) and sealed inside HPLC vials. Low molecular weight phenolics (LMWP) were analyzed using a HPLC system (Agilent Series 1200, Agilent Technologies, Waldbronn, Germany) with a G1312A binary pump, a G1329A autosampler, a G1316A thermoregulated column heater, and a G1315D diode array detector. As the stationary phase a Thermo Scientific column type was used (Thermo Fisher Scientific Inc., Waltham, USA), with a 50 × 4.6 mm internal diameter and filled with ODS Hypersil (3µm) particles. The mobile phase consisted of two solvents that eluted the samples by way of a gradient as in (Julkunen-Tiitto and Sorsa, 2001). The injection volume was 20 µl. The phenolic compounds were identified using a UHPLC quadrupole time-of-flight liquid chromatograph/mass spectrometer (UHPLC/Q-TOF MS) (6540 series, Agilent) (Supplementary Table S1; Virjamo et al., 2014). The eluents were 1.5% tetrahydrofuran + 0.25% acetic acid in Milli-Q ultrapure water (Eluent A) and 100% MeOH (Eluent B). The following gradient was used for eluent A: 0–1.5 min 100% A, 1.5–3 min 100–85% A, 3–6 min 85–70% A, 6–12 min 70– 50% A, 12–20 min 50% A, 20–22 min 50–0% A. The flow rate was 0.4 ml min−<sup>1</sup> and the injection volume was 0.5 µl. The Q-TOF parameters were: mass range 100–3,000 m/z; temperature of the drying gas and sheath gas 350◦C and flow rates 12 l/min and 11 l/min, respectively; nebulizer pressure 35 psi; capillary voltage 3,500 V; nozzle voltage 1,000 V; fragmentor voltage 80 V; skimmer voltage 65 V; octopole voltage 750 V. The reference m/z 922.0098 was used for accurate mass measurements. Accuracy of tentative compound identification (ppm) was calculated as (measured mass—calculated mass) × 10<sup>6</sup> /calculated mass (Supplementary Table S1). Compounds that were not identifiable using the Q-TOF MS were identified by comparing retention times and ultraviolet spectra to the literature.

The absorption spectra at 270 and 320 nm, along with respective retention times, were used to identify the chemical compounds and to calculate concentrations by comparing with commercial standards.

#### Analysis of Condensed Tannins (CT)

Concentrations of both MeOH-soluble and MeOH-insoluble CTs were identified using the acid butanol assay for proanthocyanidins described in Hagerman (2002). The HPLCvials were removed from the auto sampler maximum 48 h after analysis and from these 50–100 µl were used to determine the amounts of MeOH-soluble CTs. The amount of MeOH-insoluble CTs were analyzed from the residues left after the extraction process. The samples were put in 10 ml glass tubes along with enough MeOH to equal 0.5 ml in total (0.5 ml MeOH regardless for MeOH-insoluble tannins), then further mixed with 3 ml butyric acid (95% butanol, 5% hydrochloric acid), and 100 µl iron reagent (2 M HCL with 2% ferric ammonium sulfate). The glass tubes were properly sealed, mixed, and placed in boiling water for 50 min. Duplicate samples was prepared when extract amounts allowed. After cooling, the light absorption at 550 nm was determined using a spectrophotometer (UV-1800; Shimadzu Corp., Kyoto, Japan). The average between duplicate samples was used as one data value. Purified tannins from spruce needles were used as standards to calculate concentrations.

#### Data Analysis

For each compound, we used Split plot ANOVAs to test for the effect of treatment as main plot factor and age as sub-plot factor. The composition of phenolic compounds in the needles were visualized with the two first axes of a principle component analysis (PCA) using the R package vegan (Oksanen et al., 2016). We used the ordiellipse function (Oksanen et al., 2016) to plot the 95% confidence intervals (CI; based on SE-values) of the respective age × treatment centroid. All analyses were performed using R 3.2.5 (R Core Team, 2018).

#### RESULTS

#### Total N and C Concentrations

Fertilization increased N% in both current and previous year needles, and correspondingly decreased the C:N ratio (**Figure 1**). At the same time, N% was higher, and C: N ratio lower in current year needles compared with those from the previous year (**Figure 1**).

#### Phenolic Concentrations

Control needles from the previous year contained over three times more total low molecular phenolics than did current year needles (**Table 1**), and the composition differed greatly (**Figure 2**, **Table 1**). In previous year needles (controls), stilbenes constituted 33% of the total low molecular phenolics concentration, and the most abundant compounds were piceatannol glucoside and resveratrol glucoside. In comparison, stilbenes represented only 0.03% in current year needles. On the other hand, current year needles had more than double concentration of flavonols compared with the older needles, of which 3,6-dicoumaroyl astragalin was the highest compound in both needle types (**Table 1**). The other flavonols were present in minor concentrations, mostly differing significantly between the two needle cohorts, but in both directions. The same was true for most other flavonoids, but (+) catechin and gallocatechin had more than 2 and 10 times higher concentration in previous year needles, respectively. We identified picein only in current year needles, while the older needles had more than 30 times more its aglycone, 4 hydroxy acetophenone (**Table 1**). In addition, the needle cohorts contained two different unknown lignans and the same three procyanidins as well as B<sup>3</sup> in minor concentrations. Control needles from the previous year contained about the same

was sqrt-transformed.

amount of the two fractions of condensed tannins, while the current year ones had higher concentration of the insoluble fraction.

The composition of phenolic compounds changed upon fertilization in both the current and previous year's needles (**Figure 2**), but this was statistically significant only for the new needles. However, the pattern differed between the two age classes. The total concentration of stilbenes increased upon fertilization in previous year needles, while it decreased in current year needles. These patterns were seen for all three individual stilbenes that were present in both needles cohorts (significant A × T interactions, **Table 1**). Two of three (resveratrol aglycon and methylpiceatannol glucoside) stilbenes that were only present in the older needles increased upon fertilization. Of the flavonols, quercetin 3 galactoside, quercetin 3-glucoside, kaempferol 3-glucoside, and isorhamnetin decreased upon fertilization in both needle cohorts. Monocoumaroylatragalin and the 3,6-dicoumaroylastragalin,

#### TABLE 1 | Concentrations (mg g−<sup>1</sup> DW) (mean values ±1 SE) of and Split plot ANOVAs on phenolic compounds.


The effect of treatment is the main plot factor and age is the sub-plot factor. Numbering of compounds is corresponding to numbering in the PCA *(Figure 2)*. † log-transformed, ‡

sqrt-transformed. Numbers in bold indicate statistical significant results. Lines in italics indicate sums of chemical groups.

on the other hand, increased in current year needles and decreased in those from the previous year (significant A × T interactions, **Table 1**). The flavonoids apigenin 7-glucoside, apigenin aglycon and (+)-catechin, as well as two of the procyanidins decreased in fertilized plots in both needle types.

The soluble fraction of condensed tannins changed strongly upon fertilization in both needle cohorts, decreasing in current

year needles, while it increased in the older ones. The insoluble fraction, on the other hand, decreased slightly in the youngest needles, while it increased in those from the previous year (significant A × T interactions, **Table 1**).

#### DISCUSSION

Changes in the availability of a limiting resource such as nitrogen may have large impacts on ecosystem composition and functioning. In the present study, we showed that heavy fertilization increased the N concentrations of spruce needles, and induced changes in both concentrations and compositions of phenolic compounds.

The most striking result is the differences between the two needle cohorts. The current year needles reacted to fertilization in accordance with our first hypothesis and ecological theories on balance between C and N (e.g., Bryant et al., 1983; Jones and Hartley, 1999), as concentrations of many phenolic compounds were lower in needles from fertilized plots, compared with controls (**Table 1**). The difference was strongest for the soluble fraction of condensed tannins, which were reduced to less than half the concentration of controls. This decrease corresponds well with the decrease in (+)-catechin, the precursor of condensed tannins. In addition, many flavonoids had lowest concentrations in fertilized plots, as well as the total concentration of low molecular phenolics. The coumaroyl astragalins, however, was not affected by fertilization, suggesting high priority of these compounds in the young needles. Blodgett et al. (2005) and Edenius et al. (2012) also found reductions in total phenolics in fertilized P. resinosa and P.abies. Virjamo et al. (2014), on the other hand, found little effects of fertilization in needles of 1-year old P.abies seedlings, except that the concentration of (+)-catechin increased, also contrasting our results.

As expected, the previous year needles were less affected, but interestingly some individual stilbenes (and stilbenes in total), as well as both fractions of condensed tannins, were present in highest concentrations in needles from fertilized trees. It may be that increased N availability increased chlorophyll content and thus photosynthesis and C availability of fertilized trees. According to the carbon-nutrient balance hypothesis and the PCM (Bryant et al., 1983; Jones and Hartley, 1999, respectively), such an improved resource situation would be expected to give increased growth and not defense. Since defense was reduced in current year needles, this also challenges the established truth that new needles have priority when resources are distributed (resources are transported from old parts to young parts). Or, alternatively, it suggests that the defensive compounds of old needles cannot be transported, or dissolved/turned over and then transported, to younger needles.

As we hypothesized, previous year needles from control plots had higher concentrations of PSMs than the new needles. In addition, the composition of compounds differed strongly. The large picture was that the older needles accumulated large amounts of stilbenes and acetophenone, in addition to flavonoids, while the defensive system of new needles mostly consisted of flavonoids (**Table 1**). The levels of condensed tannins, however, did not differ much with age in needles from unfertilized trees. These results fits well with studies of seasonal changes in phenolics of P.abies. Both Slimestad (1998) and Ganthaler et al. (2017a) noted that flavonoids (and especially flavonols) were predominant in the first weeks after bud burst, after which stilbenes, picein and shikimic acid increased. Increases in stilbene concentrations toward late summer and autumn was also seen in some older studies (Kaufmann et al., 1974; Bjørnseth, 1977; Solhaug, 1990). Our results from previous year needles (sampled 1 year after they sprouted) suggest that the high protection levels achieved last autumn are kept over time. Flavonoids, and especially flavonols, are often found in the vacuoles of epidermal cells, where they play a role in light protection (Close and McArthur, 2002). Quercetins and kaempferols are often induced in response to enhanced levels of ultraviolet light (e.g., Ryan et al., 1998; Nybakken et al., 2012). In an experiment excluding ultraviolet (UVB) light from young scots pine plants (Pinus sylvestris), Turunen et al. (1999) showed that non-acetylated flavonol 3-glycosides were reduced, while the di-acetylated flavonols (dicoumaroyl astragalins) were unaffected. Earlier indoor studies, however, showed that also di-acetylated flavonols were induced by UVB radiation (Jungblut et al., 1995). It may be that light protection must be prioritized during early needle development in spring, when the epidermis is still thin and physical defenses like waxes are not so prominent. The high levels of dicoumaroyl astragalins, and their high resilience against change shown both here and by Turunen et al. (1999), suggest they are especially important in young conifer needles. Flavonoids may also play other defense roles, e.g., against herbivores (e.g., Bryant et al., 1983), and may be regarded as a general defense. Interestingly, the newly developed needles also contained high levels of the more complex flavonoids condensed tannins, which are also important herbivore defenses (Barbehenn and Constabel, 2011).

Stilbene levels, on the other hand, were very low in young needles. Both stilbenes and flavonoids are produced through the phenylpropanoid pathway from the precursor phenylalanine, but in the next step the enzymes stilbene synthase (STS) and chalcone synthase (CHS), respectively, lead to stilbene and flavonoid biosynthesis from the same cinnamoyl-CoA/p-coumaroyl-CoA precursor (Kodan et al., 2002). Based on this, we have no indications that stilbenes are more complex or costly to produce for the young needles, although this is scarcely studied. However, it may be that stilbenes serve a more specialized role, and are therefore less prioritized during early spring when resources (both C and N) are scarce. Stilbenes are located in the central part of spruce needles (Solhaug, 1990), as well as in the xylem of conifers (e.g. Harju et al., 2009), and should thus have other roles than light protection. Stilbenes are indeed expected to play roles in fungal protection, as they are frequently induced by pathogen attack (Chong et al., 2009; Jeandet et al., 2010; Ganthaler et al., 2017a), and varying concentrations have been associated with intraspecific variation of host plant susceptibility to infection (Lieutier et al., 2003). In addition to pathogens, the previous year needles most likely contained a higher amount of fungal endophytes, as they have been exposed to environmental inoculum for a longer time period. Little is known about how these microbes affect the biochemistry of plants, but they could potentially both consume, process or produce phenolics themselves (e.g., Hardoim et al., 2015; Yang et al., 2016).

To sum up, fertilization changed phenolic profiles of both current and previous year needles, but in different ways. Older needles partly increased their protection levels, both regarding what is thought to be specific fungal defense (stilbenes) as well as specific herbivore defense (condensed tannins). The new needles, on the other hand, showed reductions in the already low stilbene levels, as well as in condensed tannins and flavonoids. This may have implications for the susceptibility of spruce against attacks from fungal pathogens. Typically, needle diseases like rust attack the new needles, and a connection between phenolic levels and susceptibility to infection was found for spruce attacked by the needle bladder rust (Chrysomyxa rhododendri) (Ganthaler et al., 2017b). In sub-alpine spruce forests in central Europe, repeated infections by rust fungi has led to reduced timber yields and severe problems with regeneration of the species (Oberhuber et al., 1999; Ganthaler et al., 2014). Repeated rust (C. abietis) attacks on spruce are also observed in Scandinavia the latest years, but to our knowledge the impact of these are not yet studied. A changing climate with milder winters, as well as increased trade of living plants (Liebhold et al., 2012) may help pathogens to increase in abundance and outbreak frequencies, but also to reach new areas. The same is expected for insect pests. In a review, Throop and Lerdau (2004) concluded that N deposition has positive effects on general insect performance, probably though the effect it has on plants: increased N content and decreased PSMs. A study along latitudinal gradients in Europe showed that the load of sapfeeding insects on forest trees are driven by climatic patterns, indicating that also climate warming will drive plant losses to insect feeding (Kozlov et al., 2005). In conclusion, this may imply that nitrogen deposition and forest fertilization leave forests less robust in a time when pests may take advantages of a changing climate.

#### AUTHOR CONTRIBUTIONS

LN planned and supervised the chemical analyses, interpreted the data and wrote the MS. ML designed the experiment, supervised the sampling, and contributed to writing the MS. JA conducted the statistical analyses and contributed to writing the MS. RJ-T supervised the UHPLC-MS analyses and the interpretation of data, and contributed to writing the MS. MO designed the experiment, performed the fertilizations, and contributed to writing the MS.

#### REFERENCES


#### ACKNOWLEDGMENTS

We thank Mathilde Lorentzen for help with sampling the needles, as well as Annie Aasen and Claus Kreibich for performing the chemical analyses.

#### SUPPLEMENTARY MATERIAL

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


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

Copyright © 2018 Nybakken, Lie, Julkunen-Tiitto, Asplund and Ohlson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Flavanone-3-Hydroxylase Plays an Important Role in the Biosynthesis of Spruce Phenolic Defenses Against Bark Beetles and Their Fungal Associates

Almuth Hammerbacher<sup>1</sup> \*, Dineshkumar Kandasamy<sup>2</sup> , Chhana Ullah<sup>2</sup> , Axel Schmidt<sup>2</sup> , Louwrance P. Wright<sup>3</sup> and Jonathan Gershenzon<sup>2</sup>

<sup>1</sup> Department of Zoology and Entomology, Forestry and Agricultural Biotechnology Institute, University of Pretoria, Pretoria, South Africa, <sup>2</sup> Department of Biochemistry, Max Planck Institute for Chemical Ecology, Jena, Germany, <sup>3</sup> Zeiselhof Research Farm, Pretoria, South Africa

#### Edited by:

Bartosz Adamczyk, Natural Resources Institute Finland (Luke), Finland

#### Reviewed by:

Elizabeth H. J. Neilson, University of Copenhagen, Denmark Françoise Martz, Natural Resources Institute Finland (Luke), Finland

\*Correspondence:

Almuth Hammerbacher almuth.hammerbacher@fabi.up.ac.za

#### Specialty section:

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

Received: 03 September 2018 Accepted: 07 February 2019 Published: 25 February 2019

#### Citation:

Hammerbacher A, Kandasamy D, Ullah C, Schmidt A, Wright LP and Gershenzon J (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. doi: 10.3389/fpls.2019.00208 Conifer forests worldwide are becoming increasingly vulnerable to attacks by bark beetles and their fungal associates due to the effects of global warming. Attack by the bark beetle Ips typographus and the blue-stain fungus it vectors (Endoconidiophora polonica) on Norway spruce (Picea abies) is well known to induce increased production of terpene oleoresin and polyphenolic compounds. However, it is not clear whether specific compounds are important in resisting attack. In this study, we observed a significant increase in dihydroflavonol and flavan-3-ol content after inoculating Norway spruce with the bark beetle vectored fungus. A bioassay revealed that the dihydroflavonol taxifolin and the flavan-3-ol catechin negatively affected both I. typographus and E. polonica. The biosynthesis of flavan-3-ols is well studied in Norway spruce, but little is known about dihydroflavonol formation in this species. A flavanone-3-hydroxylase (F3H) was identified that catalyzed the conversion of eriodictyol to taxifolin and was highly expressed after E. polonica infection. Down-regulating F3H gene expression by RNA interference in transgenic Norway spruce resulted in significantly lower levels of both dihydroflavonols and flavan-3-ols. Therefore F3H plays a key role in the biosynthesis of defense compounds in Norway spruce that act against the bark beetle-fungus complex. This enzyme forms a defensive product, taxifolin, which is also a metabolic precursor of another defensive product, catechin, which in turn synergizes the toxicity of taxifolin to the bark beetle associated fungus.

Keywords: taxifolin, catechin, flavonoids, F3H, Picea abies, Endoconidiophora polonica, Ips typographus

# INTRODUCTION

Flavonoids are specialized or secondary plant metabolites that are universally distributed in all vascular plants and occur in many algae as well (Treutter, 2006; Mouradov and Spangenberg, 2014). Containing a C6–C3–C<sup>6</sup> backbone of two aromatic rings connected by a C<sup>3</sup> bridge, they are known to perform numerous internal functions in plants as well as in interactions with the environment.

**56**

Flavonoids influence the central processes of plant development and nutrition by regulating auxin transport (Kuhn et al., 2011), pollen tube growth (Ylstra et al., 1992) and nitrogen fixation (Ferreyra et al., 2012). They are also involved in negative or positive plant interactions with other organisms. They, for example, inhibit the growth of disease-causing microorganisms (Goetz et al., 1999; Ullah et al., 2017) and deter herbivore feeding (Dübeler et al., 1997; Ohse et al., 2017). On the other hand, these compounds play an important role in attracting pollinators to flowers and dispersal agents to fruits. Flavonoids are well known anti-oxidants that increase plant tolerance to soil salinity (Mahajan and Yadav, 2014), cold stress (Meng et al., 2015) as well as to drought stress (Song et al., 2016). They protect against UV radiation by anti-oxidant effects and direct screening (Ferreyra et al., 2012). The anti-oxidant properties of these molecules are also valuable in human nutrition and contribute to lower incidences of human diseases related to high dietary fat intake and stress (Weidmann, 2012).

Flavonoids are produced by sequential condensation of a phenylpropanoid acid coenzyme A (CoA) ester (C6–C3) to three malonyl CoA molecules to form chalcone with its typical C6–C3–C<sup>6</sup> flavonoid backbone (Treutter, 2006). In the early steps of the flavonoid pathway, chalcone is modified via isomerization and oxidation to a dihydroflavonol (Winkel-Shirley, 2001). In most plants, dihydroflavonols are a major branch point in the flavonoid biosynthesis pathway, and subsequent reactions can either produce dihydroflavonol glucosides via glycosylation, flavonols or anthocyanidins via oxidation and flavan-3-ols via reduction (**Figure 1**).

The biosynthesis of dihydroflavonols was identified as a key regulatory point for the formation of down-stream metabolites in the flavonoid pathway and is catalyzed by flavanone 3β-hydroxylase (F3H; EC 1.4.11.9; Forkmann and Stotz, 1981). F3H has been isolated and characterized from more than 50 plant species (Meng et al., 2015; Han et al., 2017) and is a 2-oxoglutarate-dependent dioxygenase that catalyzes the 3β-hydroxylation of 2S-flavanones to 2R,3Rdihydroflavonols using 2-oxoglutarate, dioxygen (O2), ferrous iron, and ascorbate as substrates and co-factors (Britsch et al., 1992). Amino acid sequences of these enzymes are conserved, with well-defined 2-oxoglutarate and iron binding sites (Britsch and Grisebach, 1986; Britsch et al., 1992; Turnbull et al., 2004). The enzyme is localized in the cytoplasm (Tu et al., 2016) and may form multi-enzyme complexes with other enzymes involved in flavonoid biosynthesis (Burbulis and Winkel-Shirley, 1999). These multi-enzyme complexes allow for co-regulation of multiple enzymatic steps in the pathway. For example, F3H is transcriptionally regulated together with up-stream genes in Arabidopsis thaliana [chalcone synthase (CHS) and chalcone isomerase (CHI)] (**Figure 1**) (Pelletier and Shirley, 1996; Burbulis and Winkel-Shirley, 1999). However, in many other plant species, this gene is co-regulated with down-stream genes, including dihydroflavonol reductase (DFR), anthocyanidin synthase (ANS), and leucoanthocyanidin reductase (LAR) (Beritognolo et al., 2002; Singh et al., 2008; Liu et al., 2013; Song et al., 2016; Ullah et al., 2018).

Flavanone-3-hydroxylase gene expression is positively regulated by the plant hormones salicylic acid (Sun et al., 2016), jasmonic acid (Meng et al., 2015), and abscisic acid (Song et al., 2016). Increased transcription of the F3H gene leads to tolerance to abiotic as well as biotic stress, including drought (Watkinson et al., 2006), saline conditions (Mahajan and Yadav, 2014), cold (Meng et al., 2015), UV radiation (Liu et al., 2013) as well as biotrophic and necrotrophic attackers (Sun et al., 2016). High expression of this enzyme results in accumulation of proanthocyanidins (PAs), which are major end products of the pathway (Beritognolo et al., 2002; Song et al., 2016), higher levels of antioxidants (Meng et al., 2015) and lower levels of reactive oxygen species (Mahajan and Yadav, 2014).

Tree species within the Pinaceae, most notably spruce (Picea spp.) and pine (Pinus spp.) are economically important keystone species that dominate temperate, boreal, and montane landscapes. These long-lived woody perennials are very susceptible to the effects of climate change (Hanewinkel et al., 2013). Warmer weather conditions, wind storms and unseasonal frost have recently resulted in a world-wide decline of spruce and pine forests (Allen et al., 2010; Bentz et al., 2010). A main driver of these declines are bark beetles, which initially attack stressed and wind-damaged trees. This results in beetles building up massive population sizes and switch from an endemic population state to an epidemic phase. Bark beetles in the epidemic phase attack healthy trees by pheromone-driven mass attacks and disperse rapidly over wide areas resulting in the loss of millions of hectares of forest per year (Boone et al., 2011). Bark beetle success in overcoming the resistance mechanisms of healthy host trees has been partially ascribed to simultaneous attacks by bark beetle-associated fungi (Krokene, 2015), which are thought to exhaust tree defenses, although this view is not shared by others (Six and Wingfield, 2011).

In an effort to preserve Pinaceae forests in areas most affected by global warming, research is being conducted to identify resistance traits against bark beetles and their associated fungi (Keeling and Bohlmann, 2006; Hamberger et al., 2011; Krokene, 2015). These studies focused on understanding the biosynthesis of terpenoid oleoresin. Resins entrap and intoxicate attacking beetles and inhibit the growth of their fungal associates (Keeling and Bohlmann, 2006; Schiebe et al., 2012). The Pinaceae also produce high concentrations of polyphenols, such as stilbenes and PAs (Raiber et al., 1995; Booker et al., 1996; Li et al., 2012; Hammerbacher et al., 2011, 2013). Recent studies suggested that PAs (also known as condensed tannins) appear to function in tree defense against bark beetle-fungus invasions (Hammerbacher et al., 2014, 2018). However, little is known about the role of other flavonoids in the defense of spruce and pine against bark beetles and their associated fungi, although circumstantial evidence suggests that they should also play an important role (Brignolas et al., 1995; Li et al., 2012).

We therefore investigated the biosynthesis and defensive role of flavonoids in the Pinaceae using a study system of Norway spruce (Picea abies), its most important bark beetle pest (Ips typographus) and the bark beetle vectored fungus, Endoconidiophora polonica, an ascomycete that stains the wood with a blue color. We artificially simulated bark beetle attack by

wounding trees and inoculating these wounds with E. polonica. Among the flavonoids, dihydroflavanol and flavan-3-ol content increased most significantly in the cambium of fungus-infected trees compared to wounded control trees. Bioassays of these compounds revealed that both classes were toxic to bark beetles as well as their fungal associates. Flavan-3-ol biosynthesis in Norway spruce was elucidated recently (Hammerbacher et al., 2014, 2018), but little is known about dihydroflavonol biosynthesis in this tree species. We identified a single gene encoding a functional F3H enzyme that produces dihydroflavonols in Norway spruce, which is transcriptionally up-regulated during simulated bark beetle attack. Suppressing expression of this gene by RNA interference (RNAi) in transgenic Norway spruce saplings resulted in lower concentrations of dihydroflavonols and flavan-3-ols, revealing its important role in regulating down-stream flavonoid accumulation in this species.

#### MATERIALS AND METHODS

#### Simulated Bark Beetle Attack: Inoculation of Norway Spruce With Bark Beetle-Associated Fungus or Sterile Agar

Seven-year old Norway spruce saplings from the clone S21K0420232, purchased from Skogforsk (Sweden), were grown in 5 L pots in an outdoor plot at the MPI-CE, Jena, Germany as described previously (Hammerbacher et al., 2018). The saplings had a stem diameter of 2–2.5 cm. E. polonica isolate K2014 (Kandasamy, isol. Ex. bark beetle gallery, Thuringian Forest, Gotha, Germany), was grown on potato dextrose agar (Difco, BD, Franklin Lakes, NJ, United States) at 25◦C for 7 days in darkness. In June 2015, a 5 mm diameter cork borer was used to remove a circular piece of bark from the lower stems (approximately 80 cm above ground level) of the spruce saplings. A similar-sized disk of the fungal culture was inserted into the wound and sealed with Parafilm. Sterile disks of potato dextrose agar were used for control treatments. In total 20 saplings were inoculated and 20 were wounded. Actively growing cambium from five trees (n = 5) from inoculated and wounded saplings were harvested after 2, 7, 14, and 28 days post inoculation (dpi) and flash frozen in liquid nitrogen. Sections from 2.5 cm above to 2.5 cm below the inoculation point were harvested from all treatments at 2 and 7 dpi as well as from the sterile agar-inoculated treatments at 14 and 28 dpi. The fungus-inoculated lesions from the 14 and 28 dpi treatments were separated into two samples, comprising (1) a section from 2.5 cm above to 2.5 cm below the point of inoculation (inner lesion) and (2) sections from 2.5 to 4 cm both above and below the point of inoculation (outer lesion).

#### Flavonoid Analysis

Samples from fungus-inoculated treatments and sterile agar-inoculated controls as well as stems of transgenic spruce carrying the F3H RNAi construct were finely ground in liquid

nitrogen using a mortar and pestle. A subsample of the resulting wood powder was lyophilized at 0.34 mbar pressure using an Alpha 1-4 LD plus freeze dryer (Martin Christ GmbH, Osterode, Germany). Approximately 20 mg of dried spruce tissue powder was extracted twice for 4 h with 800 µl analytical grade methanol containing 10 µg ml−<sup>1</sup> internal standard, apigenin-7-glucoside (Carl Roth GmbH, Karlsruhe, Germany). Flavonoids were analyzed 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) using the protocols described in Hammerbacher et al. (2014) for the following analytes: quercetin-glucoside, catechin, gallocatechin, procyanidin B1, and catechin-gallocatechin dimers. Dihydromyricetin was measured using the protocols from Hammerbacher et al. (2018). Additional multiple reaction monitoring (MRM) was used in this study to measure analyte precursor ion → fragment ions for additional flavonoid metabolites as follows: m/z 271.0 → 151.0 [collision energy (CE), −28 V; declustering potential (DP), −55 V] for naringenin; m/z 287.0 → 151.0 (CE, −20 V; DP, −75 V) for eriodictyol; m/z 303.0 → 125.0 (CE, −28 V; DP, −40 V) for taxifolin; m/z 465.0 → 285.0 (CE, −30 V; DP, −70 V) for taxifolin glucoside. To develop the MRMs, standards were purchased from Sigma-Aldrich (naringenin, eriodictyol, and taxifolin), or from TransMIT GmbH Giessen, Germany (taxifolin-glucoside, dihydromyricetin). Products of in vivo conversion of naringenin and eriodictyol by PaF3H in Escherichia coli were analyzed using a Bruker Daltronics ion trap mass spectrometer following the protocols of Hammerbacher et al. (2014). External calibration curves from 1 µg ml−<sup>1</sup> to 1 mg ml−<sup>1</sup> were used for quantification of all metabolites reported in this study, except for taxifolin glucoside for which a calibration curve of quercetin rhamnoside (TransMIT) was used with a response factor of 1.

# Identification of PaF3H Candidate Sequences and Phylogenetic Analysis

Protein sequences of previously characterized F3H genes were downloaded from NCBI and used in BLASTP searches against all Picea transcript assemblies in the conifer genome database<sup>1</sup> and on NCBI. Transcripts with high identity to known F3H sequences from P. abies, P. glauca, and P. sitchensis were assembled into full-length contigs using CLC Workbench (Qiagen, Aarhus, Denmark). Primers were designed for the 5 0 and 3<sup>0</sup> ends of the full length transcript of P. abies (F3HF GGGGACAAGTTTGTACAAAAAAGCAGGCTCAATGGCGC CCGCCGCAGTCGTGG and F3HR GGGGACCACTTT GTACAAGAAAGCTGGGTATCACGACTTTTCCTGCTCAGC AAC). RNA was extracted from approximately 50 mg of fresh, frozen (−80◦C) spruce tissue as described in Hammerbacher et al. (2018) using the Stratec Plant RNA mini kit (Stratec, Birkenfeld, Germany). RNA quality was assessed spectrophotometrically and by agarose gel electrophoresis. RNA samples with a concentration higher than 100 ng µl −1 and an A260/<sup>280</sup> higher than 1.75 was used for reverse transcription. Approximately 500 ng total RNA was reverse-transcribed to cDNA using Superscript v. II (Invitrogen) following the manufacturer's protocols. Full length PaF3H was then amplified from cDNA using Phusion Taq (New England Biolabs, Ipswich, MA, United States) and cloned into the Gateway compatible vector pDONR 207 (Invitrogen, Carlsbad, CA, United States) using BP clonase II (Invitrogen). To verify the PaF3H sequence, the entry clones were sequenced using pDONR F and R primers (Invitrogen) using Sanger chain termination sequencing and the BigDye v3.1 cycle sequencing kit (Thermo Fisher Scientific) and an ABI prism 7000 capillary electrophoresis sequencing instrument (Thermo Fisher Scientific). Protein sequences of characterized F3H enzymes from different angiosperm species were aligned with the translated nucleotide sequences of F3H from different gymnosperm species (**Supplementary Figure S1**) using the MAFFT web interface<sup>2</sup> and imported into MEGA v.6 (Kumar et al., 2008; Center for Evolutionary Medicine and Informatics, Tempe, AZ, United States). The evolutionary history was inferred by using the Maximum Likelihood method based on the JTT matrix-based model (Jones et al., 1992). The tree with the highest log likelihood (−10349.13) was constructed (500 boot strap replicates). The percentage of trees in which the associated taxa clustered together is shown next to the branches. Bootstrap values lower than 50% are not shown. The initial tree for the heuristic search was obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then selecting the topology with a superior log likelihood value. The tree is drawn to scale, with branch lengths measured as the number of substitutions per site. The analysis involved 22 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 262 positions in the final dataset (GenBank accession numbers for all sequences are included in **Supplementary Table S1**).

#### qPCR

RNAs were extracted and reverse transcribed as described above from the same samples from which the metabolites were extracted. Approximately 100 ng cDNA were used per PCR reaction. qPCR was performed with the Kappa SYBR Fast qPCR kit (Kappa Biosystems, Boston, MA, United States) and a CFX96 – Real-Time System (Bio-Rad, Hercules, CA, United States) following the manufacturer's specifications. PaF3H was amplified using PaF3H-F (GCAGAGCGTGCACAG) and PaF3H-R (GTGAGTTGAGTTCTGTGGAG) primers with an initial denaturation step at 95◦C for 2 min. followed by 40 cycles of 95◦C denaturation and 60◦C extension. The PCR was normalized using PaUBI (Schmidt et al., 2010; GB:EF681766.1) and calibrated to the 2 days sterile agar-inoculation treatment. The means of three technical replicates from each of five biological replicates per treatment group were used to calculate manually the relative transcript abundance according to Pfaffl (2001). A dilution series was made to determine the primer efficiencies, which were close to 100%. Non-template controls and non-reverse transcribed RNA controls were included. Each

<sup>1</sup>http://congenie.org/blast

<sup>2</sup>http://mafft.cbrc.jp/alignment/server/

PCR product was analyzed using a melting curve to verify single product amplification.

#### Functional Characterization of PaF3H by Expression in E. coli and RNAi in Transgenic Norway Spruce

The in vivo enzyme activity of PaF3H was assessed by transforming E. coli BL21 with PaF3H cloned into the destination vector pDEST15 (Invitrogen) using LR clonase (Invitrogen) following the manufacturer's protocols. BL21 (DE3) bacterial cells were grown on Luria-Bertani (LB) plates containing 100 µg ml−<sup>1</sup> ampicillin. For protein expression single colonies were inoculated into 5 ml LB broth with 100 µg ml−<sup>1</sup> ampicillin, and grown for 12 h at 30◦C. The 5 ml starter cultures were used to inoculate 100 ml LB medium supplemented with 100 µg ml−<sup>1</sup> ampicillin. Bacterial cultures were grown for 9 h at 18◦C (220 rpm) and protein expression was induced with 0.5 mM IPTG (isopropyl β-D-1-thiogalactopyranoside). Two hours after inducing expression with IPTG the cultures were split and naringenin or eriodictyol in DMSO was added to the culture medium to a final concentration of 200 µg ml−<sup>1</sup> . Cultures were harvested 12 h after addition of the substrates. Bacteria were removed from the culture medium by centrifugation. The medium was acidified with 1% (v v−<sup>1</sup> ) of 0.1 N HCl and extracted with three volumes of ethyl acetate. The ethyl acetate extracts were evaporated using a rotary evaporator and re-dissolved in 500 µl methanol for LC-MS analysis following the protocols described above.

In order to assess the in planta activity of PaF3H, transgenic saplings were made carrying an F3H RNAi construct. For that a 330-bp region of the PaF3H gene was PCR-amplified by using the oligonucleotides PaF3HRNAi-for (tgctctagagcagggtttggcatgcaatatgttggcg) and PaF3HRNAi-rev (cgggatcccgggtccgcgctcttgaactt) and cloned in sense and antisense orientations into the multiple cloning sites of the pTRAIN vector on either side of an intron as described by Levée et al. (2009). After HindIII digestion, the excised RNAi-cassette including also an upstream maize ubiquitin promotor was ligated into the multiple cloning site of the pCAMBIA 1305.2 vector<sup>3</sup> . Agrobacterium tumefaciens-mediated stable transformation of P. abies embryogenic tissue (Pa186/3c) was performed as described in detail by Schmidt et al. (2010). Generation of somatic transgenic seedlings was based on a protocol originally reported for white spruce from Klimaszewska et al. (2005). Seedlings of three independent kanamycin-resistant transgenic lines plus an empty vector control line were characterized by qPCR, using four plants per line. Transgenic seedlings were grown for 2 years before the stem tissue was harvested, ground to a fine powder under liquid nitrogen and analyzed using both qPCR and LC-MS (**Supplementary Figure S2**).

# Fungal Growth on Catechin and Taxifolin

Growth of E. polonica on potato dextrose medium containing either DMSO without flavonoids (control), 2 mg ml−<sup>1</sup> catechin, 2 mg ml−<sup>1</sup> taxifolin, or 2 mg ml−<sup>1</sup> of a physiologically relevant 2:1 mixture of catechin and taxifolin was determined in Petri dishes using the protocols described by Hammerbacher et al. (2014). In brief, the compounds were added to the medium in DMSO (50 mg ml−<sup>1</sup> ) just before the medium was dispersed into Petri dishes. Fungal growth of four replicates per compound was measured daily until the fungus grew to the edges of the Petri dish. Fungal growth was plotted for each Petri dish and a linear curve was fitted. The growth rate was determined by the gradient of the fitted curve.

#### Bark Beetle Feeding Assays

Adult bark beetle feeding assays were performed in 12 cm long pneumatic tubes, outer diameter 6 mm, inner diameter 4 mm (Sang-A Pneumatic, Daegu, Korea). Each assay tube was sealed at one end with Parafilm and filled with spruce bark diet prepared as follows: 7% (w/v) finely milled spruce inner bark powder was mixed with 1% fibrous cellulose (Sigma), 4% Bactoagar (Roth) in water and autoclaved for 20 min at 121◦C. In order to facilitate uniform packing of the diet in each tube, the diet mix was vortexed thoroughly and drawn up in a disposable sterile syringe. The end of the syringe was fitted onto the open end of the assay tubes to discharge the diet into them. The open ends of the assay tubes were then sealed and the medium was allowed to dry overnight at room temperature. Approximately 2 cm of diet was manually removed from one end of the assay tube to facilitate insertion of bark beetles. Stock solutions of catechin and taxifolin were prepared in DMSO and added to the spruce bark agar before filling the mixtures into the tubes. Only DMSO was added in the control diet without amended compounds. The final concentrations of the compounds in the medium were 1 mg ml−<sup>1</sup> catechin, 1 mg ml−<sup>1</sup> taxifolin, or 0.5 mg ml−<sup>1</sup> of both catechin and taxifolin (1:1). Adult beetles used in this assay were reared in the laboratory using the protocols of Anderbrant et al. (1985). The initial insects for the laboratory colony were collected from naturally infested trees in Jena, Germany. The weight of each beetle was measured to the nearest 0.01 mg, and it was placed carefully in the assay tube with its head facing towards the diet. The length of the feeding tunnel in each tube was measured after 6 and 48 h and the weight of the beetle was then measured after 48 h to determine the weight change resulting from different treatments. Previous studies have reported 6 h to be sufficient to see the difference in tunneling length and 48 h for weight change in adult bark beetles (Faccoli and Schlyter, 2007, Cardoza et al., 2008). Each treatment (control, catechin, taxifolin, or the mixture of both compounds) was replicated 15 times. The assay beetles included both males and females.

#### Statistical Analysis

All data were subjected to a Shapiro–Wilk test for normality and non-normal data were transformed with the log or (1 + log) function. Transformed or non-transformed data were analyzed using two- or one-way ANOVAs or repeated measures ANOVA (specifically for the tree inoculation trial). Following this, Tukey's post hoc test was used on data showing significant differences using a 95% confidence interval. For bark beetle feeding

<sup>3</sup>www.cambia.org

assays, a paired sample t-test was used to compare the mean weight of beetles before and after feeding. All analyses were conducted using R<sup>4</sup> .

#### RESULTS

#### Flavonoids Accumulate in Norway Spruce Upon Infection With the Bark Beetle–Associated Fungus, E. polonica

Norway spruce saplings were wounded and inoculated with the bark beetle-associated, sapwood-staining fungus E. polonica or inoculated with sterile agar as a control. Samples were harvested from both infected and control trees for analysis over a 28 days time course. Fungal lesions at 14 and 28 days post-inoculation (dpi) were separated into two sections: The inner lesion section contained fungus that established itself at the onset of the experiment and the outer lesion section was made up of newly infected tissue. Control samples equivalent to the inner lesion region were harvested at the same time points. Flavanones (naringenin and eriodictyol), dihydroflavonols (taxifolin, taxifolin-glucoside, and dihydromyricetin), flavonols (quercetin, myricetin, laricitrin, and their glycosides) as well as flavan-3-ols (catechin, gallocatechin, and dimeric proanthcyanidins) were analyzed in the cambial tissues using high-performance liquid chromatography coupled to mass spectrometry (LC-MS). Among these naringenin, eriodictyol, taxifolin, taxifolin glucoside, catechin, and procyanidin B1 could be detected above the quantification limit (**Supplementary Figure S3** and **Supplementary Table S2**).

The flavanone naringenin, as well as its 3<sup>0</sup> -hydroxylated derivative eriodictyol (**Figure 1**) accumulated to significantly higher concentrations in the inner lesions compared to the wounded control (p < 0.001) at 7 dpi and reached a maximum at 14 dpi after which the content of both flavanones remained constant (**Figures 2A,B**). In the newly infected tissue of the outer lesion, the concentration of naringenin and eriodictyol only increased significantly above control concentrations at 28 dpi (p < 0.001). A slight increase in the concentrations of both flavanones was also observed in the controls inoculated with sterile agar over the time course of 28 days, which was significant between 7 and 14 dpi (p < 0.01).

The dihydroflavonol taxifolin (2R,3R-dihydroquercetin) and its glycosylated derivative, taxifolin-3-glucoside were detected at much higher concentrations than the flavanones. The glucoside accumulated to levels as high as 5 mg g−<sup>1</sup> dry weight after 28 dpi in fungus-infected tissue (**Figures 2C,D**). The aglucone accumulated more rapidly in fungus-inoculated tissue than the glucoside and reached significantly higher levels than in the controls as early as 7 dpi (p < 0.01). Taxifolin reached a maximum in both inner as well as outer lesions at 14 dpi and remained constant until 28 dpi. Compared to the aglucone, the glucoside accumulated at a slower rate, and reached a maximum at 28 dpi, which was significantly higher than the levels in the control (p < 0.05) in both the inner as well as in the expanding outer lesion. An increase of taxifolin-3-glucoside was also recorded in the controls inoculated with sterile agar over the time course of this experiment, where treatments at 2 dpi differed significantly from treatments at 28 dpi (p < 0.05).

The last two flavonoids detected in our analysis were both flavan-3-ols that are precursors for the biosynthesis of condensed tannins (**Figures 2E,F**). Monomeric catechin could be detected at higher levels than its dimerized derivative procyanidin B1 (PAB1) in both uninfected and infected spruce stems. In the inner lesion, catechin content increased significantly above the control by 7 dpi (p < 0.001) reaching a maximum at 14 dpi, similar to the flavanones and taxifolin, but leveled off close to control levels at 28 dpi. In the expanding outer lesion, on the other hand, significantly higher levels of catechin were detected at 14 and 28 dpi compared to the control. Catechin concentrations in the sterile agar control increased during the 28 days of this experiment and were significantly higher at 14 dpi compared to 2 dpi (p < 0.01). Dimeric PAB1 content in E. polonica inoculated tissue did not increase as dramatically as catechin content, but also reached significantly higher levels compared to the control (p < 0.05). However no significant increases were observed in the control over time.

Taken together, these data show that the flavonoid content of spruce stems increases in response to fungal infection. Monomeric aglucone concentrations increased more rapidly and to higher levels in fungus-infected tissue relative to the sterile agar-inoculated controls. Monomers and aglucones also increased more rapidly compared to their glycosylated and dimeric derivatives. Among the aglucones, the content of the dihydroflavonol taxifolin and the flavan-3-ol catechin were significantly higher compared to that of the flavanones, naringenin and eriodictyol.

#### A Physiologically Relevant Mixture of the Flavonoids Taxifolin and Catechin Is More Toxic to E. polonica Than the Single Compounds

Due to their high and rapid accumulation in fungus-infected spruce tissue, the anti-fungal activity of catechin, and taxifolin was assessed in an in vitro assay. Potato dextrose agar, allowing for optimal growth of E. polonica (negative control), was amended with either 2 mg ml−<sup>1</sup> catechin, 2 mg ml−<sup>1</sup> taxifolin, or 2 mg ml−<sup>1</sup> of a 2:1 mixture of catechin and taxifolin (similar to ratios observed in the fungus-inoculated treatments, **Figure 2**). Fungal growth was measured daily and the ultimate growth rate was calculated from the gradient of the growth curve (**Figure 3A**).

Both catechin and taxifolin significantly inhibited fungal growth. However, taxifolin was more toxic to the fungus than catechin. Interestingly, a mixture of both compounds in amounts equivalent to what was observed in fungus-induced spruce tissue (**Figure 2**) was even more toxic to E. polonica than the individual compounds alone (p < 0.001; **Figure 3A**). Catechin and taxifolin therefore have a synergistic antifungal effect. Based on these data we hypothesize that taxifolin plays an important role in the chemical defense of Norway spruce against infection by the bark beetle associated fungus E. polonica.

(significance codes: <sup>∗</sup>p < 0.04; ∗∗p < 0.01; n = 5; error bars = SE).

# Catechin and Taxifolin Are Toxic to Bark Beetles and Inhibit Tunneling

To assess if catechin and taxifolin have an effect on I. typographus, an assay was conducted where bark beetles were forced to tunnel into medium containing either 1 mg ml−<sup>1</sup> catechin, taxifolin, or 0.5 mg ml−<sup>1</sup> of both compounds (similar to ratios observed in the wounded control treatments, **Figure 2**). After 6 h, the tunnel lengths of beetles feeding either on taxifolin or a 1:1 mixture of taxifolin and catechin were significantly shorter (p < 0.001) than those of beetles feeding on unamended control medium or on catechin (**Figure 3B**). This result suggests that taxifolin is an anti-feedant for I. typographus. Beetles were allowed to continue

FIGURE 3 | Taxifolin and catechin are defense compounds in Norway spruce that act against both bark beetles and an associated fungus. (A) E. polonica grew significantly slower on potato dextrose agar amended with 2 mg ml−<sup>1</sup> catechin, 2 mg ml−<sup>1</sup> taxifolin, or a 2 mg ml−<sup>1</sup> mixture of catechin and taxifolin (2:1) (n = 4; error bars = SE; p < 0.001). (B) Ips typographus tunneled significantly further into bark medium containing 1 mg m−<sup>1</sup> catechin or unamended medium compared to medium with a 1:1 catechin-taxifolin mixture (n = 15; error bars = SE; p < 0.001). (C) Bark beetles gained more weight feeding on unamended bark medium without flavonoids than on medium amended with 1 mg ml−<sup>1</sup> catechin, taxifolin, or both (n = 15; error bars = SE; p < 0.05). Different letters refer to statistically significant differences at p < 0.05.

tunneling and feeding for 48 h before they were weighed. Beetles that fed on non-amended control diet gained significant weight after feeding for 48 h (pairwise t-test before and after feeding for 48 h on diet: t<sup>14</sup> = −2.43, p = 0.03). Beetles feeding on either catechin or taxifolin amended diets did not gain weight (pairwise t-tests: catechin: t<sup>13</sup> = −0.84, p = 0.42, taxifolin: t<sup>14</sup> = −0.47, p = 0.65). Interestingly, insects feeding on the diet amended with the mixture of catechin and taxifolin gained weight, but the weight gain was not statistically significant (pairwise t-test: t<sup>14</sup> = −1.89, p = 0.08; **Figure 3C**). These results suggested that the added compounds had an anti-nutritive effect and that weight gain was affected mainly by the concentration of individual compounds rather than the composition of mixtures.

#### Taxifolin Is Synthesized by a Flavanone-3-Hydroxylase Enzyme, Which Is Encoded by a Single Gene in Norway Spruce and Is Highly Expressed in Response to Fungal Infection

The biosynthesis of catechin is well studied in Norway spruce (Hammerbacher et al., 2014), but little is known about the formation of taxifolin in this species. It is known that dihydroflavonols, such as taxifolin, are synthesized from flavanones by 2-oxoglutarate dependent dioxygenases known as flavanone 3-hydroxylase (F3H) enzymes. Only one putative protein sequence with high identity to other characterized F3H enzymes could be retrieved from the P. abies High Confidence Genes database<sup>5</sup> . On the other hand, Warren et al. (2015) annotated five open reading frames in the P. glauca genome as F3H-encoding genes. We manually re-annotated these genes and eliminated four of them, which encoded enzymes with high similarity to other oxoglutarate-dependent dioxygenase enzymes, including a flavonol synthase as well as an anthocyanidin reductase enzyme. Furthermore, none of the inferred protein sequences from Warren et al. (2015) corresponded to the transcript sequences deposited on NCBI for P. sitchensis (Ralph et al., 2008), or P. abies<sup>5</sup> , both close relatives of P. glauca. BLAST searches of available transcriptomes of P. abies, P. glauca, P. sitchensis, Pinus radiata, and P. taeda using characterized F3H protein sequences from tea (Singh et al., 2008), apple (Halbwirth et al., 2002), and petunia (Britsch et al., 1992) F3H protein sequences resulted in the retrieval of a single high-confidence F3H transcript for each species.

Phylogenetic analysis of the putative P. abies F3H protein and other F3H proteins from both angiosperms and gymnosperms revealed a clear separation of F3H protein sequences from the Pinaceae. F3H from the Pinaceae formed a clade with Gingko biloba, which is also a gymnosperm. P. abies F3H was 99% similar to P. sitchensis and 96% similar to the P. glauca F3H enzyme, but only 89% similar to the F3H enzymes from the two Pinus species included in the analysis and 84% similar to that from G. biloba. Angiosperm F3H sequences were clearly separated from the gymnosperm sequences and had on average 66–70% amino acid identity to the F3H sequence from P. abies (**Figure 4A**).

The expression of the candidate PaF3H gene was analyzed in fungus-inoculated and sterile agar inoculated-control tissue using quantitative real-time PCR (**Figure 4B**). Compared to the saplings inoculated with sterile agar, a significant increase in transcript abundance was observed in fungus-inoculated tissue 7 and 14 dpi in the inner lesion (p < 0.01). The outer lesion showed an increase in transcription at 14 dpi, which was equivalent to that observed for the inner lesion. The PaF3H gene expression correlated well with taxifolin accumulation in the same tissues (**Figure 2C**), which increased at 7 dpi in the inner lesion above control levels and reached a maximum at 14 dpi in both inner and outer lesions. This supported our hypothesis that the PaF3H gene is involved in the biosynthesis of taxifolin in defense-induced Norway spruce tissue.

<sup>5</sup>http://congenie.org

FIGURE 4 | The PaF3H gene encodes the enzyme for the production of taxifolin in Norway spruce. (A) Phylogenetic analysis revealed separate clades for F3H enzymes from gymnosperms and angiosperms (500 boot-strap replicates). Scale bar shows number of substitutions per site. (B) F3H gene expression increased significantly 7 and 14 days post inoculation compared to the wounded control (significance code: ∗∗p < 0.01; n = 5; error bars = SE). Transcript abundance was measured by qRT-PCR normalized against PaUBI and calibrated against abundance in a sterile agar-inoculated control sample. Different letters refer to statistically significant differences at p < 0.05.

### PaF3H Encodes a Functional Flavanone 3-Hydroxylase Enzyme and Is Involved in Taxifolin Biosynthesis in Norway Spruce

The open reading frame of PaF3H was cloned and the encoded protein was heterologously expressed in E. coli. The substrates naringenin or eriodictyol were added to the medium to a final concentration of 0.2 mg ml−<sup>1</sup> 2 h after induction of protein expression using IPTG. 12 h after adding the substrate to the growth medium, it was extracted with acidic ethyl acetate and analyzed using LC-MS. Addition of naringenin to the medium resulted in the accumulation of dihydrokaempferol (data not shown) and addition of eriodictyol resulted in taxifolin accumulation (**Figure 5A**). Medium-only as well as non-transformed E. coli controls were included in the assays. Peak identities and retention times were confirmed with pure standards. PaF3H therefore encodes an enzyme, which can accept both naringenin and eriodictyol as substrates (**Supplementary Table S3**).

To study the in planta role of PaF3H, an RNA-interference (RNAi) construct was made for this gene and Norway spruce callus was transformed with it. Transgenic calli were selected and induced to differentiate into saplings, which were analyzed for their flavonoid content. Three transgenic lines were obtained with active expression of the RNAi construct. Due to the RNAi constructs, transcript abundance of the PaF3H gene was around 4- to 8-fold lower in the transgenic lines compared to the trees transformed with the vector only (V-control; **Figure 5B**). Consequently the accumulation of taxifolin (**Figure 5C**), taxifolin-3-glucoside (**Figure 5D**) and catechin (**Figure 5E**) was significantly lower in the transgenic F3H-RNAi lines than in the vector control line (p < 0.01). Taken together, these data show that PaF3H encodes a functional F3H, which contributes significantly to the accumulation of taxifolin as well as down-stream products of the pathway, including taxifolin-3-glucoside and catechin.

#### DISCUSSION

Due to climate change, conifers adapted to cool, wet conditions, such as Norway spruce, are currently experiencing increasing heat and drought stress (Netherer et al., 2015). This renders them susceptible to attack by bark beetles and their fungal associates, which can cause extensive forest losses under warm and dry conditions (Krokene, 2015). To offset climate-driven tree mortality, efforts are underway to identify resistance mechanisms in conifers that may render them partially tolerant to insect and pathogen attacks (Keeling and Bohlmann, 2006; Schiebe et al., 2012). In our study, we focused on Norway spruce, which is a host for the bark beetle I. typographus and its fungal associates. This tree species produces terpenoid oleoresins and polyphenols in response to bark beetle-fungus attacks (Schiebe et al., 2012), but which compounds are important in warding off such attacks is still unclear. To address this question systematically, we focused in this study on flavonoid biosynthesis in Norway spruce and the bioactivities of the most abundant flavonoids produced by the tree during simulated bark beetle attack.

#### Two Classes of Flavonoids, Dihydroflavonols and Flavan-3-ols, Accumulate in the Cambial Zone of Norway Spruce After Fungal Infection

In addition to acetophenones, hydroxycinamic acids, coumarins, lignans, neolignans, and stilbenes, Norway spruce produces a diverse range of flavonoids, especially in foliage. Strack et al. (1989) showed, using nuclear magnetic resonance spectroscopy, that spruce needles contained nine different flavonoids, including the flavan-3-ols catechin and gallocatechin, the glucosides of kaempferol, quercetin, and isorhamnetin as well as the dihydroflavonol taxifolin. In addition, Slimestad and Hostettmann (1996) analyzed spruce needles with LC-MS and found the glucosides of myricetin, syringetin,

Hammerbacher et al. Flavonoid Defenses in Norway Spruce

mass spectrometry.

laricitrin and taxifolin, rutinosides of kaempferol, quercetin, myricetin and syringetin as well as acylated glucosides of kaempferol, quercetin, laricitrin, syringetin and isorhamnetin. However, the diversity of flavonoids in bark and wood seems more limited. Only the flavonols quercetin, kaempferol, laricitrin, myricetin and their glucosides, the dihydroflavonol taxifolin and its glucoside as well as a diverse range of flavan-3-ols have thus far been identified in the bark of spruce trees in significant concentrations (Brignolas et al., 1995; Li et al., 2012; Schiebe et al., 2012; Hammerbacher et al., 2014, 2018). In this

study we specifically analyzed the flavonoids in developing xylem and phloem and detected high levels of taxifolin, its glucoside, catechin and its dimeric derivative PAB1 as well as their metabolic precursors, naringenin and eriodictyol. Interestingly, concentrations of all these compounds increased significantly after simulated bark beetle attack compared to the wounded control.

The flavanones naringenin and eriodictyol are metabolic precursors for down-stream metabolites in the flavonoid pathway (**Figure 1**) and thus only accumulated at concentrations that

abundance was measured by qRT-PCR normalized against PaUBI and calibrated against the abundance of RNAi line 1. Metabolites were analyzed by LC-tandem

were three to fifteen times lower than the other flavonoids identified in our study. These flavanones could possibly play an important role in Norway spruce defense against bark beetles and their associated fungi, for example as precursors. However, due to their low concentrations and due to the fact that bark beetles were reported to carry bacterial symbionts that can degrade them (Cheng et al., 2018) we did not study their direct defensive properties against I. typographus and its associated fungus E. polonica.

Two of the major end products of the flavonoid pathway in Norway spruce are flavan-3-ols, the monomer catechin and the catechin-epicatechin dimer, PAB1 (Hammerbacher et al., 2014). We found that these compounds accumulated in response to inoculation with the bark beetle associated fungus E. polonica, which is in agreement with previous studies (Brignolas et al., 1995; Hammerbacher et al., 2014). Interestingly, a strong and rapid accumulation of these compounds in Norway spruce correlated with resistance to both bark beetle attack (Brignolas et al., 1995) as well as spruce rust (Chrysomyxa rhododendri) infection (Ganthaler et al., 2014). Furthermore, an allele of leucoanthocyanidin reductase 3 (PaLAR3), one of the four genes encoding the enzyme catalyzing the conversion of leucoanthocyanidin into catechin, was identified as a quantitative trait locus in Norway spruce accounting for 27% lower growth of Heterobasidion root rot fungus in sapwood (Nemesio-Gorriz et al., 2016). All this evidence therefore suggests that elevated catechin biosynthesis is a general anti-fungal defense in Norway spruce.

Another class of flavonoids, dihydroflavonols, is also induced by fungal attack on Norway spruce. We showed that high concentrations of taxifolin and its glucoside accumulated in response to infection with a bark beetle-associated fungus, as would happen during a bark beetle attack. Norway spruce was also shown to accumulate taxifolin in response to low density E. polonica infection (Evensen et al., 2000; Krajnc et al., 2014) as well as during infection by the spruce rust, C. rhododendri (Ganthaler et al., 2014). Similarly, other conifer species were also shown to accumulate taxifolin in response to fungal infection. For example, in Sitka spruce, increased taxifolin content was recorded after Heterobasidion annosum infection (Deflorio et al., 2011), and Pinus nigra responded to inoculation with Diplodia sapinea by producing both taxifolin and catechin (Bonello and Blodgett, 2003). The constitutive levels of taxifolin may also be important in spruce-fungus interactions since high constitutive levels of taxifolin glucoside content correlated positively with resistance to spruce rust (Ganthaler et al., 2014), and the absence of this compound was thought to partially explain spruce susceptibility to bark beetle attack (Brignolas et al., 1995).

# Flavanone-3-Hydroxylase Expression Is Correlated With Taxifolin and Catechin Accumulation

In contrast to flavan-3-ol biosynthesis in Norway spruce, which was elucidated in previous work (Hammerbacher et al., 2014, 2018), very little is known about the synthesis of the dihydroflavonol taxifolin. We found a single high confidence transcript encoding a F3H enzyme in publically available gene sequence resources, and a single copy of this gene was also recorded in angiosperms such as petunia (Britsch et al., 1992), Arabidopsis (Pelletier and Shirley, 1996) and tea (Mahajan and Yadav, 2014). The existence of only a single copy of F3H in many plant species might be explained by the fact that other oxoglutarate-dependent dioxygenases in the flavonoid pathway such as flavonol synthase and anthocyanidin synthase can partially complement F3H activity (Buer et al., 2010).

In this study we showed that in Norway spruce F3H gene expression increased after infection with E. polonica and this correlated closely with taxifolin and catechin accumulation. On the other hand, in our transgenic lines in which F3H transcription was silenced by RNAi, lower levels of both metabolites were detected. This is in agreement with many other studies. For example, overexpression of F3H from plants such as Lycium chinense and tea in tobacco resulted in increased expression levels of late flavonoid biosynthesis genes and a significant increase in flavan-3-ols (Mahajan and Yadav, 2014; Song et al., 2016). F3H gene expression also correlates with plant developmental stages in which flavan-3-ols accumulate. This was shown in young and old tea leaves (Singh et al., 2008) as well as during heartwood development in black walnut (Beritognolo et al., 2002). All this evidence suggests that, like in other plant species, F3H plays an important role in the accumulation of not only taxifolin, but also down-stream metabolites in the flavonoid pathway of Norway spruce such as catechin and PAs.

#### Taxifolin and Catechin Are Toxic to Both Bark Beetles and Their Fungal Associates

Accumulation of down-stream metabolites in the flavonoid pathway has often been shown to be a defense response to pathogen infection or herbivory. For example, the high accumulation of the flavan-3-ol catechin in poplar results in increased resistance to rust infection (Ullah et al., 2017). In this study we also demonstrated the in vitro anti-fungal activity of catechin against E. polonica at concentrations approximating the maximum amounts of compound produced by the saplings in this study. Interestingly, taxifolin was found to be even more toxic than catechin and a physiologically relevant mixture of catechin and taxifolin at the 2:1 ratio observed in pathogen-challenged cambium had an even greater effect than either of the compounds alone. Although taxifolin is known as a potent anti-oxidant (Weidmann, 2012) its function as an anti-microbial compound has not been fully elucidated yet. It may act by blocking the biotransformation or degradation of catechin. In previous work, the virulence of E. polonica was linked to its ability to degrade stilbene glucosides, flavan-3-ols such as catechin and procyanidin B1 and other spruce phenolics by oxidation, deglucosylation and by oxidative dimerization (Hammerbacher et al., 2013; Wadke et al., 2016). However, taxifolin was not a substrate for the fungal catechol dioxygenase enzymes involved in the oxidative cleavage of host phenolics and so might act as a competitive inhibitor (Wadke et al., 2016). Similarly, in studies on the grape

pathogen Botryotinia fuckeliana, taxifolin-rhamnoside, and other phenolics were found to be strong inhibitors of the laccases (stilbene oxidases) that are responsible for oxidative stilbene dimerization (Goetz et al., 1999). The synergistic anti-fungal effect of the taxifolin-catechin mixture might therefore result from taxifolin-mediated inhibition of enzymes necessary for the detoxification of catechin. Synergy may also act at the level of sensory perception or by preventing compensatory herbivore feeding (Gershenzon et al., 2012).

Our results suggest that both catechin and taxifolin inhibit bark beetle tunneling and that the ingestion of these compounds had a visible effect on weight gain on this species. However, there were no differences in activity between catechin and taxifolin, and a mixture of the two compounds was no more toxic than the single compounds. The bioactivity of phenolics against insect herbivory is still debated (see Barbehenn and Constabel, 2011 for detailed review), but seems to depend on insect species. For example, catechin was identified as a feeding stimulant for willow leaf beetles (Doskotch et al., 1970; Kuokkanen et al., 2003) in contrast to other studies showing that this compound is a defense against herbivore feeding in oak trees (Feeny, 1970; Moctezuma et al., 2014). However, in the case of the spruce bark beetle, a previous study using a different bioassay set-up also showed that catechin and taxifolin inhibited tunneling behavior in I. typographus (Faccoli and Schlyter, 2007), supporting our hypothesis that high concentrations of these compounds could result in increased resistance to bark beetle attack in Norway spruce. Sensitivity to phenolics might also partially explain the close association of I. typographus with blue stain fungi known to degrade these compounds (Zhao et al., 2018).

Although our bioassays showed that both taxifolin and catechin play a role in tree defense against I. typographus and its associated fungi, other compounds might play a similar or even more important role in tree resistance. Concentrations of taxifolin-glucoside as well as PAB1 also increased significantly in response to fungal infection in this study. Experimental evidence is available showing that PAs have similar anti-fungal activity as catechin. For example, Ullah et al. (2017) showed that growth inhibition of the poplar rust fungus was equivalent when challenged with catechin and PAs. Similar results were also observed in growth trials of E. polonica with different concentrations of flavan-3-ols (Hammerbacher et al., 2014). However, the bioactivities of catechin and condensed tannins against insects are still not well understood and it is unknown if monomers are more active against insect feeding than polymers. It is also not known if taxifolin and its glucoside have similar modes of action. In a study by Auger et al. (1994) the biological activities of taxifolin and its glucoside were studied in pine against Dipirion pini, the common pine sawfly. The authors reported a toxic effect of both compounds on the sawfly but could not resolve the degree of toxicity of each. It is thought that glycosides are usually less toxic than the aglycone and that glycosylated compounds are activated during the plant defense response by deglycosylation (Morant et al., 2008). It was shown that E. polonica deglucosylates stilbenes (Hammerbacher et al., 2013). It is therefore possible that taxifolin glucoside can also be deglucosylated by the fungus, producing a more toxic compound that might impact negatively on its bark beetle vector. Clearly more research must be conducted to resolve the degree of toxicity of different flavonoids and other phenolic compounds in bark beetles and their associated fungi. Transgenic trees, such as the F3H RNAi lines produced in this study would be extremely valuable in addressing these questions, as in vitro experiments alone cannot reflect the full complexity of a living plant.

#### CONCLUSION

In this study we identified two down-stream products of the flavonoid pathway in Norway spruce with activity against the bark beetle I. typographus as well as its fungal associate E. polonica. We also elucidated an unknown step in the Norway spruce flavonoid pathway for the biosynthesis of these compounds by characterizing a novel F3H enzyme. In plants, high F3H gene expression is not only known to increase resistance against biotrophic and necotrophic attackers (Mahajan and Yadav, 2014; Sun et al., 2016), but is also involved in increasing tolerance to cold, drought and heat stress (Mahajan and Yadav, 2014; Meng et al., 2015; Song et al., 2016). Selecting conifer genotypes with high F3H gene expression might therefore partially counteract the detrimental effects of climate change, as these genotypes might not only be more resistant to bark beetle attack, but also more resistant to abiotic stress which is the indirect trigger of large bark beetle outbreaks.

#### AUTHOR CONTRIBUTIONS

All authors wrote the manuscript. AH and JG funded this project out of grants received from the Max Planck Society and the University of Pretoria. AH, DK, AS, CU, and LPW performed the experiments and data analysis.

# FUNDING

This work was funded by the Max Planck Institute for Chemical Ecology and the University of Pretoria RDP program.

# ACKNOWLEDGMENTS

We would like to thank Marion Staeger and the horticulturalists of the MPI-CE for assistance in maintaining the plant material and Friz Schoemberg for grinding the samples. We also would like to thank Emily Ann Puckett for assisting with the bark beetle feeding assays.

# SUPPLEMENTARY MATERIAL

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

# REFERENCES

fpls-10-00208 February 21, 2019 Time: 17:42 # 13



dependence of the enzymes. FEBS Lett. 361, 299–302. doi: 10.1016/0014- 5793(95)00199-J


show differences in photosynthetic recovery after drought stress as reflected in gene expression profiles. Plant Sci. 171, 745–758. doi: 10.1016/j.plantsci.2006. 07.010


influence bark beetle tunneling behavior. Fungal Ecol. (in press). doi: 10.1016/j. funeco.2018.06.003

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

Copyright © 2019 Hammerbacher, Kandasamy, Ullah, Schmidt, Wright and Gershenzon. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Inducibility of Plant Secondary Metabolites in the Stem Predicts Genetic Variation in Resistance Against a Key Insect Herbivore in Maritime Pine

Xosé López-Goldar1,2,3,4 \*, Caterina Villari2,5, Pierluigi Bonello<sup>2</sup> , Anna Karin Borg-Karlson<sup>3</sup> , Delphine Grivet4,6, Rafael Zas<sup>1</sup> and Luís Sampedro<sup>1</sup>

<sup>1</sup> Misión Biológica de Galicia, Consejo Superior de Investigaciones Científicas, Pontevedra, Spain, <sup>2</sup> Department of Plant Pathology, The Ohio State University, Columbus, OH, United States, <sup>3</sup> Ecological Chemistry Group, Department of Chemistry, Royal Institute of Technology, Stockholm, Sweden, <sup>4</sup> Department of Forest Ecology and Genetics, Forest Research Centre, INIA, Madrid, Spain, <sup>5</sup> Daniel B. Warnell School of Forestry and Natural Resources, University of Georgia, Athens, GA, United States, <sup>6</sup> Sustainable Forest Management Research Institute, INIA-University of Valladolid, Palencia, Spain

Edited by:

Judy Simon, Universität Konstanz, Germany

#### Reviewed by:

Axel Schmidt, Max-Planck-Institut für Chemische Ökologie, Germany Stefanie Wienkoop, Universität Wien, Austria

> \*Correspondence: Xosé López-Goldar xlgoldar@mbg.csic.es

#### Specialty section:

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

Received: 07 August 2018 Accepted: 24 October 2018 Published: 21 November 2018

#### Citation:

López-Goldar X, Villari C, Bonello P, Borg-Karlson AK, Grivet D, Zas R and Sampedro L (2018) Inducibility of Plant Secondary Metabolites in the Stem Predicts Genetic Variation in Resistance Against a Key Insect Herbivore in Maritime Pine. Front. Plant Sci. 9:1651. doi: 10.3389/fpls.2018.01651 Resistance to herbivores and pathogens is considered a key plant trait with strong adaptive value in trees, usually involving high concentrations of a diverse array of plant secondary metabolites (PSM). Intraspecific genetic variation and plasticity of PSM are widely known. However, their ecology and evolution are unclear, and even the implication of PSM as traits that provide direct effective resistance against herbivores is currently questioned. We used control and methyl jasmonate (MJ) induced clonal copies of genotypes within families from ten populations of the main distribution range of maritime pine to exhaustively characterize the constitutive and induced profile and concentration of PSM in the stem phloem, and to measure insect herbivory damage as a proxy of resistance. Then, we explored whether genetic variation in resistance to herbivory may be predicted by the constitutive concentration of PSM, and the role of its inducibility to predict the increase in resistance once the plant is induced. We found large and structured genetic variation among populations but not between families within populations in resistance to herbivory. The MJ-induction treatment strongly increased resistance to the weevil in the species, and the genetic variation in the inducibility of resistance was significantly structured among populations, with greater inducibility in the Atlantic populations. Genetic variation in resistance was largely explained by the multivariate concentration and profile of PSM at the genotypic level, rather than by bivariate correlations with individual PSM, after accounting for genetic relatedness among genotypes. While the constitutive concentration of the PSM blend did not show a clear pattern of resistance to herbivory, specific changes in the chemical profile and the increase in concentration of the PSM blend after MJ induction were related to increased resistance. To date, this is the first example of a comprehensive and rigorous approach in which inducibility of PSM in trees and its implication in resistance was

**71**

analyzed excluding spurious associations due to genetic relatedness, often overlooked in intraspecific studies. Here we provide evidences that multivariate analyses of PSM, rather than bivariate correlations, provide more realistic information about the potentially causal relationships between PSM and resistance to herbivory in pine trees.

Keywords: genetic variation, herbivory, inducibility, maritime pine, plant secondary metabolites (PSM), phenolics, resistance, terpenes

#### INTRODUCTION

Resistance to herbivores and pathogens has been widely recognized as a key plant trait with strong adaptive value (Futuyma and Agrawal, 2009), particularly at early stages of development (Goodger et al., 2013). In long-lived plants with high apparency, like trees, resistance against pests and pathogens usually relies on high concentrations of a diverse array of plant secondary metabolites (PSM) (Feeny, 1976; Wiggins et al., 2016). Such chemical cocktail generally produces dose dependent effects, where higher concentrations in specific plant tissues targeted by herbivores results in reduced enemy performance and/or plant damage (Zhao et al., 2011a). These chemical defenses are expensive to produce and maintain (Gershenzon, 1994; Sampedro et al., 2011a). All these features are certainly present in conifers, whose higher concentrations of terpenes and phenolics are commonly associated with increased resistance against insect herbivores (Phillips and Croteau, 1999; Witzell and Martín, 2008; Schiebe et al., 2012).

Within conifer species, variation in defensive PSM is usually difficult to explain based exclusively on an arms race between trees and their biotic aggressors because many other factors often affect plant–organism interactions (Moore et al., 2013). For example, heterogeneity in the biotic and abiotic environment (Lombardero et al., 2000; Sampedro et al., 2010), the large, repetitive, and idiosyncratic structure of conifer genomes (Theis and Lerdau, 2003), trade-offs with other life functions (Sampedro et al., 2011a; Moreira et al., 2014) and other evolutionary forces, such as gene flow and genetic drift (Kremer et al., 2012), can all be involved in maintaining large intraspecific variation in PSM as well as in resistance. Phylogenetic analyses has been widely used in recent years (Rasmann and Agrawal, 2011; Pellissier et al., 2012; Moreira et al., 2014; Carrillo-Gavilán et al., 2015) to control for the statistical dependence among related species in macroecology (Freckleton et al., 2002), leading to the identification of imprinted historical events in contemporary observed traits (Pagel, 1999). At the intraspecific level, demographic processes, habitat fragmentation, gene flow and local adaptation can also greatly influence the variation in a species' traits among and within populations (Savolainen et al., 2007; De-Lucas et al., 2009), making necessary to account for the genetic structure of populations and the relative kinship among genotypes within families when exploring trait-trait or trait-environment associations (Yu et al., 2006).

Evidence is accumulating that PSM are highly plastic, as most plant species respond with large changes in the concentration and profile of PSM to the perception of a biotic threat, such as a pathogen or a herbivore (Karban et al., 1999; Heil, 2009). This sort of biochemical plasticity, referred to as inducibility, may be a very important ecological and evolutionary plant trait, as it is a source of functional phenotypic variation in resistance and thus, potentially, a trait under selection (Agrawal et al., 2015). Considering a single plant genotype, inducibility is usually estimated as the difference between the induced vs. the constitutive concentration of PSM. In large plants, estimating inducibility of a PSM may require a repeated measures design within the same subject, i.e., before and after challenge with real insect herbivory (Morris et al., 2006). However, a common undesired side-effect is that sampling plant tissues may itself be an inducer of plant defenses, which would alter the estimation of inducibility (see Moreira et al., 2012, for example). Even more challenging can be the establishment of any relationship between inducibility of PSM and the effective resistance against real herbivory. This requires knowing the concentration of PSM before and after exposure of both constitutive and induced phenotypes to real herbivory when measuring plant damage or insect performance. For small plants and to evaluate plant traits that require destructive sampling, this problem can only be solved using clonal replicates of the same genotype.

At the intraspecific level, while genetic variation and plasticity of PSM are widely recognized, the ecology and evolution of PSM are far from being clear, and the role itself of PSM as traits providing direct effective resistance against herbivores is currently under debate (Carmona et al., 2011; Cárdenas et al., 2014). Despite the common finding of positive associations between concentration of PSM and resistance to specific herbivores (Borg-Karlson et al., 2006; Erbilgin et al., 2006; Kannaste et al., 2009; Zhao et al., 2011a), cause-effect studies demonstrating the direct role of particular PSM or blends of PSM on resistance are much less frequent (Agrawal, 2005; Iason et al., 2011).

In this study we tested the hypothesis that intraspecific variation in effective resistance of maritime pine (Pinus pinaster) to a key insect herbivore, the pine weevil (Hylobius abietis), may be predicted by variation in inducibility of PSM. We used clonally replicated genotypes with known family structure from 10 populations representative of the main distribution range of maritime pine. We performed an exhaustive characterization of constitutive concentrations and profiles of PSM in the stem phloem and their associated induced changes after application of methyl jasmonate (MJ), a plant hormone involved in signaling of chewing herbivore damage (Sampedro et al., 2011b; Schiebe et al., 2012; Zas et al., 2014). Both constitutive and MJ-induced clonal copies were exposed to insect herbivory to measure plant damage as a proxy of resistance (i.e., the greater the damage, the lower the resistance). Then we explored to what extent resistance

to herbivory may be predicted by the constitutive concentration of PSM in the stem, and the role of its inducibility in predicting the increase in resistance once the plant is induced. Since the non-independence due to relatedness among genotypes, families and populations may affect the relationships between PSM and resistance, we accounted for the population structure and the relative kinship among individuals within populations.

## MATERIALS AND METHODS

#### Study System

We used P. pinaster as a model tree system, a Mediterranean conifer species with a fragmented distribution displaying large and structured among-population genetic variation (Vendramin et al., 1998; Bucci et al., 2007; Jaramillo-Correa et al., 2015). Populations largely differ in life history traits (Tapias et al., 2004; Grivet et al., 2013), including those putatively related to resistance to biotic stress and PSM (Sampedro et al., 2011b; Zas et al., 2015).

The pine weevil H. abietis is a bark-chewing insect herbivore with Palaearctic distribution, and is specialized in young conifers, feeding on several species (Leather et al., 1999). It harbors low neutral genetic variation across its large distribution range (Conord et al., 2006). The pine weevil is a harmful forest pest that causes severe damage on young conifer seedlings, and has a strong impact on plant fitness and early survival (Langström and Day, 2004; Zas et al., 2006, 2014). Host resistance has been showed to be related to plant survival and to be heritable (Zas et al., 2005). The pine weevil became a model insect for studies of resistance in pine trees (see for instance Nordlander et al., 2003; Wainhouse et al., 2005; Björklund, 2007; Nordlander et al., 2011; Sampedro et al., 2011b; Zas et al., 2014, 2017). The success in the regeneration of conifer forests is highly contingent on this insect (Langström and Day, 2004), which causes important economic losses in large part of Europe (Leather et al., 1999). Larval stages feed on stump roots of recently cut trees (Nordenhem and Nordlander, 1994), and newly emerged adults migrate in massive numbers to fresh clear-cuts in spring (Solbreck and Gyldberg, 1979; Tan et al., 2011). Generation times may vary from 1 to 3 years depending on the location in Europe (Day et al., 2004; Wainhouse et al., 2014). The longer generation times observed at higher latitudes, and the continuous masses of conifer forests, particularly in North Europe, strongly facilitate pine weevil dissemination, establishment, and subsequent damage among areas with clear-cuts for at least two up to four seasons after harvesting (Örlander et al., 1997). Apart from control and management measures to reduce the impact of the insect on plant fitness and survival (Leather et al., 1999; Nordlander et al., 2011), it is important to know the defensive potential in conifers to deal with this pernicious pest.

#### Plant and Herbivore Material

Ten maritime pine populations belonging to the collection 'CLONAPIN Bank 1' were selected along a latitudinal gradient across the natural distribution range of the species (**Figure 1**), from the Atlantic France to North Africa. Briefly, open-pollinated seeds were collected from five mother trees randomly selected in

natural stands within each of the ten populations, resulting in fifty half-sib families (i.e., known 'mother' and unknown 'father'). A separation of at least 50 m between selected mother trees was performed to minimize inbreeding (González-Martínez et al., 2006). Five seeds from each mother tree were sown and cultured in the facilities of SERIDA (Asturias, Spain), producing the complete clonal ortet collection, which comprised 250 genotypes (10 populations × 5 families × 5 half-sib genotypes). These genotypes were then clonally propagated by recurrent early cutting as described in Majada et al. (2011). A total of 205 out of the 250 original genotypes (belonging to 48 families out of the original 50) were eventually included due to low germination rate and/or rooting efficiency.

The cuttings were grown at SERIDA Agricultural Station (Grado, Asturias, Spain) in 200 mL containers containing peat and perlite (70:30 v/v), then transported to Norfor Forest Nursery (Figueirido, Pontevedra, Spain) and transplanted to 2 L pots. In September 2012, when pine cuttings were around 2-years old, the plants were transferred to an environmentally controlled greenhouse in the MBG-CSIC (Pontevedra, Spain). Plants were fertilized with a slow release fertilizer (Granum, Soaga SL, Vilanova de Arousa, Spain, NPK 11-22-9), supplemented with a foliar fertilizer (Fertimón, Progando SL, Galicia, Spain), and watered twice a week. Pines were allowed to acclimatize until the start of the experiment by December 2012. During this period, mean temperatures and relative humidity were recorded during the day (25.6 ± 0.06◦C and 57.2 ± 0.08%, respectively) and during the night (18.2 ± 0.05◦C and 69.0 ± 0.05%, respectively).

Adult weevils were collected at a clear-felled pine forest in Monte Castrove (Pontevedra, Spain, 42◦ 2702300N, 8 ◦ 4301900W) using Nordlander traps (Nordlander, 1987) and maintained in culture chambers at 10◦C several weeks until the start of the herbivory bioassays (see below). During that period, the pine weevils were fed weekly with fresh twigs of maritime pine.

#### Experimental Design

fpls-09-01651 November 20, 2018 Time: 16:5 # 4

Plants were assigned to treatments following a factorial split-plot design replicated in 5 blocks, with population as the whole plot factor (ten populations) and the factorial combination of family (3–5 families per population) and the MJ induction treatment (2 levels: control and induced) as the split factor. Two clonal replicates of each genotype from each family were included in each of the two MJ-induction treatments. One of these copies was used for chemical analyses and the other for evaluation of the effective resistance against the pine weevil. The experiment comprised 814 plants from 10 populations, with 3–5 families per population, 5 half-sibs (genotypes) per family, and 4 clonal replicates per genotype (two of them, control and MJ, used for chemical analysis; and two, control and MJ, used in the weevil bioassays). A detailed graphical scheme of the experimental design can be found in the **Supplementary Figure SM1**.

#### Induction Treatment

Plasticity of PSM in response to biotic stress has been extensively studied in maritime pine. Large changes in the concentration and chemical profile of total or specific PSM using either the exogenous application of MJ (Moreira et al., 2009), real herbivory (Moreira et al., 2013a), or both (Sampedro et al., 2011b; Moreira et al., 2012) have been reported in this and other pine species (Moreira et al., 2014; Carrillo-Gavilán et al., 2015). In December 2012, plants allotted to induction were transferred to a separate chamber prior to treatment application to avoid crosscontamination with control plants. Then, MJ-induced plants were sprayed with a suspension of 25 mM of 95% MJ (Sigma– Aldrich, #39270-7) in deionized water with 2.5% ethanol (v/v). Control plants were treated just with the carrier solution. Both solutions were sprayed over the aboveground part of each plant to runoff (∼10 mL per plant). MJ concentration and time elapsed between MJ application and sampling (1 month, see below) were based on previous studies (Moreira et al., 2013b; Zas et al., 2014). This lapse of time allowed plants to develop induced and longlasting anatomical structures (i.e., resin ducts and polyphenolic parenchyma cells) that are involved in the production of induced defensive PSM (Moreira et al., 2015).

#### Plant Sampling, Measurements and Analysis of Plant Secondary Metabolites

One month after MJ-application, two of the clonal copies of each genotype (one treated with MJ and one control) were destructively harvested by blocks, cutting the stem aboveground. Needles were separated from the stem and two 2-cm-long stem subsamples were immediately collected from its middle part, flash-frozen in liquid nitrogen, and stored at −80◦C. We then analyzed the composition and concentration of PSM in 3 out of the 5 blocks (260 plants). Briefly, terpenoids were extracted from one of the frozen stem subsamples with HPLC-grade hexane (HiPersolv Chromanorm #83992.320) and separated, identified and quantified by gas chromatography-mass spectrometry (GC-MS) and GC-flame ionization detection (GC-FID). Phenolics were extracted from the other stem subsample with HPLC-grade methanol (HiPersolv Chromanorm #152506X) and separated, identified and quantified by ultra high performance liquid chromatography-mass spectrometry (UHPLC-MS) and UHPLCdiode array detection (UHPLC-DAD). A detailed description of the analytical procedures and data processing prior to statistical analyses can be found in the **Supplementary Methods SM1**).

We identified 118 PSM that were classified into eight chemical groups [monoterpenes, sesquiterpenes, diterpenes, flavonoids, hydroxycinnamic acids (HCAs), lignans, eugenols, and fatty acids]. For the purpose of the current study, we only considered those compounds present in at least 5 individuals of the control or MJ-induced plants for each population. Ninety-eight chemical compounds were finally included in the eight chemical groups (**Supplementary Tables SM1**, **SM2**).

#### Herbivory Bioassays

At the time plants not subjected to herbivory were harvested for chemical analyses, the remaining two clonal replicates (one control and one treated with MJ) were used for evaluating the resistance of each genotype to the pine weevil, following López-Goldar et al. (2016), with modifications. Briefly, a transparent plastic tube was fitted around each plant and a randomly selected pair of weevils (one male and one female) were put on each plant and allowed to feed along the entire plant for 48 h. A fine-mesh was fitted at the top of the tube to prevent weevil escape. Weevils were starved in a Petri dish with moist filter paper (48 h at 15◦C in dark) and pre-weighed prior the herbivory treatment. After 48 h of exposure to herbivory along the entire plant, weevils were removed and the damage caused by the weevils in the stem of each plant was measured at the nearest mm<sup>2</sup> using calibrated area templates.

#### Genetic Variation in Resistance to Herbivory

Variation in the damage caused by H. abietis in the stem was analyzed by fitting a mixed model using the PROC MIXED procedure in SAS v9.4 (Littell et al., 2006), with the induction treatment (MJ), population (P), family within population [F(P)] and MJ × P and MJ × F(P) interactions as fixed factors. The MJ × P and MJ × F(P) interactions are informative of the variation in the responses to the induction treatment among populations and families within populations, respectively. Block (B), the B × P interaction and the genotype [MJ × F(P) × B] were considered random factors. The B × P interaction was considered a random factor in order to analyze the main effect of P with the appropriate error term in the split-plot design (Littell et al., 2006). The genotype was also included as a random factor to account for the dependence between the two clonal replicates subjected to herbivory within each block. Damage was square root-transformed to achieve normality. Heterogeneous residual variance models across P were fitted. Weevil weight, as a proxy of weevil size, and stem diameter of the plant, both known to affect feeding rate (Wainhouse et al., 2004; Wainhouse et al., 2005), were included as covariates because they improved model fit. The same model was fit to explore the genetic variation among genetic groups of populations (GP) and within GP (families within each genetic group). Genetic groups of populations were

identified and characterized elsewhere using molecular markers (Jaramillo-Correa et al., 2015; Serra-Varela et al., 2015), and corresponded to: Atlantic France (PLEU and MIMI), Atlantic Iberia (PTOV, CDVO, ARMY and SCRI), Central Iberia (COCA, ASPE), Southwestern Iberia (ORIA) and Morocco (TAMR) (**Figure 1**).

In order to quantify the relative contribution of each factor to the genetic variation among populations and families within populations in resistance to the pine weevil, the previous mixed model was run again assuming all factors as random factors. Variance components of each factor were estimated by restricted maximum likelihood.

#### Relationships Between PSM and Resistance to Herbivory

To explore the causal relationships between resistance and PSM in the stem of P. pinaster, we performed pairwise correlations and multiple regression analyses between the weevil damage and PSM concentration at the genotypic level. Exploring the implication of PSM in resistance to herbivory in individuals within populations across the main distribution range will allow us to predict which PSM or blend of PSM within a genotype might be the optimal to determine the highest resistance to herbivory in the entire species. The clonal structure of our pine collection allowed for this genotypic-based analysis.

Because genotypes belong to a hierarchical collection of families within populations and to different populations in maritime pine, we accounted for the non-independence among genotypes within families and populations to avoid genetic effects due to relatedness in the correlation analyses (see **Supplementary Methods SM2**). For that, we included the population structure (Q) and kinship (K) matrices at the genotypic level using a mixed-model approach as described in Yu et al. (2006). Briefly, Q and K matrices were constructed based on 126 SNPs: for the Q matrix, we performed a Bayesian cluster analysis in STRUCTURE v 2.3.4 (Pritchard et al., 2000), whereas for the K matrix we used SPAGeDi (Hardy and Vekemans, 2002) based upon Loiselle et al. (1995) kinship coefficients. The Q and K matrices were then incorporated to the mixed models in order to account for the multiple level of relatedness among individuals (see **Supplementary Methods SM2**).

In addition, we fitted a mixed model for constitutive (C) and for the inducibility (MJ – C) of each total and individual PSM and weevil damage. For each defensive mode, population (P) and family within population [F(P)] were considered as fixed factors. Block (B) and the B × P interaction were considered random factors as above. Q and K matrices were incorporated into each mixed model to account for the genetic relatedness among genotypes. The predicted values of the variables from each plant in each mixed model were then used for the pairwise correlations and multiple regression analyses. Because PSM data were available in 3 out of 5 blocks, pairwise correlations and multiple regression analyses with weevil damage were performed only in 102 genotypes. Given that we were not able to extract the predicted values for all 102 genotypes after implementing Q and K matrices, and also not all clonal copies of the same genotype were available, sample sizes were slightly lower than the original 102 genotypes for each defensive state.

Two types of evaluations were performed to explore the relationships between PSM and weevil damage in the stem of P. pinaster. First, to observe the basal resistance without interference of the induction treatment, we explored the relationships between constitutive (C) PSM and damage by H. abietis. Second, we explored the relationships between the estimated inducibility of each total and individual PSM (MJ – C) with the inducibility of resistance to the pine weevil after MJ application. For each type of analysis, Pearson correlations between the concentration of each total and individual PSM and weevil damage were performed using PROC CORR in SAS v9.4. To estimate all the pairwise correlations without inflating Type I error due to multiple testing, p-value adjustments were performed using the False Discovery Rate (FDR) for p < 0.05 (Benjamini and Hochberg, 1995). Also, for the same two evaluations (constitutive and inducibility), multiple regression analyses were performed with PROC REG in SAS v9.4 to explore to what extent the PSM blend predicted herbivory resistance using the stepwise model selection method on the 98 individual PSM. We applied a more stringent threshold than default parameters (α = 0.05; default, α = 0.15) when including and excluding variables during the stepwise method to reduce type I error due to the large number of PSM used. Predictor variables were previously standardized (mean = 0, variance = 1) and influential observations (measured with Cook's D) that did not meet the requirements were excluded from the analysis. Normality of the residuals was achieved using squareroot transformation of the response variable (weevil damage) when necessary.

# RESULTS

#### Sources of Variation in Resistance to Herbivory by Hylobius abietis in Pinus pinaster

Resistance to weevil herbivory largely and significantly varied among populations but not between families within populations (companion table in **Figure 2**). The highest constitutive resistance was observed for the Spanish Atlantic population PTOV, whereas the most constitutively susceptible was the French Atlantic population MIMI (**Figure 3**). Large genetic variation was also found among genetic groups of populations (F4,<sup>140</sup> = 4.73, P = 0.001), revealing similar patterns of constitutive resistance among populations of the same genetic group. The highest constitutive resistance was found for the Atlantic Spanish group (61.5 ± 7.3 mm<sup>2</sup> of damage), whereas the Atlantic French group was the most susceptible (94.9 ± 12.3 mm<sup>2</sup> of damage).

The MJ-induction treatment strongly increased resistance to herbivory in the stem of all populations (companion table in **Figure 2**), reducing the overall damage caused by the pine weevil to less than half compared to non-induced plants (**Figure 3**, inset panel). The French population PLEU showed the highest induced

FIGURE 2 | Summary of the mixed model testing the effects of the induction treatment with methyl-jasmonate (MJ), population (P), family within population [F (P)] and the MJ × P and MJ × F(P) interactions on the early resistance to herbivory by the pine weevil on plants from ten maritime pine populations. Variance components for each effect and the residual are shown in the companion pie chart. Weevil weight and plant diameter were included as covariates. Significant p-values (p < 0.05) are highlighted in bold. A total sample size of 205 genotypes was used.

resistance, whereas COCA emerged as the less resistant after MJ application (**Figure 3**). The inducibility of resistance (i.e., the plasticity to MJ application) was similar for all populations and families [no significant MJ × P and MJ × F(P) interactions, companion table in **Figure 2**], but differed among genetic groups (MJ × GP interaction, F4,<sup>140</sup> = 3.36, P = 0.012). Interestingly, greater inducibility of resistance after induction was found in both Atlantic genetic groups (71 and 67% in Atlantic France and Atlantic Spain, respectively) than in the other genetic groups (Central Spain, 37%; Southeastern Spain, 34%; and Morocco, 45%).

panel). A total sample size of 205 genotypes was used.

The main effects in our experimental design, namely MJinduction, population, family and their interactions explained 36% of the total variance in resistance. MJ-induction was by far the most relevant factor, absorbing most of the explained variation in resistance among the studied effects (>80% of the explained variance), with the contribution of all remaining factors being comparatively much lower (**Figure 2**).

#### Concentration of PSM in the Stem as Predictor of Damage by the Pine Weevil at the Genotypic Level

We performed pairwise correlations in order to explore the relationship between concentration of total and individual PSM in the stem with weevil damage, while accounting for population structure and genetic relatedness among individuals. Total and individual PSM did not show any significant negative relationship with weevil damage at the constitutive level (**Table 1** and **Supplementary Table SM3**). When we explored the inducibility, total flavonoids and total lignans showed significant negative relationships with weevil damage (**Table 1**). Consistently, when considering individual PSM, significant negative correlations were found for the inducibility of two lignan hexoside derivatives, and marginally significant for three lignan xyloside derivatives and taxifolin derivative 1, after FDR correction (**Supplementary Table SM3**). As an exception, eugenol was the only PSM that was positively related with weevil damage at the constitutive level (**Supplementary Table SM3**).

We used multiple regression analyses looking for the relative predictive value of individual PSM in the relationship between PSM and damage by H. abietis in the stem, while accounting for population structure and relative kinship among genotypes in P. pinaster. At the constitutive level, 18 individual PSM (10 terpenes and 8 phenolics; **Table 2**) significantly predicted 70% of the total variation in constitutive weevil damage (F14,<sup>66</sup> = 18.97, P < 0.001). Six PSM (4 terpenes and

TABLE 1 | Pearson correlations between damage by the pine weevil, as a proxy of resistance, and the concentration of total plant secondary metabolites (PSM) in the stem phloem of 102 genotypes from 10 populations, representing the main distribution range of maritime pine.



Analyses were performed separately for constitutive total PSM (C), and for the inducibility of total PSM (MJ – C) after methyl-jasmonate induction. Significant correlations are highlighted in bold after adjustment for multiple testing correction using false discovery rate (FDR) for p ≤ 0.05 (Benjamini and Hochberg, 1995). Data used for correlations were obtained from mixed models accounting for population structure and relative kinship among individuals. Genotypic sample size of each defensive state is indicated in brackets. Due to the fact that some clonal replicates (control or MJ-induced) were not available for several genotypes, data points in the inducibility dataset were slightly lower than the original number of genotypes used.

2 phenolics) were negatively related with weevil damage, as shown by their negative standardized regression coefficients (**Table 2**). With regard to inducibility, 12 PSM (8 terpenes, 3 phenolics and 1 fatty acid; **Table 2**) significantly predicted 68% of the total variance of weevil damage (F11,<sup>83</sup> = 17.32, P < 0.001). Eight PSM (6 terpenes, 2 phenolics) showed standardized regression coefficients that were negatively related to the variation in weevil damage after MJ-induction (**Table 2**).

#### DISCUSSION

In our study, large and structured genetic variation in resistance to herbivory by H. abietis in the stem was found among P. pinaster populations, as well as significant differences among genetic groups of populations in resistance after MJ-induction. In addition, the resistance patterns were strongly predicted by the multivariate concentration and profile of PSM of genotypes after accounting for the genetic relatedness, rather than by bivariate correlations with weevil damage. Furthermore, differences in the resistance patterns between the constitutive and inducibility strategies were observed. While a high proportion of the variance for weevil damage was explained by the constitutive PSM, the directionality of the pattern (i.e., resistance or susceptibility) was unclear. Conversely, inducibility of PSM, besides explaining also a large proportion of the variance in weevil damage, predicted greater inducibility of resistance mainly due to relevant plastic and specific changes in the profile and increases in the concentration of the chemical compounds after induction.

#### Inter- but Not Intrapopulation Genetic Variation and Increased Resistance to Herbivory After Induction

Our results show that P. pinaster populations are variable in resistance against the bark chewing herbivore H. abietis, but we found no evidence of variation between families within populations. Using the same clonal collection, Elvira-Recuenco et al. (2014) reported similar results for intraspecific variation in constitutive resistance against the fungal pathogen Fusarium circinatum. Atlantic populations showed greater constitutive resistance (with MIMI population as the exception in our study) than those from the Mediterranean Basin in both studies. This was, perhaps, an unexpected outcome, given the wide taxonomic dissimilarity of the two enemies. F. circinatum is considered an alien invasive pathogen in Europe (Wingfield et al., 2008) where, to date, it has been detected in nature in only one adult P. pinaster plantation in the northern Iberian Peninsula (Iturritxa et al., 2013). These results suggest that similar functional defensive traits may be effective against both biotic enemies, and that populations from growthprone environments are better defended against these biotic agents at early stages. Similarly, we found extensive genetic variation among genetic groups of populations. This suggests that the differences in resistance levels among individual populations might be due to greater similarity among genetically closely related populations than distantly related ones. We did not find, however, genetic variation among families within populations for resistance, which is contrasting with previous results using the same pine species (Zas et al., 2005; Elvira-Recuenco et al., 2014). This finding might be due to the strong effect of the induction treatment, which may have overwhelmed most of the variation between families, or likely to sample size, that might be too low for detecting family variation.

Boosting MJ signaling pathways by exogenous application of MJ to produce a phenotype with increased expression of chemical and physical defenses is a well-known experimental tool in conifers (Heijari et al., 2005; Zhao et al., 2011a; Moreira et al., 2013b; Zas et al., 2014) and was confirmed in this study, with MJ being by far the most relevant factor (>80% of the total explained variance in resistance). The lack of a significant MJ by population interaction suggests that the patterns of resistance to herbivory elicited by MJ are highly regulated and conserved in P. pinaster, and probably derived from common ancestral defensive mechanisms shared in the conifer lineage (Hudgins et al., 2004). Results from previous studies that explored the genetic variation in the inducibility of resistance to insect herbivores after MJ application reported contrasting results. MJ induction reduced herbivory damage in P. radiata (Moreira et al., 2013b) and P. pinaster (Sampedro et al., 2011b), but family genetic variation in the inducibility


TABLE 2 | Summary of the multiple regression analysis explaining the resistance to the pine weevil using the concentration of individual PSM as predictor variables in constitutive (upper part of the table) and in inducibility (lower part of the table) defensive modes in 102 genotypes from 10 maritime pine populations.

The regression coefficients (β) and partial R<sup>2</sup> of PSM included in the model after stepwise selection method are shown. PSM associated with resistance to weevil damage (negative β) are typed in bold. Linear regression model at each defensive mode is significant at p < 0.05. Genotypic sample size of each defensive state is indicated in brackets. Due to the fact that some clonal replicates (control or MJ-induced) were not available for several genotypes, sample size for the inducibility dataset was slightly lower than the original number of genotypes used. Unk P#, unidentified phenolic compound.

of resistance was only detected for the folivore Thaumetopoea pityocampa and not for H. abietis (Sampedro et al., 2011b; Moreira et al., 2013b), suggesting that resistance mechanisms would be differentially regulated within populations depending on the herbivore identity and/or the target tissue. Despite that we did not find a MJ × population interaction, we did find it out among genetic groups of populations, indicating that the inducibility of resistance after MJ application is genetically structured in the species, i.e., resistance plasticity is more similar among genetically closely related populations than among distantly related ones. Our results highlight the evolutionary relevance and adaptive value of induced defenses and inducibility as putative key traits explaining variation in resistance. In addition, our findings are in agreement with the predictions regarding the intraspecific patterns of variation in plant defense (Hahn and Maron, 2016), although more profound analyses are required to appropriately address the postulates of the theory, which are out of scope in this study.

### Accounting for Biases in Relationships Between PSM and Resistance Across and Within Species

While associations between concentration of PSM and resistance have been widely reported in conifers (Borg-Karlson et al., 2006; Erbilgin et al., 2006; Kannaste et al., 2009), cause-effect studies directly exploring the particular role of individual PSM or PSM blends with resistance have been little explored (Agrawal, 2005; Iason et al., 2011). Moreover, recent studies addressing the implication of PSM in resistance to herbivory across and within species show contrasting results (Moles et al., 2011; Rasmann and Agrawal, 2011; Abdala-Roberts et al., 2016; Moreira et al., 2018), putting the defensive role of PSM under strong criticism (Carmona et al., 2011; Cárdenas et al., 2014). Such controversy might have arisen due to, among other factors, a lack of control for genetic relatedness among species, populations and individuals, and the consideration of broad groups of PSM (e.g., all terpenoids or all phenolic) instead

of small chemical groups or individual PSM. First, in order to explore the relationships between PSM and resistance, the genetic relatedness across or within species must be accounted for to avoid spurious trait-trait or trait-environment relationships. In this sense, phylogenetic correction allows to account for the variation in traits due to common ancestry across related species (Freckleton et al., 2002), whereas the implementation of population structure and relative kinship among individuals allow to account for the variation due genetic relatedness among populations and among individuals within families in a single species (Yu et al., 2006). Second, the role of PSM can be better and precisely described when explored in small chemical groups or at the individual scale, because exploring relationships only among broad classes of compounds with resistance might overstate the implication of some individual PSM from those large groups that are not adaptively configured for defense (Hamilton et al., 2001). Furthermore, previous studies that criticize the adaptive value of PSM against resistance to herbivory (Carmona et al., 2011; Cárdenas et al., 2014) explored not only the effect of a small fraction of the whole diversity of chemical compounds (mostly few types and broad groups of phenolics), but also overlooked the inducibility of PSM as a relevant trait associated with resistance. In this sense, calling into question the role of PSM in resistance to herbivory based only on the effect of concrete and broad groups of PSM and when their plasticity was not measured, may obscure the adaptive value of PSM due to biased overgeneralizations. An exhaustive characterization of the constitutive chemical profile and its inducibility of the individuals would be required to maximize the range of metabolites putatively related to resistance and to identify those that might accurately predict the variation in resistance across and within species.

#### PSM as Resistance Markers Against Herbivory

Plant secondary metabolites have been found to play an essential role in plant resistance to herbivory in many and phylogenetically distant plant species (for instance, see Agrawal et al., 2002; Zhao et al., 2010; Rasmann et al., 2014; Mason et al., 2016); thus, pointing to their adaptive value (Futuyma and Agrawal, 2009; Agrawal, 2010). Developing novel and heritable defensive traits that allow to resist herbivore attacks is of particular relevance in long lived plants like trees (Futuyma and Agrawal, 2009).

In our study, i) we maximized the range of phenotypic variation in the concentration of PSM potentially related to resistance against herbivores by exploring the constitutive allocation to PSM and their inducibility in the stem, ii) minimized possible genetic biases by using clonal replicates, iii) and accounted for relatedness among genotypes within families and populations of P. pinaster during the analyses. Our results showed that inducibility of flavonoids and lignans were strongly related to the inducibility of resistance to the pine weevil when exploring pairwise relationships. Previous studies have suggested that flavonoids and lignans are involved in plant defense (see reviews in Treutter, 2006; Witzell and Martín, 2008), acting as direct defenses against biotic agents (Wallis et al., 2010; Villari et al., 2016) as precursors in cell wall lignification for physical defense (Ralph et al., 2006), and being involved in protection against oxidative stress (Treutter, 2006). Only the concentration of eugenol at the constitutive level was found to be positively related with H. abietis damage (i.e., negatively related to resistance), which contrasts with previous findings where it was found to be an antifeedant for the pine weevil (Borg-Karlson et al., 2006; Eriksson et al., 2008). It is noteworthy that those previous studies were performed under laboratory conditions using purified extracts of eugenol, thus overlooking the interactive effects with other chemical constituents in the PSM blend that would affect insect feeding behavior. An alternative explanation for this finding is that the effect of eugenol was observed to be short-lasting, with decreased anti-feedant properties after several hours of exposure (Borg-Karlson et al., 2006), maybe due to the high volatility of this PSM. It is possible that the very low constitutive concentration of eugenol had an attractant effect rather than repellent in higher concentrations. Surprisingly, we did not detect any relationship between the constitutive concentration of terpene subgroups or individual terpenes, and their inducibility, with resistance to the pine weevil, despite previous results supporting this link (Danielsson et al., 2008; Kannaste et al., 2009; Schiebe et al., 2012; Lundborg et al., 2016). Noteworthy, however, when we explored the whole phenotypic range, including both constitutive and MJ-induced clonal replicates, all total terpene subgroups, total eugenols and total lignans, and many individual terpenes and phenolics were strongly and positively related to resistance against herbivory in the stem (data not shown) as previously found in multiple studies (Porth et al., 2011; Sampedro et al., 2011b; Schiebe et al., 2012; Zas et al., 2014; Lundborg et al., 2016). This result, together with the lack of significant pairwise associations between constitutive PSM (and inducibility) and resistance, suggest that MJ may have artificially forced the pairwise relationships when exploring the whole phenotypic variation. Previous results reporting PSMresistance relationships upon control and induced plants should thus be interpreted with care. Based on our findings, we suggest that the strong effect of MJ in the reduction of weevil damage is product of a large multivariate induction rather than smaller and punctual individual effects on specific PSM.

In fact, when we explored the predictive power of the whole PSM blend in the stem, comprising 98 chemical compounds of different chemical nature, on resistance to H. abietis damage, we found strong evidence that resistance is largely mediated by the concentration of a specific blend of PSM and also by changes in the chemical profile after induction. On the one hand, constitutive resistance was explained at the multivariate level by terpene volatiles (mono- and sesquiterpenes) and many simple phenolics or precursors of complex molecules (for example, eugenols; Dudareva et al., 2013). This indicates that undamaged plants might initially invest more in simple carbon-based compounds in the stem as a cost-saving strategy when the biotic risk of attack is expected to be low (Karban, 2011). Among these PSM, β–phellandrene + limonene and lignan hexoside 1 have been showed to be strong and effective herbivore repellent and antifeedant, respectively, in conifers and other tree species, with increased resistance found in plants with greater concentrations of these PSM (Schiebe et al., 2012;

Villari et al., 2016). While β–phellandrene and limonene coeluted during the chromatographic run in our study, limonene was previously observed to be significantly induced by MJ, whereas β–phellandrene was not affected (Zhao et al., 2010) and seedlings with greater concentrations of limonene in the stem were less attacked by H. abietis (Danielsson et al., 2008). On the other hand, plasticity of PSM revealed an important allocation toward more complex PSM, and changes in the chemical profile of the PSM blend after MJ application. Induction by MJ strongly increased the investment in larger molecules in the stem, such as diterpene resin acids, which are known to be toxic and function as deterrents (Larsson et al., 1986; Tomlin et al., 1996), but also contribute to the sealing of wounded tissue and mechanically impede feeding activity (Ericsson et al., 1988; Phillips and Croteau, 1999). Two phenolics previously found in the pairwise correlations, Unk P5 and lignan hexoside derivative 2, predicted almost half of the variance in increased resistance relative to the whole inducible PSM blend, supporting prior work (Ralph et al., 2006; Villari et al., 2016), and suggesting that these compounds are probably key in resistance, needing thus further investigation.

Despite their general ecological function, the role of fatty acids in conifer resistance is not well understood. We found that the inducibility of oleic acid is associated with susceptibility to pine weevil in maritime pine, confirming that fatty acid composition may affect host suitability for insect pests (Ishangulyyeva et al., 2016). However, due to the important investment in resources that induced defenses usually imply (Cipollini and Heil, 2010), it is possible that the production of induced PSM comes, in part, from the utilization of fatty acids after remobilization and digestion of lipid bodies in the plastids (Bernard-Dagan, 1988).

Interestingly, abietic acid, one of the most well-known diterpene resin acids associated with resistance against insect herbivores (Wagner et al., 1983; Tomlin et al., 1996; Powell and Raffa, 1999), showed an opposite pattern in our study. Recent studies found, however, that gut microbiota of stem chewers and borers can help overcome host defenses, particularly diterpene resin acids, when they are present at relatively low concentrations (Boone et al., 2013; Berasategui et al., 2017). Indeed, such microbiota-mediated processes may confer advantages to such insect guilds (Fäldt, 2000 and references therein), a phenomenon that may well apply to H. abietis and thus needs to be examined further.

#### CONCLUSION

Our study supports a large body of evidence accumulated over decades that host chemical responses are involved in resistance to herbivory at the intraspecific level. However, our study provides the first example of a comprehensive approach in which PSM inducibility in trees, and its role in resistance, are analyzed in a rigorous genetic framework that allows the exclusion of spurious associations due to among-subject relatedness. It follows that multivariate analyses of PSM, conducted in this way, provide more realistic information about the potentially causal relationships between PSM and resistance to herbivory than evaluating single pairwise relationships or using broad groups of compounds.

#### DATA ACCESSIBILITY

The SNP dataset analyzed in this study is publicly available in the Zenodo repository (doi: 10.5281/zenodo.1445313).

# AUTHOR CONTRIBUTIONS

LS and RZ designed the experiment, provided the reagents, performed the sampling, helped with the statistical and GC analyses and the interpretation of the results and improved the different versions of the manuscript. PB provided the infrastructure and know-how for the UHPLC analyses of phenolic compounds and helped with the interpretation of the associated data. CV helped with the UHPLC analyses of phenolic compounds and with the interpretation of the associated data. AB-K provided the infrastructure and background for GC-MS analyses of terpenoid compounds and helped with the interpretation of the associated data. DG developed the SNPs database, constructed the population structure and kinship matrices, and provided expertise on controlling for genetic relatedness. XL-G performed the sampling, all the chemical analyses, most of the statistical analyses, produced the results, wrote the first draft with RZ and LS, and lead on the improvement of subsequent versions of the manuscript. PB, CV, AB-K, and DG contributed equally to the final version of the manuscript.

# FUNDING

This research was supported by the Spanish Government via the MINECO/FEDER grants FENOPIN (AGL2012-40151-C03- 01), FUTURPIN (AGL2015-68274-C03-02-R) and AdapCon (CGL2011-30182-C02-01/02), by the European Union's Horizon 2020 Programme for Research & Innovation (grant id 773383), and by the Xunta de Galicia – GAIN grant IN607A-2016/13. XL-G received financial support from Fundación Pedro Barrié de la Maza, and from the FPI and Short-Term Stay Grant programs (MINECO). We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

#### ACKNOWLEDGMENTS

We thank Dr. Juan Majada (SERIDA, Asturias) and Ricardo Alía (CIFOR-INIA, Madrid) and their corresponding research groups for conceiving, growing and providing the plant material; José Carlos Rodríguez and Santiago Rodal (Viveros Figueirido, NORFOR, Pontevedra, Spain) for plant maintenance and culture, César Cendán and the field personnel of Misión Biológica de Galicia (MBG-CSIC) for their help during plant sampling, and María Luz Pato and Helena Pazó for the help with the plant

sampling, measurements and chemical analyses. Also thanks to Jean-Christophe Cocuron for the support provided with the UHPLC-DAD-MS analyses, and Santiago C. González-Martínez, Carmen García Barriga, and Diana Barba Egido for the help with SNP selection, DNA extraction, and SNP editing, respectively.

#### REFERENCES


#### SUPPLEMENTARY MATERIAL

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





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

Copyright © 2018 López-Goldar, Villari, Bonello, Borg-Karlson, Grivet, Zas and Sampedro. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Gender Dimorphism Does Not Affect Secondary Compound Composition in Juniperus communis After Shoot Cutting in Northern Boreal Forests

Sari Stark1,2 \* and Françoise Martz<sup>2</sup>

<sup>1</sup> Arctic Centre, University of Lapland, Rovaniemi, Finland, <sup>2</sup> Production System Unit, Natural Resources Institute Finland (Luke), Rovaniemi, Finland

Due to a difference in plant resource allocation to reproduction, the males of dioecious plants may be more growth-orientated, whereas females may allocate more resources for synthesizing secondary compounds. This mechanism is considered to cause gender-specific differences in the plant responses to the loss of plant biomass. Here, we tested gender dimorphism in the responses of common juniper (Juniperus communis) to shoot cutting in four juniper populations located in northern boreal forests in Finland. We collected shoots from uncut junipers and from junipers subjected to shoot cutting in the previous year, and analyzed them for their shoot growth as well as phenolic and terpenoid concentrations. There were no differences in foliar phenolic or terpenoid concentrations between the males and the females. Shoot cutting increased phenolic but not terpenoid concentrations, similarly, in both males and females. Our study reveals that the nature of gender dimorphism may differ among species and locations, which should be considered in theories on plant gender dimorphism. Given the similar phenolic and terpene concentrations in both genders, the different sexes in the northern juniper populations might experience equal levels of herbivory. This lack of gender dimorphism in biotic interactions could result from the high need of plant secondary metabolites (PSM) against abiotic stresses, which is typical for juniper at high latitudes.

Keywords: Juniperus communis, phenolics, terpenoids, compensatory growth, secondary metabolism, boreal forest

#### INTRODUCTION

It has long been considered that in dioecious plant species, the growth rates and chemical defenses may differ between female and male individuals (Ågren et al., 1999; Cornelissen and Stiling, 2005; Avila-Sakar and Romanow, 2012). The gender-specific differences in plant carbon allocation likely result from a higher resource investment to reproduction in females, which may lead to a lower growth rate and higher concentrations of plant secondary metabolites (PSM) in females compared with males. Gender dimorphism is also suggested to influence plant responses to the loss of biomass in response to herbivory (Cornelissen and Stiling, 2005). In slow-growing woody plants, the loss of biomass often induces an increase in PSMs, which constitutes an important part of the plant resistance against invertebrate herbivory (Haukioja and Koricheva, 2000) and

#### Edited by:

Bartosz Adamczyk, Natural Resources Institute Finland (Luke), Finland

#### Reviewed by:

Rubén Retuerto, University of Santiago de Compostela, Spain Bruce Kimball, Animal and Plant Health Inspection Service (USDA), United States

> \*Correspondence: Sari Stark sari.stark@ulapland.fi

#### Specialty section:

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

Received: 08 October 2018 Accepted: 10 December 2018 Published: 21 December 2018

#### Citation:

Stark S and Martz F (2018) Gender Dimorphism Does Not Affect Secondary Compound Composition in Juniperus communis After Shoot Cutting in Northern Boreal Forests. Front. Plant Sci. 9:1910. doi: 10.3389/fpls.2018.01910

mammalian browsing (Bryant et al., 1991, 2014). Whether this is a defense reaction or an indirect result of environmental constraints that regulate the trade-off between plant growth and synthesis of secondary phenolic compounds has remained unclear (Tuomi et al., 1990; Jones and Hartley, 1999; Mattson et al., 2005). The capacity of plants to compensate for the lost biomass through increased growth rates (i.e., plant compensatory growth) constitutes another important means of plant tolerance to herbivory (Lehtilä, 2000; Tiffin, 2000; Peinetti et al., 2001; Cromsigt and Kuijper, 2011). Owing to gender dimorphism, males may more commonly respond to the loss of plant biomass by compensatory growth, whereas females respond by increasing plant secondary compounds (Cornelissen and Stiling, 2005).

The secondary compounds in plants exhibit a diverse spectrum of biological functions, and several factors regulate their synthesis (Koricheva et al., 1998; Stamp, 2003; Theis and Lerdau, 2003). Plant secondary metabolites can roughly be gathered in three classes of chemical compounds, namely alkaloids, phenolic compounds and terpenes, with 1000s of compounds in each class (Chomel et al., 2016). PSMs have important roles in plant development, plant–plant and plant–microbe/insect/herbivore interactions; they also govern the mechanisms of allelopathy, influencing intraand interspecific competition between plants. Despite wellestablished importance of gender dimorphism, the majority of studies investigating the effects of environmental stresses on PSMs has focused on plant phenolics (e.g., Cornelissen and Stiling, 2005) while less is known about gender dimorphism of terpenoids. For understanding the ecological consequences of gender-specific differences on plant growth and survival, more studies are needed on the role of sexual dimorphism that would consider the multiple functions of the different classes of PSMs.

Common juniper (Juniperus communis) of the Cupressaceae family is an evergreen dioecious gymnosperm shrub that has one of the widest global distributions of any gymnosperm (Thomas et al., 2007). In Europe, the distribution of juniper ranges from the Mediterranean to the Arctic and within each area, is found in a very wide range of habitats, such as old pastures, forests, and peatlands. Juniper foliage is rich in secondary substances, particularly terpenoids and phenolics (Martz et al., 2009). Previous studies have demonstrated strong gender dimorphism in juniper. For example, juniper populations may be biased toward males, because females have higher mortality rates in resource-limited conditions (Ortiz et al., 2002). The males of Juniperus thurifera start flowering younger (Gauquelin et al., 2002). Different juniper genders also exhibit differing physiological responses to shading from neighboring plants (Verdú et al., 2004). Some studies on junipers, however, have not found gender-specific differences in growth or survival (Marion and Houle, 1996; Ortiz et al., 2002), or these differences have been highly site-specific (DeSoto et al., 2016). Northern boreal juniper populations provide an excellent opportunity to investigate gender dimorphism in plants, because these populations grow in stressful environments and have high concentrations of both phenolics and terpenes in their biomass (Martz et al., 2009).

FIGURE 1 | Junipers in northern boreal forest sites (A: Savukoski, B: Salla) in Finland.

Here, we analyzed gender dimorphism in juniper by analysing the secondary compounds and responses to shoot cutting in northern boreal forests in Finland. In northern boreal forest, junipers commonly form polycormic tall shrubs instead of trees (**Figure 1**). We assumed that if there are differences in the reproductive effort between the genders, male and female individuals should differ in growth rates and the concentrations of PSM as well as in their responses to a loss of shoot biomass. We predicted that (1) concentrations of secondary metabolites should be higher whereas growth rates lower in females than males. As woody plants respond to biomass loss both by compensatory growth and production of secondary metabolites, we further predicted that (2) shoot cutting should induce a stronger increase in the concentrations of PSMs in females than males, whereas the compensatory growth after biomass loss should be higher in males than females.

# MATERIALS AND METHODS

# Study Sites and Selecting Junipers for Study

For our investigation, we selected four sites in northern Finland with dense juniper population in a uniform area (**Table 1**). All sites are located in the northern boreal vegetation zone and can be classified as mid-successional boreal forests with Scots pine (Pinus sylvestris) as the dominant tree species. We compared the role of gender and the loss of biomass in these populations, because the local communities had used these sites for shoot collection for natural products. Common juniper extract is used in some pharmaceutical and technical preparations, cosmetic products, and as a food additive (Mäkitalo et al., 2006; Stark et al., 2010). As part of an applied research project to investigate the recovery rate of junipers from shoot gathering and for creating recommendations for sustainable gathering (see Stark et al., 2010), we searched for several sites with a history of commercial shoot gathering during years 2002–2004. In each site, we randomly chose male and female junipers in a collected area and in an adjacent uncollected area (i.e., ranging from a few 100 m to 5 km; **Table 1**). Individuals with shoot cut in 2004 – generally implemented by cutting with knives or by striping – were identified in collaboration with collectors. In the sites used



The total number of samples includes all individuals (males and females in the uncut and cut areas within each site).

for shoot gathering, the collectors generally gather shoots from each juniper as they go through the area. Consequently, the different juniper individuals are subjected to shoot cutting in a randomized way that does not depend on juniper characteristics such as size or growth rate. When selecting juniper individuals for our study, we also ensured that each of them was cut by a visual identification of cutting marks.

#### Sampling

The chemical quality of needles and shoot growth rates were collected and analyzed in July 2005, 1 year after the previous shoot cutting. In the field, morphological data were recorded on the selected junipers [N = 6 per treatment in Keminmaa (S1), N = 10 in Rovaniemi (S2), N = 8 in Salla (S3), and Savukoski (S4)]. For each shrub, we recorded: (1) height of tree, (2) diameter of the crown of the shrub from two directions, which, together with height, was later used for calculating an approximation of juniper shrub volume, (3) dry weight of the current year shoots, (4) needle coverage (%), (5) percentage of top dead shoots. We collected a sample of 15–40 (depending on the size of the juniper) current-year shoots from each individual. Current-year shoots (including stem) were stored in paper bags and transported to the laboratory within 1–2 days where they were immediately air dried (60◦C 1 day), milled and stored in sealed plastic bag at +4 ◦C in the dark until analysis. We analyzed the mean fresh and dry weight (drying at 60◦C for 48 h) of the shoot samples, indicative of their growth rate, by calculating the number of shoots and analysing the total weight of the dried samples.

#### Chemical Analyses

PSM were extracted from dry powder and analyzed as previously described (Martz et al., 2009). Briefly a one-step extraction protocol was developed to extract terpenoid in hexane and soluble phenolics in methanol 75%. Milled juniper needles (0.5 g) were extracted with 4 ml of methanol:H2O (3:1, v/v) + 4 ml n-hexane by shaking 2 h in the dark at room temperature. Isoborneol (200 µg) was added before shaking as internal standard for terpenoids.

After centrifugation, the upper organic fraction containing terpenoids was removed, concentrated and analyzed by gas chromatography with a HP-5 (30 m × 320 µm × 0.25 µm) column (Agilent Technologies, Santa Clara, CA, United States) using the following conditions: injector 200◦C, flame ionization detector 280◦C, helium as carrier gas (1 ml/min) and the following temperature programme: 50–80◦C at 15◦C/min, 80–100◦C at 3◦C/min, 100–160◦C at 10◦C/min, 160–200◦C at 3 ◦C/min, 200–250◦C at 15◦C/min. Terpenoids were identified and quantified as described previously (Martz et al., 2009). The internal standard isoborneol was used for calculation of the efficiency of recovery of each sample during the whole extraction process. Limonene and α-pinene were used to draw calibration curves for all monoterpenes (average curve used) and β-caryophyllene was used, similarly, for all sesquiterpenoids. Only monoterpenoids and sesquiterpenoids were quantified in this study, and their sum labeled as "terpenoid content."

Soluble phenolics present in the lower aqueous fraction were analyzed by HPLC (Waters, Milford, MA, United States) with a Spherisorb ODS II column (4.6 × 250 mm, particle size 5 µm) column (Waters, Milford, MA, United States) using a binary solvent system (solvent A: 1% ammonium formiate, 10% formic acid in water; solvent B: 1% ammonium formiate, 10% formic acid in methanol). The elution programme was as follow: 0–5 min: 0% B, 5–45 min: 0–100% B, 45–86 min: 100% B; 86–90 min: 100–0% B, 90–120 min: 0% B, at 35◦C and at 1 ml/min. Detection and quantification were made at 280 nm using a UV/visible diode-array detector (Waters PDA 996). The following compounds were used to draw calibration curves for quantification: catechin (for proanthocyanidins and unknown compounds), rutin (for all flavonols) and apigetrin (for all flavones).

#### Calculation and Statistics

The volume of the individual juniper plants was computed using the height and the averaged width (2 values). The raw data have been used in **Figures 1**, **2**; values are the mean of the four sites ± SE of the mean. For statistical analysis, data were transformed when required to meet the assumption of normality. Log10 transformation was applied to shoot biomass, volume of the shrub, total phenolics, % of monoterpenes, Grp3, Grp4, Grp5, U1, α-pinene. Square root transformation was used for total terpenoids and the % of U2. The experiment followed a blockdesign, with site used as blocks. The Linear Mixed Model was used with gender, cutting and their interaction as fixed factors and site as a random factor. Variance component was used as a covariance structure. All statistical testing was conducted with the IBM SPSS Statistics Software (Version 25.0, SPSS Inc., Chicago, IL, United States).

#### RESULTS

#### Composition of Phenolics and Terpenoids in the Common Juniper Shoots

The chemical composition of juniper shoots was similar to that previously reported from samples in Finland (Martz et al., 2009). Thirty-six peaks were identified by HPLC analysis. According to their UV spectra and comparison with authentic standards, the peaks were gathered into seven groups. In order of decreasing abundance, we identified: proanthocyanidins (PAs) (including catechin; Group 1), flavones (mainly apigenin

derivatives; Group 2), flavonols (including quercetin derivatives; Group 3), and unidentified compounds ("others," Group 4). Although unidentified as well, a few compounds with similar spectra were detected and thus were not included in group 4: group 5 (two compounds with λmax of 284 and 311 nm) and group 6 (two compounds with λmax of 280, 296, and 305 or 277, 305, and 328 nm). In addition, two unidentified individual compounds, which were not included in group 4 because of their significant abundance, were detected as follows: U1 (λmax = 266 nm with a shoulder at 300 nm, possibly a neolignan) and U2 (λmax = 273 nm) (Martz et al., 2009; **Table 2** and **Supplementary Table S1**). The terpenoid composition appeared much more variable and the most abundant compounds were of germacrene types (B, D-ol, D) (sesquiterpenoids), α-pinene and several unknown compounds.

#### Effect of Cutting on Growth Rate

A significant increase in shoot biomass of juniper shoots was measured due to previous-year cutting (**Figure 2** and **Table 3**). In the uncut plots, male and female did not show any difference in growth. However, females tended to have a higher growth rate after cutting, although the interaction between cutting and gender was only marginally significant (P = 0.085; **Table 3** and **Supplementary Table S2**). Previous-year cutting did not affect the needle coverage nor the proportion of top dead shoots (**Table 3**). The site had no significant effect on modeling the shoot biomass, volume of the shrub, needle coverage or top dead shoots (Wald Z-test: p = 0.241, 0.256, 0.468, and 0.857, respectively).

### Effect of Cutting on PSM

Previous-year cutting significantly increased the total phenolic content in juniper shoots in both male and female individuals (**Figure 3** and **Table 3** and **Supplementary Table S2**). The phenolic composition was as well affected by cutting with significant increases in flavonols and compounds of Group 6 as well as decreases in abundances of compounds in Groups 4 and 5 (**Table 2**). Flavonols represent a major group of phenolics in juniper shoots, and a detailed analysis of the individual compounds showed that 3 compounds of the four flavonols detected increased in response to cutting (**Supplementary Table S1**). Only hyperin, the most abundant flavonol, was not affected by cutting. Although the abundance of all flavones (Group 2) did not significantly increase in response to cutting, two specific compounds (apigenin derivative 1: λmax 272, 332 nm and U31: λmax 268, 335) were more abundant in shoots from collected junipers. Several compounds in Group 4 showed the same decreasing abundance after cutting. Two other unknown compounds (U11 in Group 5: λmax 291, 316 nm and U27 in Group 6: λmax 277, 305, 328 nm) were as well significantly affected due to cutting (**Supplementary Table S1**). Neither the terpenoid content nor its composition was affected by cutting (**Tables 2**, **3**). The site had no significant effect on modeling the phenolic or terpenoid concentrations (Wald Z-test: p = 0.252 and 0.241, respectively).

#### Effect of Gender Dimorphism on PSMs

Generally only little significant difference was observed between males and females in uncut areas and response to cutting (**Tables 2**, **3** and **Supplementary Table S2**). Only one phenolic (U18: λmax 280, 367 nm, in Group 4) was marginally more abundant in males than females shoots (**Supplementary Table S1**). A gender-specific difference was observed in the abundance of germacrenes with significantly higher abundance in females compared to males, especially in uncut areas (**Table 3**).

# DISCUSSION

Studies on dioecious plants have documented gender-related differences in growth rates and concentrations of secondary metabolites with higher growth but lower PSM concentrations in males compared with females (e.g., Ågren et al., 1999; Cornelissen and Stiling, 2005). Consequently, in dioecious plant species, the intensity of herbivory may be biased toward males (Cornelissen and Stiling, 2005). Studies have found exceptions to these generalizations (Avila-Sakar and Romanow, 2012), and among the species belonging to the dioecious Juniperus genus, previous evidence on gender dimorphism has also been mixed. Contrasting the generalization, Massei et al. (2006) found higher concentrations of terpenes and phenolics in male than female individuals in Mediterranean J. oxycedrus macrocarpa, and McGowan et al. (2004) that the females of prostrate juniper J. communis ssp. nana in northern Scotland

TABLE 2 | Phenolic and terpenoid composition in juniper shoots and statistical significance of cutting, gender, or cutting by gender interaction effects calculated using the Mixed Linear Model.


Concentrations are expressed as % of total soluble phenolic or terpenoid content, respectively (values are mean ± SE, n = 31). Statistical significant effects are indicated by ∗∗ (p < 0.01) or <sup>∗</sup> (p < 0.05). An arrow indicates the direction of the change: CUT compared to UNCUT, and FEMALE compared to MALE.

were more frequently subjected to herbivory than males. No gender-specific difference in juniper growth was observed along an altitudinal gradient in Sierra Nevada despite decreasing reproductive success along this gradient (Ortiz et al., 2002). Here, except for a few individual compounds, we found no effects of gender in phenolic and terpenoid concentrations in J. communis in northern boreal forests. Although, as predicted, shoot cutting significantly increased phenolic concentrations, the level of PSMs increased, similarly, in both males and females. Further, theories on gender dimorphism state that males should be more growth-orientated than females (e.g., Ågren et al., 1999; Cornelissen and Stiling, 2005), but we found marginally (P < 0.10) higher growth after shoot cutting in females compared with males.

The lack of clear gender dimorphism in PSMs and growth in northern juniper populations compared with e.g., boreal deciduous trees (Nissinen et al., 2018; Ruuhola et al., 2018; Zhang et al., 2018) is noteworthy, because junipers in these systems seem to have a particularly high level of chemical defense. Along a large geographical gradient from the southern to northern boreal forests, Martz et al. (2009) found that secondary compound concentrations in the common juniper significantly increased with latitude. Interestingly, these trends were similar in both terpenoids and phenolics, although these compound groups exert largely differing ecological functions. Terpenoids may have a key role in the defense against mammalian herbivory (Mutikainen et al., 2000), and the low palatability of juniper to herbivores is commonly derived from oils found in the needles, cones and wood, dominated by monoterpenes (Thomas et al., 2007). Despite the low palatability, junipers are commonly subjected to herbivory in different habitats (Livingston, 1972; Miller et al., 1982; Fuhlendorf et al., 1997), and in the boreal forests, herbivory on juniper is dominated by shoot consumption by moose (Alces alces) during the winter. Contrasting with terpenoids, phenolics may have a higher importance in the plant protection against invertebrate herbivory (Koricheva et al., 1998), and photooxidative stress (Close and McArthur, 2002). The need for antioxidative compounds in relation to light intensity may increase in conditions of both nutrient deficiency and low temperature (Close and McArthur, 2002), which could explain why their concentrations increase with latitude (Stark et al., 2008; Martz et al., 2009, 2010). Traditional theories also suggested that greater concentrations of phenolics may result from allocation of extra photosynthesized carbon to secondary metabolites as a result of nutrient deficiency (Bryant et al., 1983; Mattson et al., 2005).

Although both phenolic and terpenoid concentrations in juniper seem to show a similar gradient with increasing latitude (Martz et al., 2009), here, we found that the responses of PSM concentrations to shoot cutting varied greatly among the different compounds and compound groups. More specifically, cutting lead to increased concentrations of soluble phenolics and a higher abundance of flavonoids (specific flavones and flavonol glycosides). As we analyzed samples 1 year after the shoot cutting, our analyses depict delayed inducible reactions that take place after a time lag from the loss of biomass (Tuomi et al., 1990). Although previous studies have indicated that delayed inducible



Statistical significant effects are indicated by ∗∗ (p < 0.01) or <sup>∗</sup> (p < 0.05).

reactions in conifers is relatively uncommon when compared with the deciduous trees (Nykänen and Koricheva, 2004), increasing phenolic concentrations in response to shoot cutting agree with previous findings on e.g., Pinus species (Honkanen et al., 1999; Roitto et al., 2008). Flavonols in juniper are mainly quercetin derivatives that exert strong antioxidant activity due to their chemical features (Rice-Evans et al., 1997); thus, shoot cutting led to higher abundances of compounds with a high antioxidant capacities that are more efficient in relation to the carbon cost to their synthesis (sensu Close and McArthur, 2002).

Contrasting with phenolics, we found no overall shoot cutting effect on terpenoids. In line with our investigation, the previous study by Honkanen et al. (1999) also found increasing phenolic concentrations in response to defoliation in boreal coniferous tree Pinus sylvestris L. but no effects on terpenoids, despite they are considered the main class of antiherbivore defensive compounds. Contrasting with this idea, a meta-analysis by Nykänen and Koricheva (2004) concluded that the damage in woody plants quite commonly reduces the concentrations of terpenes. This may be true also in the case of juniper, as reduced yields of essential oils in junipers after a severe browsing damage have previously been found (Markó et al., 2008). Condensed tannins and flavonoids, among many other phenolic compounds are derived from phenylalanine via the phenylpropanoid pathway (Vogt, 2011) and terpenoids are synthesized via acetyl-CoA, pyruvate and glyceraldehyde-3-phosphate in parallel cytosolic of plastid metabolic pathways (Singh and Sharma, 2015). This shows that disturbance of the

general carbon metabolism due to, for example biotic/abiotic stress or compensatory growth will have consequences on the PSM content.

Although the mechanisms underlying our findings remain uncertain, our study adds to previous evidence showing unclear or ambiguous effects of gender on growth, secondary compound concentrations and reproduction in Juniperus species (e.g., Marion and Houle, 1996; Ortiz et al., 2002; Verdú et al., 2004; Massei et al., 2006; DeSoto et al., 2016). For example, aged juniper populations are male-biased, but this bias does not seem to be easily explained by the genderspecific differences in the cost of reproduction (Gauquelin et al., 2002; Ortiz et al., 2002). Earlier reviews have already concluded that there is a need to revise theories predicting how gender dimorphism affects plant performance and responses to environmental stresses (Avila-Sakar and Romanow, 2012; Vega-Frutis et al., 2013). For woody plants, the tolerance of herbivory is a major component of plant resistance, because the probability of herbivory is high due to large size and long life span. However, the recovery potential of these species could be driven by the type of herbivory they commonly experience (Haukioja and Koricheva, 2000). Noteworthy, García et al. (2001) found that frugivory at Juniperus

communis in the Mediterranean mountains depended largely on population characteristics rather than on individual attributes. If herbivory, such as browsing by moose, is commonly centered on locations with numerous and dense juniper populations to provide large food quantity, the chemical quality of the different individuals might not exert a primary role in the food selection of herbivores. Under these conditions, the capacity for regrowth could outweigh the importance of PSMs in herbivory tolerance. Further, the lack of clear gender dimorphism in northern boreal juniper populations could also result from the high need of PSMs at high latitude to protect from abiotic stresses (Martz et al., 2009). Although herbivory might even constitute one of the driving forces behind the evolution of dioecy in plants (Bawa, 1980; Cornelissen and Stiling, 2005; Avila-Sakar and Romanow, 2012; Vega-Frutis et al., 2013), studies on juniper populations have demonstrated that genderrelated differences in growth and resource storage may be a consequence of local adaptation to environmental conditions (DeSoto et al., 2016). Protection against abiotic stresses through the PSMs could have such major significance for plant success under northern conditions that it might override any genderspecific differences in the carbon allocation for synthesizing PSMs.

Ecological factors that limit plant success in each specific conditions could also explain why we detected a marginally higher capacity for compensatory growth in females despite a presumably higher resource cost for reproduction (sensu Cornelissen and Stiling, 2005). Over a large geographical gradient from the Mediterranean to the sub-Arctic vegetation zone across Europe, García et al. (2000) concluded that the juniper population viability in the north may in fact be under less pressure compared with juniper populations at mid-latitudes, because these populations are free from seed predation. Experimental evidence on dioecious plants has suggested that increasing resources may weaken gender-specific differences within plant species, and consequently, the gender-specific differences in phenolic concentrations could be pronounced under high resource limitation (Palumbo et al., 2007). Further, woody plants seem to recover from the loss of biomass better under low than high resource availability, possibly because under these conditions, plants generally grow below their potential

#### REFERENCES


maximum growth rate (Hawkes and Sullivan, 2001). As northern conditions with low nutrient availability, low temperatures and high light likely require high defense through PSMs (Martz et al., 2009), the combination of high need for chemical defense and low pressure on reproduction (García et al., 2000) might direct gender dimorphism from gender-specific differences in defense toward the importance of compensation after the loss of biomass.

## AUTHOR CONTRIBUTIONS

SS designed the experiments and performed field sampling and measurements. FM ran the laboratory analyses and conducted statistical tests. SS and FM jointly wrote the manuscript.

# FUNDING

This work was supported by the European Commission, Regional Development Fund (Projects "Mette", 70025/05, 2005–2006, "Lappi luo", 2105, 533/3560-2008, 2008–2011).

#### ACKNOWLEDGMENTS

We are greatly indebt to the juniper collectors for the information on the shoot collection sites around Finnish Lapland. Furthermore, we sincerely thank Riitta Nielsen for helping with the juniper shoot extractions, and Aarno Niva and Tarja Posio for invaluable input in the juniper growth measurements. We acknowledge Juha Hyvönen and Pasi Rautio for their help in statistical analyses.

#### SUPPLEMENTARY MATERIAL

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




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

The handling Editor declared a shared affiliation, though no other collaboration, with one of the authors FM at the time of review.

Copyright © 2018 Stark and Martz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

,

# Quantification and Localization of Formylated Phloroglucinol Compounds (FPCs) in Eucalyptus Species

Bruna Marques dos Santos 1,2, Juliane F. S. Zibrandtsen1,2†, Disan Gunbilig1,2

# Edited by:

Judy Simon, Universität Konstanz, Germany

#### Reviewed by:

Manoj Gajanan Kulkarni, University of KwaZulu-Natal, South Africa Bartosz Adamczyk, Natural Resources Institute Finland (Luke), Finland Pinarosa Avato, Università degli Studi di Bari, Italy

\*Correspondence:

Elizabeth Heather Jakobsen Neilson en@plen.ku.dk

#### †Present address:

Juliane F. S. Zibrandtsen, Syngenta Ltd, Manchester, United Kingdom Federico Cozzi, BIOMIN Research Center Technopark, Tulln an der Donau, Austria

#### Specialty section:

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

Received: 28 September 2018 Accepted: 05 February 2019 Published: 26 February 2019

#### Citation:

Santos BM, Zibrandtsen JFS, Gunbilig D, Sørensen M, Cozzi F, Boughton BA, Heskes AM and Neilson EHJ (2019) Quantification and Localization of Formylated Phloroglucinol Compounds (FPCs) in Eucalyptus Species. Front. Plant Sci. 10:186. doi: 10.3389/fpls.2019.00186 Elizabeth Heather Jakobsen Neilson1,2,6 \* <sup>1</sup> Section for Plant Biochemistry, Department of Plant and Environmental Sciences, University of Copenhagen, Copenhagen, Denmark, <sup>2</sup> VILLUM Center for Plant Plasticity, Department of Plant and Environmental Sciences, University of Copenhagen,

Mette Sørensen1,2, Federico Cozzi 3†, Berin A. Boughton4,5, Allison Maree Heskes 1,2,6 and

Copenhagen, Denmark, <sup>3</sup> Section for Molecular Plant Biology, Department of Plant and Environmental Sciences, University of Copenhagen, Copenhagen, Denmark, <sup>4</sup> School of BioSciences, University of Melbourne, Parkville, VIC, Australia, <sup>5</sup> Metabolomics Australia, School of BioSciences, University of Melbourne, Parkville, VIC, Australia, <sup>6</sup> Center for Synthetic Biology 'bioSYNergy', Department of Plant and Environmental Sciences, University of Copenhagen, Copenhagen, Denmark

The Eucalyptus genus is a hyper-diverse group of long-lived trees from the Myrtaceae family, consisting of more than 700 species. Eucalyptus are widely distributed across their native Australian landscape and are the most widely planted hardwood forest trees in the world. The ecological and economic success of Eucalyptus trees is due, in part, to their ability to produce a plethora of specialized metabolites, which moderate abiotic and biotic interactions. Formylated phloroglucinol compounds (FPCs) are an important class of specialized metabolites in the Myrtaceae family, particularly abundant in Eucalyptus. FPCs are mono- to tetra-formylated phloroglucinol based derivatives, often with an attached terpene moiety. These compounds provide chemical defense against herbivory and display various bioactivities of pharmaceutical relevance. Despite their ecological and economic importance, and continued improvements into analytical techniques, FPCs have proved challenging to study. Here we present a simple and reliable method for FPCs extraction, identification and quantification by UHPLC-DAD-ESI-Q-TOF-MS/MS. The method was applied to leaf, flower bud, and flower samples of nine different eucalypt species, using a small amount of plant material. Authentic analytical standards were used to provide high resolution mass spectra and fragmentation patterns. A robust method provides opportunities for future investigations into the identification and quantification of FPCs in complex biological samples with high confidence. Furthermore, we present for the first time the tissue-based localization of FPCs in stem, leaf, and flower bud of Eucalyptus species measured by mass spectrometry imaging, providing important information for biosynthetic pathway discovery studies and for understanding the role of those compounds in planta.

Keywords: Corymbia, Eucalyptus, formylated phloroglucinol compounds, macrocarpal, MALDI-mass spectrometry imaging, sideroxylonal, specialized metabolites

**94**

# INTRODUCTION

The Eucalyptus genus (Myrtaceae) is composed of long-lived trees that dominate a vast array of climatic regions across Australia. Eucalyptus trees possess several traits which have made some species economically important, such as fast growth, good wood quality as well as disease and insect resistance (Grattapaglia et al., 2012). They are now grown all over the world across a diverse range of climates, providing renewable resources to produce essential oil, paper, pulp, timber, and other biomaterials. The biogeographical diversity and success of those plants are partly due to their ability to produce a plethora of specialized metabolites, such as terpenes, flavonoids, cyanogenic glucosides, and phloroglucinols. These compounds play an important role in moderating interactions with the environment and combating biotic and abiotic stresses.

Phloroglucinols are an important class of specialized metabolites widely distributed in different natural sources such as plants, marine organisms and microorganisms. In recent years, phloroglucinol derivatives, especially formylated phloroglucinol compounds (FPCs) have been a hot research topic due to their structurally interesting skeletons and important bioactivities including antimicrobial (Faqueti et al., 2015), anticancer (Qin et al., 2016) and antimalarial effects (Bharate et al., 2006). Accordingly, large efforts have been made to characterize new FPCs structures from plants (e.g., Shang et al., 2016; Liu et al., 2018; Qin et al., 2018). FPCs are mono to tetra-formylated phloroglucinol based derivatives often with an attached terpene moiety that occur in the Myrtaceae family, primarily in Eucalyptus species (Eschler et al., 2000). Due to the terpenoid moiety, FPCs have also been called formyl phloroglucinol meroterpenoids (FPMs) (Shang et al., 2016), but herein we refer to them as FPCs.

The simplest FPCs are fully substituted, formylated acylphloroglucinols, such as jensenone (**Figure 1**). The units of jensenone form the basis of dimeric acylphloroglucinols, such as sideroxylonals and grandinal. The formylated acylphloroglucinols can also form adducts with monoand sesqui-terpenes, such as euglobals and macrocarpals (Eschler et al., 2000; Moore et al., 2004b). The first FPC characterized was grandinol, isolated from Eucalyptus grandis and described as a root inhibitor (Crow et al., 1977). Since this first discovery, the Eucalyptus genus has proven to be a rich source of FPCs, with more than 70 compounds characterized in 39 species, predominantly in the subgenus Symphyomyrtus (**Supplementary Table S1**), with macrocarpals and sideroxylonals being the most common groups of FPCs reported in this genus (Moore et al., 2004b).

Sideroxylonals are compounds with a 2-phenylchromane skeleton. The typical structural characteristics of those compounds are the four formyl groups located in the aromatic rings at the positions C-3, C-5, C-3', and C-5', an isobutyl at C-7 and the isopropyl substituent is at C-10'. The differences between individual sideroxylonals appear in the stereochemistry at C-7 and C-10' (Sidana et al., 2010). There are three characterized sideroxylonals from E. sideroxylon and E. melliodora (**Supplementary Table S1**).

The macrocarpals possess an unusual skeleton that can be divided in two domains: one domain comprising a phloroglucinol dialdehyde moiety (common to all macrocarpals) and a second terpenoid domain (Alliot et al., 2013). Macrocarpal A (Murata et al., 1990) was the first macrocarpal to be isolated from E. macrocarpa and has its structure elucidated, showing an globulol skeleton in the terpenoid moiety. Seventeen other macrocarpals have since been isolated from various Eucalyptus species (**Supplementary Table S1**).

More recently, other species from the Myrtaceae family were discovered to be an abundant source of FPCs, with guava (Psidium guajava), Rhodomyrtus spp. and Eugenia spp. possessing 34, 7, and 4 new compounds, respectively (**Supplementary Table S1**). Despite their prevalence in the Myrtaceae family, especially in the important Eucalyptus genus, little is known about the biosynthesis and role of FPCs in planta.

Whilst several biosynthetic pathways for FPCs have been proposed (Yang et al., 2007; Yin et al., 2007; Shao et al., 2012; Gao et al., 2013; Jian et al., 2015; Qin et al., 2016; Tang et al., 2017), no formal studies have been carried out. Localization techniques, such as matrix assisted laser desorption/ionization mass spectrometry imaging (MALDI-MSI), can provide important information of specific cell or tissue localization of specialized metabolites, which can aid pathway discovery and give insight into the metabolite function (Boughton et al., 2016; Andersen et al., 2017; Heskes et al., 2018). To date, the only documented ecological role of FPCs is the strong deterrent effect against marsupial herbivores (Lawler et al., 2000; Moore et al., 2005) and insects (Matsuki et al., 2011). The limited work to investigate the role FPCs in planta is likely hindered by challenges related to identification and quantification. Furthermore, chemical synthesis studies have attempted to produce FPCs (Singh et al., 2010), however this is a difficult and costly process. Consequently, few analytical standards are available in the market, and these are obtained from the isolation and purification from many kilograms of Eucalyptus leaves.

The most widely cited method for extraction and quantitative determination of FPCs was described 15 years ago using HPLC-UV detection at 275 nm (Wallis et al., 2003). This is surprising considering the fast evolution of analytical chemistry techniques and significantly increased performance over this time. Highperformance liquid chromatography, particularly when coupled with tandem mass spectrometry (HPLC-MS/MS), is the most suited method for the analysis of complex mixtures of phenolic components from vegetal origin due to high-sensitivity and specificity (Santos et al., 2013). The fragmentation pattern generated by high resolution tandem mass spectrometry is extremely valuable, because it can guide the identification or differentiation of structurally related compounds (Neilson et al., 2011; Heskes et al., 2012). Recently, HPLC-MS/MS allowed the putative identification of 13 FPCs among 70 phytoconstituents in an E. sideroxylon leaf extract (Okba et al., 2017), however the method was not developed specifically for FPCs, and may underestimate the diversity of this class of specialized metabolites. Here we present a fast, simple, and reliable method for FPCs extraction, detection, and quantification from

Number of each compound corresponds to text and other figures.

complex biological samples using ultra-high-performance liquid chromatography coupled to ultra-violet diode array detection and electrospray ionization quadrupole time-of-flight tandem mass spectrometry (UHPLC-DAD-ESI-Q-TOF-MS/MS) system, which enabled the detection of 49 FPCs in one single leaf extract. Furthermore, to our knowledge we report for the first time the tissue-based localization of FPCs in leaf, stem and flower bud of two Eucalyptus species. Finally, a literature review of all characterized FPCs is presented.

### MATERIALS AND METHODS

#### Plant Material

Samples of the following eight Eucalyptus and one Corymbia species were harvested around Melbourne, VIC, Australia: E. camphora ssp. humeana, E. camaldulensis, E. cladocalyx, E. leucoxylon, E. sideroxylon, E. viminalis, E. yarraensis, and C. ficifolia. All Eucalyptus species belong to the subgenus Symphyomyrtus. All samples were harvested in biological triplicates from adult trees, kept on ice during transport to the laboratory and stored in −80◦C until further analysis.

#### Species Details and Tree Locations

E. camaldulensis (river red gum) is a large white-flowered tree widely distributed across Australia. Samples of leaves were harvested at the University of Melbourne (37◦ 47'21.4"S 144◦ 57'24.9"E) and in Royal Park (37◦ 46'57.7"S 144◦ 56'32.7"E) on the 4th of December 2015 and the 9th of March 2016, respectively.

E. camphora ssp. humeana (mountain swamp gum; henceforth referred to as E. camphora) is a small to mediumsized white flowered tree of south-east Australia. Samples of leaves, flower buds, and flowers were harvested in Buxton at Maroondah Highway by the Igloo Road House (37◦ 25'22.8"S, 145◦ 42'32.4"E) on the 6th of January 2015.

E. cladocalyx (sugar gum) is a small to medium-sized white-flowered tree of South Australia. Samples of leaves, flower buds, and flowers were harvested at the University of Melbourne (37◦ 79'66.0'S, 144◦ 95'96.4'E), in the Royal Botanic Garden (37◦ 83'26.7'S, 144◦ 98'28.7'E) and at Creswick Campus 37◦ 42'25.8'S, 143◦ 89'92.7'E) on the 23rd of January 2015, 12th of February 2015 and the 5th of February 2015, respectively.

E. globulus (blue gum) is a large white-flowered tree endemic to south-east Australia and Tasmania. Leaves of tree individual adult trees were harvested at the University of Melbourne, Creswick Campus (37◦ 42'25.8'S, 143◦ 89'92.7'E) on the 5th of February 2015.

E. leucoxylon (yellow gum) is a small to medium-sized pink, red or yellow-flowered tree of south-east Australia, Kangaroo Island, Flinders Ranges, and Mount Lofty Range. Samples of leaves and flowers were harvested in Royal Park (a pink-flowered individual at 37◦ 47'21.4"S 144◦ 57'24.9"E) and at Monash University, Clayton Campus (a bright yellow-flowered individual at 37◦ 54'24.5"S 145◦ 08'.9"E and a red-flowered individual at 37◦ 54'28.2"S 145◦ 08'24.1"E) on the 27th and 30th of November 2015, respectively.

E. sideroxylon (red ironbark) is a medium to large-sized white, pink, red, or light yellow-flowered tree naturally distributed in east and south-east Australia. Samples of leaves and flowers were harvested in Princes Park (two creamy/white-flowered individuals at 37◦ 47'21.4"S 144◦ 57'24.9"E) and at Monash University, Clayton Campus (a pink-flowered individual at 37◦ 54'26.7"S 145◦ 08'18.0"E) on the 27th and 30th of November 2015, respectively.

E. viminalis (manna gum) is a large white-flowered tree common in New South Wales and Victoria, on Phillip Island and Tasmania. Samples of leaves and flowers were harvested on Phillip Island (38◦ 29'04.4"S 145◦ 15'56.2"E) on the 16th of February 2016.

E. yarraensis (Yarra gum) is a medium-sized white-flowered tree of south-eastern Australia. Samples of adult leaves and flower buds were harvested in Coldstream at Maroondah Highway (37◦ 73'46.2'S, 145◦ 37'66.4'E) on the 6th of January 2015.

Corymbia ficifolia (red-flowering gum) is a small to mediumsized red-flowered tree of south-west of Western Australia. Samples of leaves and flowers were harvested at the University of Melbourne (37◦ 47'46.7"S 144◦ 57'33.0"E) and at Monash University, Clayton Campus (37◦ 54'50.3"S 145◦ 07'48.3"E) on the 12th and 13th of January 2016, respectively.

#### Solvents and Chemicals

Acetonitrile and methanol of HPLC-grade with purity ≥99.9% were purchased from Sigma Aldrich (Schnelldorf, Germany). PierceTM LC-MS grade formic acid was purchased from Thermo Fisher Scientific. Sodium formate for internal mass calibration of each analysis was purchased from Sigma Aldrich (Steinheim, Germany). FPCs analytical standards were purchased from BOC Sciences (NY, USA) for macrocarpal J (1), macrocarpal N (2), macrocarpal A (3), macrocarpal L (4), and sideroxylonal A (5) shown in **Figure 1**. Care must be taken when sourcing commercial FPCs standards. In our experience, so called authentic standards provided by a prominent supplier were determined in our laboratory by standard analytical chemistry techniques (1H, <sup>13</sup>C NMR, and MS/MS) to be other unrelated chemical compounds and not suitable for the purpose.

# Extraction Procedure

Frozen plant material was ground in liquid nitrogen with mortar and pestle and ∼100 mg was weighed in a 2 mL screw-lid Eppendorf tube. The plant material was boiled in 500 µL of extraction solvent (85% methanol and 0.1% formic acid in water) for 5 min and immediately cooled on ice for 10–15 min. Subsequently, the samples were centrifuged (15,000 × g, 5 min at 4 ◦C) and the supernatant transferred into a new brown glass vial and stored at −80◦C until further analysis. The plant material was oven-dried for 24 h at 70◦C and the dry weight was recorded. Before the LC-MS analysis, samples were diluted 5 times in 0.1% formic acid in water and filtered through a membrane filter (0.45µm, Merck Millipore).

# Preparation of Standard Curve

Authentic FPCs standards were freshly dissolved in methanol (1 mg /mL) in brown glass vials and diluted in water to different concentrations, depending on the range of the calibration curve. The solubilized FPCs standards had a stable signal intensity for maximum 3 days. Therefore, all standards were freshly prepared and analyzed immediately.

# UHPLC-DAD- ESI-Q-TOF-MS/MS

Chromatographic separation was performed on a Dionex Ultimate 3000RS UHPLC (Thermo Fisher Scientific) system with DAD detector, cooling auto-sampler at 10◦C and column oven at 40◦C. For data acquisition and processing, Compass DataAnalysis software (version 4.3, Bruker Daltonics) was used. Five microliter of the extracts were separated using a Phenomenex Kinetex <sup>R</sup> column (150 × 2.1 mm) packed with 1.7µm C18 material with pore size of 100Å. Extracts were eluted under gradient conditions at constant flow rate of 0.3 mL/min as follows: 50% solvent A (0.05% formic acid in water) and 50% solvent B (0.05% formic acid in acetonitrile) linearly increasing to 100% solvent B in 20 min, hold for 13 min, and finally decreasing to initial conditions and re-equilibrating the column for 10 min.

The UHPLC system was coupled to compactTM (Bruker Daltonics) mass spectrometer with an electrospray ionization source. Eluted compounds were detected from m/z 50–1200 in negative ion mode under the following instrument settings: nebulizer gas, nitrogen at 2 bar; dry gas, nitrogen at 8 L/min and 250◦C; capillary, 4,500 V; in-source CID energy, 0 V; hexapole RF, 100 Vpp; quadrupole ion energy, 4 eV; collision gas, nitrogen; MS/MS bbCID collision energy, 20 eV; collision RF 500 Vpp; transfer time, 100 µs; prepulse storage, 5 µs; spectra rate, 1 Hz; number of precursor ions, 3; active exclusion, 3 spectra, exclusion release, 12 s. Direct infusion of 10 mM sodium formate in the mass spectrometer at the beginning of the LC-MS run using a syringe pump setup with flow rate of 3 µL/min was performed to allow internal mass calibration of each analysis. All samples and standards were run in full scan mode for accurate quantification, and subsequently ran in autoMS/MS mode to acquire fragmentation information for the correct identification of FPCs.

Raw data was processed using the software DataAnalysis 4.2 from Bruker Daltonics. Automatic internal calibration was applied for each sample using the negative ions of the sodium formate cluster, in high precision calibration mode, search range m/z of 0.05 and intensity threshold of 1,000. Extracted ion chromatograms (EIC) for specific [M-H]<sup>−</sup> ions were used to locate compounds. The identification of FPCs was based on the UV absorbance at 275 nm, measured [M-H]<sup>−</sup> with < ± 2 ppm error when compared to the accurate [M-H]<sup>−</sup> and the presence of the diagnostic fragment ions m/z 249, 207, and 181 observed in the authentic analytical standards. Calibration curves covered the range 0.5–75µM (**Supplementary Table S3**) and were used for absolute quantification of the five compounds corresponding to the analytical standards. Other macrocarpals were relatively quantified using the calibration curve for macrocarpal A (3). Other sideroxylonals were relatively quantified using the calibration curve for sideroxylonal A (5). The sum of the peak areas at m/z 485.2557, 489.2858, 471.2752, and 499.1610 were used for the relative quantification of total FPCs.

#### Statistical Analysis

Total FPCs concentration was statistically analyzed in SigmaPlot for Windows version 13.0 using One-Way ANOVA. The dataset was square root transformed to meet the normality and equal variance assumptions.

#### Sample Preparation for MALDI-MSI

To locate FPCs in the different tissues, matrix-assisted laser desorption/ionization mass spectrometry imaging (MALDI-MSI) was carried out according to the procedure described in Schmidt et al. (2018). In short, the plant material was embedded in denatured albumin (boiled egg-white) and gently frozen over a bath of liquid nitrogen. Albumin embedded tissue was cryo-sectioned (30µm) (Leica CM3050 S cryostat, Leica Microsystems, Wetzler, Germany), mounted onto double-sided carbon tape attached to glass slides and freeze-dried overnight. 2,5-Dihydroxybenzoic acid (DHB) was used as matrix and sublimed onto tissue sections using a custom-built sublimation apparatus. MS images were acquired in positive ion mode using a

Bruker SolariX 7-Tesla Hybrid ESI/MALDI-FT-ICR-MS (Bruker Daltonics, Bremen, Germany).

A mass range of 100–2,000 m/z was employed with the instrument set to broadband mode with a time domain for acquisition of 2M providing an estimated resolving power of 130,000 at m/z 400. The laser was set to 30–45% power using the minimum spot size (laser spot size ∼10 × 15µm), smartwalk selected and random raster enabled resulting in ablation spots of ∼35–40µm in diameter, a total 1,500–2,000 shots were fired per spectra (pixel) at a frequency of 2 kHz within a 40 × 40µm array. Optical images of tissue sections were acquired using an Epson Photosmart 4,480 flatbed scanner using a minimum setting of 4,800 d.p.i. Optical images are presented as an inverted image, generated in Adobe Photoshop CS2 (Adobe). The data were analyzed using Compass FlexImaging 4.1 (Bruker Daltonics), data was normalized using RMS normalization and

individual ion images were scaled as a percentage of maximum signal intensity to enhance visualization. Potassium adducts of sideroxylonals [[M + K]<sup>+</sup> m/z 539.1314, 0.01 ppm error] and other FPCs [[M + K]<sup>+</sup> m/z 493.2351, −0.04 ppm error, m/z 511.2456, 0.09 ppm error, and m/z 525.2249, 0.02 ppm error]. A full list of FPCs corresponding to these specific masses is presented in the **Supplementary Table S1**.

#### RESULTS

#### Detection and Quantification of FPCs

Formylated phloroglucinol compounds (FPCs) were successfully extracted, detected, and quantified from different Eucalyptus tissues by UHPLC-DAD-ESI-Q-TOF-MS/MS. Representative UV chromatogram and extracted ion chromatograms from E. globulus leaves are shown in **Figure 2**. The maximal absorbance UV wavelength of 275 nm is commonly used to detect FPCs (Eschler et al., 2000), therefore in **Figure 2A** it is possible to see a peak for each of the five analytical standards highlighted, but also many other peaks with the same m/z value. In general, FPCs require a high concentration of organic solvent to be eluted from C18 columns, and in this study the first FPC peak elutes after ∼7.5 min, at 68% solvent B showing m/z value of 489.2861. The last compound eluted is sideroxylonal A, at ∼23.3 min, with 100% solvent B.

Authentic analytical standards of compounds **1-5** were subjected to MS/MS and characteristic fragmentation patterns identified (**Figure 3**). For the macrocarpals, compounds **1–4**, the most common and intense fragment is the ion m/z 207, which was subsequently used as a diagnostic fragment ion for this group of FPCs. For compound **5**, the ions m/z 249 and 181 were the most intense and subsequently used as diagnostic fragments for sideroxylonals.

A total of 49 peaks were detected with characteristic FPCs features (**Table 1**). Identification was based upon the combination of UV absorbance at 275 nm, comparison of TABLE 1 | FPCs detected in methanol extracts of Eucalyptus samples (leaf, flower bud, and flower) analyzed by UHPLC-DAD-ESI-Q-TOF-MS/MS.


Eca, E. camaldulensis; Ec, E. camphora; Eg, E. globulus; El, E. leucoxylon; Es, E. sideroxylon; Ev, E. viminalis; Ey, E. yarraensis; n.d., not detected; L, leaf; FB, flower bud; F, flower. N.B. Only leaves were analyzed for E. camaldulensis and E. globulus.

measured [M-H]<sup>−</sup> with theoretical [M-H]<sup>−</sup> within < ± 2 ppm error and the presence of the diagnostic fragments for FPCs. Among all peaks detected, 32 were phloroglucinolsesquiterpene-coupled compounds, corresponding to the m/z value of 471 such as compounds **1** and **2**, m/z 489 such as compound **3**, and m/z 485 such as compound **4** (**Figure 1**). Four peaks corresponded to phloroglucinol dimers with an m/z value of 499, such as compound **5**, jensenal, grandinal, and sideroxylonal B and C. The peak with m/z 451.2497 was putatively identified as eucalyptal D (C28H36O5), as there are no other FPCs described that would correspond to this m/z value. Finally, nine peaks with m/z 453 and three peaks with m/z 467, also corresponded to phloroglucinol-sesquiterpene-coupled compounds (**Supplementary Table S1).**

From all species analyzed, E. camphora and E. globulus had the highest concentration of total FPCs in leaves, with 65 and 41 mg g <sup>−</sup><sup>1</sup> DW, respectively (**Figure 4**, **Supplementary Table S2**). Eucalyptus camphora also had high concentration of FPCs in flower buds and flowers, with 13 and 12 mg g−<sup>1</sup> DW, respectively. Interestingly, three Eucalyptus species showed a tendency to accumulate more FPCs in flowers compared to the leaves. Eucalyptus leucoxylon, E. sideroxylon, and E. viminalis contained ∼40, 5, and 3 times more total FPCs in the flowers compared to leaves, respectively **Figure 4**, **Supplementary Table S2**. Eucalyptus yarraensis presented very low amounts of FPCs in leaves and flower buds, and it is the only species that does not contain any sideroxylonals. Eucalyptus cladocalyx and C. ficifolia did not show any traces of this class of specialized metabolites in the tissues analyzed.

#### Localization of FPCs in Eucalyptus Tissues by MALDI-MSI

The spatial distribution of FPCs in flower bud, leaf and seedling stem from two Eucalyptus species were investigated by matrixassisted laser desorption ionization-mass spectrometry imaging (MALDI-MSI; **Figures 5**, **6**). MALDI-MSI revealed that FPCs were associated with the subdermal embedded glands in all tissues analyzed. In addition, FPCs were located in the epidermal layer adjacent to some embedded glands in the leaf and present within the stamens of the E. camphora flower bud. Ions corresponding to sideroxylonals were not detected.

#### Detailed Review of FPC Prevalence and Characterization

Due to their apparent prevalence and importance in Eucalyptus and other Myrtaceae species, we conducted a detailed review of all characterized FPCs (**Supplementary Table S1**). Since the characterization of grandinol from E. grandis in 1977 (Crow et al., 1977), an impressive number of 162 FPCs have been described with more than 80 compounds characterized in 39 Eucalyptus species, predominantly in the subgenus Symphyomyrtus. More recently, other plants from the Myrtaceae family where discovered to be abundant sources of FPCs, including 34 new FPCs identified from guava (Psidium guajava) (for example psiguajadials A-K by (Tang et al., 2017); for more references see **Supplementary Table S1**), seven identified from Rhodomyrtus spp. and four identified from Eugenia spp. The structural differences between the most common FPCs found in Eucalyptus species vs. the new compounds identified in other genus inside the Myrtaceae family are illustrated in the **Supplementary Figure S4**.

## DISCUSSION

#### An Improved Method to Detect and Quantify FPCs

The most cited method for extraction and quantitative determination of FPCs was described 15 years ago (Wallis et al., 2003) using only HPLC-UV detection at 275 nm, with some modifications for sideroxylonals (Wallis and Foley, 2005). Near infrared reflectance spectroscopy (NIRS) has also been used to quickly predict sideroxylonal concentration in leaves but it is non-specific. In this method, the light reflected when a sample is exposed to light in the near-infrared spectrum depicts the chemical bonds in the sample. These spectra have to be calibrated against reference values obtained by analyzing a portion of the samples using traditional analytical chemistry methods such as HPLC (Foley et al., 1998; Wallis et al., 2002). It is evident that analytical limitations have restricted previous investigations of FPCs to the quantification of sideroxylonals only. Therefore, the qualitative and quantitative variation in this very diverse group of specialized metabolites remains largely unknown and the total concentration of those compounds might have been underestimated.

Here we present a fast, simple, and reliable method for FPCs extraction, detection, and quantification from complex biological samples using an UHPLC-DAD-ESI-Q-TOF-MS/MS system. Unlike the previous methods, the sample extraction presented here is much simpler and faster and is suited for small amounts of plant material, such as a single leaf. This allows future studies to investigate the qualitative and quantitative variation of FPCs e.g., in leaves of an individual tree and during seedling development.

In this study, we show that FPCs require high amounts of organic solvent to be eluted from the C18 column, in agreement with previous reports (Eyles et al., 2003; Moore et al., 2004b). Sideroxylonal A, for example, can only be eluted after 3 min at 100% solvent B. This suggests that most of the typical methods used in untargeted metabolomics studies using reverse-phase columns would miss this compound. The chromatograms also show many peaks corresponding to the calculated accurate mass of FPCs, and the identification of those peaks as FPCs was possible due to the typical UV absorbance at 275 nm and the comparison with authentic analytical standards accurate mass and fragmentation pattern. This shows the high diversity of FPCs in complex Eucalyptus extracts.

A specific fragmentation pattern for macrocarpals and sideroxylonals, the most abundant groups of FPCs reported in Eucalyptus samples (Moore et al., 2004b), has been found in this study. The characteristic fragment ion at m/z 249 from compound **5** results from the retro-Diels-Alder cleavage of the molecular ion (Soliman et al., 2014). It is the most abundant

and thus the diagnostic fragment ion for sideroxylonals, followed by a less abundant peak of m/z 181, suggested to be the diformyl phloroglucinol moiety by Chenavas et al. (2015). For macrocarpals, the diagnostic fragment ion detected in the compounds **1-4** is m/z 207, in agreement with previous reports (Eyles et al., 2003; Okba et al., 2017), and we here suggest it to be the isopentyl diformyl phloroglucinol moiety + C2H2. The fragmentation pattern and the high-resolution mass spectra shown here have great value for future studies allowing the confident identification of FPCs, which have probably been overlooked in many other plant species due to the challenges of detecting them. Our method provides a detection limit as low as 0.5µM of FPCs in a very complex plant extract (standard curves available in excel file—**Supplementary Material**).

#### Qualitative and Quantitative FPCs Variation in Eucalyptus

Here we present the comprehensive list of FPCs identified in Eucalyptus tissue. Our results show that FPCs are more abundant and diverse than previously reported. For example, 49 FPCs peaks were identified in a single E. globulus leaf extract, undoubtedly representing some new uncharacterized compounds. Previously, Okba et al. (2017) putatively identified 13 FPCs in a bulk leaf extract of E. sideroxylon using HPLC-MS/MS. Using the current method, we report 34 FPCs identified in leaves of the same species. This significant increase in the number of FPCs detected reflect how method optimization increases detection sensitivity. The phloroglucinol-sesquiterpene compounds with m/z 471 and 485 are the most frequently detected peaks, and potentially corresponding to the m/z value of many different isomeric FPCs, such as macrocarpals.

The presence and abundance of FPCs is a highly variable trait. There is considerable intraspecific variation in their concentration as reported by Lawler et al. (2000). In the present study, E. globulus had one of the highest concentrations of total FPCs, dominated by macrocarpals (**Figure 4**), while in the study by Eschler et al. (2000) macrocarpals and sideroxylonals were not detected in this species. Even among the three biological replicates of E. sideroxylon used in the present study, we identified an individual that had no FPCs, while the two others possessed significant concentrations of sideroxylonals and macrocarpals. A striking example of variation in FPC abundance is demonstrated by the mosaic trees E. melliodora and E. sideroxylon. These individual trees had higher FPC

stem were prepared for matrix-assisted laser desorption ionization-mass spectrometry imaging (MALDI-MSI). Corresponding ion maps of different FPCs are shown: yellow; m/z 493.2351; [M+K]+, green m/z 511.2449; [M+K]+, pink m/z 525.2249; [M+K]+, and a combination of all selected ions overlaid. FPCs are associated to the embedded glands, adjacent to the epidermis of the flower bud, to the embedded glands beneath the nectary, and to the stamens. For the stems, FPCs are localized to embedded glands within the cortex. All images are root mean square (RMS) normalized, with internal scaling of 50% for the ions m/z 493.2351 and 511.2449, and 25% for m/z 525.2249. Scale bar = 500µM. C, cortex; P, pith; V, vascular tissue.

concentration in one single branch, conferring resistance to insect herbivory (Padovan et al., 2012).

Eucalyptus camphora and E. globulus presented high concentrations of total FPCs in expanded leaves, with 65 and 41 mg g−<sup>1</sup> DW, respectively. These concentrations are in a similar range to previous reports. For example, the concentration of sideroxylonals have been reported to reach up to 52 mg g−<sup>1</sup> DW in E. melliodora (Wallis et al., 2002) and up to 100 mg g−<sup>1</sup> DW in E. loxophleba ssp. lissophloia (Wallis and Foley, 2005).

Prior to this study, FPCs had only been detected in the reproductive tissue on one other occasion, with sideroxylonal C identified from flowers of E. albens (Neve et al., 1999). Here we show the presence of FPCs in flowers of four additional Eucalyptus species, and for three of them (E. leucoxylon, E. viminalis, and E. sideroxylon), the concentration of FPCs shows a tendency to be higher in the flowers compared to mature leaves. Interestingly, tissue-specific qualitative variation was observed, with some FPCs only detected in the leaves, or only in the reproductive tissue (**Table 1**). These observations could suggest tissue-specific roles for different FPCs.

The absence of FPCs in E. cladocalyx and C. ficifolia may suggest that these species do not synthesize these compounds. It may also be possible, however, that the FPCs are synthesized in other tissues, or at a different ontogenetic stage. For example, FPCs have been detected in the wood of E. globulus and E. nitens, following a wounding event (Eyles et al., 2003) and in fine roots of 12 species from the subgenera Symphyomyrtus and Eucalyptus (Senior et al., 2016). Qualitative polymorphism has been recorded for FPCs (Padovan et al., 2012) and other specialized metabolites classes in Eucalyptus e.g. cyanogenic glucosides (Neilson et al., 2011). Interestingly, E. cladocalyx synthesizes high levels of the cyanogenic glucoside prunasin (Hansen et al., 2018), which could suggest an alternative chemical defense strategy for this species. Corymbia spp. have been reported to contain FPCs previously (Eschler et al., 2000). The absence of FPCs from C. ficifolia leaves may be due to the lack of embedded glands (Brooker and Nicolle, 2013). Given that we have shown that FPCs are localized to the embedded glands, it could suggest that C. ficifolia does not possess the structural cells to biosynthesize or store these metabolites.

#### Localization and Possible Roles of FPCs in Eucalyptus

To provide further insights into FPC metabolism, we conducted MALDI-MSI on the two highest producing species, E. camphora and E. globulus. MALDI-MSI can greatly aid in the understanding of physiological roles and provide valuable information for biosynthetic pathway elucidation (Boughton et al., 2016; Heskes et al., 2018). Here we show that FPCs are associated with the subdermal embedded glands in all the tissues analyzed. Localization of other Eucalyptus specialized metabolites such as mono- and sesqui-terpenes, flavones and monoterpene glucose esters have been reported to be present in embedded glands (Goodger et al., 2010, 2016). Localization of specialized metabolites to storage structures is widely prevalent across the plant kingdom, largely as a way to prevent autotoxicity (Knudsen et al., 2018). Upon tissue disruption, such as that caused by herbivory, the compounds are released to provide a layer of chemical defense. Accordingly, the localization of FPCs to the embedded gland supports their role as chemical defense compounds. In the case of Eucalyptus, the content of embedded glands was shown to be in two phases leading to the suggestion that monoterpene glucose esters are spatially separated from the mono- and sesqui-terpenes within the lumen (Heskes et al., 2012). In the results presented here, spatial organization of different FPCs within embedded glands could also be suggested, as the ion m/z 525.2249 typically locates to the center of the gland lumen, whilst the ions m/z 493.2351 and m/z 511.2449 locate to the outer perimeter of the glands. Higher resolution imaging of these embedded glands could provide further insight into more complex spatial separation and localization of FPCs. Interestingly, FPCs also co-localized to the leaf epidermal layer of several glands. It could be speculated that this may provide a metabolite highway by which the lipophilic mono- and sesqui-terpenes could be emitted from the embedded gland to the atmosphere.

The presence of FPCs in the stamens of flower buds could suggest a possible role in defense against florivores, or other roles such as attracting specific pollinators. Given the presence of oil glands and FPCs directly below the nectary, it could be speculated that the FPCs may also be present in Eucalyptus nectar to either attract pollinators or deter nectar robbers. To our knowledge, no study has investigated the concentration of different specialized metabolites in Eucalyptus nectar, and this is an intriguing avenue of research to pursue.

Overall, the role FPCs play in Eucalyptus leaves has been strongly linked to chemical defense against marsupial folivores. For example, the concentration of total FPCs was the most important variable determining feeding by marsupial folivores on Eucalyptus species (Lawler et al., 2000; Moore et al., 2005; Jensen et al., 2015), and were implicated in habitat patchiness in Australian forests (Lawler et al., 2000). Interestingly, the highest concentrations of total FPCs in this study were found in four koala-preferred species: E. camphora, E. globulus, E. camaldulensis, and E. viminalis (Moore et al., 2004a; Higgins et al., 2011). This could reflect the ongoing "arms race" between plant and herbivore, where koalas can tolerate significant levels of these compounds. Indeed, koalas can tolerate up to 50 mg of sideroxylonal per gram of dry matter (equivalent to 5%) (Moore et al., 2004a) and show specific behavioral traits when dealing with eucalypt diets of varying toxicity (Marsh et al., 2007). The metabolism and detoxification of FPCs by the koala is currently unknown. FPCs can also determine patterns of damage by Christmas beetles in Eucalyptus plantations (Steinbauer and Matsuki, 2004; Andrew et al., 2007; Matsuki et al., 2011). Andrew et al. (2007) suggests that the potential for evolution by natural selection of sideroxylonal concentrations is not strongly constrained by growth costs and that both growth and defense traits can be successfully incorporated into breeding programs for plantation trees. Therefore, FPCs are important not only in natural forest ecology, but also for commercial Eucalyptus plantations.

# FPC Prevalence, Characterization, and Future Research Directions

In recent years, FPCs have become a hot research topic due to their structurally interesting skeletons and important bioactivities including antimicrobial (Sidana et al., 2010; Faqueti et al., 2015), anticancer (Qin et al., 2016), and antimalarial effects (Bharate et al., 2006). The wide reported bioactivities of FPCs were reviewed in depth by Brezáni and Šmejkal (2013).

FPCs were once believed to be exclusive to Eucalyptus species, but according to our literature review, 91 new compounds were discovered during the past 10 years, with 42 FPCs isolated from other species in the Myrtaceae family (**Supplementary Table S1**). Their structural differences and similarities are illustrated in the **Supplementary Figure S4**. The FPCs found in Eugenia umbelliflora contain, in the phloroglucinol part, both an aldehyde and a butyroyl chain (e.g., Eugenial C), whereas the analogous compound Macrocarpal G widely occurring in the Eucalyptus species contain two aldehyde substituents (Faqueti et al., 2015). Most FPCs found in Rhodomyrtus plants have the skeleton of the phloroglucinol coupled to an eudesmane (sesquiterpene) moiety (Shou et al., 2014) (e.g., Rhodomyrtal A), showing similarity to grandinol, the very first FPC characterized which was isolated from Eucalyptus grandis. Psidium meroterpenoids are a subgroup of FPCs, which are exclusively reported from the species Psidium guajava. Structurally, they are characterized by the presences of 3,5-diformyl-benzyl phloroglucinol moiety and a dihydropyran ring junction (Tang et al., 2017), e.g., Psiguajadial C, F, and K in the **Supplementary Figure S4**.

Plant species from the Myrtaceae family have been used as medicinal plants by native populations in Brazil, South Africa, Australia, and China for centuries. This could be partially explained by the recent discovery of FPCs in Psidium and Eugenia species and the broad range of bioactivity demonstrated by those compounds (Shao et al., 2010; Faqueti et al., 2012, 2015). Considering the challenges to detect and quantify those compounds, we hypothesize that FPCs are more widespread in the plant kingdom and have been overlooked so far. The improved method presented here will be invaluable for the continued exploration of FPCs prevalence.

It is evident that FPCs have important ecological functions based on their role in deterring herbivores. However, many aspects of FPC prevalence and their potential multifunctional roles across different ecosystems remains elusive. A simple, reliable method for FPC detection and quantification will be of great benefit for impending studies. Future work into FPC biosynthesis will also aid the understanding of the role and regulation of these compounds across the Eucalyptus genus and beyond. Furthermore, with atmospheric CO<sup>2</sup> and temperatures increasing at unprecedented rates, it is important to investigate how different climatic factors will influence FPC concentrations, and the impact this may have on ecological and commercial systems. Previous studies have shown that environmental stresses related to climate change—such as elevated atmospheric CO<sup>2</sup> and ozone—can increase specialized metabolite concentration in Eucalyptus (Gleadow et al., 1998; Kanagendran et al., 2018). Our improved method to detect and quantify FPCs will allow future studies to investigate how those compounds respond to environmental changes. An increase in the concentration of FPCs in Eucalyptus leaves will have significant implications in terms of the palatability of foliage and defense against herbivores, directly affecting the food chain.

# CONCLUDING REMARKS

Here we demonstrate a fast, simple, and reliable extraction method suitable for FPCs in combination with UHPLC, DAD and MS/MS analysis. For the first time, authentic analytical standards were used to provide high resolution mass spectra and fragmentation patterns, which have great value for the correct identification of those compounds in complex biological samples. This class of specialized metabolites has been overlooked likely due to the challenges related to their identification and quantification. Therefore, we believe that the results presented here will allow future studies to identify FPCs with high accuracy, which is essential for better understanding the role of those compounds in planta, particularly through ontogenetic development and in response to biotic and abiotic stresses. The impressively large amount of FPCs detected in flowers opens the debate for different roles of these compounds. FPCs are also considered very bioactive molecules with potential applications in the pharmaceutical industry. The tissue-based localization presented here provides important information on the spatial distribution of FPCs in Eucalyptus, and, can directly contribute to pathway discovery studies by providing target tissues for gene expression analyses such as transcriptomics.

# DATA AVAILABILITY

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

# AUTHOR CONTRIBUTIONS

EHJN conceived the study and research plans. BMS, DG, MS, JFSZ, and FC performed the experiments. EHJN, AMH, and BAB supervised experimental work. BAB provided technical assistance. BMS and EHJK analyzed the data. BMS wrote the manuscript with contributions from MS and EHJN. All authors have reviewed and approved the final manuscript.

# FUNDING

This work was financially supported by a Young Investigator Program fellowship from the VILLUM Foundation (Project No. 13167), by a Danish Independent Research Council Sapere Aude Research Talent Post-Doctoral Stipend (Grant No. 6111- 00379B) awarded to EHJN and by the VILLUM Research Center for Plant Plasticity awarded to Prof. Birger Lindberg Møller (Project number 7523).

#### ACKNOWLEDGMENTS

We would like to thank David Pattison and Eleni Lazardi for their excellent technical help in the PLEN Core Metabolomics Platform. We thank Mohammed Saddik Motawia for checking the chemical structures of FPCs and conducting the NMR on the incorrect standards. We thank David Vernon for his generous support for wildlife and early career researchers. We thank Gustavo Avelar Molina for the design of the plots on **Figure 4**. We thank Metabolomics Australia (School of BioSciences, University of Melbourne), an NCRIS initiative

#### REFERENCES


under Bioplatforms Australia, for conducting MALDI-MSI. Finally, we would like to acknowledge all the researchers that have contributed to the characterization of new formylated phloroglucinol compounds, and we apologize if any reference was unintentionally left out.

#### SUPPLEMENTARY MATERIAL

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


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

Copyright © 2019 Santos, Zibrandtsen, Gunbilig, Sørensen, Cozzi, Boughton, Heskes and Neilson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.