Phylogenetic analysis of secondary metabolites in a plant community provides evidence for trade-offs between biotic and abiotic stress tolerance

Plants' responses to conflicting stresses may result in physiological trade-offs due to the inter-dependent and costly nature of physiological investments. Physiological tradeoffs have been proved within species, but to what extent these trade-offs are the result of phylogenetic constraints remains poorly known. Environmental stresses can vary widely in different biomes, and therefore assessing physiological tradeoffs across species must account for this variation. One way of doing so is to assess it within a community, where the co-occurring species have faced a shared combination of filters to establish. Considering a representative sample of species in a single community, we use a macroevolutionary approach to test the hypothesis that plant physiological trade-offs are evolutionarily conserved within this community (i.e., closely-related species tend to solve the trade-offs similarly). We analyze the content of five metabolites in thirty co-occurring plant species, capturing their range of contrasting exposures to abiotic and biotic stresses (growing solitary and in vegetation patches). Our results support that species investment in response to abiotic stress (i.e., proline and abscisic acid content) is traded off against their investment to face biotic stress (i.e., jasmonic acid and salicylic acid), shown by the contrasting loadings of these two groups of metabolites in the first axes of a principal component analysis (PCA). In addition, the metabolic strategies observed in this community are evolutionarily conserved, as closely related species tend to have similar scores in this PCA, and thus resemble each other in their balance. This is shown by a significant phylogenetic signal in the species’ scores along the first axes of the PCA. Incorporating the evolutionary history of plant species into physiological studies can help to understand the response of plants to multiple stresses currently acting in ecological communities.


