Sheep feeding preference as a tool to control pine invasion in Patagonia: influence of foliar toughness, terpenoids and resin content

Herbivores modulate the structure and composition of plant communities, including plant invasions. This is conditioned by plant palatability which can be reduced by its chemical or physical traits. The effects that ungulates browsing has on pine invasions are variable and the empirical evidence on the causes of this variability is scarce. We experimentally explored how sheep browsing preference varies between seedlings of pine species with different invasiveness; Pinus contorta (high invasiveness), P. ponderosa (medium invasiveness), P. radiata (medium invasiveness) and P. jeffreyi (low invasiveness). Secondly, we quantified anti-herbivory chemical compounds and physical traits of these species and related them with sheep preference observed. The browsing incidence of P. contorta was 68%, P. ponderosa 58%, P. radiata 29%, and P. jeffreyi 84%. Among anti-herbivory traits analyzed, α-pinene concentration had a negative effect on the probability of a terminal bud being browsed and on browsing intensity. Meanwhile, foliar toughness was negatively related to browsing intensity and water concentration was positively related to browsing intensity. Also, the most invasive species, P. contorta, was highly damaged. Thus, sheep herbivory could be slowing pine invasion rate; suggesting that could be considered a tool to control early invasions, especially for this particular species.

Sheep browsing preference. To evaluate the differential preference of sheep (Ovis aries) among nonnative pines with different invasive capacity, we exposed seedlings of four species to Merino wethers herbivory (castrated male sheep) under a recommended stocking rate on that site. Pine species were P. contorta (high invasiveness), P. ponderosa and P. radiata (moderate invasiveness) and P. jeffreyi (low invasiveness). Within the study area, we installed five experimental enclosures of 25 × 25 m with a 1.50 m tall fence (each enclosure constituted a replicate). Enclosures had hardware cloth 0.6 m tall buried in the soil to prevent access by other herbivores, such as European hare (L. europaeus). Each enclosure was stocked with sheep at a density typically recommended for the study area (0.2-0.3 wethers/ha/year for grassland) 61,62 . This density was based on the assessment of the forage condition of the sites and the area of the enclosure, according to this, the density estimated corresponded to two wethers in each enclosure for three days. More details about the estimation of the stocking rate can be found as Supplementary Methods S1 online. The experiment was conducted in October 2015, during the austral spring, the season with the maximum availability and quality of forage in the community, minimizing the risk of overgrazing. Also, we established the enclosures adjacent to each other, and no sheep were alone, so as not to alter their gregarious behavior. A similar design was used successfully in previous studies in the region 49,60 . Environmental conditions and forage availability between enclosures corresponding to a replicate were similar. In each enclosure (a replicate) we planted 20 seedlings of each pine species (100 individuals of each species in total), alternating species so that there were not two consecutive individuals of the same species.
The number of browsed seedlings, maximum seedling height before and after the treatment, number of terminal buds damaged, number of browsed branches, the number of defoliated seedlings (i.e. whether the seedling was only the needles damaged, not branches) and the probability of seedling survival immediately after the treatment were recorded for all individuals (100 individuals of each of the four species). The probability of seedling survival was estimated based on the conservative assumption that only seedlings with all branches browsed would die. All seedlings were removed at the end of the experiment.
Because our research was conducted on non-regulated animals in our country, we did not require any ethical approval. However is important to mention that all animals were treated under general animal welfare criteria following the guidelines of National Institutes of Animal production 63,64 , with constant clean water and food availability (besides the pine seedlings presented during the experiment).

