Tritrophic interactions follow phylogenetic escalation and climatic adaptation

One major goal in plant evolutionary ecology is to address how and why tritrophic interactions mediated by phytochemical plant defences vary across species, space, and time. In this study, we tested three classical hypotheses about plant defences: (i) the resource-availability hypothesis, (ii) the altitudinal/elevational gradient hypothesis and (iii) the defence escalation hypothesis. For this purpose, predatory soil nematodes were challenged to hunt for root herbivores based on volatile cues from damaged or intact roots of 18 Alpine Festuca grass species adapted to distinct climatic niches spanning 2000 meters of elevation. We found that adaptation into harsh, nutrient-limited alpine environments coincided with the production of specific blends of volatiles, highly attractive for nematodes. We also found that recently-diverged taxa exposed to herbivores released higher amounts of volatiles than ancestrally-diverged species. Therefore, our model provides evidence that belowground indirect plant defences associated with tritrophic interactions have evolved under two classical hypotheses in plant ecology. While phylogenetic drivers of volatile emissions point to the defence-escalation hypothesis, plant local adaptation of indirect defences is in line with the resource availability hypothesis.


Results
Herbivore-induced VOCs facilitate predator recruitment. To address our aims, we studied VOCs production and EPN attraction on plants with and without root herbivore attack across 18 species of Festuca that grow in the Alps (Figs. 1, S1). On average, across all species, we found that Melolontha herbivory decreased root biomass by about 40% (see supplementary material for description of the statistical methods; Fig. 2a [31.81, 160.60], pMCMC(868) = 0.002). Predator recruitment is likely mediated by the abundance and identity of herbivore-induced root VOCs released in rhizosphere, which show strong species-level variations ( Fig. S3; PERMANOVA; species effect: F 17,148 = 3.16, p = 0.001). However, we found no effect of the herbivore treatment on volatile profiles (herbivore treatment effect; F 1,148 = 0.81, p = 0.60; root biomass effect; F 1,148 = 2.08, p = 0.04, and species by treatment interaction; F 17,148 = 0.86, p = 0.88), indicating that, although larvae ate abundantly on the roots, root VOCs induction by M. melolontha larvae, across the multivariate VOCs matrix, was weak. Herbivore effect on VOCs induction was however strong when assessing the total amount of VOCs produced (Table S4), as well as when assessing individual compounds such as the sesquiterpenoids E-β-caryophyllene and isosativene (Table S4, Fig. S4). Finally, we observed a positive association between the total constitutive production of VOCs and their inducibility ( Fig. S5; correlation corrected for spurious correlations: r = 0.83, p < 0.001; PGLS, F 1,16 = 35.82, p < 0.001). Plant species adapted to high elevation, cold and humid habitats, recruit more predatory nematodes. The climate of each species' home site from along the elevational gradient was an important factor determining EPN recruitment. Root herbivore induced-species originating from higher elevation, colder and more humid habitats (Fig. S7) produced almost twice more VOCs and attracted about twice as many nematodes than their lowland congeners from warmer and drier habitats ( Fig. 3; see PCA1 axis ranging from negative to positive values, and indicating variation from cold and humid high elevation habitats to warm and dry low elevation habitats, Table 1). Multivariate phylogenetic-corrected regression models of the different climatic variables indicated that the number of precipitation days per growing season (r = 0.7, F 1,16 = 15.54, p = 0.001), humidity (r = 0.60, F 1,16 = 9.21, p = 0.01) and number of frost days(r = 0.52, F 1,16 = 5.85, p = 0.03), which are all intimately associated with elevation, were the best predictors of nematode attraction. Furthermore, we explored the relationship between VOCs production and EPNs attraction with soil variables (Fig. S8, Table S3) and found a significant relationship between herbivore-induced VOCs and soil variables (Fig. S9, distance-based redundancy analysis (dbRDA): F 8,9 = 1.67, p = 0.01). The soil effect was mainly driven by soil relative humidity (HR), which was also positively associated with the terpenes β-pinene and E-β-caryophyllene, while soils richer in organic matter (high CEC, pH and CN values) were negatively associated with the production of the terpenes aromadendrene, β-vatirenene and isosativene. On the other hand, none of the soil variables could explain EPN recruitment measured in the bioassays (overall linear model; r = 0.36, F 8,9 = 0.18 and p = 0.99).
Phylogenetic escalation for VOCs production and EPNs recruitment. The radiation of Festuca along elevation gradients and across climatic niches carries phylogenetic signal that is expected under a Brownian motion model of trait evolution (K = 0.5, p = 0.04; K = 0.4, p = 0.06 respectively; Figs. S10 and S11, Table S5). Furthermore, our data was remarkably consistent with classic theory: we found phylogenetic escalation in the total amount and diversity of inducible VOCs (Table 1, Fig. 3). Specifically, after 7-9 divergence events from the split that separated the predominantly high altitudinal and latitudinal fine-leaved Festuca species from their broad-leaved congeners, between 13 and 16 million years ago, the plants' inducible VOCs production and their diversity increased by an average of 4 folds when induced by herbivory. Additionally, we found a weak negative correlation (r = 0.43, estimate = −52.12, F 1,16 = 3.79, p = 0.07) between branch length and the first axis of the principal component analysis of the VOCs, representing the axis of maximal scent variation across species (Fig. S10). This indicates that more ancestrally-diverged species tend to have similar chemical composition characterized by monoterpenes and low-weight hydrocarbon chains (VOCs 4,5,6,9,10,13,15), while more recently-diverged species tend to produce additional compounds, which, most importantly, included sesquiterpenes.