Introduction
Plants simultaneously deal with stresses of multiple origins, which might require a differential prioritization of some stress-responses against others to maximize fitness. The optimal solution to these trade-offs might be constrained by the evolutionary history of the species, in a similar way that plant responses to a single stress, such as herbivory, are macroevolutionary constrained (Pearse et al. 2017). The balance in the investments to face multiple stresses has received considerable attention within species (Eränen et al. 2009;Ariga et al. 2017;Berens et al. 2019;Pihain et al. 2019), for instance showing genetic correlations between alleles related to osmotic stress tolerance, and resistance to aphids and caterpillars (Thoen et al. 2017). However, as far as we know, it remains unknown whether the response of plants to multiple stresses is phylogenetically constrained, resulting in speciesspecific combinations of metabolic profiles that tend to resemble between closely related species. As the diversity and magnitude of different stresses vary across biomes, potential evolutionary constraints might only emerge when differences across environments are taken into account. Therefore, a first step would be to approach these evolutionary constraints within a community, where the co-occurring species have shared a combination of filters to establish.
Plants can respond to abiotic and biotic stress through different physiological mechanisms. Some plant physiological responses can be simultaneously associated with both biotic and abiotic stress, like an increase in indole acetic acid (IAA) (Epstein and Ludwig-Müller 1993), which can induce root development and growth in response to either water stress or root damage due to pathogens or herbivores. However, other physiological responses are likely related to either biotic or abiotic stresses. It is widely accepted that abscisic acid (ABA) is the main hormone controlling the stomata opening under water stress conditions (Sah et al. 2016;Gómez-Cadenas et al. 2015;Munemasa et al. 2015) and it is key for the adaptation to dry environments. Similarly, a thoroughly described plant response to limited water availability is to accumulate proline to adjust osmotic potentials (e.g., Dobra et al. 2010). On the other hand, biotic stresses caused by arthropod herbivores and necrotrophic pathogens induce jasmonic acid (JA) accumulation, as shown elsewhere (Glazebrook 2005;Bosch et al. 2014;Erb et al. 2012), which is involved in plant resistance to these pathogens. It is also well-known that salicylic acid (SA) is vital for the establishment of the systemic acquired resistance (SAR), an immune response of plants that provides long-lasting protection to infection by a broad spectrum of phytopathogens (Hartmann and Zeier 2019;Klessig et al. 2018). Nevertheless, recently it has been shown that phytohormone signaling to abiotic and biotic stresses (i.e., ABA and SA) can be compromised, for instance showing a differential prioritization along time (Berens et al. 2019). Trade-offs between plant responses to biotic/abiotic stress may also result from investing in the production of one metabolite at the expense of other defenses, such as reducing proline biosynthesis may induce the increase of lignin content (Guan et al. 2019).
Plant communities shaped by facilitation provide a convenient scenario to capture species variation in response to multiple stresses for two reasons. On the one hand, plant facilitation is prevalent in arid ecosystems where the abiotic stress is severe. Nurse plants alleviate abiotic stress by providing shade, moisture, and nutrients to the shared rhizosphere, facilitating the establishment of other species (Foronda et al. 2019). Although facilitated species find a milder abiotic microhabitat under nurse plants, they also cope with new neighbors that might impose increasing biotic stress (Ristok et al. 2019). Therefore, an increase of biomass in vegetation patches, compared to plants growing solitary, can trigger attraction to pathogens with low host specificity, such as fungi (Gilbert and Webb 2007;Kluger et al. 2008;Hersh et al. 2012;Spear and Mordecai 2018), or herbivores (Novotny and Basset 2005). In this sense, clearly defined microhabitats in ecosystems governed by facilitation (i.e. growing solitary vs. in a vegetation patch) provide a wide range of ecological contexts from the abiotic stressful bare ground to the biotic stressful environments driven by plant neighbors and their associated pathogens and herbivores. Thus, conspecific individuals growing in these contrasting conditions can be considered to provide a range of variability of plant responses to biotic and abiotic stresses. On the other hand, environmental stresses can vary widely in different biomes, and therefore physiological trade-offs across all the angiosperms might be blurred by these differences. Plant communities, like those governed by facilitation, provide a set of species that have been co-occurring in the same community for long periods, facing a shared combination of filters to establish. This scenario prevents the mixture of contrasting macro-environmental conditions, resulting in a convenient framework to assess potential evolutionary constraints in plant responses to multiple stresses.
Here, we profited from a Mediterranean semiarid environment where facilitative interactions generate contrasting ecological conditions by structuring the plant community in vegetation patches (Delalandre and Montesinos-Navarro 2018). In this system, we test whether plant species show a trade-off in their investment in metabolites related to abiotic and biotic stress tolerance. Specifically, we expect that if that trade-off exists, those species with a higher investment in proline and ABA (associated with responses to abiotic stress) will show lower levels of JA and SA (associated with biotic stress) and vice versa. In case there is evidence supporting our first hypothesis, we will also test whether this trade-off is evolutionarily conserved across the plant species that co-occur in the community. Understanding potential trade-offs in the metabolic responses across coexisting plant species can contribute to link evolutionary and physiological processes taking place at the community level.