Anti-herbivores chemical compounds and physical traits. To examine whether seedlings of the four
Pinus species studied (Pinus contorta, P. ponderosa, P. radiata, and P. jeffreyi) differ in their chemical defenses, physical characteristics and plant quality; we measure toughness and monoterpene content in needles of 20 seedlings per species, and 15 out of those 20 seedlings were used to measured non-volatile resin and water content of the needles and branches. It was not possible to perform the sheep preference experiment and these measurements on the same individuals for two reasons. Firstly, since we were studying seedlings (not juveniles or adults), the amount of biomass left after sheep browsed was not enough to measure the chemical defenses and physical traits in the same individual. Secondly, it is known that pine trees have induced defenses 65,66 , so if we measured the chemistry of the seedling after the sheep fed on them, the chemistry would have changed in response to the damage, and we wanted to relate the browsing preference of sheep to the chemistry and physical traits for the four pine species to tissues that were not previously damaged.
To measure needle toughness, 60 fresh needles from 20 seedlings per Pinus species were removed. Needle toughness (g/cm 2 ) was measured with a force gauge penetrometer (type 516; Chatillon, Largo, FL, USA) and three measurements were averaged per seedling. The remaining plant material was placed individually in plastic bags and stored in a − 80 ºC freezer for later analyses. For monoterpene extraction, needles were ground with liquid nitrogen in a chilled mortar and pestle to minimize the possible loss of monoterpene. For each sample, the frozen needle powder was put in a separate borosilicate glass vial, and weighted exact weights (between 0.5 and 0.6 g). After recording the weights, 2 mL of GCgrade n-hexane (Merck, Germany) with 0.1 μL mL-1 ( +)-fenchone (Sigma-Aldrich) as an internal standard was added to each vial. Vials were immediately closed with bakelite caps and parafilm to avoid evaporation and mixed with a vortex. The ground needles were left to soak for 7 days at ambient temperature. After a week, 2 mL aliquots of each sample were placed into a clear chromatography autosampler vials for GC analysis. A 1 μL subsample of each aliquot was then injected on a PerkinElmer Clarus 680 Gas Chromatograph with PerkinElmer Clarus 600 T Mass Spectrometer (GC-MS) fitted with an Elite-5 MS column (60 m × 0.25 mm × 0.25 μm; Perkin Elmer) and coupled with a mass spectrometer (PerkinElmer Clarus 600 T Mass Spectrometer) used as detector. The carrier gas used was Helium at a flow rate of 1 mL/min for 20 mL/min split mode. The injector temperature was set at 250 °C to volatilize the liquid sample and the time-temperature history of the oven enclosing the column consisted of three steps: an initial heating at 60 °C during 1 min, a temperature ramp of at 5 °C/min from 60 to 130 °C, and a final holding at 130 °C during 15 min. We could identify four monoterpenes (α-pinene, β-pinene, 3-carene and β-phellandrene) by comparing retention times of known standards and mass spectra using the NIST Mass Spectral Search Program for the NIST/EPA/ NIH Mass Spectral Library (version 2.0f) and TurboMass V 5.4.2 software to calculate the area under the absorbance curve for each peak. Since we were unable to obtain a pure standard of the monoterpene β-phellandrene, this monoterpene was measured using the α-phellandrene pure standard being α and β-phellandrenes similar in retention times and mass spectra. Concentrations for each compound were calculated using 9-point calibration curves with injections of known amounts of pure standards and fenchone (the internal standard). All standards were purchased from Sigma Aldrich (Saint Louis, MO). Scientific RepoRtS | (2020) 10:12113 | https://doi.org/10.1038/s41598-020-68748-y www.nature.com/scientificreports/ In order to measure the non-volatile resin content, needles and branches were analyzed together. We cut 15 cm long branches with its needles in pieces, weighted to the nearest 0.0001 g and placed in labeled borosilicate glass tubes. Then each tube was filled with hexane until the material was covered, mixed with a vortex and placed in an ultrasonic bath at 45 °C for 20 min. The tubes were closed and the plant material soaked for 24 h at room temperature. The entire extraction was repeated twice to maximize extraction outputs. The solution was then filtered with paper filters in pre-weighted borosilicate glass tubes and left to evaporate under the fume hood. Each tube has a residual content that corresponded to the non-volatile resin and was weighted to the nearest 0.0001 g. Since non-volatile resin content was obtained from wet weights, a separate subset of 15 cm branches from each sample was dried (oven dry 60 °C until constant weight) and weighed to obtain dry mass conversion factors. The dry weight (DW) of the sample was calculated using the wet weights of the samples and multiplied by each conversion factor. The resin content of the sample was expressed as mg of non-volatile resin /g of stem DW. To assess plant quality, percent water content was calculated as [(wet weight − dry weight)/wet weight] × 100, from the same 15 cm long branches used to measure total non-volatile resin content. Data analysis. Sheep browsing preference. We used generalized linear mixed models (GLMM) 67,68 to test the browsing preference of sheep for four species of pines. The models were fitted with the enclosure as a random effect to control for the possible variability between replicates in adjacent enclosures. The predictor variable was pine species, a factor with four levels. We modeled the following six response variables: (a) the probability of a seedling being browsed (browsing incidence: number of browsed seedlings/total number of seedlings), for this variable, we consider that a seedling was browsed if the terminal bud, and/or any branches, and/or the needles were damaged, (b) the relative reduction in height (initial height − final height/initial height), (c) the probability of a terminal bud being browsed, (d) the proportion of browsed branches per seedling (browsing intensity: number of browsed branches/total number of branches), (e) the probability of a seedling being defoliated (whether only the needles were damaged and it was not damaged the stem or branches), and (f) the probability of seedling survival. Models a, c, d, e, and f were fitted with a logit link function and binomial error distribution and model b was fitted with a logit link function and beta error distribution.
Because in the preference trial the four species were presented to sheep at the same time, we performed posthoc comparisons that compared one species with the other three in turn and not a paired comparison. For that analysis, the lsmeans function from "lsmeans" R package 69 was applied and the "Scheffe" adjust method was used.
The function glmer from the R package "lme4" 70 was used to fit the browsing incidence, the probability of a terminal bud being browsed, the browsing intensity, the probability of a seedling being defoliated, and the probability of seedling survival models. The function glmmadmb from the "glmmADMB" package 71 was used to fit the relative reduction in height (initial height − final height/initial height) model.