Discussion
We tested three major predictions of classic plant defence theory regarding plant direct defences aboveground for the potential to explain a plants' ability to recruit predators (i.e. indirect defences through tritrophic interactions) belowground. On one hand, findings showing that Festuca species adapted to high elevation cold and humid habitats had high predator recruitment after herbivore attack and chemical defence induction, provided evidence for the resource availability theory. On the other hand, we showed that increasingly-derived Festuca species were associated with highest emissions and diversity of volatile organic compounds (VOCs), thereby supporting the defence-escalation theory.
First, we found that root herbivory induced plants to be more attractive to soil-dwelling entomopathogenic nematodes (EPNs), in general, across 18 grass species of Festuca, when compared to undamaged plants. Aboveground herbivore-induced VOCs were shown to attract predatory arthropods (e.g., wasps and mites) in several dozens of systems and in a variety of habitats 2 . Later, root herbivores were also shown to induce root VOCs for attracting soil-dwelling EPNs 35 , but the effect of VOCs induction on predator recruitment were never tested across several plant species simultaneously. In summarizing work on chemically-mediated belowground tritrophic interactions, terpene-based compounds seem to play a major role in root defence 36 . For instance, corn plants genetically modified for enhanced E-β-caryophyllene production were more protected by EPNs against Diabrotica virgifera corn root larvae in the field 37 . Nonetheless, several other VOC classes can also induce attraction -as well as repulsion -of EPNs 35 . Because we found no direct correlation between total VOC production and EPN attraction, the unique identity of the VOC blends should be the main driver of EPN behaviour. Accordingly, previous work on Festuca rubra showed that EPN recruitment relies on the balance between the emissions of . Phylogenetic and climatic effect on root volatile organic compound (VOCs) production and entomopathogenic nematodes (EPNs) attraction across 18 species of Festuca growing in the Alps. Each dot represents one species, and dots are color-coded based on their optimal climatic distribution. Green colours represent colder and more humid climates, while red colours represent warmer and drier habitats. Patristic distance represents the total branch length from the last common ancestor, and climate represents the first axis of a principal component analysis ( Fig. S6; the axis ranges from negative to positive values, and indicates variation from cold and humid high elevation habitats to warm and dry low elevation habitats). VOCs C are the volatiles emitted by undamaged roots, VOCs I are volatiles emitted by herbivore damaged roots and VOCs I-C is the difference between the two and represents each species root volatile inducibility potential. VOCs H represents the Shannon diversity values. VOCs PCA1 represents the first axis of the PCA ordination of all VOCs emitted at the induced state (Fig. S10). EPNs C indicates the number of EPNs attracted to undamaged roots. EPNs I indicates the number of EPNs attracted to herbivore-damaged roots and EPNs I-C is the difference between the two. The lines show significant correlations (linear model (lm) and phylogenetic corrected model (pgls) for phylogeny and climate, respectively; p < 0.05 for filled lines and p < 0.1 for dotted lines, see Table 1). three attractants (α-curcumene, 1-undecene, nonadecatriene) and one deterrent (decane) compound 38 . Contrary to theory based on resource allocation within plants 39 , we observed a positive association between the total constitutive production of VOCs and their inducibility. Such positive association is not unprecedented 21,[40][41][42][43] . For instance, positive correlations between constitutive and inducible defences were observed across >50 species of milkweeds (Asclepias) and across 21 wild cotton (Gossypium) species. Therefore, we consider our findings as a common pattern resulting from repeated evolution independently of the phylogenetic history, the geography of diversification, and of the cost of induction across species.
Second, our findings show that species growing at high elevation colder and more humid habitats, which are generally less productive 44 , harbour enhanced emission of root-herbivore induced VOC and better predatory recruitment power. These observations point to the prominent role of the resource availability hypothesis and refute the latitudinal/elevational gradient hypothesis. Our most strongly-supported climatic and soil predictors of induced VOCs production and EPN recruitment were related to temperature and humidity, and organic matter rich soils, respectively. Speciation into more climatically-stressed and resource-poor high elevations has thus been followed by a generalized adaptation to produce more attractive VOC blends. In the case of Festuca, we observed a relatively strong phylogenetic signal for elevation and climatic niche, indicating that more closely-related species tend to share more similar altitudinal zones (e.g., species of the F. pratensis clade in the phylogeny; Fig. S12) 45 . Nonetheless, we inferred several phylogenetically-independent transitions to high elevation niches across the phylogeny, as shown by the observation that representative alpine species span multiple clades in the phylogeny (Amphigenes: F. pulchella; Eskia-Dimorpha: F. quadrifolia; Exaratae: F. violacea; Festuca: F. alpina, F. halleri; Fig. 2, Figs. S12-14), thus having independently colonized the alpine zone regardless of phylogenetic history. Altogether, our results indicate that fescue radiation into alpine niches from different ancestors has resulted in optimal adaptation towards the production of highly attractive VOCs blends. This suggests that the resource-availability hypothesis is associated with a convergent evolutionary pattern. However, parallel to the plant adaptation perspective, we cannot exclude that insect adaptation may have also driven differences in VOC production along elevation. At low elevation, where plants and root herbivores co-occur frequently, selective pressures may have promoted insect adaptations enabling the control of plant VOC emissions. Along this line, the composition of oral secretions of adapted herbivores has been shown to limit VOC induction 46 . At the opposite end of the elevation gradient, a relaxation in plant-herbivore interactions may have hampered such insect adaptation process, thereby increasing VOC emissions at high elevation.
Third, we provide the first evidence of defence escalation for plant traits mediating tritrophic interactions between plants herbivores and predators. The continuous model of evolution we used indicates that defence escalation has occurred gradually via continuous character displacements in VOCs, but also according to subsequent speciation events; indeed we found strong correlation between VOCs (I, I-C, H) and the patristic distance and the number of intervening nodes from common ancestor (Table 1). Festuca ancestors likely originated in the www.nature.com/scientificreports www.nature.com/scientificreports/ lowland southern Palaearctic region and then diversified preferentially upwards into high mountains and also into high latitudinal Northern and Southern Hemisphere regions and colder climates [47][48][49] . Nonetheless, even when accounting for this directional radiation, we still found consistently stronger inducibility in species that inhabit higher altitudes. To date, other tests of the defence escalation hypothesis include interactions between plants and herbivores (i.e. direct defences). For instance previous work highlighted an increase in more toxic angular furanocoumarins in wild parsnip species 50 , an increase in tolerance (i. e., the ability to regrow) in milkweeds 51 or an increase in alkaloid diversity in Ficus species 52 . While this pattern might apply to a wide range of defence traits, several studies failed to find correlation between phylogenetic branch lengths or speciation events and defence traits. For instance, the total amount of cardenolides decreases with increasing phylogenetic node numbers in milkweeds 51 . In the genus Streptanthus, while the total amount of glucosinolates varies independent of the phylogeny, the production of aliphatic glucosinolates increases over evolutionary time, while glucosinolate diversity declines 53 . Thus, our findings align with previous studies in that not all traits show consistent patterns of defence escalation. Indeed, while induction of VOCs in Alpine Festuca species correlates with phylogenetic branch length, the constitutive production of VOCs varies independently of the phylogeny.
Altogether, we have shown that for high-elevation Festuca species, herbivores induce consistently more and more diverse VOC emissions and attract more EPNs than their low-elevation counterparts. Moreover, our results suggest that higher EPN recruitment in harsh conditions relies on specific VOC profiles rather than on total amount of VOC. This is in line with the assumption that the cost of herbivory depends on the plant's opportunities to acquire resources. In this case, the short growing season of alpine environments makes the period of resource acquisition brief and intermittent during the growing season. While high-elevation soils contain lower diversity and abundance of arthropods and predators 6 we cannot completely exclude moderate amounts of root herbivory by different guilds, and thus we expect selection of the induction of optimally-attractive VOC blends upon unexpected herbivore attack. Finally, with high altitude conditions, we might also expect pleiotropy of plant defence traits, in which the production of specific secondary metabolites can ameliorate plant fitness in harsh environments. Indeed, a recent review has indicated that VOC emission rates for phototrophs in extreme environments are frequently higher than that of organisms from less stressful environments 54 . Therefore, the reason for widespread VOC emissions from different extreme environments and its association with herbivore pressure and phylogenetic history deserves further attention, as these compounds may have important roles in stress resistance and adaptation to extreme habitats.