Study site and field sampling design
We performed the study in gypsum outcrops located within 20 km 2 in Eastern Spain (38° 29′ N, 0° 44′ O; elevation: 568 m). The climate is semi-arid, with an annual mean rainfall of 414 mm, and a variation of 55 mm between the driest and the wettest months. Mean daily maximum and minimum temperatures range from 3.3 to 13.3 °C in January, and from 18.4 to 30.6 °C in August. The plant community is a scrubland with abundant chamaephytes such as Helianthemum squamatum (L.) Dum. Cours., Teucrium libanitis Schreb. and Helianthemum syriacum (Jacq.) Dum.Cours. (Delalandre and Montesinos-Navarro 2018).
Plant sampling took place in November 2018. We selected the most representative thirty plant species in the gypsum outcrops. This design includes all the dominant species (Delalandre and Montesinos-Navarro 2018), and around 80% of the species of the community. The number of species used (30 species) responds to the ecological constraint of the species richness of the community and, although it represents a very small proportion of all the angiosperms, it properly characterizes almost all the species present in the studied community. Phylogenetic comparative methods can be used for different purposes. Therefore, while the assessment of a phylogenetic signal can explore evolutionary patterns across broad taxonomic levels (Gómez et al. 2010;Werner et al. 2014), it has also been widely applied to assess the phylogenetic basis of niche differentiation at the community level (Webb et al. 2002) and explore the ecological and evolutionary mechanisms underlying this differentiation (Ackerly 2003), progressively improving its correct implementation (Münkemüller et al. 2015). In order to contextualize the magnitude of our sample size, we can compare it with previous studies that have also tested for a phylogenetic signal in other plant traits within a community. Focusing on plant traits related to plant-soil feedback, phenology, flammability, plant-mycorrhizal interactions or germination, most of the studies considering a particular community used similar sample sizes as our study, ranging from 17 to 60 plant species (Anacker et al. 2014;CaraDona and Inouye 2015;Simpson et al. 2016;Li et al. 2017;Chen et al. 2017;Salazar et al. 2018). Only in two studies, the sample size was above 100 species which corresponded to megadiverse ecosystems such as the Brazilian Atlantic forest or subtropical rainforest in Taiwan (Chang-Yang et al. 2016;Gastauer et al. 2017). Whether the phylogenetic signal was significant or not depended on the specific trait analyzed rather than on the sample size, as some studies with only 17 plant species found a significant signal while others with more than 30 species found nonsignificant phylogenetic signals (Montesinos-Navarro et al. 2012; Nagahama and Yahara 2019). As our design covers most of the species in the community studied, we consider that the 30 species analyzed can properly characterize the constraints among those species that have shared the same environmental filters present in the selected ecosystem.
For each species, we sampled a pair of individuals of a similar size, each of them in: (a) a harsh abiotic environment with low biotic stress (i.e., on the bare ground), and (b) a milder abiotic microhabitat with increasing biotic stress imposed by neighbors (i.e., in a vegetation patch). This sampling pursues to maximize the intraspecific variation captured for each species. We collected 10 g. of fresh leaf tissue for each individual within four hours on a single day avoiding water loss by refrigerating them during the whole process, and froze them at − 20 ºC until we analyzed the samples. Every plant sampled seemed healthy and showed no sign of pathogen infection or herbivory.

Proline concentration
We performed proline analysis as described by Bates et al. (1973) with some modifications. We extracted 50 mg of ground frozen leaf tissue by sonication for 30 min in 5 ml of 3% sulfosalicylic acid (Panreac, Barcelona, Spain). After centrifuging at 4000×g at 4 °C for 20 min, the supernatant was mixed with glacial acetic acid and ninhydrin reagent (Panreac) in a 1:1:1 proportion (v:v:v). We incubated the reaction mixture at 100 °C for 1 h. After centrifugation at 2000×g at 4 °C for 5 min, proline concentration was spectrophotometrically determined at 520 nm. We performed a standard curve with pure proline (Sigma-Aldrich, St. Louis, MO, USA).

Phytohormone analysis
We carried out hormone extraction and analysis as described in Durgbanshi et al. (2005) with few modifications. Briefly, 100 mg of ground frozen leaf tissue was extracted in 2 mL of ultrapure water using a mill ball equipment (MillMix20, Domel, Železniki, Slovenija) after spiking with 50 ng of [ 2 H 6 ]-ABA, [ 13 C 6 ]-SA, dihydro jasmonic acid (DHJA) and [ 2 H 2 ]-IAA. Samples were centrifuged at 4000×g for 10 min at 4 °C after the extraction and supernatants were recovered and pH adjusted to 2.8-3.2 with acetic acid. We partitioned the extracts twice against 2 mL of diethyl-ether and the organic layer evaporated under vacuum in a centrifuge concentrator (Speed Vac, Jouan, Saint Herblain Cedex, France). The residue was resuspended in 0.5 mL of methanol: water 10:90 and filtered through 0.22 µm polytetrafluoroethylene membrane syringe filters (Albet S.A., Barcelona, Spain) and directly injected into an ultra-performance liquid chromatography system (Acquity SDS, Waters Corp., Milford, MA, USA). Chromatographic separations were carried out on a reversed-phase C18 column (Gravity, 50 × 2.1 mm 1.8-µm particle size, Macherey-Nagel GmbH, Germany) as stationary phase using a methanol:water (both supplemented with 0.1% acetic acid) gradient at a flow rate of 300 µL min −1 as mobile phase. We quantified hormones with a triple quadrupole mass spectrometer (Micromass, Manchester, UK) connected online to the output of the column through an orthogonal Z-spray electrospray ion source.