Anti-herbivores chemical compounds and physical traits.
To identify differences in total monoterpene between species (that follow a normal distribution) we performed an ANOVA and pairwise post-hoc comparisons of the coefficients were tested using the function glht from the "multcomp" package 72 . The variables α-pinene, β-pinene, 3-carene, β-phellandrene, resins, toughness, and water content did not follow a normal distribution (Shapiro-Wilks Normality Test) so, to test differences in these variables between species, the non-parametric Kruskal-Wallis test was applied, and pairwise post-hoc comparisons of the coefficients were tested using the function kruskalmc from the "pgirmess" R package 73 .

Relation between sheep preference and pines' anti-herbivores chemical compounds and physical traits.
With the aim to have an overview of the relation between sheep herbivory and the anti-herbivores compounds assessed we explore two aspects of the experiments detailed before. First, to understand how anti-herbivore traits drive the first approach of sheep and pines (i.e. as browsing repellents), we assessed the relation between the concentrations of chemical compounds and physical traits with the probability of a sheep browse only the terminal bud of the seedlings. On the other hand, to understand how anti-herbivore traits drive sheep browsing intensity on pines (i.e. sheep continues to eat the pine after the first approach), we assessed the relationship between these anti-herbivore defenses and browsing intensity. To identify the relation between the chemical and physical variables (predictor variables) with the probability of a terminal bud being browsed or the browsing intensity (response variables) we performed three generalized linear mixed models (GLMM) for each response variable: (1) total monoterpenes and total resins as predictor variables; (2) monoterpenes separately, i.e. α-pinene, β-pinene, 3-carene and β-phellandrene as predictor variables and (3) toughness and water content as predictor variables. The models were fitted with a logit link function and binomial error distribution. Before performing these regressions we checked for collinearity and, following Dormann et al. 74 , we only include variables with a correlation <|0.7|. Only β-pinene and 3-carene were highly correlated (0.77), so for model 2, we do not include β-pinene. Because the number of individuals tested for sheep preferences (20 seedlings of each species in five enclosures, 100 seedlings of each species in total) was greater than individuals tested for chemical and physical content (15 seedlings of each species to test resins and water, and 20 seedlings of each species to test monoterpenes and foliar toughness), we performed the regressions considering the median of each predictor variable in relation to the individuals tested in the preference experiments. We are aware that this analysis has limitations, but it is appropriate to give a general idea of the principal relations between the anti-herbivory compounds and the sheep preference observed. Also, this approach is used in other biochemical studies due to the high cost and difficulties to assess these traits in larger sampler sizes 75 . Predictor variables were standardized. The models were nested, fitting species as random effects to consider the possible variability between them. In the cases where there was no variation between species, we fitted a generalized linear model (GLM) without random effects. The function glmer from the R package "lme4" 70