Material and Methods
Study system. With more than 500 species described, the genus Festuca represents a major evolutionary line of the tribe Poeae. These tufted grasses are herbaceous perennials with a height range of 10-200 cm and a cosmopolitan distribution, growing on every continent except Antarctica 47,55 . Within the entire European Alpine range, there are about 50-70 species depending on their classification 56 . From this range, we collected different taxa at 31 sites (Fig. S1, Table S1), ultimately selecting 18 confirmed species with which we were able to perform belowground olfactometer behavioral assays (see below). The plants were collected from pristine sites that represent the optimal niche of each species; www.infoflora.ch) (Fig. S11). At each site, we excavated 20-30 individuals with their own soil and potted them in plastic pots (20 cm diameter, 25 cm high) by completing the native soil with standard soil mixture composed of 2/3 regular potting soil (5 parts compost, 4 parts peat substitute based on wood fiber and 1 part topsoil, pH in water = 7.45, conductivity = 2 mS cm −1 , Ricoter, Aarberg, Switzerland) and 1/3 sand in garden plastic pots (25 cm diameter, 30 cm high) in an outdoor common garden at the Botanical Garden of Neuchâtel, Switzerland. For phylogenetic reconstruction of the Festuca species see details in supplementary methods, Tables S1-S2, and Figs. S12-S14.
At each collection site, we also collected soil samples to characterize edaphic variables associated with each species. We sampled about 1 Kg of soil from 10-30 soil cores (5 cm wide and 20 cm deep). Soils samples were homogenized, dried at 40 °C for 48 hours, sieved at 2 mm, and ground using agate mortars for subsequent physicochemical analyses. Specifically, we measured (1) relative humidity (HR) by weighting soils before and after drying at 105 °C for 3 days; (2) organic matter content (OR) through loss on ignition by weighting before and after burning 10 gr of soil at 450 °C for 2 hrs; (3) pH in 2.5 soil volume of deionized water; (4) cation exchange capacity (CEC) following the cobaltihexamine chloride method; (5) carbon to nitrogen ratio (CN) using an elemental analyser (FLASH2000, Thermo Fisher Scientific, Waltham, Massachusetts, United States); (6) total carbonates (CaCO3) using a calcimeter (Bernard type), and (7) bioavailable phosphorous (P_bio) using the Kjeldahl digestion method 57 . Additionally, for each species we extracted ecological indicator values for soil fertility (Soil_N) based on Landolt et al. 58 . To measure EPNs infectivity potential of each site we ran an infectivity test shortly after soil collection (maximum one week) to avoid loss of potential EPNs present in soil. Three soil subsamples (200 g) were added with five larvae of the greater wax moth (Galleria mellonella) and stored in the dark 59 . After six days we checked for dead larvae, which were placed in White traps for nine days and monitored for EPNs emergence. The experiment was run twice on the same soils to corroborate initial results. The average values of each soil variable are summarized in Table S3.