Statistical analyses
To assess the combination of the five metabolites contents that best explain the variation observed across individuals, we performed a phylogenetic principal component analysis. Non-phylogenetic principal component analyses assume that the sample consists of independent data points, which is likely violated by phylogenetic data from related species. If phylogeny is ignored, type I error of the statistical estimators and hypothesis tests can be increased, which can be avoided using phylogenetic principal component analyses (Revell 2009). However, these procedures provide residuals and scores in the original, phylogenetically dependent, species space, instead of in a phylogenetically independent space, so that the scores from these analyses can still be analyzed to explore phylogenetic patterns (Revell 2009). For three out of the thirty species, the values of SA were one order of magnitude larger (Santolina chamaecyparissus: 5386.6; Coris monspeliensis: 1076.9) or lower (Sedum sediforme: 15.7) than in any other species (the median of the rest of species: 151.8). We excluded these species from the analyses (Table S1, in Appendix S1), and the rest of the values were log-transformed to fulfill the assumptions of parametric statistics. The phylogenetic principal component analysis was performed using the functions "phyl. pca" and "scores" to obtain the PC scores for individuals based on a rotation computed for the species means, both functions in the "phytools" package (Revell 2012) in R version 3.6.1 (R Core Team 2018). Hereafter, for simplicity, we refer to the phylogenetic principal component analysis as PCA.
We assembled the phylogenetic relationships among the plant species in the community with the R function "S.PhyloMaker" (Qian and Jin 2016), which uses an expanded version of the time-calibrated angiosperm species-level phylogeny (Zanne et al. 2014). Most of the studies and the statistical methods to assess the evolutionary conservatism of traits use a single trait value for each species, ignoring within-species variation (Ives et al. 2007). However, we used Ives et al. (2007) approach to take into account variation within species. Furthermore, we tried to maximize the intraspecific variation captured for each species by sampling the most contrasting situations representing abiotic versus biotic stress (i.e., plants growing alone vs. growing in vegetation patches, in a semi-arid community). Thus, we used the scores in the PC1 of the paired individuals (growing alone and in a vegetation patch) to calculate the mean and standard error (SE) for each species, and computed 1 3 the phylogenetic signal based on the statistics Blomberg's K (Blomberg et al. 2003) following Ives et al. (2007). This test demonstrates correct Type I error rate and good power for simulated datasets with 20 or more species, and therefore it is appropriate despite our low sample size. This approach constructs a phylogenetic tree that combines the phylogenetic covariance among species with the intra-specific variation. It includes the later by elongating the terminal branches of the tree according to the value of the intra-specific variation relative to the variance of the evolutionary process. For three species (Thymus vulgaris, Plantago albicans, Brachypodium phoenicoides), some metabolites could not be quantified in some of the individuals (Table S1, in Appendix S1), preventing the calculation of the SE and their inclusion in this analysis. Statistical significance was assessed with 999 randomized repetitions using the function "phylosig" in the package "phytools" for R (Revell 2012). K values range from 0 to ∞, with values K > 0 indicating the existence of a phylogenetic signal. Values between 0 and 1 indicate low trait conservatism and values higher than 1 are indicative of high (i.e., more than the expected under Brownian motion) phylogenetic signal.