Results
Sheep browsing preference. In general terms, the browsing incidence on P. contorta, P. ponderosa, P. radiata, and P. jeffreyi seedlings, as estimated by the model, was 68%, 58%, 29%, and 84% respectively (Table 1). Pinus radiata seedlings were significantly less browsed than the seedlings of the other three pine species, while P. jeffreyi seedlings were significantly more browsed than the seedlings of the other species (Fig. 1a, results of post-hoc comparisons can be found as Supplementary Table S2.1 online). The relative reduction in P. contorta, P. ponderosa, P. radiata, and P. jeffreyi seedling height estimated by the model was 27%, 24%, 23%, and 18% respectively (Table 1) with the reduction in height of P. contorta seedlings significantly higher than the reductions in height of P. radiata and P. jeffreyi seedlings; and reductions in height of P. ponderosa seedlings were significantly greater than the reduction in height of P. jeffreyi seedlings ( Fig. 1b; results of post-hoc comparisons can be found as Supplementary Table S2.1 online).
The probability that a terminal bud of a P. contorta, a P. ponderosa, a P. radiata and a P. jeffreyi seedling being damaged, as estimated by the model, was 44%, 29%, 31%, and 19% respectively (Table 1); with the proportion of terminal bud browsed of P. contorta seedlings significantly higher than the proportion of terminal bud browsed of P. jeffreyi (Fig. 1c, results of post-hoc comparisons can be found as Supplementary Table S2.1 online).
The browsing intensity of P. contorta, P. ponderosa, P. radiata and P. jeffreyi seedlings estimated by the model was 28%, 29%, 27%, and 18%, respectively (Table 1); with the browsing intensity of P. jeffreyi significantly lower than the browsing intensity of P. contorta and P. ponderosa (Fig. 1d, results of post-hoc comparisons can be found as Supplementary Table S2.1 online).
In terms of the probability of a seedling having only the needles damaged, i.e., being defoliated, the probability estimated by the model for P. contorta, P. ponderosa, P. radiata and P. jeffreyi seedling was 10%, 30%, 1.7%, and 56% respectively (Table 1). These percentages differed significantly for the four species (Fig. 1e, results of posthoc comparisons can be found as Supplementary Table S2.1 online).
The probability of survival of a P. contorta, a P. ponderosa, a P. radiata and a P. jeffreyi seedling immediately after the treatment was 93%, 90%, 91% and 99% respectively (Table 1), with the probability of survival of P. jeffreyi Table 1. Results of the models performed to analyze sheep herbivory preference between four invasive pine species. We present the estimated value (Logit scale), standard error, z-value and their associated p value, and probability and lower and upper limits of the 95% confidence interval. Significant p values are highlighted in bold. www.nature.com/scientificreports/ significantly higher than the probability of survival of P. ponderosa and P. radiata seedlings ( Fig. 1f; results of post-hoc comparisons can be found as Supplementary

Anti-herbivores chemical compounds and physical traits.
The mean concentration of α-pinene in seedlings of P. contorta, P. ponderosa, P. radiata, and P. jeffreyi was 0.44, 0.98, 0.40 and 3.02 μL/g respectively; with the concentration of P jeffreyi seedlings significantly higher than the other three species, and the concentration of P. ponderosa seedlings significantly higher than the concentration of P. contorta and P. radiata seedlings (Kruskal-Wallis chi-squared = 61.177, df = 3, p value = 3.294e−13; Fig. 2a, results of post-hoc comparisons can be found as Supplementary Table S3.2 online). The mean foliar toughness of P. contorta, P. ponderosa, P. radiata and P. jeffreyi seedlings was 98.25, 173.25, 47.75 and 209 g/cm 2 respectively; with the toughness of P. ponderosa and P. jeffreyi significantly higher than the toughness of P. contorta and P. radiata seedlings (Kruskal-Wallis chi-squared = 67.654, df = 3, p value = 1.357e−14; Fig. 2b, results of post-hoc comparisons can be found as Supplementary Table S3.2 online).
The percent water content of P. contorta, P. ponderosa, P. radiata and P. jeffreyi seedlings was 0.56, 0.61, 0.58 and 0.46% respectively; with the percent water content of P. ponderosa significantly higher than the percent water content of P. contorta and P. jeffreyi and the percent water content of P. radiata significantly higher than the percent water content of P. jeffreyi (Kruskal-Wallis chi-squared = 54.388, df = 3, p value = 9.274e−12; Fig. 2c, results of post-hoc comparisons can be found as Supplementary Table S3.2 online).
The concentrations of total monoterpenes, β-pinene, 3-carene, β-phellandrene, and non-volatile resins differed between pine species (Fig. 2d, e, f, g, and h, respectively). Mean values and results of statistical analyses are detailed in Tables   Black point refers to the mean and black line to the median, whiskers correspond to SE. Different letters above boxplots indicate significant differences. Data of total monoterpenes, α-pinene β-pinene, β-phellandrene, 3-carene, and toughness was taken from the needles and data of non-volatile resins content and water content was taken from branches and needles together.

Relation between sheep preference and pines' anti-herbivores chemical compounds and physical traits.
Here, we report the general inferences made from parameters with p < 0.05. Details about final models, parameter estimates, significances, and confidence intervals can be found in Tables 4 and 5. Out of all anti-herbivores defenses analyzed, only α-pinene had a significant negative effect on the probability of a terminal bud being browsed ( Table 4). The concentrations of α-pinene and foliar toughness were negatively related to browsing intensity (Table 5). Conversely, water concentration was positively related to browsing intensity (Table 5).

Discussion
We show that all pine species were browsed by sheep; however, the magnitude and the pattern of damage differed. Regarding the browsing incidence, P. jeffreyi was the most preferred species, followed by P. contorta, and P. ponderosa, and lastly P. radiata. That is, sheep mostly preferred the least invasive and the most invasive species. But, regarding browsing intensity, sheep browsed more intensely the most invasive pine species. Since the four species were damaged, sheep browsing at this stocking rate could be affecting the growth rate and retarding seedlings from growing into adults 78 , slowing down the invasion rate, although not stopping it.
Regarding the chemical compounds, we observed that the total concentration of monoterpenes does not influence the attractiveness of pines to sheep, nor the subsequent damage as we did not find a significant effect on the probability of browse on terminal buds not browsing intensity. But looking specifically across monoterpenes, the concentration of α-pinene was the major deterrent to browsing. Therefore, in general terms, the concentration of total monoterpenes does not prevent the browsing of pines by sheep but the presence of certain types of terpenes could be conditioning their pattern of herbivory. The negative effects of α-pinene concentration on sheep browsing observed in our study were also observed in another study 79 . Estell et al. 79 performed an experiment assessing intake by sheep of alfalfa pellets treated with six chemical compounds (camphor, limonene, cis-jasmone, β-caryophyllene, borneol, or α-pinene), each one at five different concentrations. They observed a negative effect of chemical concentration on the intake of camphor and α-pinene. The general lack of anti-herbivore effect by monoterpenes concentration (except α-pinene) observed in our study could be due to the sheep domestic condition, probably individuals with high reaction to adverse stimuli were not selected during domestication [80][81][82] . Table 3. Mean concentration of chemical anti-herbivory compounds in four invasive species of pine and results of Kruskal-Wallis test.  Table 4. Results of models that assess the probability of a terminal bud being browsed in relation to chemical compounds and physical traits of pine seedlings. We present the estimated value (Logit scale), standard error, z-value and their associated p value, and the lower and upper limits of the 95% confidence interval. Significant p values are highlighted in bold. M1: Results of the model that assess the probability of a terminal bud being browsed in relation to total monoterpenes and total resins concentration. α is the baseline probability, β 1 is the regression coefficient that represents the effects of totals monoterpenes; β 2 is the regression coefficient that represents the effects of totals resin content. M2: Results of the model that assess the probability of a terminal bud being browsed in relation to α-pinene concentration. α is the baseline probability, β 1 is the regression coefficient that represents the effects of α-pinene concentration. M3: Results of the model that assess the probability of a terminal bud being browsed in relation to foliar toughness and water content. α is the baseline probability, β 1 is the regression coefficient that represents the effects of foliar toughness; β 2 is the regression coefficient that represents the effects of water content. www.nature.com/scientificreports/ We found that foliar toughness negatively influenced browsing intensity. The differential browsing observed between the four species (terminal bud, lateral branches or needle damage) may be due to the differences in morphology and/or degree of lignification of each species (all seedlings were younger than two years old, thus these differences were not due to age, but rather intrinsic differences of the species). Pinus jeffreyi and P. ponderosa seedlings presented the highest foliar toughness and also, they were the most defoliated. Cingolani et al. 56 measured plant functional traits and estimated sheep selectivity and plant grazing response, they observed that characteristics that indicate high structural quality for herbivores, between them low toughness, were in fact avoided suggesting that chemical defenses could be playing a more important role for this herbivore. In fact, Villalba and Provenza 83 posed the hypothesis that food structure and nutritional composition interact to determine the lambs' preferences, so, when a particular nutrient is needed, the biochemical composition is more important than the structure in determining the preference. This pattern was also observed by Evju et al. 84 in a study of sheep selectivity, they found that sheep select large, late-flowering herbs that have higher nutritional qualities (low leaf C/N ratio). In our study, P. jeffreyi and P. ponderosa, species with high browsing incidence, could have higher nutritional qualities for sheep, even when the physical structure was not the ideal for a grazer herbivore (high toughness). Regarding nutritional qualities, in our study, we only assessed seedlings' water content. We observed that this variable had a significant positive influence on browsing intensity. It would be important to complement this work with studies of the influence of other nutritional compounds on sheep's preference, like carbon, nitrogen, and phosphorus content.

Mean concentration per species
It is important to clarify that due to the well-known induced response of Pinaceae to herbivory 65,66 , the preference experiments and the anti-herbivory traits were measured in different plant material. Measuring these different variables in different plants might create differences that are not possible to account for. Even though all the seedlings came from the same source and were from the same pool, it is possible that other factors influenced the preference observed besides the traits assessed in our work, like differential microsite growth condition of each individual or morphological differences that were not possible to includes in the analyses, like plant height. However, we think that the performed analyses were the best possible approach to further our understanding of the relations between the anti-herbivory compounds and the sheep preference observed.
Some studies suggest that plant invasions could be controlled with mammal herbivory in other systems 28-30, 85, 86 . In Patagonia, Sarasola et al. 27 suggested that domestic ungulates could be one of the factors that are controlling pines invasions. However, in that study, the authors did not detail the species of ungulate, and for Patagonia, no empirical evidence about sheep herbivory on pines has been recorded to date. We observed that sheep damaged pine seedlings according to this decreasing gradient of preference: Pinus jeffreyi, Pinus contorta, Pinus ponderosa, and Pinus radiata. This suggests that the most invasive species, P. contorta, could potentially be controlled with sheep herbivory. In another study in New Zealand, Crozier and Ledgard 50 found that Pinus radiata and Pinus contorta were the most consumed species by sheep, followed by Pinus ponderosa. Our results correspond partially with this gradient of consumption since P. radiata was the least preferred species. In another study done in Puerto Patriada, Argentina, it was found that P. radiata individuals from different cohorts growing outside of plantations had varying non-volatile resin content of needles and branches (Dimarco personal communication). The difference in the sheep preference gradient found in our study and the one found by Crozier and Ledgard in New Zealand could be explained by the distinctive levels in chemical defenses that different genotypes of the same plant species could possess 46,87 , and also by the growing conditions of the pine seedlings used by Crozier and Ledgard and in our study. Table 5. Results of models that assess the browsing intensity in relation to chemical compounds and physical traits of pine seedlings. We present the estimated value (Logit scale), standard error, z-value and their associated p value, and the lower and upper limits of the 95% confidence interval. Significant p values are highlighted in bold. M1: Results of the model that assess the probability of browsing intensity in relation to total monoterpenes and total resins concentration. α is the baseline probability, β 1 is the regression coefficient that represents the effects of totals monoterpenes; β 2 is the regression coefficient that represents the effects of totals resin content. M2: Results of the model that assess the probability of browsing intensity in relation to α-pinene concentration. α is the baseline probability, β 1 is the regression coefficient that represents the effects of α-pinene concentration. M3: Results of the model that assess the probability of browsing intensity in relation to foliar toughness and water content. α is the baseline probability, β 1 is the regression coefficient that represents the effects of foliar toughness; β 2 is the regression coefficient that represents the effects of water content. www.nature.com/scientificreports/ Even though P. contorta was preferred by sheep in the present study, we observed in a previous study that the stocking rate implemented will determine the effectiveness of the invasion control 49 . The stocking rate must be estimated taking into account the palatability, the pine age, the regeneration rate, and the degree of invasion of the focal species. If the species to control are unpalatable, it could be necessary stock management concentrating high stocking rates in key seasons (i.e. spring) 88 . Regarding the age, previous studies observed that the concentration of anti-herbivores compounds changes with the ontogeny 89 . So, we suggest performing an herbivory invasion control over seedling younger than 2 years old 50 to guarantee high browse intensity of the seedling in order to avoid branch resprouting and ensure mortality. Lastly, Ledgard 48 , suggested that P. contorta is among the planted pine species with the highest regeneration rate. Thus, taking into account all the variables mentioned, for an effective herbivory control of P. contorta invasion (one of the most invasive species in the Southern Hemisphere) we suggest apply a high stocking rate in news fronts of invasion during short periods but with an annual frequency.
Our study suggests that sheep browsing could be a practical tool to control pine spread beyond the plantations. We observed that sheep consumed seedlings across our pine species at different intensities. Herbivory pine control could be more effective for pine species with low α-pinene concentration, low foliar toughness, and high nutritional qualities. Forestry practices accompanied by appropriate stock management plans could be an efficient method to control early pine invasions.

Data availability
The datasets generated during the current study are available from the corresponding author upon request.