Root herbivores and predatory nematodes.
In order to study Festuca interspecific variation in belowground tritrophic interactions along environmental gradients, we used the dietary generalist herbivore, Melolontha melolontha (Coleoptera: Scarabaeidae). The larvae develop in soils by feeding on roots for up to three years, causing severe damage to plants, particularly to Poaceae in low-elevation grasslands, where herbivore population densities can reach 100 individuals/m 2 , but the beetle larvae practically disappear in the Alpine zone 60 . M. melolontha larvae were collected in pastures around Switzerland and stored at 10 °C in potting soil. Although M. melolontha is generally found between 300 and 600 m above sea level (asl), Scarabaeidae root herbivores can easily cause damages in mountainous zones 61 . While we used a single herbivore species in order to compare (2020) 10:2074 | https://doi.org/10.1038/s41598-020-59068-2 www.nature.com/scientificreports www.nature.com/scientificreports/ herbivore effects on different plant species, we assumed that the generalist diet of M. melolontha enables to extrapolate herbivore effects of closely-related Scarabaeidae larva damaging mountainous zones. As soil predators, we used Heterorhabditis megidis (Rhabditida: Heterorhabditidae) entomopathogenic nematodes (EPNs) (Andermatt Biocontrol, Andermatt, Switzerland). EPNs mode of action can be summarized in three steps: i) penetration into the arthropod host by the third instar resistant stage juveniles, ii) release of symbiotic bacteria in the hemolymph of the host and iii) EPNs feeding on the newly-formed bacterial soup for 2-3 generations inside the host 62 . Species of EPNs in the genera Steinernema and Heterorhabditis, including H. megidis, are found in most habitats, including the Swiss Alps. In a previous survey of nematodes along elevation gradients, we collected H. megidis up to 2000 m asl, and we observed that EPNs, in general, range from 500 to 2750 m. asl 63 . Most EPNs actively seek their hosts in the soil and, while CO 2 can serve as a general cue for host searching, more specific plant-associated volatile compounds have been shown to inform EPNs of the presence of a root insect herbivore in the surroundings 4,64 . H. megidis relies on host-plant cues for several systems such as coniferous plants, grasses and other angiosperms 36 . Additionally, nematodes belonging to Heterorhabditis genera cause relatively fast larval mortality (i.e., after 7 days) for different Scarabaeidae species 65 . Belowground olfactometer bioassays. After a period of about 10 months (from summer 2016 to spring 2017), including an overwintering period of about 6 months in which all aboveground parts of the plants were cut to the ground, 5 plants per species were randomly selected from the initial pool and divided into two equally-sized tufts for performing olfactometer bioassays. For three species with wide distribution ranges (F. heterophylla, F. ovina, F. rubra), we doubled the sample size in order to include two distinct collection sites per species. Individual tufts of each species were separated in half to obtain unique clonal pairs of plants, which were potted in custom-made two-arm belowground olfactometer pots as described in Kergunteuil et al. 6 , one pair of clones per olfactometer arm and filled with the same soil mixture as previously described. One month later, one plant per olfactometer was infested with three M. melolontha L2 larvae that were previously starved for three days. Two days after the plant infestation, the central pieces of the olfactometers were filled with 10% moist sand (Neogard, Gontenschwil, Switzerland) and assembled with each pair of plants. Three days after the initial damage, 2000 six-day-old H. megidis EPNs were released in the central chamber of the olfactometers and allowed to move toward either uninfested or infested plants for 24 hours. The behavioural choices of nematodes were recorded by counting the individuals collected in the arms connected to 'uninfested' and 'infested' treatments after the Baermann funnel extraction 4 .
Chemical analyses of VOCs. Immediately after the olfactometer bioassays described above, all the plants were removed from the pots, and roots were carefully washed and ground to a fine powder in liquid nitrogen. The root material was then stored at −80 °C before chemical analyses of volatile organic compounds. All analyses were done with a gas chromatograph (GC) (Agilent 7890A) coupled with a mass spectrometer detector (MSD) (Agilent 5975C). The samples were prepared and introduced into the GC with the use of robotic multipurpose samplers (MPS, Gerstel GmBH). The compounds were separated on Agilent HP-5MS columns (30 m length × 0.25 mm i.d., and 0.25 µm film thickness). In all cases, the MSD transfer line temperature was set at 280 °C and the ion source and quadrupole temperatures were set at 230 °C and 150 °C, respectively. The electron impact (EI) mode was used with a scanning over the mass range of 33-250 m/z. The chemical compounds were trapped through the use of solid-phase microextraction (SPME) headspace technique. The samples were incubated in a 20 ml glass vial for 3 minutes at 35 °C before inserting a 100 µm polydimethylsiloxane (PDMS) coated fiber (Supelco, Bellefonte, USA) into the headspace for 20 minutes. Afterward, the compounds were thermally desorbed from the fiber for 210 seconds (splitless mode, 250 °C, 6.5 psi pressure, 210 mL/min purge flow, helium carrier gas) before injection onto the GC column. The initial column temperature of 50 °C was held for 1 min, then it was ramped at 6 °C/ min until 250 °C (hold time 1 min), and finally a 3 min post run at 260 °C. The helium flow rate was 0.9 mL/min (constant flow mode). The analysis of root-emitted volatiles was based on 26 compounds that were not detected in blank samples. While 3 volatiles were identified using pure standards, the additional volatiles were tentatively identified by comparing the mass spectra with the NIST05 mass spectra library (Table S4). The peak area for each detected compound was divided by the peak area of the internal standard (Tetralin ® , Sigma-Aldrich, St. Gallen, Switzerland), and the root emissions of volatiles were finally given as tetralin equivalent nanograms of compound released by gram of root biomass collected 66 . For subsequent statistical analyses, we calculated the total amount of VOCs produced per species, and the Shannon diversity index H (function diversity from the package vegan 67 ). We did not include the number of individual peaks (i.e. peak richness), since it is highly correlated with both the total amount (t 1,158 = 2.92, p < 0.001), and the Shannon diversity (t 1,158 = 6.30, p < 0.001).