Results
The content of the five metabolites was highly variable across the 30 studied species (Tables S1-S3, in Appendix S1). The first and second principal components of the PCA using the five metabolites explained a 40% and 24% of the total variance respectively (Fig. 1, scores of the first principal component in Table 1). The loadings of the five metabolites in the first principal component represent a combination concordant with different responses to biotic versus abiotic stress (Fig. 1). On one pole of the axis, the first principal component segregated the metabolites related with abiotic-stress response (loadings proline = 0.6, and ABA = 0.3; Table S1 in Appendix S1). On the other pole of the axis were segregated the metabolites usually related to biotic stress tolerance (JA = − 0.9 and SA = − 0.5; Table S3 in Appendix S1). Finally, a metabolite related to both types of stresses was located between the previously described extremes (IAA = 0.2, Table S2) (Fig. 2). The distribution of samples scores along the PC1 was independent of the time of the day of each collection (i.e., hours since the first sample was collected: F = 0.3, df = 1, p-value = 0.6).
When the mean and SE of the PC1 scores per species were considered, the species scores (i.e., strategy defined by each species) was phylogenetically conserved (i.e., significant phylogenetic signal, Blomberg's K = 0.2, p-value = 0.03), indicating that close relatives tended to show a similar score, and thus a similar strategy regarding their investment in the 5 metabolites analyzed (Fig. 2).

Discussion
Plant responses to stress combinations result in emergent metabolic profiles that are unlikely predicted by the effect of individual stresses (De Vos et al. 2006;Makumburage et al. 2013;Zandalinas et al. 2018). This might result in constrained physiological investments in stress tolerance when the plants face multiple stresses simultaneously. Our results show that in the studied community, plant species tradeoff their investment in abiotic-stress tolerance against their investment in biotic stress tolerance defining two metabolic strategies: some plant species tend to invest in abiotic stress-related metabolites (proline, ABA) while others tend to invest in biotic stress-related metabolites (JA and SA). Interestingly, these strategies seem to be evolutionarily conserved among the species that co-occur in this community, with close relatives tending to resemble each other (Fig. 2). The magnitude of such phylogenetic signal, although significant, was not very high, suggesting that despite being to some extent evolutionarily constrained, the metabolic strategy also had some flexibility. In addition, it is worth to highlight that this study is community focused (with a limited number of species), and therefore does not intend to draw a general conclusion about a phylogenetic constraint through the whole phylogeny of angiosperms. On the contrary, it shows that when we consider the set of species that had overcome a shared combination of filters to establish in a given environment, they show evolutionarily conserved strategies to deal with biotic/abiotic stresses.
Plant species might deal with physiological constraints to respond to both abiotic and biotic stresses, as for instance, responses involving ABA pathway (associated with abiotic stress) can be influenced by the SA or JA pathways (associated with biotic stress) (Anderson et al. 2004;Fujita et al. 2006;De Torres Zabala et al. 2009;Kissoudis et al. 2015;Ximénez-Embún et al. 2018). Potential tradeoffs between the plants' investment in their responses to abiotic versus biotic stresses have been mainly approached at the individual (or within species) level (Eränen et al. 2009;Ariga et al. 2017;Berens et al. 2019), but as far as we know, this is the first attempt to consider it at a macroevolutionary scale, testing  Table S1 for an evolutionary conservatism of the tradeoff across co-occurring species. Our analyses take into account intra-specific variation by sampling two individuals in contrasting environments, which may not be sufficient to capture the extent of intraspecific variation in metabolic responses fully, but may be a broad approximation of the maximum withinspecies variation. A deeper characterization of intraspecific variation will allow assessing potential phylogenetic patterns in species' plasticity when multiple individuals of the same species are exposed to a wider diversity of ecological contexts varying in their balance of biotic/abiotic stresses. Some studies have explored the evolutionary conservatism of plastic responses for instance considering phenology shifts (Rafferty and Nabity 2017), but this is a far unexplored topic regarding plant metabolic responses to multiple stresses. Despite the intraspecific variation in metabolite contents resulting from our sampling, under contrasting environmental conditions, our results show that the sampled individuals of the same species tend to have similar metabolic strategies and close relatives resemble each other in their metabolic strategy. Similar phylogenetic patterns have been found showing shared metabolic responses to folivory between species within the same genus (Rivas-Ubach et al. 2016b;, or unique and conserved metabolic responses to drought stress levels within species belonging to the same family (Sanchez et al. 2012). Further research exploring whether the intraspecific trade-off in the investment to respond to biotic/abiotic stresses is also evolutionarily conserved, may contribute to a further understanding of these evolutionary processes. The application of profiling techniques to plant physiology and ecology has revealed considerable plasticity of metabolomes under different ecological contexts (Sardans et al. 2011;Rivas-Ubach et al. 2012;2016a). Understanding the phylogenetic contributions of the metabolome is critical in the current context of global change (Edwards et al. 2007;Kuntner et al. 2014;González-Orozco et al. 2016), as the underlying mechanisms can be important to predict the potential of plants to respond to different ecological scenarios.