International Journal of Molecular Sciences Article We Are What We Eat: A Stoichiometric and Ecometabolomic Study of Caterpillars Feeding on Two Pine Subspecies of Pinus sylvestris Albert Rivas-Ubach 1,2, * , Josep Peñuelas 2,3 , José Antonio Hódar 4 , Michal Oravec 5 , Ljiljana Paša-Tolić 1 , Otmar Urban 5 and Jordi Sardans 2,3 1 Environmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, Richland, WA 99354, USA; Ljiljana.PasaTolic@pnnl.gov 2 CREAF, Center for Ecological and Forestry Applications, Cerdanyola del Vallès, 08913 Catalonia, Spain; josep.penuelas@uab.cat (J.P.); j.sardans@creaf.uab.cat (J.S.) 3 CSIC, Global Ecology Unit CREAF-CEAB-CSIC-UAB, Cerdanyola del Vallès, 08913 Catalonia, Spain 4 Grupo de Ecología Terrestre, Departamento de Biología Animal y Ecología, Facultad de Ciencias, Universidad de Granada, 18071 Granada, Spain; jhodar@ugr.es 5 Global Change Research Institute, Czech Academy of Sciences, Bĕlidla 4a, CZ-603 00 Brno, Czech Republic; oravec.m@czechglobe.cz (M.O.); urban.o@czechglobe.cz (O.U.) * Correspondence: albert.rivas.ubach@gmail.com  Received: 19 November 2018; Accepted: 20 December 2018; Published: 24 December 2018  Abstract: Many studies have addressed several plant-insect interaction topics at nutritional, molecular, physiological, and evolutionary levels. However, it is still unknown how flexible the metabolism and the nutritional content of specialist insect herbivores feeding on different closely related plants can be. We performed elemental, stoichiometric, and metabolomics analyses on leaves of two coexisting Pinus sylvestris subspecies and on their main insect herbivore; the caterpillar of the processionary moth (Thaumetopoea pityocampa). Caterpillars feeding on different pine subspecies had distinct overall metabolome structure, accounting for over 10% of the total variability. Although plants and insects have very divergent metabolomes, caterpillars showed certain resemblance to their plant-host metabolome. In addition, few plant-related secondary metabolites were found accumulated in caterpillar tissues which could potentially be used for self-defense. Caterpillars feeding on N and P richer needles had lower N and P tissue concentration and higher C:N and C:P ratios, suggesting that nutrient transfer is not necessarily linear through trophic levels and other plant-metabolic factors could be interfering. This exploratory study showed that little chemical differences between plant food sources can impact the overall metabolome of specialist insect herbivores. Significant nutritional shifts in herbivore tissues could lead to larger changes of the trophic web structure. Keywords: plant-insect; metabolomics; stoichiometry; processionary moth; scots pine; secondary metabolites; herbivory 1. Introduction Plants contain all the necessary nutrients for herbivore insects although the absolute amount can vary even within the same individual of a same plant species [1]. Most herbivore insects are specialized to feed on specific plant species or families whereas only a few are generalists, being able to feed on diverse plant species [2]. In this regard, several studies have proven the co-dependency between host plants and herbivores demonstrating a tight interdependent plant-insect co-evolution [3–5]. A balanced diet for an insect herbivore is critical for metabolic homeostasis, performance, and fitness [6–8]. Although insects present diverse mechanisms to regulate nutrient intake [1], the availability and Int. J. Mol. Sci. 2019, 20, 59; doi:10.3390/ijms20010059 www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2019, 20, 59 2 of 19 nutritional content of food during larval development can directly influence the phenotype of adults [9]. A classic experiment demonstrated that caterpillars of Heliothis zea fed with diets with carbohydrates and proteins in a ratio of 1:4 showed the highest performance [10]. Later, Simpson et al. [11] showed that increases in insect performance is more related with the total abundance of carbohydrates and proteins in food rather than their proportions. In natural ecosystems, insects are often constrained in environments with non-optimal quality of nutrients [12] and natural selection could favor those herbivore insects able to regulate the nutrient intake at the needed proportions and concentrations [1]. The selection of foliage by insect herbivores has been widely discussed during the last decades. Commonly, foliar N and P concentrations have been identified as two of the main factors in foliar selection by insect herbivores; being more selected those plant individuals with higher nutrient abundance [13–17]. However, although plant selection based on specific nutritional requirements may occur in several insect species [18,19], it is still unclear if this is a general behavior [20–23]. Insect herbivores can regulate food intake to maintain a balanced diet for optimal fitness. Larval diet has proven to lead changes not only on larval performance but also on adult reproductive traits and body composition [24]. N concentration and digestibility of food are critical factors determining survival of Lepidoptera larvae, especially during early stages of larval development [18]. In fact, significant changes in larval performance and adult female fitness have been already reported in larvae feeding on different plant varieties of the same species [25]. It is known that insect-predators may include other nitrogen-rich predators in their diets, or even promote cannibalism [26], under situations of N limitation [27]. However, it is still unknown whether insect-prey selection by predators is actually driven by nutritional concentrations. Changes in plant nutritional concentrations could lead to cascade effects through trophic webs due to alteration of population dynamics (food abundance for predators) and the imbalance in nutritional intake by predators [16]. It is thus necessary to determine whether herbivores can maintain their elemental and metabolic homeostasis under changes in the plant source intake to further understand any potential impact on the structure of trophic webs. Ecometabolomics, the study of the ecosystem structure and function through metabolomics analyses [28–33], has proven to be useful to understand the metabolic changes of organisms under the pressure of biotic and/or abiotic stressors [22,31,34–38]. The metabolome is defined as the total set of metabolites present in the organism at a particular abundance and moment [39]. The metabolome responds quickly to environmental fluctuations and stressors [28], and can be thus considered as the chemical phenotype of the organism [39–41]. Metabolomics analyses have demonstrated to be sensitive enough as to detect specific metabolic features upregulated or downregulated between different plant genotypes of the same species under herbivore attack [42], in understanding the success of an insect species feeding in different organs of a same plant [43,44] or even the impact of plant chemical defenses on insect folivores [45]. Metabolomics should thus allow the detection of any alleged metabolic shift in insect herbivores feeding on different plant species or varieties. This metabolomic information can, moreover, provide valuable information about the overall nutritional status of herbivores as most elements (i.e., C, N, P, . . . ) do not act as themselves but as molecular compounds [46]. All the larvae stages of the pine processionary moth (hereafter; PPM), Thaumetopoea pityocampa (Denis & Schiffermüller), remain on the same host, being thus an excellent subject to study plant-insect relationships. Caterpillars of PPM are considered conifer-specialist herbivores since they are able to feed on several species of conifers [20,47,48]. The different stages of caterpillars of PPM elapse from end summer to beginning spring with a peak of folivory in winter [49]. The caterpillar of PPM is considered a pest and constitute a severe problem for many pine populations in the Mediterranean region [20]. Sierra Nevada Natural Park (South-East Spain) is a unique natural study site to examine whether specialist insects vary their overall metabolome structure and chemical composition when feeding on two different plant hosts. Two sympatric subspecies of P. sylvestris, the autochthonous P. sylvestris subsp. nevadensis (hereafter nevadensis) and the introduced P. sylvestris subsp. iberica (hereafter iberica), coexist together in some areas in Sierra Nevada Natural Park [50] since the mid-twentieth century and both subspecies are recently threatened by the PPM [18]. The main aim of this study is to examine Int. J. Mol. Sci. 2019, 20, 59 3 of 19 whether caterpillars of PPM feeding on different plant hosts (subspecies of P. sylvestris (nevadensis and iberica)) present distinct metabolomes and nutritional properties. Although nevadensis and iberica belong to the same species, their overall metabolome structure and nutrient concentrations have proven to be significantly different [23,36]. We hypothesize that metabolomic differences between pine subspecies can lead to changes of insect herbivore metabolomes despite having physiological and metabolic mechanisms to regulate nutrient intake for maintaining body homeostasis. At the same time, we expect finding a direct relationship between overall nutrient concentration (C, N, P, etc.) in caterpillars and their plant hosts. We conducted liquid chromatography coupled to mass spectrometry (LC-MS)-based metabolomics, elemental, and stoichiometric analyses of caterpillars of PPM feeding on nevadensis (hereafter C-nevadensis) and on iberica (hereafter C-iberica), and leaves from attacked and non-attacked pines of each subspecies in Sierra Nevada Natural Park. 2. Results The list of metabolites identified in caterpillars and pines are shown in Table 1. We identified 12 metabolites in pines which were not found in caterpillars; caryophyllene, sabinene, apigenin, catechin, catechol, epicatechin, epigallocatechin, kaempferol, luteolin, taxifolin, δ-tocopherol, and gibberellic acid 3. Table 1. Identified molecular compounds in caterpillar tissues and pine needles by liquid chromatography coupled to mass spectrometry (LC-MS). Averaged Relative Abundance Compound Class Compound Name (Deconvoluted Ion Chromatogram Peak Area) Caterpillars Pine Needles Hexoses 37,234,258 75,422,781 (Glucose, Fructose, Mannose, Galactose, . . . ) Pentoses (Arabinose, Ribose, Xylose, . . . ) 9,983,894 12,228,327 Carbohydrates Deoxy-Hexoses (Deoxy-glucose, 592,288.1 52,9697.9 Deoxy-galactose, Fucose, . . . ) Hexoses-OH (Sorbitol, Mannitol, . . . ) 3,303,206 69,191.52 Pentoses-OH (Xylitol, Arabitol, . . . ) 64,255,464 29,779,498 Alanine 2.55 × 108 14,910,060 Arginine 73,078,590 10,791,686 Aspartic acid 1,060,014 603,470.4 Glutamic acid 21,812,174 9,902,287 Glutamine 18,562,019 35,001,363 Glycine 9,044,821 3,829,374 Isoleucine 4.38 × 109 4.98 × 108 Leucine 55,187,379 1,556,441 Amino acids Lysine 3.12 × 108 5,959,588 Methionine 1.76 × 108 65,552,885 Phenylalanine 2.24 × 109 2 × 108 Proline 1.5 × 109 53,596,908 Serine 10,026,643 7,239,244 Threonine 18,231,142 8,067,142 Tyrosine 2.14 × 108 24,445,586 Tryptophan 2.03 × 109 4.8 × 108 Valine 4.39 × 108 47,193,674 Adenine 1,422,770 41,186,385 Guanine 428,695.9 281,113.5 Nucleobases Thymine 2,393,343 19,099,857 Uracil 5,393,886 3,712,866 Nucleosides Adenosine 1,733,495 1,297,581 Citric acid 4,984,835 76,318,569 Malic acid 6,556,566 93,520,587 Krebs cycle related Pyruvic ac 2,293,721 8,504,367 organic acids Succinic ac 1.68 × 108 740,971.5 α-Ketoglutaric acid 37,118.36 235,423.5 Caryophyllene 13,334,699 Terpenes Sabinene 2,857,821 Int. J. Mol. Sci. 2019, 20, 59 4 of 19 Table 1. Cont. Averaged Relative Abundance Compound Class Compound Name (Deconvoluted Ion Chromatogram Peak Area) Caterpillars Pine Needles Apigenin 21,878,949 Catechin 1.54 × 108 Catechol 39,448,216 D -Pinitol 34,467,901 30,806,094 Epicatechin 219,812.1 Epigallocatechin 376,194.4 Ferulic acid 754,250.3 972,992.9 Phenolics Kaempferol 1,199,232 Luteolin 442,780.3 Phenil-phenol 883,635.4 110,703.6 Quercetin 1,089,097 805,365.7 Quinic acid 52,693,061 4.26 × 108 Taxifolin 8,041,666 Vanillic acid 917,952.3 457,177.3 Abscisic acid 502,407.6 3,853,610 Choline 2.44 × 108 67,441,330 δ-tocopherol 3,376,661 Other Gibberellic acid 3 1,050,073 Vitamin B5 (Panthotenic acid) 38,952,234 1,464,870 Vitamin B2 (Riboflavin) 6,437,087 94,010.85 Shikimic acid 26,206,952 1.27 × 108 Averaged relative abundance of the identified metabolites in caterpillars and/or pines correspond to the average of the peak area of the deconvoluted ion chromatogram of each compound for each group of samples (caterpillars and pines) (see Section 4.7 of Material and Methods for more detailed information). Permutational multivariate ANOVA (PERMANOVA) on caterpillar metabolomes, nutrient, and stoichiometry data showed that herbivores feeding on different P. sylvestris subspecies have significantly distinct overall elemental, stoichiometric, and metabolome composition (p < 0.05) (Table 2). PERMANOVA accounted for about 10.3% of the overall variation between caterpillar groups (C-iberica and C-nevadensis) (Table 2). One-way ANOVAs were performed for each individual variable to test for statistical significance between caterpillar groups. We found that 34.9% of the variables and a 32.1% of identified variables (elemental, stoichiometric and metabolomic features) changed significantly between caterpillar groups at p < 0.05. These percentages of significantly changing variables between caterpillar groups were 18.2% (all variables) and 9.6% (identified variables) for a significance level of p < 0.01. The list of metabolites showing marginal significant changes between caterpillars (p < 0.1) raised to 46.4% for the entire set of variables and 36.5% for the known variables (see Table S1 for univariate analyses of known variables). Table 2. Permutational multivariate analysis of variance (PERMANOVA) of the complete dataset of caterpillars (Caterpillar-Dataset) including the metabolomics and stoichiometric data with caterpillar group (C-iberica and C-nevadensis) as categorical independent factor. Source of Degrees of Sums of Mean Pseudo-F R2 p Value Variation Freedom Squares Square Caterpillar group 1 1.30 × 1019 1.30 × 1019 2.524 0.103 0.0223 Residuals 22 1.13 × 1020 5.14 × 1018 0.897 Total 23 1.26 × 1020 1 Principal component analysis (PCA) performed with caterpillars alone showed clear separation between C-nevadensis and C-iberica along the principal component (PC) 1 (Figure 1). In general, the variable plot of the PCA showed that C-nevadensis had significantly higher N concentration and marginally higher P concentration than C-iberica. C-nevadensis had significantly lower C:N and Int. J. Mol. Sci. 2018, 19, x 5 of 19 Principal component analysis (PCA) performed with caterpillars alone showed clear separation Int. J. Mol. Sci. 2019, 20, 59 5 of 19 between C-nevadensis and C-iberica along the principal component (PC) 1 (Figure 1). In general, the variable plot of the PCA showed that C-nevadensis had significantly higher N concentration and marginally marginally higher P lower significant concentration than C:P ratios C-iberica. than C-nevadensis C-iberica. The overallhad significantly carbohydrate lower abundance relative C:N and marginally significant lower C:P ratios than C-iberica. The overall carbohydrate relative abundance was higher in C-iberica but only hexoses alcohol (hexoses-OH) showed statistical significance. Relative was higher in C-iberica but only hexoses alcohol (hexoses-OH) showed statistical significance. abundances of most amino acids did not change significantly between caterpillar groups with the Relative abundances of most amino acids did not change significantly between caterpillar groups exception of glutamine and valine that showed higher relative abundance in C-nevadensis and with the exception of glutamine and valine that showed higher relative abundance in C-nevadensis isoleucine, leucine, and tyrosine with higher abundance in C-iberica. The relative abundances of all and isoleucine, leucine, and tyrosine with higher abundance in C-iberica. The relative abundances of identified nucleobases were significantly different between caterpillar groups; relative abundances of all identified nucleobases were significantly different between caterpillar groups; relative adenine, thymine, abundances and uracil of adenine, were higher thymine, in C-iberica and uracil whileinguanine were higher C-ibericawas in relative while guaninehigher was inabundance relative in higher C-nevadensis. abundance in C-nevadensis. Succinic acid and α-ketoglutaric acid had significantly higher in Succinic acid and α-ketoglutaric acid had significantly higher relative abundance C-iberica. relative We found several abundance plant related in C-iberica. compounds We found severalin plant caterpillar metabolomes; related compounds D-pinitol, ferulic in caterpillar acid, phenil-phenol, quinic acid, quercetin, vanillic acid, abscisic acid, and shikimic acid. Phenil-phenol metabolomes; D-pinitol, ferulic acid, phenil-phenol, quinic acid, quercetin, vanillic acid, abscisic acid, and vanillic and acid shikimic in higher acid. relative and Phenil-phenol abundance in C-nevadensis. vanillic acid in higher relative abundance in C-nevadensis. Figure Figure Variable 1. 1. Variableand andcasecaseplot plotof ofthe the principal component(PC) principal component (PC)11ofofthe theprincipal principal component component analysis analysis (PCA) (PCA)including includingallallelemental, elemental,stoichiometric, stoichiometric, andand metabolomic metabolomicdata dataofofthe thedataset datasetofofcaterpillars caterpillars (Caterpillar-Dataset). (Caterpillar-Dataset). Only Onlyelemental, elemental,stoichiometric, and identified stoichiometric, and identifiedmetabolomic metabolomicvariables variables areare shown shown in in thethe variable variableplot. plot.Each Eachbarbarrepresents represents the the loading value for loading value for each eachvariable variableononPC1. PC1.Elemental Elemental andand stoichiometric stoichiometric variables variablesare areshown shownin inred. red. Metabolite groupsare Metabolite groups areshown showninindifferent different colors: colors: dark dark blue, blue, saccharides saccharides (carbohydrates);green, (carbohydrates); green,amino aminoacids; acids; light light blue, nucleobases/nucleosides; nucleobases/nucleosides; orange, orange, organic organic acids acids related related toto tricarboxylicacid tricarboxylic acidcycle; cycle;violet, violet, phenolic phenolic compounds; compounds;and anddark darkred, other. red, Caterpillars other. Caterpillars feeding on Pinus sylvestris ssp. iberica (C-iberica) are shown in blue triangles and caterpillars feeding on Pinus sylvestris ssp. nevadensis (C-nevadensis) are shown in red circles. Asterisks and crosses next to elemental, stoichiometric and metabolomic variables denote statistical significance (* p < 0.05) and marginal significance († p < 0.1), respectively. Int. J. Mol. Sci. 2018, 19, x 6 of 19 feeding on Pinus sylvestris ssp. iberica (C-iberica) are shown in blue triangles and caterpillars feeding on Pinus sylvestris ssp. nevadensis (C-nevadensis) are shown in red circles. Asterisks and crosses next Int. J. to elemental, Mol. stoichiometric Sci. 2019, 20, 59 and metabolomic variables denote statistical significance (* P < 0.05) and6 of 19 marginal significance († P < 0.1), respectively. The relative The relative differences differences(%)(%)ofofelemental elementaland andstoichiometric stoichiometric variables variables of of C-iberica C-iberica respect respect to to C- C-nevadensis were +1.4% for C (not significant; hereafter N.S.), − 3.5% for N (p = 0.023), nevadensis were +1.4% for C (not significant; hereafter N.S.), −3.5% for N (P = 0.023), −9.3% for P (P = − 9.3% for P (p = 0.072), 0.072), +9.05% +9.05% for C:Nfor(P C:N (p = +12.43% = 0.039), 0.039), +12.43% for C:P for (P =C:P (p = +3% 0.076), 0.076), for+3% N:P for N:P−1.87% (N.S.), −1.87% (N.S.), for N:K for N:K (N.S.), and +5.01% for K:P (Figure 2). From the plant side, not-attacked (N.S.), and +5.01% for K:P (Figure 2). From the plant side, not-attacked trees (hereafter NATs) of trees (hereafter iberica of ibericasignificantly NATs)showed showed significantly higher P and higher lowerP C:P, and and lower C:P, and marginally marginally significant significant lower lower N:P than N:P NATs than of NATs of nevadensis. nevadensis. Elements and Elements and stoichiometric stoichiometric variables didvariables did not not differ differ significantly significantly between between attacked attacked branches of attached trees (hereafter AT.ABs) of both subspecies while branches of attached trees (hereafter AT.ABs) of both subspecies while not-attacked branches ofnot-attacked branches of attacked attacked trees trees (hereafter (hereafter AT.NABs) AT.NABs) ibericahad ofofiberica hadlower lowerN:PN:Pand andmarginally marginallylower lower C:P C:P ratios ratios and and higher P concentration than AT.NABs of higher P concentration than AT.NABs of nevadensis.nevadensis. Figure 2. 2. Relative Relative difference difference (%) (%) of of elemental elemental and and stoichiometric stoichiometric variables variables between between caterpillar caterpillar groups groups (C-iberica relative to C-nevadensis) and between between pine pine subspecies subspecies for for each each folivory folivory level level(iberica (ibericaNATs NATs relative to relative nevadensis NATs; to nevadensis iberica AT.ABs NATs; iberica nevadensis AT.ABs; AT.ABs relative to nevadensis iberica AT.NABs AT.ABs; iberica AT.NABs relative to nevadensis AT.NABs). Positive values indicate higher elemental concentration nevadensis concentration or stoichiometric stoichiometric ratio value in C-iberica or iberica trees while negative values indicate higher elemental concentration value concentration or or stoichiometric ratio stoichiometric ratio value value in C-nevadensis or nevadensis trees. Asterisks Asterisks and and crosses denote statistical p < 0.05) and marginal significance (†† P significance (* P p < 0.1), respectively. Relationships between Relationships betweentotal totalcarbohydrate carbohydrate abundance abundance and and amino amino acidacid abundance abundance vs.N,C, vs. C, N, C:N, C:N, and, C:P, C:P, and, N:P ratios N:P ratios in caterpillars in caterpillars were were performed performed (Figure(Figure 3). C 3). wasC correlated was correlated solelysolely to to total total carbohydrate abundance and N was correlated solely to total amino carbohydrate abundance and N was correlated solely to total amino acid abundance. Our results acid abundance. Our results showed showed different different tendencies tendencies between between caterpillar caterpillar groups; groups; correlation correlation carbohydrates carbohydrates vs C,vs. C, C:N, C:N, C:P, C:P, and and N:P N:P werewere positive positive andand significant significant (marginally (marginally significantfor significant forC:N) C:N)when when considering considering both both caterpillar groups caterpillar groupstogether. together.For ForC-iberica C-iberica alone, wewe alone, found foundsignificant positive significant correlations positive between correlations total between total carbohydrates and C, C:N, and C:P and marginally positive correlation with N:P. However, not carbohydrates and C, C:N, and C:P and marginally positive correlation with N:P. However, we did we findnot did significant correlation find significant betweenbetween correlation total carbohydrates and C, C:N, total carbohydrates andand C:P in C, C:N, andC-nevadensis, only N:P C:P in C-nevadensis, was marginally only positivelypositively N:P was marginally correlatedcorrelated with total carbohydrates (Figure 3). (Figure with total carbohydrates On the other 3). Onhand, correlation the other hand, between total correlation amino acids between and N was total amino acidsnegative and N when considering was negative whencaterpillar groups considering together groups caterpillar and for C-nevadensis. together Total and for amino acids vs. C-nevadensis. N correlation Total amino acids wasvs.marginally positively N correlation was correlated marginallyinpositively C-iberica. Total amino acids in C-iberica were negatively correlated with C:P (p < 0.1) and N:P and correlated in C-iberica. Total amino acids in C-iberica were negatively correlated with C:P (P < 0.1) no significance wereN:P and found andwith C:N. For C-nevadensis, no significance were found we withfound C:N. positive correlationwe For C-nevadensis, between found total amino positive acids vs. correlation C:N and C:P but no significance was found against N:P ratio (Figure 3). between total amino acids vs C:N and C:P but no significance was found against N:P ratio (Figure 3). Int. J. Mol. Sci. 2019, 20, 59 7 of 19 Int. J. Mol. Sci. 2018, 19, x 7 of 19 Figure Pearson lineal Figure3.3. Pearson linealcorrelations correlationsbetween between total total carbohydrate carbohydrate abundance abundance versus versus C, C:N, C, C:N, C:P, C:P, and and N:PN:P ratios, ratios, andand total total amino amino acidacid abundance abundance versus versus N,N, C:N, C:N, C:P,and C:P, andN:P N:Pfor forboth bothcaterpillar caterpillargroups groups (C-nevadensis (C-nevadensisand andC-iberica) C-iberica)separately separately (n (n == 12) 12) and and together together (n = 24). Coefficient Coefficient ofof determination determination 2 ), (R(R 2),statistic statisticFisher Fishervalue value (F), (F), and and pP value value (p) (P) are are shown. shown. DataData corresponding corresponding to to caterpillars caterpillars feeding feeding ononPinus Pinussylvestris ssp.iberica sylvestrisssp. iberica(C-iberica) (C-iberica) are are shown shown in blue triangles and caterpillars caterpillars feeding feedingon Pinus onPinus sylvestris sylvestrisssp. nevadensis (C-nevadensis) ssp. nevadensis (C-nevadensis)are areshown shown in in redred circles. circles. Lineal Lineal correlations correlations for C-iberica for C-iberica and and C-nevadensis C-nevadensis are are shownshown in blue in blue andand red,red, respectively. respectively. Lineal Lineal correlations correlations for both for both C-iberica C-iberica andand C- C-nevadensis nevadensis are areshown shownininblack. black.Correlations Correlations with P >P 0.1 with > 0.1 areare represented represented byby dashed dashed lines. lines. The ThePCA PCAperformed performedincluding including caterpillar caterpillar and and pines pines together showed showed clearclear clustering clustering between between caterpillars caterpillarsand andpines pinesalong alongthe thePC1 PC1(Figure (Figure 4b). 4b). In addition, the PC2 of of the the PCA PCA clustered clusteredboth bothpine pine subspecies (nevadensis and iberica) and caterpillar subspecies (nevadensis and iberica) and caterpillar groups (C-nevadensis and C-iberica). Interestingly, (C-nevadensis and C-iberica). Interestingly, C-nevadensis C-nevadensis cases cases clustered clustered closer to nevadensis closer to nevadensis pines than to ibericaiberica pines, pines, andand C-iberica C-iberica cases cases clustered clusteredcloser iberica pines closer to iberica pinesthan nevadensis thantotonevadensis along along thethe PC2PC2 (accounting (accounting for8.63% for the the 8.63% of theof the total total variance) variance) showing showing thusthus a certain a certain degree degree of ofmetabolic metabolicresemblance resemblancebetweenbetweenplant planthosts hosts and and herbivores.One-way herbivores. One-way ANOVA ANOVA andand Tukey’s Tukey’s HSDHSD post-hocpost-hoc test detected test detected significant significant differences differences between between all analyzed all groups analyzed ofgroups samples of along samplesthealong PC2, the beingPC2, being nevadensis nevadensis and ibericaandtheiberica most the most distant distant groups of samples (Figure 4b). Caterpillars were also significantly clustered along the PC2 (F =68.49; groups of samples (Figure 4b). Caterpillars were also significantly clustered along the PC2 (F = 68.49; p P< <0.0001). 0.0001).InInaddition, addition,one-way one-way ANOVA ANOVA on on the the metabolomic metabolomic distance between between caterpillars caterpillars and and pinescalculated pines calculatedalong alongthethePC2 PC2of ofthe thePCA PCAshowed showed the the “C-iberica “C-iberica vs. iberica” iberica” distance distance smaller smallerthan than the“C-iberica the “C-iberica vs. vs. nevadensis” nevadensis”distance distance(Figure (Figure4c). The 4c). same The pattern same waswas pattern found for C-nevadensis; found for C-nevadensis;“C- nevadensis vs.vs. “C-nevadensis nevadensis” nevadensis” distance distancewaswassmaller smallerthan than“C-nevadensis “C-nevadensisvs. vs.iberica” iberica”distance. distance.InIngeneral, general, caterpillartissues caterpillar tissueshad hadhigher higher concentration concentration of of N, N, P, and K, P, and K, together together with with higher higherN:K N:Kandandlower lowerC:N, C:N, C:P, and K:P ratios than pines. Few of the plant-related metabolites present in C:P, and K:P ratios than pines. Few of the plant-related metabolites present in caterpillars (vanillic caterpillars (vanillic acid,and acid, andphenyl-phenol) phenyl-phenol)were were found found atat higher higher relative relative abundance abundance thanthan inin plant plant tissues. tissues.InIngeneral, general, caterpillars had higher relative abundance of amino acids and lower caterpillars had higher relative abundance of amino acids and lower abundance of phenolic compoundsabundance of phenolic compounds than pines. than pines. Int. J. Mol. Sci. 2019, 20, 59 8 of 19 Int. J. Mol. Sci. 2018, 19, x 8 of 19 Figure4.4. Variable Figure Variable(a)(a)and andcase plotplot case (b) of (b)the ofprincipal component the principal (PC) 1 vs. component the1PC2 (PC) vs.ofthe the PC2 principal of the component analysis (PCA) including all elemental, stoichiometric and metabolomic principal component analysis (PCA) including all elemental, stoichiometric and metabolomic data data of the of dataset of caterpillars and pines (Pine-Caterpillar-Dataset). Only elemental, stoichiometric the dataset of caterpillars and pines (Pine-Caterpillar-Dataset). Only elemental, stoichiometric and and identified identified metabolomic metabolomic variables variables areare shown shown in variable in the the variable plotThe plot (a). (a).location The location of a variable of a variable (element, (element, stoichiometric stoichiometric ratio or metabolite) ratio or metabolite) on theplot on the variable variable planeplot plane represents represents the value the loading loadingonvalue PC1 vs. on for PC2 PC1such vs PC2 for such specific specific variable. variable. and Elemental Elemental and stoichiometric stoichiometric variables variables are shownare inshown in red. red. Metabolite Metabolite groups groups are shown in are showncolors: different in different colors: dark blue, dark(carbohydrates); sugars blue, sugars (carbohydrates); green,light green, amino acids; aminoblue, acids; light blue, nucleobases/nucleosides; orange, organic acids related to tricarboxylic nucleobases/nucleosides; orange, organic acids related to tricarboxylic acid cycle; violet, phenolic acid cycle; violet, phenolic compounds; dark compounds; darkplot red, other. Case red,(b) other. shows Case plot (b) shows caterpillars feedingcaterpillars feeding on on Pinus sylvestris ssp.Pinus iberica sylvestris ssp. iberica (C-iberica) in blue circles and caterpillars feeding on Pinus (C-iberica) in blue circles and caterpillars feeding on Pinus sylvestris ssp. nevadensis (C-nevadensis)sylvestris ssp. in nevadensis (C-nevadensis) in red circles. Needle samples of attacked branches (AT.ABs) and non- red circles. Needle samples of attacked branches (AT.ABs) and non-attacked branches from attacked attacked branches from attacked trees (AT.NABs) are shown in dark and light blue for Pinus sylvestris trees (AT.NABs) are shown in dark and light blue for Pinus sylvestris ssp. Iberica, and dark and light ssp. Iberica, and dark and light red for Pinus sylvestris ssp. nevadensis. Different letters next to each red for Pinus sylvestris ssp. nevadensis. Different letters next to each colored arrow denote significant colored arrow denote significant difference between groups of samples after Tukey’s HSD post-hoc difference between groups of samples after Tukey’s HSD post-hoc test (p < 0.05). One-way ANOVA test (P < 0.05). One-way ANOVA plot (c) contrasting the metabolomic distances along the PC2 of the plot (c) contrasting the metabolomic distances along the PC2 of the PCA between caterpillars and pines. PCA between caterpillars and pines. Different letters indicate significant difference after Tukey’s HSD Different letters indicate significant difference after Tukey’s HSD post-hoc test (p < 0.05). Not attacked post-hoc test (P < 0.05). Not attacked trees are not represented in the PCA. Detailed discussion of the trees are not represented in the PCA. Detailed discussion of the effects of the caterpillar on pines effects of the caterpillar on pines metabolism has been already published elsewhere [23,36]. metabolism has been already published elsewhere [23,36]. Int. J. Mol. Sci. 2019, 20, 59 9 of 19 3. Discussion Our analyses found evidence supporting the hypothesis that caterpillars have different metabolomes and elemental composition when eating in different subspecies of P. sylvestris (Figure 1) and such difference accounted for over 10% of the total variability between caterpillar groups (Table 1). As previously shown [23,36], pine subspecies and different folivory levels (i.e., AT.ABs vs. AT.NABs) varied in their overall metabolomic structure. At the elemental level, we found C-iberica nutritionally poorer compared to C-nevadensis having higher C:N and C:P and lower total N and P concentrations (Figures 1 and 2). However, the elemental and stoichiometric differences between caterpillar groups were not directly related to the nutritional concentrations of pines they fed on (Figure 2). In fact, NATs of iberica trees had slightly higher total concentration of P and lower C:P ratio compared to nevadensis suggesting that iberica trees could be more attractive than nevadensis from a nutritional point of view for oviposition by adult PPM females although unselective oviposition for PPM moths have been previously reported [18]. Nevertheless, once trees are attacked by PPM, the differences between subspecies in foliar P concentration and C:P ratios lost significance (p < 0.1 for AT.NABs and N.S. for AT.ATBs) (Figure 2) suggesting thus a significant adjustment of nutrients in iberica in response to folivory attack. This result clearly showed that nutrient transfer through trophic webs is not necessarily directly related with the nutritional concentration of food source and other factors such as the plant metabolome could play critical roles regulating nutrient transfer across trophic web levels. For example, folivory rates can change depending on nutrient concentration in food source, being those herbivores feeding on sources with high C:N and C:P ratios the most active to cope with nutrient demand [51]. The overall metabolomic differences between pine subspecies, as previously detailed and discussed elsewhere [23], could potentially cause nutritional imbalances in herbivores, forcing them to change their folivory rates to keep their nutritional requirements. Several lines of evidence proved that organismal and environmental C:N, N:P, and C:P ratios can play critical roles determining the ecosystem structure and function [14,52,53]. Therefore, significant differences in nutrient concentration between caterpillar groups (C-nevadensis with higher N, P concentrations and lower C:N and C:P ratios than C-iberica; Figure 2) could lead to significant changes in the population dynamics of the PPM, a species known for having episodic outbreaks depending on different factors such as food quality. The number of nevadensis trees in Sierra Nevada is recently much lower compared to iberica due the intensive felling during the twentieth century [20]. The replacement of autochthonous trees by iberica may cause significant changes in overall nutrient composition of PPM populations which could lead to shifts in population dynamics and consequently impact the trophic web structure [53–55]. Potassium (K) is an essential element for insects due to its high concentration in the hemolymph [56] and we did not find significant changes between caterpillar groups (Figures 1 and 2). In plants, K is an important element to maintain a good performance of the physiological processes [57,58], especially at low water availability [57,59]. We found higher K concentration in caterpillars compared to pine needles (Figure 4) suggesting an accumulation of K in insect tissues, as expected with certain nutrients through trophic web levels. However, we found higher K:P and lower N:K ratios in pines compared to caterpillars (Figure 4) suggesting that the requirement of K in pines in proportion to other critical nutrients (N and P) is significantly higher than in caterpillars. This high difference is likely by the distant physiological functional traits between plants and insects. For example, K in trees play a critical function maintaining a large osmolite gradient from roots to leaves in order to sustain the cellular and physiological homeostasis [57]. Contrasting the elemental stoichiometry of plants vs. insects as shown here highlight the important role that K plays in plant physiology and should be considered in plant stoichiometric studies, especially under the pressure of stress conditions [29,57,59]. Among all the identified carbohydrates, only the group of hexoses-alcohol changed significantly between caterpillars (Figure 1). However, C-iberica had higher overall carbohydrate abundance than C-nevadensis (Figure 1) which could partly explain their higher tissue C:N and C:P ratios. Carbohydrates are C-rich compounds and, with exception of amino-sugar species, do not contain Int. J. Mol. Sci. 2019, 20, 59 10 of 19 heteroatoms and are solely composed of C, H, and O [46]. Relationships of overall carbohydrate abundance vs. C concentration, C:N, and C:P ratios in caterpillars were significant when both groups were analyzed together (Figure 3) showing a significant body stoichiometry-metabolome relationship in insects, as previously reported in plants [29]. Such relationships were also significant for C-iberica considered alone but not for C-nevadensis, suggesting higher C allocation into carbohydrate metabolism in C-iberica; metabolic pathways generally related to rapid growth [29]. A general tendency in carbohydrate increases in relation to other molecular compounds could generate a nutritional imbalance in C-iberica and impact the trophic web structure by altering the population dynamics of PPM and the nutrient intake of their predators [26,27]. In addition, the overall increase of carbohydrates in C-iberica could be directly linked with the higher relative abundance of tricarboxylic acid cycle (TCA) organic acids (Figure 1). Higher abundances of TCA-related compounds and carbohydrates suggest higher respiratory metabolism in C-iberica with probably higher energy generation through the acetate oxidation directly derived from carbohydrates. However, our metabolomics results alone cannot definitely confirm this mechanism and additional analyses, i.e., proteomics, would help to decipher how those two groups of caterpillars assign C from carbohydrates to specific metabolic pathways. Interestingly, N:P body ratio did not significantly change between caterpillar groups (Figures 1 and 2) but it was positively correlated with total carbohydrate abundance (Figure 3) even though the carbohydrates identified in our study do not contain heteroatoms. This result would be consistent with the expected in the frame of growth rate hypothesis when lower N and P concentrations and higher C:N and C:P ratios is expected to be related with higher N:P ratios [14,16]. Only a few amino acids changed significantly between caterpillar groups (Figure 1). However, we found that total abundance of amino acids was negatively correlated with N and positively with C:N body ratios of caterpillars when analyzing both caterpillar groups together and in C-nevadensis (Figure 3) suggesting that a large part of N concentration in C-nevadensis may be part of non-soluble structural proteins or other N-rich compounds such nucleotides as guanine and cytosine with five and three N atoms, respectively. Guanine is paired with cytosine and adenine with thymine/uracil and our analyses showed that C-nevadensis had higher relative abundance of N-rich nucleotide pairs (Guanine-Cytosine with eight atoms of N per pair) and lower relative abundance of adenine-thymine/uracil pairs (seven atoms of N per pair) than C-iberica. This result could partially explain the higher N concentration in C-nevadensis and suggests that the differences between nevadensis and iberica pines at both elemental and metabolomic levels, although similar still significantly different between them [36], could change the genome expression and protein synthesis in insect herbivores. From an evolutionary perspective, those changes could be selected trough several generations favoring different DNA pairs depending on the environmental (or food source in this case) N:P ratio as previously shown in plants under poor N availability [60]. Due to practical limitations in ecometabolomics research in general, and the complex logistics of this study, caterpillars were processed with the digestive content included. The dry weight of Lepidoptera larvae gut content can reach up to 20% of total caterpillar biomass when full [61]. We acknowledge this methodological limitation in our study which could lead to a bias in our results. However, the effects of temperature on feeding activity of the caterpillar of the PPM is well known and it only occurs when night temperatures are >0 ◦ C and the temperature of the previous day was >6 ◦ C [49,62,63]. Actually, some predictive models have been already developed to predict the feeding activity of the caterpillars of the PPM. For example, Battisti et al. (2005) [49] proposed a model to understand the feeding time of PPM caterpillars in areas of expanding range by using a threshold temperature of 0 ◦ C combined with a more conservative activation temperature threshold of 9 ◦ C, instead of 6 ◦ C. According to the meteorological station at La Cortijuela (~1600 m a.s.l), night temperatures were below −5 ◦ C for at least three consecutive days before caterpillar sampling and day temperatures did not reach +5 ◦ C (see Material and Methods section). Therefore, caterpillars probably did not go out of their winter tents to feed for at least few days before they were sampled (4 March 2011). We even expect lower temperatures in the study site at Collado the Matasverdes as it Int. J. Mol. Sci. 2019, 20, 59 11 of 19 is located over 300 m. a.s.l. from La Cortijuela. Given the non-optimal meteorological conditions for PPM caterpillars at the moment of sampling, we assume that caterpillars were already in a status of semi-starvation and the majority of the analyzed biomass corresponded to the insect and not to the content of guts. Furthermore, we also expect that most molecular components in the gut would have been substantially modified through the different digestion processes if not already assimilated or excreted by the caterpillars. We detected several plant-related compounds in the metabolome of caterpillars such as polyphenolics (Figure 1 and Table 1) which could be directly detected from the gut content, however, the relative abundance of compounds such as phenil-phenol and vanillic acid was higher in caterpillar tissues than in pines (Figure 4). If caterpillars of the PPM digested such compounds or had specific mechanisms to excrete them [64], the detected signal in larvae tissues would never be larger than in pine needles, the original source, thus suggesting bioaccumulation in insects [65]. It is known that insects have higher levels of certain vitamins than plants [66,67]. We found higher relative abundance of some vitamin-b complex compounds (choline, vitamin B5 (panthothenic acid), and vitamin B2 (riboflavin)) in caterpillars compared to pine needles (Figure 4). The requirement of vitamins for insects is well known, especially the B-complex group [68] which play important roles in cell metabolism. The relative abundance of riboflavin was higher in C-nevadensis compared to C-iberica, and in addition to its important role in the carbohydrate, lipid, and amino acid metabolism, it is also involved in the activation of defensive responses in plants [69,70] and animals [71] against biotic stresses (fungal, bacterial, and viral pathogens). Nevertheless, those responses are commonly species-specific and may vary depending on the infectious agent [70]. The potential sequestration of plant secondary metabolites (i.e., iridoids) by folivores to be used against parasitoids has been already suggested [72]. Glucosinolates are an example of defensive metabolites present in abundance in Brassicaceae plants that are accumulated and used by herbivores against parasitoids [73]. We also found significant higher accumulation of riboflavin and other compounds such as vanillic acid and phenil-phenol in C-nevadensis (Figure 1) that could be used by caterpillars for their own defense against parasitoids. Even so, further focused research is necessary to understand whether caterpillars of PPM accumulate and use specific plant secondary compounds against biotic stressors. In summary, although animals have complex physiological mechanisms to regulate their metabolism and maintain homeostasis and fitness, caterpillars of the PPM feeding on distinct pine subspecies presented significantly distinct overall stoichiometry and metabolome structure accounting for over 10% of their total variability (Figure 1, Table 2). In addition, although plants and insects have very distinct metabolome structure, metabolomes of caterpillars presented a certain resemblance to the metabolome of the pine subspecies where they fed on (Figure 4b,c), suggesting thus that the overall insect herbivore composition is linked to the metabolome of their hosts. Nutrient transfer (i.e., N and P) trough trophic levels is not necessarily a net sum from the abundance in food sources. The metabolites where N and P are assigned to can play critical roles determining the elemental composition of higher trophic levels. From the ecological point of view, the stoichiometric and metabolomic shifts of caterpillars feeding on iberica pines could significantly impact the nutrient transfer along trophic webs by altering the PPM population dynamics and nutritional composition of larvae, causing structural and functional changes at higher trophic levels due to the larger abundance of iberica trees compared to the autochthonous nevadensis. This exploratory study using elemental, stoichiometric and metabolomic data to understand the relationship between Pinus syslvestris and caterpillars of the processionary moth raised some questions in the field of plant-insect interactions research that warrant an accurate and focused investigation: 1. Do insect herbivores have the ability to accumulate specific plant secondary metabolites to repeal parasitoids? If this ability is found in multiple insect groups, is it consequence of common ancestry or is consequence of multiple independent acquisitions (phenotypic convergence)? 2. Why nutrient concentration in insect herbivore tissues is not necessarily correlated with the nutrient abundance of their food source? Are there molecular compounds in the food source that Int. J. Mol. Sci. 2019, 20, 59 12 of 19 can potentially impact the nutrient transfer across different trophic levels by forcing changes in herbivore physiology? 3. How flexible are specialist insects feeding in diverse plant species? Which are the key plant molecular factors constraining the food diversity of specialist herbivores? 4. Material and Methods 4.1. Study Site Sampling was conducted in early-March 2011 (late-winter) in Collado de Matasverdes (37.05◦ N, 3.27◦ W; 1900 m a.s.l.), one of the sites where the subspecies of P. sylvestris nevadensis and iberica coexist [46] in Sierra Nevada National Park (Granada, SE Spain). The climate is Mediterranean, with cold winters and hot summers with usually severe drought. The mean annual precipitation is 945 mm and the mean annual temperature is 9.8 ◦ C. January is the coldest month with an average monthly minimum temperature of −0.1 ◦ C and July is the warmest month with an average maximum temperature of 30.1 ◦ C. Rainfall is concentrated mainly in spring and autumn [74]. The maximum temperatures recorded in the meteorological station of La Cortijuela (~1600 m a.s.l; 37.08◦ N, −3.47◦ W) for the sampling period (1–4 March 2011) were 13.9, 2.0, 2.9, and 2.9 ◦ C, respectively. The minimum temperatures for the sampling period were −7.3, −7.3, −6.7, and −5.9 ◦ C, respectively. Daily temperatures after 5 pm were always below 0 ◦ C during the sampling period. Pines and caterpillars were collected on 2nd and 4th March, respectively. 4.2. Experiment Design and Sampling of Needles and Caterpillars Twenty-four adult iberica and nevadensis trees, >5 m in height and >45 years old, were semi- randomly selected as study cases (total n = 48). Twelve of those trees of each subspecies had no signs of caterpillar attack (NATs), and the other 12 trees had caterpillars of PPM on the canopy, easily located by their winter tents (2–4 per tree). From each affected tree (iberica and nevadensis), one of the winter tents of PPM were removed with a pole, opened and 30–50 caterpillars were collected, pooled, packed in paper bags and quickly frozen in liquid nitrogen. According to previous studies, no difference in fitness has been observed between caterpillars feeding on iberica and nevadensis [18]. All sampled caterpillars corresponded to their final instar (5th instar) of their biological cycle (~5–6 cm. long). In addition, 60–100 needles were collected from the NATs, from the attacked branches of the attacked trees (AT.ABs; branches with caterpillar tents) and from the non-attacked branches of the attacked trees (AT.NABs; branches without caterpillar tents). All plant samples were also packed and rapidly frozen in liquid nitrogen. To be consistent and due to the significant physiological differences between sunlit and shade leaves [75], all needles collected from NATs, AT.ABs and AT.NABs were all oriented to the South (sunlit). Winter tents were also south-faced. All samples were collected during a narrow window of time (11:00 am to 3:00 pm) under constant environmental conditions to avoid differences caused by plant diurnal rhythms. 4.3. Caterpillar and Foliar Processing for Elemental and Metabolomics Analyses Briefly, frozen caterpillars and pine needles were lyophilized and subsequently ground with a ball mill operating at 1600 rpm (~30 Hz) for 8 min. The fine homogeneous powder produced from caterpillars and pine needles was stored at −80 ◦ C until the extraction of the metabolites for liquid chromatography coupled to mass spectrometry (LC-MS) analyses. 4.4. Elemental Analysis C and N concentrations from caterpillar and pine needle samples were determined with a CHNS-O Elemental Analyzer (EuroVector, Milan, Italy). For each analysis, 1.4 mg of dried sample powder was used. P and K concentrations from caterpillars and pine needles were determined by acid digestion extraction in a microwave reaction system under high temperature and pressure [76] and analyzed Int. J. Mol. Sci. 2019, 20, 59 13 of 19 by an ICP-OES (Optic Emission Spectrometry with Inductively Coupled Plasma) (Perkin-Elmer Corporation, Norwalk, CT, USA). Briefly, 250 mg of dried sample powder were added into a Teflon tube. Each tube subsequently received 5 mL of nitric acid and 2 mL of H2 O2 and samples were digested in a MARSXpress microwave reaction system (CEM, Mattheus, NC, USA). The digestions were transferred into 50-mL glass flasks, resuspended in nanopure water to a final volume of 50 mL and analyzed by ICP-OES (Perkin-Elmer Corporation). 4.5. Extraction of Metabolites for Liquid Chromatography-Mass Spectrometry (LC-MS) Analyses Polar and semi-polar metabolites from caterpillars and pine needles were extracted as described elsewhere [77] with minor modifications. Briefly, two sets of 2 mL centrifuge tubes were labeled: set A was to perform the extraction and set B was to keep the extracts from set A. For each sample, 30 mg of caterpillar powder and 100 mg of pine powder were weighted into their corresponding tube of set A. Subsequently, one mL of MeOH/H2 O (80:20) was added to each tube, vortexed for 15 min, sonicated at 24 ◦ C for 5 min and centrifuged at 15,000× g for 5 min. Thus, 0.6 mL of the supernatant (extract) from each tube of set A was transferred to the corresponding tubes of set B. This procedure was performed twice to the tubes of set A to perform two extractions on the same sample. The tubes of set B were centrifuged at 15,000× g for 5 min and each extract was collected by crystal syringes, filtered through 0.22 µm pore syringe microfilters and transferred to a labeled HPLC vial. Vials with extracts were stored at −80 ◦ C until the LC-MS analysis. 4.6. LC-MS Analysis Separation of metabolites was performed by LC with a reversed-phase C18 Hypersil gold column (150 × 2.1 mm, 3 µ particle size; Thermo Scientific, Waltham, MA, USA) and an Ultimate 3000 HPLC system (Thermo Fisher Scientific/Dionex RSLC, Dionex, Waltham, Massachusetts, USA). Chromatography operated at constant temperature of 30 ◦ C at a flow rate of 0.3 mL per minute. For each sample, 5 µL were injected. We used water (0.1% acetic acid) (A) and acetonitrile (B) as mobile phases. Both A and B were previously filtered and degassed for 10 min in an ultrasonic bath. The elution gradient was set to begin at 90% A (10% B), maintained for 5 min and then changed to 10% A (90% B) during the next 15 min, and held for 5 minutes. The elution gradient was thus gradually recovered linearly to the initial conditions (90% A; 10% B) over the next 5 min. The chromatographic column was washed and stabilized for 10 min before injecting the next sample. LC was coupled to an LTQ Orbitrap XL high-resolution mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA) equipped with an HESI II (heated electrospray ionization, Thermo Fisher Scientific) source for mass spectrometry analyses. All samples were injected twice; once with the HESI operating in positive ionization mode (+H) and once in negative ionization mode (–H). The mass spectrometer operated in FTMS (Fourier Transform Mass Spectrometry) full-scan mode with high-mass resolution (60,000) and mass range acquisition of 50–1000 m/z. A caffeine standard was injected every 15 samples to monitor the resolution and sensitivity of the spectrometer. Experimental blank samples were analyzed every 10 samples for sample background determination. 4.7. Processing of LC-MS Chromatograms Chromatograms of both positive and negative ionization modes were analyzed separately. RAW files from caterpillars and attacked trees were processed together. Metabolomes of NATs were not included in this study as have been already discussed in previous studies [23,36]. Solely the elemental composition and the C:N:P:K stoichiometry was used from NATs in this study. MZmine 2.14.2 was used to process the raw data files obtained from the mass spectrometer [78]. All chromatograms were baseline corrected, deconvoluted, retention time normalized, aligned, and metabolic features were assigned to specific metabolites (See Table S2 for full parameter details). Metabolites were automatically annotated to specific metabolites by matching the detected features with the total exact mass and retention time of our in-house library generated from the measurements of hundreds standards in Int. J. Mol. Sci. 2019, 20, 59 14 of 19 positive and negative ionization modes. This metabolite assignation method corresponds, thus, to a second level of identification confidence [79]. Numerical datasets were exported and the different identified variables corresponding to the same molecular compound were summed to end up with a single variable per identified metabolite (see Rivas-Ubach et al. [23] for details). The numerical values of the variables of the datasets correspond to the deconvoluted peak areas of the chromatograms detected by the spectrometer which are directly proportional to the abundance of the metabolic feature. For this reason, we used the term relative abundance hereafter when referring to differences in the amount of metabolites between groups of samples. 4.8. Statistical Analyses After processing the metabolomics RAW files with caterpillars and pines together, we coupled the nutrient, stoichiometry and metabolomics data into a single dataset and it was subsequently filtered through the following steps: (1) All zero values of the dataset were replaced for missing data (NA). (2) For each variable (metabolite feature, element and stoichiometric ratio), outlier values were determined for each cell individually. A cell corresponds to each group of samples defined by the combination of each factor and level. Therefore, our study is mainly composed by 6 cells: Caterpillars feeding on iberica (C-iberica), caterpillars feeding on nevadensis (C-nevadensis), AT.ABs of iberica, AT.NABs of iberica, AT.ABs of nevadensis and AT.NABs of nevadensis. Since only stoichiometric and elemental data was used from NATs of both pine subspecies, this dataset was filtered separately. Detected outliers were replaced for NA and were defined as: U pper Outliers → value > Q75 + 2 × IQR Lower Outliers → value < Q25 − 2 × IQR (1) where Q75 represents the third quartile, Q25 represents the first quartile and IQR is the interquartile range (IQR = Q75–Q25) of each variable and cell. (3) Variables with less than 60% of data within all cells were removed from the dataset. (4) Variables with signal to noise lower than 15 (determined by the signal from blanks analyzed during the sequence) were removed from the dataset. After the initial dataset filtering, two definitive datasets were generated for statistical analyses: Caterpillar-Dataset and Pine-Caterpillar-Dataset. The Caterpillar-Dataset included the elemental, stoichiometric and metabolomics data of caterpillars alone (C-iberica and C-nevadensis). This dataset consisted in a total of 9494 variables; 9 of them corresponded to elemental and stoichiometric variables (C, N, P, K, C:N, N:P, C:P, N:K, and K:P) and the remaining 9485 variables were metabolic features which 44 of them were identified by our in-house metabolite library (Table 1). Additionally, for the Caterpillar-Dataset, we generated two additional variables by summing all the peak areas corresponding to the identified sugars (total carbohydrate abundance) and amino acids (total amino acid abundance), separately. Those two variables were posteriorly correlated with elemental and stoichiometric variables (Figure 3). The Pine-Caterpillar-Dataset included the data from the attacked-pines and caterpillars together consisting in the same 9 elemental and stoichiometric variables and a total of 9643 metabolomic features, 55 of them identified (Table 1). For this dataset, the values for each metabolomic feature of each sample of the Pine-Caterpillar-Dataset was scaled by the total intensity of its chromatogram to allow comparisons of the relative abundance of compounds between pines and caterpillars. This study mainly focuses on the relationship between plant-insect stoichiometry and metabolome structure, therefore, the metabolomics datasets were processed considering exclusively the original individual food source PPM: the attacked trees. Therefore, the metabolomics data of NATs were not analyzed in this study and only elemental and stoichiometric data were slightly discussed. The results focused on the comparison of the metabolomic changes between pines attacked by the caterpillars of the PPM were already discussed in previous publications [23,36,80]. Int. J. Mol. Sci. 2019, 20, 59 15 of 19 To test for significant differences between metabolomes and stoichiometry of caterpillars feeding on different pine subspecies, a PERMANOVA test was conducted on the Caterpillar-Dataset using the Euclidian distance with the caterpillar group as fixed factor and setting the number of permutations at 10,000 (Table 2). To represent the variability between C-iberica vs. C-nevadensis, the Caterpillar-Dataset was subjected to principal component analysis (PCA) (Figure 1), one of the most used ordination analyses in metabolomics studies [81–83]. Additionally, all elemental, stoichiometric and identified metabolomic feature were individually submitted to one-way ANOVAs to test for differences between caterpillar groups (Table S1). An additional PCA with the Pine-Caterpillar-Dataset was performed (Figure 4) to plot together caterpillars and pines in the same multidimensional space and detect trends between cases of caterpillars and pines along the first two PCs. The coordinates of pines and caterpillars along the PC2 were subjected to one-way ANOVA to test for significant clustering between groups of samples along the axis. Additionally, the metabolomics distances of each caterpillar sample (C-nevadensis and C-iberica) with each pine subspecies sample (iberica and nevadensis) along the PC2 were determined (“C-nevadensis—nevadensis” distance; “C-nevadensis—iberica” distance; “C-iberica—nevadensis” distance; “C-iberica—iberica” distance). We used the linear segment between the PC2 coordinate of each caterpillar and pine sample as the distance unit representing the metabolomics distances between caterpillars and pines for such axis. In total, 1152 distances along the PC2 were determined (12 caterpillar samples × 2 caterpillar groups × 24 pine samples (AT.ABs and AT.NABs) × 2 pine subspecies). All distances were subsequently submitted to one-way ANOVA, followed with a HSD Tukey’s post-hoc test to test for differences between them. All statistical analyses were performed in R [84]. One-way ANOVAs were performed with the function aov from the package “stats” [84]. HSD Tukey’s tests were performed with the HSD.test function of the package “agricolae” [85]. Linear correlations were performed using rcorr function from the package Hmisc [86]. The PERMANOVA analysis was conducted with adonis function in the package “vegan” [87]. PCAs were plotted using the PCA function from “FactoMineR” package [88] with the missing data from the dataset imputed using imputePCA function from package “missMDA” [89]. All graphs were first performed in R [84] and subsequently graphically treated by Adobe Illustrator CS6 (San José, CA, USA). Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/1/59/s1. Author Contributions: Conceptualization of the research, A.R.-U., J.A.H., J.P., and J.S.; Field methodology, A.R.-U., J.A.H., and J.S.; Laboratory methodology, A.R-U., M.O., O.U.; and LC-MS file Acquisition, M.O. and O.U; Dataset generation, A.R.-U.; Statistical Analyses and visualization of the data, A.R.-U. and J.S.; Writing—Original Draft Preparation, A.R.-U., J.P., J.S, and J.A.H.; Writing—Review & Editing, A.R.-U., J.P., J.S, J.A.H., M.O., L.P.-T., and O.U; Supervision, J.P., J.S., J.A.H., O.U., and L.P.-T.; Project Administration, J.A.H., J.P., and O.U.; Funding Acquisition, J.P., J.A.H., L.P.-T., and O.U. Funding: This research was funded by the research fellowship (JAE) from the CSIC (A.R.-U), the European Research Council Synergy grant SyG-2013-610028 IMBALANCE-P (J.P., J.S.), the Spanish Government projects CGL2016-48074-P and OAPN 022/2008 (PROPINOL) (J.P., J.S.), the Catalan Government project SGR 2014-274 (J.P., J.S.), DOE Office of Science User Facility sponsored by the Office of Biological and Environmental Research at the Pacific Northwest National Laboratory (A.R.-U), and by the Ministry of Education, Youth and Sports of the Czech Republic (grant No. CZ.02.1.01/0.0/0.0/16_013/0001609, and No. LO1415) (O.U., M.O). Acknowledgments: The authors thank Laia Mateu-Castell, Sara Férez, Ander Achotegui-Castells and Joan Llusià for their field and laboratory support. Conflicts of Interest: The authors declare no conflict of interest. References 1. Behmer, S.T. Insect Herbivore Nutrient Regulation. Annu. Rev. Entomol. 2009, 54, 165–187. [CrossRef] [PubMed] 2. Bernays, E.A. Evolution of Feeding Behavior in Insect Herbivores. Bioscience 1998, 48, 35–44. [CrossRef] 3. Pauw, A.; Stofberg, J.; Waterman, R.J. Flies and Flowers in Darwin’s Race. Evolution 2009, 63, 268–279. [CrossRef] [PubMed] Int. J. Mol. Sci. 2019, 20, 59 16 of 19 4. Agrawal, A.A.; Petschenka, G.; Bingham, R.A.; Weber, M.G.; Rasmann, S. Toxic cardenolides: Chemical ecology and coevolution of specialized plant-herbivore interactions. New Phytol. 2012, 194, 28–45. [CrossRef] [PubMed] 5. Toju, H. Weevils and camellias in a Darwin’s race: Model system for the study of eco-evolutionary interactions between species. Ecol. Res. 2011, 26, 239–251. [CrossRef] 6. Heimpel, G.E.; Rosenheim, J.A.; Kattari, D. Adult feeding and lifetime reproductive success in the parasitoid Aphytis melinus. Entomol. Exp. Appl. 1997, 83, 305–315. [CrossRef] 7. Andersen, L.H.; Kristensen, T.N.; Loeschcke, V.; Toft, S.; Mayntz, D. Protein and carbohydrate composition of larval food affects tolerance to thermal stress and desiccation in adult Drosophila melanogaster. J. Insect Physiol. 2010, 56, 336–340. [CrossRef] 8. Bauerfeind, S.S.; Fischer, K. Testing the plant stress hypothesis: Stressed plants offer better food to an insect herbivore. Entomol. Exp. Appl. 2013, 149, 148–158. [CrossRef] 9. Kasumovic, M.M. The multidimensional consequences of the juvenile environment: Towards an integrative view of the adult phenotype. Anim. Behav. 2013, 85, 1049–1059. [CrossRef] 10. Waldbauer, G.P.; Cohen, R.W.; Friedman, S. Self-Selection of an Optimal Nutrient Mix from Defined Diets by Larvae of the Corn Earworm, Heliothis zea (Boddie). Physiol. Zool. 1984, 57, 590–597. [CrossRef] 11. Simpson, S.J.; Sibly, R.M.; Lee, K.P.; Behmer, S.T.; Raubenheimer, D. Optimal foraging when regulating intake of multiple nutrients. Anim. Behav. 2004, 68, 1299–1311. [CrossRef] 12. Le Gall, M.; Behmer, S.T. Effects of Protein and Carbohydrate on an Insect Herbivore: The Vista from a Fitness Landscape. Integr. Comp. Biol. 2014, 54, 942–954. [CrossRef] [PubMed] 13. Mattson, W.J. Herbivory in Relation to Plant Nitrogen Content. Annu. Rev. Ecol. Syst. 1980, 11, 119–161. [CrossRef] 14. Sterner, R.; Elser, J. Ecological Stoichiometry: The Biology of Elements from Molecules to the Biosphere; Princetion University Press: Princeton, NJ, USA, 2002. 15. Ngai, J.T.; Jefferies, R.L. Nutrient limitation of plant growth and forage quality in Arctic coastal marshes. J. Ecol. 2004, 92, 1001–1010. [CrossRef] 16. Sardans, J.; Rivas-Ubach, A.; Peñuelas, J. The elemental stoichiometry of aquatic and terrestrial ecosystems and its relationships with organismic lifestyle and ecosystem structure and function: A review and perspectives. Biogeochemistry 2012, 111, 1–39. [CrossRef] 17. Chen, Y.; Coleman, T.W.; Jones, M.I.; Flint, M.L.; Seybold, S.J. Foliar nutrients explain goldspotted oak borer, Agrilus auroguttatus, adult feeding preference among four California oak species. Entomol. Exp. Appl. 2013, 149, 57–66. [CrossRef] 18. Hódar, J.A.; Zamora, R.; Castro, J. Host utilisation by moth and larval survival of pine processionary caterpillar Thaumetopoea pityocampa in relation to food quality in three Pinus species. Ecol. Entomol. 2002, 27, 292–301. [CrossRef] 19. Metspalu, L.; Kruus, E.; Jögar, K.; Kuusik, A.; Williams, I.H.; Veromann, E.; Luik, A.; Ploomi, A.; Hiiesaar, K.; Kivimägi, I.; Mänd, M. Larval food plants can regulate the cabbage moth, Mamestra brassicae population. Bull. Insectol. 2013, 66, 93–101. 20. Hódar, J.A.; Castro, J.; Zamora, R. Pine processionary caterpillar Thaumetopoea pityocampa as a new threat for relict Mediterranean Scots pine forests under climatic warming. Biol. Conserv. 2003, 110, 123–129. [CrossRef] 21. Loaiza, V.; Jonas, J.L.; Joern, A. Grasshoppers (Orthoptera: Acrididae) select vegetation patches in local-scale responses to foliar nitrogen but not phosphorus in native grassland. Insect Sci. 2011, 18, 533–540. [CrossRef] 22. Rivas-Ubach, A.; Gargallo-Garriga, A.; Sardans, J.; Oravec, M.; Mateu-Castell, L.; Pérez-Trujillo, M.; Parella, T.; Ogaya, R.; Urban, O.; Peñuelas, J. Drought enhances folivory by shifting foliar metabolomes in Quercus ilex trees. New Phytol. 2014, 202, 874–885. [CrossRef] [PubMed] 23. Rivas-Ubach, A.; Sardans, J.; Hódar, J.A.; Garcia-Porta, J.; Guenther, A.; Oravec, M.; Urban, O.; Peñuelas, J. Similar local but different systemic metabolomic responses of closely related pine subspecies to folivory by caterpillars of the processionary moth. Plant Biol. 2016, 18, 484–494. [CrossRef] [PubMed] 24. Colasurdo, N.; Gélinas, Y.; Despland, E. Larval nutrition affects life history traits in a capital breeding moth. J. Exp. Biol. 2009, 212, 1794–1800. [CrossRef] [PubMed] 25. Moreau, J.; Benrey, B.; Thiery, D. Assessing larval food quality for phytophagous insects: Are the facts as simple as they appear? Funct. Ecol. 2006, 20, 92–600. [CrossRef] Int. J. Mol. Sci. 2019, 20, 59 17 of 19 26. Fagan, W.F.; Siemann, E.; Mitter, C.; Denno, R.F.; Huberty, A.F.; Woods, H.A.; Elser, J.J. Nitrogen in insects: Implications for trophic complexity and species diversification. Am. Nat. 2002, 160, 784–802. [CrossRef] 27. Denno, R.F.; Fagan, W.F. Might nitrogen limitation promote omnivory among carnivorous arthropods? Ecology 2003, 84, 2522–2531. [CrossRef] 28. Peñuelas, J.; Sardans, J. Ecological metabolomics. Chem. Ecol. 2009, 25, 305–309. [CrossRef] 29. Rivas-Ubach, A.; Sardans, J.; Peŕez-Trujillo, M.; Estiarte, M.; Penũelas, J. Strong relationship between elemental stoichiometry and metabolome in plants. Proc. Natl. Acad. Sci. USA 2012, 109, 4181–4186. [CrossRef] 30. Brunetti, C.; George, R.M.; Tattini, M.; Field, K.; Davey, M.P. Metabolomics in plant environmental physiology. J. Exp. Bot. 2013, 64, 4011–4020. [CrossRef] 31. Sardans, J.; Peñuelas, J.; Rivas-Ubach, A. Ecological metabolomics: Overview of current developments and future challenges. Chemoecology 2011, 21, 191–225. [CrossRef] 32. Bundy, J.G.; Davey, M.P.; Viant, M.R. Environmental metabolomics: A critical review and future perspectives. Metabolomics 2008, 5, 3–21. [CrossRef] 33. Macel, M.; Van Dam, N.M.; Keurentjes, J.J.B. Metabolomics: The chemistry between ecology and genetics. Mol. Ecol. Resour. 2010, 10, 583–593. [CrossRef] [PubMed] 34. Shulaev, V.; Cortes, D.; Miller, G.; Mittler, R. Metabolomics for plant stress response. Physiol. Plant. 2008, 132, 199–208. [CrossRef] [PubMed] 35. Rivas-Ubach, A.; Barbeta, A.; Sardans, J.; Guenther, A.; Ogaya, R.; Oravec, M.; Urban, O.; Peñuelas, J. Topsoil depth substantially influences the responses to drought of the foliar metabolomes of Mediterranean forests. Perspect. Plant Ecol. Evol. Syst. 2016, 21, 41–54. [CrossRef] 36. Rivas-Ubach, A.; Sardans, J.; Hódar, J.A.; Garcia-Porta, J.; Guenther, A.; Paša-Tolić, L.; Oravec, M.; Urban, O.; Peñuelas, J. Close and distant: Contrasting the metabolism of two closely related subspecies of Scots pine under the effects of folivory and summer drought. Ecol. Evol. 2017, 7, 8976–8988. [CrossRef] 37. Gargallo-Garriga, A.; Sardans, J.; Pérez-Trujillo, M.; Oravec, M.; Urban, O.; Jentsch, A.; Kreyling, J.; Beierkuhnlein, C.; Parella, T.; Peñuelas, J. Warming differentially influences the effects of drought on stoichiometry and metabolomics in shoots and roots. New Phytol. 2015, 207, 591–603. [CrossRef] 38. Gargallo-Garriga, A.; Sardans, J.; Pérez-Trujillo, M.; Rivas-Ubach, A.; Oravec, M.; Vecerova, K.; Urban, O.; Jentsch, A.; Kreyling, J.; Beierkuhnlein, C.; et al. Opposite metabolic responses of shoots and roots to drought. Sci. Rep. 2014, 4, 6829. [CrossRef] [PubMed] 39. Fiehn, O. Metabolomics—The link between genotypes and phenotypes. Plant Mol. Biol. 2002, 48, 155–171. [CrossRef] 40. Patti, G.J.; Yanes, O.; Siuzdak, G. Metabolomics: The apogee of the omics trilogy. Nat. Rev. Mol. Cell Biol. 2012, 13, 263–269. [CrossRef] 41. Johnson, C.H.; Ivanisevic, J.; Siuzdak, G. Metabolomics: Beyond biomarkers and towards mechanisms. Nat. Rev. Mol. Cell Biol. 2016, 17, 451–459. [CrossRef] 42. Mirnezhad, M.; Romero-González, R.R.; Leiss, K.A.; Choi, Y.H.; Verpoorte, R.; Klinkhamer, P.G.L. Metabolomic analysis of host plant resistance to thrips in wild and cultivated tomatoes. Phytochem. Anal. 2010, 21, 110–117. [CrossRef] [PubMed] 43. Jansen, J.J.; Allwood, J.W.; Marsden-Edwards, E.; van der Putten, W.H.; Goodacre, R.; van Dam, N.M. Metabolomic analysis of the interaction between plants and herbivores. Metabolomics 2009, 5, 150–161. [CrossRef] 44. Gómez, S.; Steinbrenner, A.D.; Osorio, S.; Schueller, M.; Ferrieri, R.A.; Fernie, A.R.; Orians, C.M. From shoots to roots: Transport and metabolic changes in tomato after simulated feeding by a specialist lepidopteran. Entomol. Exp. Appl. 2012, 144, 101–111. [CrossRef] 45. Glauser, G.; Marti, G.; Villard, N.; Doyen, G.A.; Wolfender, J.-L.; Turlings, T.C.J.; Erb, M. Induction and detoxification of maize 1,4-benzoxazin-3-ones by insect herbivores. Plant J. 2011, 68, 901–911. [CrossRef] [PubMed] 46. Rivas-Ubach, A.; Liu, Y.; Bianchi, T.S.; Toli, N.; Jansson, C.; Paša-Tolić, L. Moving beyond the van Krevelen diagram: A new stoichiometric approach for compound classification in organisms. Anal. Chem. 2018, 90, 6152–6160. [CrossRef] [PubMed] 47. Avtzis, N. Development of Thaumetopoea pityocampa schiff.(Lepidoptera: Thaumetopoeidae) in relation to food consumption. For. Ecol. Manag. 1986, 15, 65–68. [CrossRef] Int. J. Mol. Sci. 2019, 20, 59 18 of 19 48. Battisti, A. Host-plant relationships and population dynamics of the Pine Processionary Caterpillar Thaumetopoea pityocampa (Denis & Schiffermuller). J. Appl. Entomol. 1988, 105, 393–402. [CrossRef] 49. Battisti, A.; Stastny, M.; Netherer, S.; Robinet, C.; Schopf, A.; Roques, A.; Larsson, S. Expansion of geographic range in the pine processionary moth caused by increased winter temperatures. Ecol. Appl. 2005, 15, 2084–2096. [CrossRef] 50. Robledo-Arnuncio, J.J.; Navascués, M.; González-Martínez, S.C.; Gil, L. Estimating gametic introgression rates in a risk assessment context: A case study with Scots pine relicts. Heredity 2009, 103, 385–393. [CrossRef] 51. Johnson, S.N.; McNicol, J.W. Elevated CO2 and aboveground-belowground herbivory by the clover root weevil. Oecologia 2010, 162, 209–216. [CrossRef] 52. Elser, J.J.; Bracken, M.E.S.; Cleland, E.E.; Gruner, D.S.; Harpole, W.S.; Hillebrand, H.; Ngai, J.T.; Seabloom, E.W.; Shurin, J.B.; Smith, J.E. Global analysis of nitrogen and phosphorus limitation of primary producers in freshwater, marine and terrestrial ecosystems. Ecol. Lett. 2007, 10, 1135–1142. [CrossRef] [PubMed] 53. Sardans, J.; Rivas-Ubach, A.; Peñuelas, J. The C:N:P stoichiometry of organisms and ecosystems in a changing world: A review and perspectives. Perspect. Plant Ecol. Evol. Syst. 2012, 14, 33–47. [CrossRef] 54. Andersen, T.; Hessen, D.O. Carbon, nitrogen, and phosphorus content of freshwater zooplankton. Limnol. Oceanogr. 1991, 36, 807–814. [CrossRef] 55. Mayntz, D.; Toft, S. Nutrient composition of the prey’s diet affects growth and survivorship of a generalist predator. Oecologia 2001, 127, 207–213. [CrossRef] [PubMed] 56. Ashhurst, D.E. The connective tissues of insects. Annu. Rev. Entomol. 1968, 13, 45–74. [CrossRef] 57. Sardans, J.; Peñuelas, J. Potassium: A neglected nutrient in global change. Glob. Ecol. Biogeogr. 2015, 24, 261–275. [CrossRef] 58. Clarkson, D.T.; Hanson, J.B. The Mineral Nutrition of Higher Plants. Annu. Rev. Plant Physiol. 1980, 31, 239–298. [CrossRef] 59. Sardans, J.; Peñuelas, J.; Coll, M.; Vayreda, J.; Rivas-Ubach, A. Stoichiometry of potassium is largely determined by water availability and growth in Catalonian forests. Funct. Ecol. 2012, 26, 1077–1089. [CrossRef] 60. Acquisti, C.; Elser, J.J.; Kumar, S. Ecological Nitrogen Limitation Shapes the DNA Composition of Plant Genomes. Mol. Biol. Evol. 2009, 26, 953–956. [CrossRef] 61. Mathavan, S.; Muthukrishnan, J. Use of larval weight and length as indices of gut content weight in some Lepidopterous larvae. Oecologia 1980, 44, 317–318. [CrossRef] 62. Hoch, G.; Toffolo, E.P.; Netherer, S.; Battisti, A.; Schopf, A. Survival at low temperature of larvae of the pine processionary moth Thaumetopoea pityocampa from an area of range expansion. Agric. For. Entomol. 2009, 11, 313–320. [CrossRef] 63. Robinet, C.; Baier, P.; Pennerstorfer, J.; Schopf, A.; Roques, A. Modelling the effects of climate change on the potential feeding activity of Thaumetopoea pityocampa (Den. & Schiff.) (Lep., Notodontidae) in France. Glob. Ecol. Biol. 2007, 16, 460–471. 64. Després, L.; David, J.-P.; Gallet, C. The evolutionary ecology of insect resistance to plant chemicals. Trends Ecol. Evol. 2007, 22, 298–307. [CrossRef] [PubMed] 65. Nishida, R. Sequestration of Defensive Substances from Plants by Lepidoptera. Annu. Rev. Entomol. 2002, 47, 57–92. [CrossRef] [PubMed] 66. Paoletti, M.G. Ecological Implications of Minilivestock; CRC Press: Boca Raton, FL, USA, 2005; ISBN 9781482294439. 67. Banjo, A.; Lawal, O.; Songonuga, E. The nutritional value of fourteen species of edible insects in southwestern Nigeria. Afr. J. Biotechnol. 2006, 5, 298–301. 68. Sulanc, M.; Emre, I. Effects of B group vitamins and choline chloride on the development and protein synthesis in the male larvae of Pimpla turionellae L. (Hym., Ichneumonidae). J. Appl. Entomol. 2000, 124, 151–153. [CrossRef] 69. Liu, F.; Wei, F.; Wang, L.; Liu, H.; Zhu, X.; Liang, Y. Riboflavin activates defense responses in tobacco and induces resistance against Phytophthora parasitica and Ralstonia solanacearum. Physiol. Mol. Plant Pathol. 2010, 74, 330–336. [CrossRef] 70. Taheri, P.; Tarighi, S. Riboflavin induces resistance in rice against Rhizoctonia solani via jasmonate-mediated priming of phenylpropanoid pathway. J. Plant Physiol. 2010, 167, 201–208. [CrossRef] [PubMed] Int. J. Mol. Sci. 2019, 20, 59 19 of 19 71. Verdrengh, M.; Tarkowski, A. Riboflavin in innate and acquired immune responses. Inflamm. Res. 2005, 54, 390–393. [CrossRef] 72. Peñuelas, J.; Sardans, J.; Stefanescu, C.; Parella, T.; Filella, I. Lonicera Implexa Leaves Bearing Naturally Laid Eggs of the Specialist Herbivore Euphydryas Aurinia have Dramatically Greater Concentrations of Iridoid Glycosides than other Leaves. J. Chem. Ecol. 2006, 32, 1925–1933. [CrossRef] 73. Gols, R.; Harvey, J.A. Plant-mediated effects in the Brassicaceae on the performance and behaviour of parasitoids. Phytochem. Rev. 2009, 8, 187–206. [CrossRef] 74. Achotegui-Castells, A.; Llusia, J.; Hódar, J.A.; Peñuelas, J. Needle terpene concentrations and emissions of two coexisting subspecies of Scots pine attacked by the pine processionary moth (Thaumetopoea pityocampa). Acta Physiol. Plant. 2013, 35, 3047–3058. [CrossRef] 75. Rajsnerová, P.; Klem, K.; Holub, P.; Novotná, K.; Večeřová, K.; Kozáčiková, M.; Rivas-Ubach, A.; Sardans, J.; Marek, M.V.; Peñuelas, J.; et al. Morphological, biochemical and physiological traits of upper and lower canopy leaves of European beech tend to converge with increasing altitude. Tree Phys. 2015, 1, 47–60. [CrossRef] [PubMed] 76. Sardans, J.; Montes, F.; Peñuelas, J. Determination of As, Cd, Cu, Hg and Pb in biological samples by modern electrothermal atomic absorption spectrometry. Spectrochim. Acta Part B At. Spectrosc. 2010, 65, 97–112. [CrossRef] 77. t’Kindt, R.; De Veylder, L.; Storme, M.; Deforce, D.; Van Bocxlaer, J. LC-MS metabolic profiling of Arabidopsis thaliana plant leaves and cell cultures: Optimization of pre-LC-MS procedure parameters. J. Chromatogr. B. Analyt. Technol. Biomed. Life Sci. 2008, 871, 37–43. [CrossRef] [PubMed] 78. Pluskal, T.; Castillo, S.; Villar-Briones, A.; Orešič, M. MZmine 2: Modular framework for processing, visualizing, and analyzing mass spectrometry-based molecular profile data. BMC Bioinform. 2010, 11, 395. [CrossRef] [PubMed] 79. Sumner, L.W.; Amberg, A.; Barrett, D.; Beale, M.H.; Beger, R.; Daykin, C.A.; Fan, T.W.-M.; Fiehn, O.; Goodacre, R.; Griffin, J.L.; et al. Proposed minimum reporting standards for chemical analysis. Metabolomics 2007, 3, 211–221. [CrossRef] 80. Rivas-Ubach, A.; Hódar, J.A.; Sardans, J.; Kyle, J.E.; Kim, Y.-M.; Oravec, M.; Urban, O.; Guenther, A.; Peñuelas, J. Are the metabolomic responses to folivory of closely related plant species linked to macroevolutionary and plant-folivore coevolutionary processes? Ecol. Evol. 2016, 6, 4372–4386. [CrossRef] 81. Kim, H.K.; Choi, Y.H.; Verpoorte, R. NMR-based metabolomic analysis of plants. Nat. Protoc. 2010, 5, 536–549. [CrossRef] 82. Van den Berg, R.A.; Hoefsloot, H.C.; Westerhuis, J.A.; Smilde, A.K.; van der Werf, M.J. Centering, scaling, and transformations: Improving the biological information content of metabolomics data. BMC Genom. 2006, 7, 142. [CrossRef] 83. Rivas-Ubach, A.; Pérez-Trujillo, M.; Sardans, J.; Gargallo-Garriga, A.; Parella, T.; Peñuelas, J. Ecometabolomics: Optimized NMR-based method. Methods Ecol. Evol. 2013, 4, 464–473. [CrossRef] 84. R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria, 2013. Available online: http://www.R-project.org (accessed on 22 December 2018). 85. De Mendiburu, F. Agricolae: Statistical Procedures for Agricultural Research. R package Version 1.2-3. 2015. Available online: https://CRAN.R-project.org/package=agricolae (accessed on 22 December 2018). 86. Harrell, F.E. Hmisc: Harrell Miscellaneous. R Package Version 4.1-1. 2018. Available online: https://CRAN. R-project.org/package=Hmisc (accessed on 22 December 2018). 87. Oksanen, J.; Guillaume-Blanchet, F.; Fierndly, M.; Kindt, R.; Legendre, P.; McGlinnd, D.; Minchin, P.; O’Hara, R.; Simpson, G.; Solymos, P.; et al. Vegan: Community Ecology Package. R Package Version 2.1-1. 2018. Available online: https://CRAN.R-project.org/package=vegan (accessed on 22 December 2018). 88. Husson, F.; Josse, J.; Le, S.; Mazet, J. FactoMineR: Multivariate Exploratory Data Analysis and Data Mining. R Package Version 1.32. 2016. Available online: https://CRAN.R-project.org/package=FactoMineR (accessed on 22 December 2018). 89. Husson, F.; Josse, J. missMDA: Handling Missing Values with Multivariate Data Analysis. R Package Version 1.9. 2015. Available online: https://CRAN.R-project.org/package=missMDA (accessed on 22 December 2018). © 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).