Statistical analyses. Effect of herbivory on root biomass, nematode recruitment and VOCs production.
To address the impact of the herbivore treatment, we first calculated standardized effect sizes for root biomass and EPNs recruited between herbivore-infested plants minus the control plants using Hedge's g statistics (function cohen.d in package effsize in R 68 ). We next used Bayesian phylogenetic mixed models (BPMMs), as implemented in the R package MCMCglmm 69 , to test the effect of the herbivore treatment on average root biomass, EPNs recruitment, total and all individual VOCs emitted across all species, while taking into account phylogenetic relationships among Festuca species. Differences of VOCs production across species and herbivore treatment were also assessed with a permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis distance matrices (function adonis from the package vegan 67 ), and average species and treatment differences were visualized using non-metric multidimensional scaling (NMDS) (function metaMDS in vegan) (Fig. S3). All VOCs data were log-transformed to normalize their distribution (function decostand in vegan). To test for potential negative correlations between constitutive VOCs (species means) and their inducibility (i.e., the difference in mean VOCs values for each species between control and damaged plants), we conducted tests to control Scientific RepoRtS | (2020) 10:2074 | https://doi.org/10.1038/s41598-020-59068-2 www.nature.com/scientificreports www.nature.com/scientificreports/ for potential spurious correlations 70 in Matlab (version 7.5.0.342 -R2007b, MathWorks Inc., USA). We also performed phylogenetically controlled analyses to tests for correlations between mean constitutive VOCs levels and their inducibility.
Phylogenetic signals and tests of the defence escalation hypothesis. We first estimated the phylogenetic signal (Blomberg's K) for all individual VOCs, their sum and their Shannon diversity, and the climatic niche (PCA1) using species-level averages with the function phylosig of the package phytools in R 71 . A K value of 1 indicates that trait values are consistent with the tree topology and a random walk model (i.e., trait similarity is directly proportional to the extent of shared evolutionary history). A K value close to 0 indicates no influence of shared ancestry on trait values (i.e., phylogenetic independence). When K > 1, the traits have more phylogenetic signal than expected from a random walk model on the given phylogeny 45 . Next, to test for the defence escalation hypothesis, we calculated the patristic distance from root to tip (i.e. the number of apomorphic step changes separating two taxa on a cladogram), for determining the amount of divergence of each species from their common ancestor. We then performed correlation analyses between the patristic distance and the constitutive, induced and inducibility of VOCs and nematode recruited trait values of the 18 species of Festuca. We also regressed the first axis of a principal component analysis (PCA) of all VOCs (function dudi.pca, packages ade4 72 ) against patristic distance. Finally, in addition to using branch length, we also regressed traits values against the number of intervening nodes from root to tip for each species (e.g. as in Agrawal and Fishbein 51 ), although the two measures are highly correlated; r = 0.88, t 1,16 = 6.89, p < 0.001).
Climate effect on VOCs production and nematode recruitment. To measure the effect of climate on constitutive, induced and inducibility of VOCs production and EPN attraction, we first constructed the climatic niche optima for each species based on occurrence data, and the first axis of the ordination of the climatic variables (Figs. S6 and S7). We then regressed this first PCA axis with each response variable using phylogenetic corrected analyses (PGLS) (function pgls, package caper 73 ). In all cases, to reduce Type I errors, we optimized branch length thorough maximum likelihood 74 . We next dissected the importance of each individual climatic variable by running a multivariate PGLS regression model.
Soil effect on VOCs production and nematode recruitment. To address potential correlations between the soil physico-chemical variables and the whole VOCs matrix we mapped the soil variables (Table S3, Fig. S8) on the VOCs ordination based on a distance-based redundancy analysis (dbRDA) (function cca in the package vegan). We performed the analysis using both Euclidean and Bray-Curtis distances, and obtained very similar results, and thus used the Euclidean distance results to show. Finally, we performed a multivariate regression analysis for testing the effect of all soil variables on EPNs attraction.

Data availability
Data underlying this article can be accessed on Dryad Digital Repository at http://dx.doi.org/XXX/doi.goes.here, and used under the Creative Commons Attribution licence.